Numerical simulation of the coffee-ring effect inside containers with time-dependent evaporation rate

In this work, using a mathematical model and numerical simulation, we investigate the effect of time-dependent evaporation rates on stripe formation inside containers. This pattern formation is driven by the coffee-ring effect. The coffee particles inside a container move according to random walk and under the gravitational force. Because of the time-dependent evaporation rate, we can observe stripe formation inside a container after evaporation of the coffee particle-laden liquid. Various numerical experiments are performed to demonstrate the proposed model can simulate the stripe formation in a container.


Introduction
When a drop of coffee evaporates on a solid substrate, the coffee particles form a ring-like stain at the contact line. This phenomenon is a consequence of faster evaporation at the contact line, driving the coffee particles to convectively migrate toward the contact lines, which is known as the coffee-ring effect [1,2]. Figure 1a, b shows coffee droplets on a transparent film on a grid paper before and after evaporation, respectively. After evaporation of coffee droplet, we can observe the coffee-ring pattern as shown in Fig. 1b.
The physical phenomenon related to the coffee-ring effect has been researched in several experiments with various colloidal suspensions containing nanoparticles [3], air bubbles [4], biological fluids [5], and bacteria [6]. Poulichet et al. [7] showed that the coffee-ring effect can be observed without evaporation when the solvent of a drop transfers into another liquid. There are many physical applications using the principle of coffee-ring effect such as ink-jet printing [8], making crystal morphologies [9], central spot formation in droplets [10], uniform droplet printing of graphene micro-rings [11], quantitative immunoassays [12], coffee-ring effect in polymer-nanocrystal mixtures [13], and micro-technologies in biological field [14]. For more detailed review on coffee-ring effect, see [15]. The coffee-ring effect has also been investigated in various numerical experiments based on the Monte Carlo method [16][17][18], finite element method [19], and finite difference method [20]. Figure 2a, b, c shows photographs of coffee solution contained in a conical flask, stripe pattern after evaporation, and magnified view of a part in (b), respectively. As we can see in Fig. 2b, c, there is a stripe pattern, i.e., the coffee-ring effect is alternatively strong (thick ring) and weak (thin ring) in the vertical direction. This phenomenon happens because of time-dependent evaporation rate in room temperature. During 24 h in a day, evaporation rate keeps changing due to temperature, humidity, wind speed, and surface area, for example.
Communicated by Peter Duck.  Most of the previous numerical studies were performed to investigate the coffee-ring effect of evaporating drops on the substrate domain. Motivated by the stripe pattern formation phenomenon with time-dependent evaporation, the objective of the present study is to present a mathematical model which can simulate the stripe pattern formation driven by coffee-ring effect and time-dependent evaporation rate. This paper is organized as follows. In Sect. 2, we present the governing equation and its numerical solution algorithm. The numerical results are presented in Sect. 3. Finally, conclusions are drawn in Sect. 4.

Governing equation and numerical solution algorithm
Let Ω t be a time-dependent liquid domain at time t in a solid container. The time-dependent liquid domain Ω t is defined by a level-set function (a signed distance function) φ(x, t), i.e., Ω t = {x ∈ R 3 |φ(x, t) ≥ 0}.
Let Γ t = x | x = (x, y, z * ) ∈ Ω t and z * = max (x,y,z)∈Ω t z be the liquid surface at time t contacting the air.
The liquid domain contains the coffee particles which are defined as the Lagrangian points moving under Brownian motion and gravitational force. Figure 3a, b shows a time-dependent liquid domain Ω t including coffee particles and a top-level surface of liquid domain Γ t with the scalar field φ at time t = 0 and t > 0, respectively.
For the sake of simplicity of modeling, we assume coffee particles have only mass and no volume and the fluid surface moves under the time-dependent evaporation rates. Assuming that the force on the particle is where X is Lagrangian variable for a particle, t is time, g = (0, 0, −g) is the scaled gravitational force. Here, ψ(X) = (ρ cos(θ ) sin(ϕ), ρ sin(θ ) sin(ϕ), ρ cos(ϕ)) represents random force of Brownian particle, where ρ is a random variable with the probability density function f (ρ) = e −0.5ρ 2 / √ 2π and between [−3, 3] which is 99.998% confidence interval. Therefore, we can safely restrict ρ within the interval [−3, 3]. The random variables θ and ϕ are uniformly distributed on the interval [0, π]. Using the explicit Euler's method, we discretize Eq. (1) as where X n k is the k-th particle position at time t = nΔt, Δt is the time step, and N k is the total number of particles. Equation (2) can be rewritten as , Z n+1 k ) be the n + 1 time-level and k-th particle position computed by using the evolution equation (3). A particle becomes pinned if it satisfies the following three conditions: If a particle moves outside the domain, then it goes back to fluid boundary, see the schematic representation of this process in Fig. 4. The pinning position of a pinned particle is defined as the closest point in Γ t from X n+1 k as illustrated with a diamond symbol. If a particle without becoming pinned goes outside the liquid domain (φ < 0), then the particle moves to the closest point of φ = 0, see the square and triangle symbols.

Numerical experiments
Unless otherwise stated, z = S(t) is the position of the liquid surface at time t, the total number of particles N k is set to 5 × 10 5 , N p is the number of pinned particles, the gravitational force is fixed to g = (0, 0, −1), and the size of bins in the histograms is set to 0.03. We first perform the numerical simulation for the motion of particles without evaporation in a cylindrical container. The static liquid domain is defined as follows: Then, we can define a signed distance function as where ∂Ω t is the boundary of the domain Ω t and dist(x, ∂Ω t ) = min y∈∂Ω t |x − y| is the distance function between x and ∂Ω t . The temporal step size Δt = 0.001 and α = 0.5 are used. The total number of particles is N k = 10 5 . Because there is no evaporation, we assume there is no pinning of particles. The evolution of particles in a cylindrical container without evaporation can be seen in Fig. 5a. The number of particles with respect to height is shown in Fig. 5b. As can be seen from the histogram in Fig. 5b, particles move downward due to the gravitational force.

