Average cluster size inside sediment left after droplet desiccation

In this work, we continue to study the formation of particle chains (clusters) inside the annular sediment during the drying of a colloidal droplet on a substrate. The average value of the cluster size was determined after processing experimental data from other authors. We performed a series of calculations and found the value of the model parameter allowed to get numerical results agreed with the experiment. Also, a modification of the previously proposed algorithm is analyzed here.


Introduction
Evaporation-induced self-assembly (EISA) methods are useful in applications such as inkjet printing, production of photonic crystals, application of transparent conductive coatings, development of biosensors, and others [1]. Modeling the formation of colloidal particle precipitation structures is important. It allows us to understand the main mechanisms of such processes, identify key parameters, and study how to influence the system to obtain the required patterns. For example, the model [2] describes the self-assembly of particles on the free surface of a droplet. The capillary interaction of particles and their transport by surface flow resulting from evaporation is considered in this paper. The authors [2] performed calculations for spherical and triangle droplets using the lattice Boltzmann method (LBM). Numerical results showed the dependence of the sediment on the curvature of the free surface. Clusters of particles formed at the liquid-vapor interface appear on the substrate after complete evaporation of the liquid. Another mechanism for the formation of particle clusters as a result of their capillary attraction is their self-assembly on the substrate where the thickness of the liquid layer is less than the particle size [3]. The calculations were performed using the modification of LBM. The feature of the algorithm is the ability to simulate the dynamics of soft particles [3]. The mass transfer of hard particles in a drying film was studied in [4]. The LBM-based model takes into account the hydrodynamic interaction of particles and their diffusion. The authors [4] performed calculations for different values of the Peclet number associated with different evaporation rates. At high values of the Peclet number, a layer of particles is formed near the free surface of the film due to the rapid movement of the two-phase boundary as a result of intense evaporation. A uniform distribution of particles is observed at a small value of the Peclet number. It is explained by diffusion transfer. The periodic predominance of capillary flow or solutal Marangoni flow sometimes leads to the formation of concentric precipitation of colloidal particles [5]. The authors [5] used the lattice model and the Monte Carlo method to study the effect of the surfactant concentration on this process. The magnetic interaction of particles in a liquid and the formation of chains was studied in [6]. The model is based on the Soft Sphere Discrete Element Method. At each time step, Newton's equations were solved to predict the motion of particles. Calculations were performed for small viscosity values [6]. The model [7] takes into account the diffusion and sedimentation of particles. These particles have no volume, but they have mass. The model is based on the phase field method and the Monte Carlo method. Two shapes of the droplet surface were considered: a spherical segment and an asymmetric shape from the experiment [8]. In both cases, a large number of deposited particles per unit area near the periphery is shown (the coffee stain effect). The model predicts a significant accumulation of particles in areas of low curvature for an asymmetric droplet shape. In the experiment [8], the particles accumulated mainly in areas of high curvature. The reason for this is the capillary flow, which is not taken into account in the model [7]. The formation of particle chains (clusters) inside the annular sediment during the evaporation of a colloidal droplet on a substrate was described in [9]. The model takes into account diffusion, advection, and capillary attraction of particles. Capillary interaction of particles occurs near the boundary (fixation radius R f ) where the thickness of the liquid layer is comparable to the particle size. This process was mimicked as follows. In one time step, the particle was shifted to the nearest neighboring particle in the small neighborhood R n at the presence of the last one. Calculations were performed for only one parameter value R n . In the current work, we consider the set of parameter R n values and determine the value which leads the average cluster size is similar to the experiment [10]. Also, here we analyze the modification of the algorithm associated with multiple attempts to shift those moving particles whose location has not changed in a time step as a result of a collision with a neighboring particle.