Coffee-ring formation with constant evaporation rate
Now, let us consider coffee-ring formation with constant evaporation rate. The time-dependent liquid domain is defined as follows: where S(t) = 3 − 1.25t is the position of the liquid surface at time t and schematically illustrated in Fig. 4. Note that the evaporation rate is defined as dS(t)/dt, which is the rate of length with respect to time. Then, we can define a time-dependent signed distance function φ(x, t) as described in Eq. (5). We use Δt = 0.001 and α = 1, 1.5, 2. Figure 6a shows the evolution of particles in a cylindrical container using α = 1. The numbers of pinned particles are shown in Fig. 6b, c, d using α = 2, 1.5, and 1, respectively. We can observe the more pinned particles because the particles move farther if α is larger .

Coffee-ring formation with changing evaporation rate
Using the proposed method, we perform the numerical simulation for the stripe pattern formation in a cylindrical container with the time-dependent evaporation rate. The time-dependent liquid domain is defined as follows:  where S(t) is a given time-dependent function and S(0) = 3. To simplicity of exposition, let us consider the following time-dependent evaporation rate using the sine function: where A is the mean evaporation rate, β is the perturbation amplitude, and γ is the frequency. Equation (8) models effectively the time-dependent evaporation rate. We should note that we can use other forms of the time-dependent evaporation rate. The solution of Eq. (8) with the initial condition S(0) = 3 is given as The parameters used are Δt = 0.0001, α = 1, A = 1.25, β = 3, and γ = 6. The first (a) and second (b) rows of Fig. 7 show the evolution of particles with N k = 10 5 and N k = 5 × 10 5 , respectively. Figure 8a, b shows the histogram of the pinned particles at time t = 24000Δt with N k = 10 5 and N k = 5 × 10 5 , respectively. From this numerical experiment, we can clearly observe that stripe patterns are formed for both the cases.

Effect of container shape
To investigate the effect of container shape, we perform a numerical test for stripe pattern formation in a conical container with time-dependent evaporation rate using the proposed method. The time-dependent liquid domain is defined as where S(t) = −At + β(cos(2πγ t) − 1)/(2πγ ) + 3. The parameters used are Δt = 0.001, α = 1, A = 1.25, β = 3, and γ = 6. Figure 9a shows the temporal evolution of particles with the time-dependent evaporation rate in a conical container. Figure 9b is the histogram of particles with respect to height at final time t = 2365Δt. We can observe there is a periodic pattern of pinned particles from z = 1 to z = 3 even though the container shape is not cylindrical. The fluid surface area is getting smaller; however, there are more particles around contact lines between the fluid and container.
The gravitational force can be written as g = (g sin(θ ), 0, −g cos(θ )). Here, for computational simplicity, we apply the tilted gravitational force instead of tilting the cylinder domain. Two different gravitational parameters are used as g = 1 and 2. Figure 10a, b shows the snapshots of evolution at time t = 0, 200Δt, and 400Δt with respect to g = 1 and 2, respectively. As we can observe, a smaller gravitational parameter causes the coffee particles to settle more slowly and more particles are pinned on the side of wall boundary. On the other hand, a larger value of gravitational parameter causes the coffee particles settle faster and the number of pinned particles on both sides of wall boundary is relatively small. Around the bottom of the container, we can also see that the coffee particles are more gathered in the tilted direction under the influence of gravity.

Effect of tilt angle
Finally, we perform the numerical simulation to study the formation of coffee ring with respect to tilt angle. We consider two different tilting angles: 10 • and 55 • . The time-dependent liquid domain is as follows.
φ(x, y, z, t) = 1 if x 2 + y 2 ≤ 1 and 0 ≤ z ≤ S(t), −1 otherwise, where S(t) = −At + β(cos(2πγ t) − 1)/(2πγ ) + 2. For 10 • (θ = π/18 rad) and 55 • (θ = 11π/36 rad), the gravitational force is set to g = (g sin(θ ), 0, −g cos(θ )) and the other parameter values are the same as in   Figure 11a, b shows the temporal evolution of particles with a time-dependent evaporation rate in 10 • and 55 • tilted cylinders, respectively. Figure 12a, b shows the histograms of the numbers of pinned particles with respect to height at final time t = 400Δt in 10 • and 55 • tilted cylinders, respectively. Here, the number of bins in the histograms is set to 40. The left and right histograms are the numbers of pinned particles in the area x < 0 and x ≥ 0, respectively. As can be clearly observed in the histograms of Fig. 12a, b, as the angle of inclination increases, the number of particles pinned on the right contact wall increases. This phenomenon happens because coffee particles rest on the right contact wall boundary more than the right contact wall boundary.

Conclusions
We proposed the mathematical model for the stripe formation in a cup of coffee with a time-dependent evaporation rate. Using the Monte Carlo method, the coffee particles are defined as Lagrangian points moving under Brownian dynamics and gravity. Temporal evolution of the governing equation is solved using the explicit Euler's method. The numerical simulations were performed in the three-dimensional space. From the numerical results, we investigated that the distribution of particles and stripe pattern formations depend on the evaporation rates, the shape of container, the gravitational force, and the tilt of the cylinder. In a future work, we will use a pinning boundary condition [21] to extend the current study to investigate the hydrodynamic effect on the coffee-ring formation.