Problem statement and model description
Here we consider an axisymmetric sessile droplet on a substrate in the mode of a pinned threephase boundary. Colloidal particles are suspended in the liquid. They are subject to transfer by flow, diffusion, and capillary attraction. At the initial time, the contact angle θ 0 = π/18. The contact radius of the droplet with the substrate R = 60 µm. The radius of spherical particles r p = 0.35 µm. As the liquid dries, these particles precipitate. The deposit is a monolayer of particles since the droplet is thin (its height h R). Time of full droplet evaporation t max = 10 s. A more detailed description of the problem statement is given in [9].
The mathematical model [9] is semi-discrete. Each particle is considered separately, but the liquid is described with a continuum approach. It is assumed that the particle velocity is equal to the fluid flow velocity, which is calculated by the formulā It was derived from the mass conservation law (references and formula derivation are given in [9]). The model considers the radial flow averaged over the drop height. Here r and t is the radial coordinate of the particle center and the current time of the process, respectively. Coordinate of the particle at the next time step as a result of displacement by the flow where δt is the time step and τ is the time step number. Capillary forces act on the particles in the region of the fixation boundary R f (t) = R 2 − 4r p R/θ(t). The time dependence of the contact angle is written as θ(t) = θ 0 (1 − t/t max ). We mimic the capillary attraction of particles as follows. If there are particles in the vicinity of the current particle that have precipitated or are also subject to capillary attraction, then it moves close to the nearest neighboring particle. By neighborhood, we mean a circle with a radius R n . Some particles precipitate if r > R f + r p . In this case, the local thickness of the liquid layer is less than the particle size, h < 2r p .
The diffusion motion of particles is modeled by the Monte Carlo method. A random offset angle α ∈ [−π; π) is generated at each time step τ . The new position of the particle is calculated using the formulas x τ +1 = x τ + cos(α) √ 2Dδt and y τ +1 = y τ + sin(α) √ 2Dδt, where x and y are the Cartesian coordinates of the particle center. The diffusion coefficient is calculated using the Einstein formula D = kT /(6πηr p ), where k is the Boltzmann constant. The values of the temperature T and the viscosity η of the liquid are taken for water under normal room conditions (D ≈ 6 · 10 −13 m 2 /s). The value of the time step δt = 0.1 ms was selected based on a series of computational experiments in such a way that the Einstein relation for the mean square displacement l 2 = 2Dt max was fulfilled. Here and N p is the number of particles.

Algorithm
For each particle, we store coordinates and status in memory. Let's denote movable particles with green and black colors. Here green particles are subject to convective and diffusive transfer. We do not use a separate status for particles diffusing near R f to simplify the algorithm [9]. Capillary forces act on black particles. Unmovable particles that have precipitated are marked in red. The particle status changes according to the rules in Figure 1 (right).

Figure 1.
Scheme to the algorithm description (left) and rules for changing the particle status (right).
For more details about the algorithm, see Figure 1 (left) and pseudocode 1. At the beginning of the process, all particles are green by default. If a green particle touches the fixation radius, it is repainted black. The green particle cannot go beyond the fixation radius, unlike the black one. After all, the movement in that area is mainly due to capillary forces. But due to the movement of the fixation radius itself, the green particle may be behind it (after that, it will be repainted red). A black particle becomes red after it is beyond the moving boundary R f as a result of the motion of the last one or after it shifts due to capillary forces. while (not all particles are displaced) & (there is a shift of at least one particle) do 10: shuffle the array of particle numbers.

11:
for i ← 1, N p do 12: red particles are skipped. 13: if (green particle) then 14: calculate the new particle coordinates due to diffusion.

17:
end if 18: Calculate the new particle coordinates due to advection.

19:
if (no collision) then 20: move the particle. if (the particle is black) & (there is a black or red particle within R n ) then 24: Calculate the new particle coordinates due to the capillary shift.

25:
if (no collision) then 26: move the particle close to the nearest one. Write the particle coordinates and colors to a file for the current time step. 32: end for At each time step, we calculate the new value of R f . Then we go through all the non-red particles and check the conditions for changing their status. Next, we shuffle the numbers of particles so that their subsequent search is random. After that, we repeatedly iterate over the moving particles that have not yet moved at the current time step and try to move each of them in turn. Failure to shift can be caused by a collision (neighbor particle overlapping) or the exit of a green particle beyond the R f boundary as a result of a diffusive motion. Repeated iterations are completed if all the moving particles were displaced at the current time step or there were no displacements at all during a repeat of the loop. To check that there is no collision of a particle with any other, it is advisable to consider the distance only to the nearest neighbors. To do this, we have to store additionally and periodically update information about subdomains containing particles (numbers of particles included in the subdomain).

Determination of cluster sizes
To determine the value of the average cluster size, an experimental picture was taken from [10] for further processing. The data we are interested in was extracted as follows. The image was imported into the bitmap editor GIMP. A transparent layer was applied on top of it, which we worked with later. The brush size was set to match the size of the particles in the image [10]. All particles not included in the annular sediment were drawn on the transparent layer manually on top of the experimental image. At the same time, we tried to choose different colors for each cluster for the convenience of further calculation.
Algorithm 2 Algorithm for finding clusters.
1: Read the coordinates of the particles x i , y i and their cluster numbers if r > r * then 5: remove this particle from the array.
if particles i and j are not from one cluster then 33: We assign a value max(c i , c j ) to the particle with min(c i , c j ). Sometimes it was difficult to visually determine whether a particular particle belongs to a particular cluster or not. In the process, we were guided by the following rule. If a particle is located close to another particle in the cluster, or the distance between these particles is much smaller than the particle size, then we assume that the current particle also belongs to this cluster. As a result, the layer with the original image was removed. Based on the received image, we manually calculated the average cluster size.
Processing of numerical simulation results was performed using the algorithm 2 for finding clusters. To determine the contact of two particles, the condition s = d p ± εd p was used, where s is the distance between two particles, the diameter of a spherical particle d p = 2r p , and ε is the constant (ε = 0.01 is used here). This nuance was also taken into account when detecting a collision of two particles (s < d p − εd p ). The calculation of the average cluster size based on the results of ten repeated tests was performed using a script written in Python.

Results and discussions
A series of calculations were performed for different values of R n in the range from 3r p to 22r p . Then the numerical result was approximated using the least-squares method (Figure 2). The plot in Figure 2 resembles a sigmoid. The value of the average cluster size determined by us based on experimental data [10] is N c ≈ 8.2 (average number of particles in the chain). The model predicts this value for R n ≈ 7.9r p .   Figure 2. Dependence of the average cluster size on the parameter R n (the standard deviation was used to estimate the error).
For comparison, the Figure 3 shows the final sediment structures based on experimental data [10] and our simulation results. Only the particles located inside the annular sediment are shown. The subset of particles belonging to the ring was subtracted from their set using the condition r > r * . The inner radius of the ring r * ≈ 0.75R was determined based on calculation data [9] (r * corresponds to the place of a high gradient of the number of particles per unit area). The internal structure of such rings was discussed in detail earlier [9].
Here we observe that for the given parameters, the model predicts the absence of single particles (Figure 3). The visualization of experimental data shows a very small number of such particles. In both cases, there are clusters of different sizes (from two to several dozen particles in a cluster). Some clusters are similar in shape. It should be noted that based on the results of modeling, tree-shaped clusters are expected to appear. Their root element is located closer to the periphery. These "trees" grow towards the center of the area. It is explained by the nature of movement and the shape of the fixation boundary R f . Such branching structures are not observed in the experiment. This distinctive feature of the numerical results is most likely related to the approximate description of the domain geometry [9]. After all, the shape of the drop at the final stage of evaporation is more like a film, in which local ruptures can occur [1]. Thus, the action of capillary forces occurs in a larger area. The model boundary R f and its neighborhood are only a rough approximation of such a region. Besides, the mimic of the capillary interaction of particles is used here. In a good way, it is necessary to solve the n-body problem numerically. Only in this problem, instead of gravity, there are capillary forces. Visualisation of clusters from experiment [10] (left) and simulation results for R n = 8r p (right).
Analysis of the algorithm 1 showed a weak nonlinearity of the dependence of the simulation time on the number of particles N p (Figure 4). Also, it was found that the average (over time steps) rate of movable particle displacement failures for the modified algorithm (with multiple attempts to particle displacement) is 1-2% according to the results of calculations. For comparison, we tested the original algorithm (with a single attempt to shift the particles) and got a failure rate of 3-4%. Five tests were run simultaneously for each type of calculation (the error is less than the marker size on the Figure 4). Thus, for a sufficiently small value of the time step, the error associated with the failure of the particle shift is small. Repeated attempts to shift particles can slow down the program and not give a significant difference in the result.   Figure 4. Dependence of the simulation time on the particle number (the program was written in Python, Intel (R) Core (TM) i9-9900K CPU 3.6 GHz was used).

Conclusion
Numerical results [9] are qualitatively agreed with the experiment [10]. As the droplet dries, chains of particles (clusters) are formed inside the annular sediment. In the current work, we have processed experimental data [10], found out the average cluster size, and defined the value of a model parameter R n to predict a quantitatively agreed value N c . But it is also worth mentioning the qualitative difference between the simulation and experimental results, related to the shape of some clusters. The simple model under consideration often predicts tree-shaped clusters, which is not typical for the experiment [10]. Further study is needed in this direction to better understand the processes that are taking place. Some questions can be answered by using the high-speed shooting of the formation of such clusters. Additional experimental data are also needed, for example, on the shape of the free surface of the liquid at the time stage when these clusters begin to appear. Further creation of more complex and accurate models that can predict the internal structure of such precipitation is important.