#### Abstract

The deposition of colloidal silica particles during the evaporation of sessile drops on a smooth substrate has been modeled by the simultaneous solution of the Navier–Stokes equations, the convective-diffusive equation for particles, and the diffusion equation for evaporated vapor in the gas phase. Isothermal conditions were assumed. A mapping was created to show the conditions for various deposition patterns for very dilute suspensions. Based on values of the Peclet (Pe) number and Damkholer numbers (Da and Da_{−1}), the effects of adsorption and desorption were discussed according to the map. Simulations were also done for suspensions with a high particle concentration to form a solid phase during the evaporation by using a packing criterion. The simulations predicted the height and width of the ring deposit near the contact line, and the results compared favorably to experimental particle deposition patterns.

#### 1. Introduction

Using the evaporating sessile drop could be a simple way to make uniform patterns for printing sensors and pharmaceutical drugs and applying novel film coatings [1–4]. The deposition pattern of colloidal silica particles during the evaporation of sessile drops on smooth substrates has been investigated experimentally and computationally. Deegan et al. proposed that the capillary flow inside an evaporating droplet carries the solute radially to the contact line, which forms a ring-patterned deposit [5]. This ring pattern was also observed by other groups for the deposition of polystyrene latex colloids on glass substrates [6, 7]. While some studies focus on the influence of microfluidic flow on the deposition [8, 9], the deposition patterns can also vary from a ring to a uniform layer according to the substrate roughness [10], particle morphology [11], the evaporation speed [12] and mode [13], and the composition of the solvent [14]. Sommer compared the deposition pattern in an evaporating sessile drop between polystyrene latex particles and the hydroxyapatite particles [15]. While both are depositing onto a titanium disk, a ring formed in the polystyrene-laden solution while a uniform layer was obtained from the hydroxyapatite solution.

Bhardwaj et al. studied the role of the Derjaguin–Landau–Verwey–Overbeek (DLVO) interactions on the final deposition patterns [16]. The particle (titania) surface potential and the substrate (glass) surface potential were measured under different pH values ranging from 1.4 to 11.7 by adding HCL (0.1 M) or NaOH (0.01 M). The competition between the hydrodynamic radial flow and the particle/surface interaction has a significant effect on the final deposition. If the attractive force between the particles and the substrate is large enough and the Peclet number is low enough, the particles tend to stick to the surface before they are swept to the contact line by the radial flow and finally form a uniform deposition pattern.

To calculate the particle deposition growth, Deegan found that the mass of the ring grows versus time based on the power law by assuming depth-average velocity and particle concentration, which is valid only in very thin droplets [5, 17]. Fischer solved the particle deposition profiles under three types of evaporation flux distribution along the droplet surface using lubrication theory, also assuming a very thin droplet and neglecting particle diffusion [18]. Hu and Larson used Brownian dynamics simulation to calculate the particle deposition profile based on the analytical form of the velocity in both radial and vertical components [19].

The adsorption and desorption of particles to the substrate can be modeled as a first-order reaction [20]. Bhardwaj et al. solved the vapor concentration, fluid profile, and particle concentration distribution inside the droplet using the finite element method [21]. The depinning phenomena of the droplet were also considered. In their simulation, once the particle concentration reaches the maximum packing *X* = 0.7, the particle deposit starts growing and the depinning happens. Widjaja and Harris solved the vapor concentration, fluid profile, and particle concentration distribution inside the droplet simultaneously by solving the Stokes equations of fluid motion and the convective-diffusive equations using the finite element method [22]. The particle deposition profile was calculated assuming that the flux of the particle deposition follows a first-order reaction. A phase diagram of the final deposition profile under different Peclet (Pe) and Damkholer (Da) numbers was obtained.

Desorption of particles can happen since there are flows inside the drop that can move the particles that have deposited. The fact that flows facilitate desorption of particles from the fluid-solid interface has been revealed previously [23, 24]. Yet limited research has incorporated the desorption behavior of particles to the studies of drop evaporation. The lack of desorption in the model may lead to incorrect prediction of the deposition profiles.

In this paper, the vapor, fluid, and particle phases were also solved simultaneously using the finite element technique. In addition to solving the particle concentration in the drop phase, the equation for particle deposition at the fluid/substrate interface was included, and the deposition profile was calculated from the interface concentration. The addition of the interface concentration makes it possible to include not only the adsorption effect but the desorption effect. The results are interpreted in two parts. First, a parametric study is performed with different convection, adsorption, and desorption effects, i.e., final deposition profiles were obtained under different Pe, Da, and Da_{−1} (desorption Damkohler number). Second, the maximum packing criterion was applied to account for the formation of particle accumulation near the contact line. Different from Bhardwaj’s setup [21], this paper works on a drop whose contact line does not move, like a water drop with suspended silica particles on a glass substrate.

#### 2. Numerical Simulation

##### 2.1. System Description

The finite element technique is used to solve the equation governing the vapor concentration, fluid velocity, and the particle concentration inside the droplet. The droplet was assumed to be axisymmetric, and the fluid properties such as viscosity, density, and the surface tension are constant throughout the entire evaporation process. While the volume of the droplet is decreasing, the contact line stays pinned, and the simulations were stopped when the contact angle decreased to about 5°. The governing equations and boundary conditions for each of the three phases will be shown in the following section.

##### 2.2. Governing Equations and Boundary Conditions

The system contains three phases: vapor, liquid, and solid (particle) phases. The evaporation of the droplet is assumed to undergo a pseudosteady state evaporation process, and the vapor is stagnant; therefore, the governing equation simplifies to the Laplace equation as follows (dimensionless form):

The boundary conditions for the vapor phase are as follows:where *c* is the dimensionless vapor concentration, *r* and *z* are the dimensionless radial and height distance from the center of the drop base, *H* is the relative humidity, and **n** is the unit normal vector. Explanation of the boundary condition (3) is in the SI. The characteristic length for the entire system was chosen to be the drop base radius *R*_{c}; the characteristic vapor concentration is the saturated vapor concentration of the liquid in the evaporating droplet. Once the vapor concentration is solved, the dimensionless evaporation flux can be calculated according to Fick’s Law, , where **n**_{s} is the unit normal vector at the drop-air interface.

To solve the fluid flow inside the evaporating droplet, we assume the drop to be incompressible and Newtonian and solve the continuity equation (5) and the Stokes equation (6). Here, we consider the case where the droplet is small enough so that the shape of the drop is part of the spherical cap where gravity has a negligible effect on the drop shape. The fluid properties (e.g., density, viscosity, and surface tension) are assumed to be constants during the evaporation process. The contact line is assumed to stay pinned until the contact angle is about 5°. The governing equations and the boundary conditions are listed as follows (dimensionless form):where **V** is the dimensionless velocity vector, .

The boundary conditions are (dimensionless form)where is the nodal velocity along the drop surface, 2*H*_{S} is twice the local mean curvature of the surface, and Ca is the capillary number, which is defined by Ca = *μv*_{c}/*σ*, where *μ* is the fluid viscosity and *σ* is the drop-air interface surface tension. is the characteristic velocity which is defined by *J*_{c}*/ρ* where *ρ* is the liquid density and *J*_{c} is characteristic evaporation flux defined by where is the diffusivity of water/air and *c*_{0} is the saturated vapor concentration. For water diffusing into air, and the saturation concentration for water, .

The convective-diffusive equation was used to solve the particle concentration inside the droplet as follows (dimensionless form):where the Peclet number, is the diffusivity of the silica particles in water with the value of 4.37 × 10^{−12}. The characteristic particle concentration is *C*_{P,0}, which is the initial particle concentration and is assumed to be initially uniform in the drop phase. With this model, the particles are assumed small enough so that Brownian forces dominate over gravitational forces, and thus, gravitational effects on the particles are negligible.

We assumed that there are no particles flowing across the axisymmetric line. The boundary condition along the droplet surface is the particle mass conservation at the interface during the evaporation, which is

The boundary condition on the substrate is the deposition of the particles onto the substrate, in which the deposition flux *J*_{p} is the sum of the adsorption and desorption [25] as follows:

Here, is the Damkholer number which is a ratio of particle adsorption over particle diffusion; is the desorption Damkholer number which is a ratio of particle desorption over particle convection, where is the characteristic time. *k* and *k*_{−1} are the adsorption and desorption rate constants. It is assumed that the rates of adsorption/desorption of particles onto the substrate are the same as the particle adsorption/desorption onto particles that have already deposited, and they both take the form of a first-order reaction. Details on the derivation and explanation of this equation are discussed in the SI. Γ_{S} is the dimensionless surface particle concentration and is also one of the variables to be solved. The governing equation to solve along the nodes on the substrate is as follows (dimensionless form):

Equations (8)∼(11) were used to solve for *C*_{P} and Γ_{S}. In this study, we investigated the effect of the dimensionless groups Pe, Da, and Da_{−1} on particle final deposition patterns. In order to compare the final deposition patterns for different Pe numbers, the final deposition will be plotted as Γ_{S}/Pe versus position *r*.

Furthermore, the case with high particle concentrations was simulated. When particle concentration is high, packing of particles happen near the contact line, as what will be discussed in Section 3.3. Therefore, the contact line region was treated as follows.

The particle concentration for each element was checked at each time step to see if the critical concentration *C*_{P,max} is achieved. If the concentration in any element achieves the critical concentration, the convective-diffusive equation in this “solid” element will not be solved anymore, and the nodal position was “fixed.”

Once the “solid” elements are determined, the *C*_{P} values on the nodes inside the solid element will be equal to the critical value, noted as *C*_{P,max} = 130.86. The common nodes between the solid element and the neighboring fluid elements are fixed and are applied under no-flux boundary condition of particles as in equation (10) in the Method section.

We assumed that the particles in the “solid” elements are still wetted by the liquid. The fluid dynamics were still being solved, i.e., evaporation is still considered on the surface of the solid phase and drives the liquid in the drop phase to flow toward the contact line. A sketch is shown in the SI. In the solid phase, the equations of fluid dynamics keep the same as the liquid phase. A “superficial velocity” resulting from the fluid dynamics equations describes the averaged velocity in the solid phase.

##### 2.3. Code Validation

After the vapor concentration profile was calculated, the evaporation flux along the droplet surface and the total evaporation rate were computed. The total evaporation rate as a function of the contact angle computed from our FEM algorithm was compared to the analytical results from Picknett and Bexon [26] as shown in Figure 1.

Figure 1 shows that our evaporation rates overlap Picknett and Bexon’s curve over a contact angle of 5° to 77°. The mesh has to be very fine close to the drop surface for the entire evaporation period in order to reach the overlap. The largest discrepancy between the total evaporation rate from our simulation and Picknett and Bexon’s theoretical value is about 0.43%, which happens at the beginning of the evaporation process when the contact angle is the largest.

The change in the drop mass during the evaporation was calculated in two ways. The first method is the integration of the volume times the density at each time step, and the other method is to subtract the sum of the total evaporation rate at each time interval from the initial mass. The overlapping of these two mass versus time verifies the correctness of the fluid velocities and mesh movement calculation in the simulation, as shown in Figure 2. Also, in Figure 2, the results for overlap with the results for , since Ca is the ratio of viscous stress and capillary stress and thus is expected not to affect the evaporation speed. The total evaporation time is around 0.3, which is the same as Widjaja’s results [27].

The streamlines inside the droplet phase are shown in Figure 3, and the evaporation flux distribution along the drop/vapor interface is shown in Figure 4. The mesh is refined enough both in the droplet phase and the vapor phase, which is indicated by the details of the velocity profile on the top and at the contact angle, as shown in Figure 5.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

The way to verify the particle concentration calculation is to check the mass conservation of the particles inside the droplet during the evaporation. The particle mass in the drop phase was calculated according to the following equation (12) (dimensionless form):

The mass accumulation on the surface was calculated from equation (13) (dimensionless form) as follows:

The sum of these two should be a constant throughout the entire evaporation process since there is no particle loss across the drop/vapor interface.

#### 3. Simulation Results and Discussion

##### 3.1. Particle Deposition without Desorption

An evolution of how the particles deposit for the case Pe = 1 and Da = 10 but without desorption, i.e., Da_{−1} = 0, is shown in Figure 6. The deposition profile stayed almost the same after *t* = 0.276 because most of the particles had already deposited onto the substrate before this time, as is also shown in Figure 7. The final deposition pattern is the same as the results from Widjaja and Harris in 2008. Under this condition, the deposition profile grows from a uniform pattern at the earlier stage to a somewhat center deposition in the later stage.

As mentioned in Section 2.3, Figure 7 is also a code validation for the particle deposition calculation. The solid curve is the mass accumulation of the particles on the substrate, while the dashed curve is the particle mass left in the drop phase that has not deposited. The sum of these two (the top line) is a constant throughout the entire evaporation process.

Figure 8 shows the “phase diagram,” which is similar to Widjaja and Harris’s results [22] but used different initial aspect ratio and humidity condition. When Pe is small, particle transportation to the contact line is weak and particles tend to stay in the middle of the drop. When Da is large, particles deposit on the substrate quickly, and the deposition reduces the chances for the particles to be convected to the contact line. The center-like pattern happens when Pe is small (around 0.01∼5) and Da/Pe is large (around 5 d∼100). When Pe is high (>10), the edge (or near-edge) deposition happens for all the Da numbers we have used.

##### 3.2. Particle Deposition with Desorption

When the particles that have adsorbed onto the substrate can desorb from the substrate, the final deposition and the deposition rate are both affected. Figure 9(a) shows one of the data points in Figure 8, where desorption does not happen. Figures 9(b)–9(d) are under the same Pe and Da numbers but with different extent of desorption.

**(a)**

**(b)**

**(c)**

**(d)**

When Da_{−1} is small, the deposition pattern is still a center deposition, and the particle deposition rate is similar to the case where there is no desorption. With an increasing Da_{−1}, the edge of the deposition pattern becomes higher, and the particle deposition rate becomes slower. The edge of the deposition is higher (the center deposition transitions to an edge deposition). When Da_{−1} is as large as 25, there is still a significant number of particles left suspended in the drop phase as the contact angle approaches 5°, i.e., termination of the simulation.

Increasing Da_{−1} results in a change to an edge deposition from a “near-edge” deposition as well as from a center-type deposition. According to the phase diagram in Figure 8, when there is no desorption, a center-type deposition pattern is obtained for Pe = 1 and Da = 5, and a “near-edge” deposition pattern forms for Pe = 10 and Da = 50. These two sets of Pe and Da both finally change to an edge deposition, as shown in Figures 9 and 10. The particle deposition rate is slower for Da_{−1} increasing in the case of large Pe and Da as shown by the right column of Figure 10, which is similar to the case of Pe = 1 and Da = 5 (Figure 9).

**(a)**

**(b)**

**(c)**

The effect of Pe on the deposition and the deposition rate are shown in Figure 11. With the same Da number and ratio of Da_{−1}/Da, increasing Pe decreases the deposition rate, and an edge deposition pattern is formed.

Some changes happen to the phase diagram (Figure 8) when desorption exists. The new phase diagram with Da_{−1} = 10 is shown in Figure 11. The effects of desorption on the phase diagram are as follows. When Pe is as small as 0.01, Da_{−1} = 10 has little effect on the deposition pattern. When Pe goes above 0.1 and below 5, the deposit still goes from a center deposition pattern to a near-edge deposition pattern as Da/Pe decreases, but the change happens at a higher Da. Similarly, at a large Pe (10 and above), the change from a near-edge deposition pattern to an edge deposition pattern also happens at a larger Da.

##### 3.3. Particle Deposition without Adsorption

An edge deposition pattern is what happens in a sessile water droplet on a glass substrate suspended with silica particles. The final deposition pattern is shown in Figure 12. Figure 12(a) shows the pattern recorded by a digital camera. When the details close to the contact line are looked into using an SEM, as shown in Figure 12(b) , the entire deposition pattern was actually composed of a dense assembly of particles at the outer ring and many discrete random clusters in the region from the location of the center of the initial sessile drop to the ring deposit.

The Pe number in this case is about 138.5. The exact value of Da for this case is unknown and varies due to variation in the solution chemistry of the sessile drop. The silica particles and the glass substrate in the water medium both carry negative charges on their surfaces, resulting in a large particle/substrate interaction repulsion force. Using the DLVO theory and the method to estimate the deposition rate constant based on the method of Ruckenstein and Prieve in 1976, we can roughly know that both Da number and Da_{−1} are much smaller than 1 in this case. With such a small Da number, Da_{−1} number, and large Pe number, the particles tend to stay in the drop phase rather than being adsorbed on the substrate, and the deposition rate is very slow.

The results with the previous simulation method did not match the experimental deposition patterns. While an edge pattern was observed in Figure 13, it was not a very sharp edge pattern as what was observed in the previous simulations for Pe > 10 and Da < 1. Then, a simulation was done with parameters that is specifically for the case of silica particles in a water drop on glass, with Pe = 138.5 and Da much smaller than 1. For this simulation, particles were swept to the contact line very quickly and accumulated to a very large concentration value close to the contact line, causing difficulties in the simulation.

The large concentration emerges because the previous simulations are for very dilute suspensions of colloidal particles. The particles in this case do not approach the critical concentration or maximum packing volume fraction for solid spheres, but this is not a reasonable assumption for the experimental case in Figure 12. Therefore, a critical concentration of particles at the maximum packing density is considered in the simulation [21]. Once the average particle concentration in an element reaches the maximum packing density, the suspension in that element converts into a “solid” phase in which no additional particles can be added to the element.

The initial particle concentration is 0.65 wt.% in our experiments, so the dimensionless concentration at which the solid phase forms is *C*_{P,max} = 130.86.

Since the droplet remained pinned at the contact line for most of the evaporation process in our experiment and the liquid phase readily wets the surface of the silica colloidal particles, the depinning of the droplet was not considered in our simulation when the solid phase is growing. Furthermore, we compute a “superficial velocity” in the regions of the solid phase since computing the detailed flow fields in the void spaces of the solid phase would be intractable. The simulation results of the particle deposition and mass accumulation with the maximum packing criterion are shown in Figure 14 and Figure 13.

Since Da and Da_{−1} are both zero in this case, Γ_{S} (the surface mass concentration of the dried deposit) in Figure 14 was computed in a different way from the case with a nonzero Da. Because Da = 0, all particles are assumed to be in the bulk phase. We assumed that all the particles deposit vertically onto the substrate at the termination of the simulation. Therefore, the “surface concentration” is calculated according to the particle distribution at the end of the simulation. Since there is still a little liquid at the termination, the drop phase is still a spherical-cap shape. Then, for the projected areas on the substrate, there are associated annular volume regions in the drop phase. Therefore, Γ_{S} is the mass of SiO_{2} particles in the regions divided by the projected areas.

The height of the solid deposit on the substrate is a function of the radial position. The height was computed by dividing Γ_{S} by *C*_{P,max} (the critical particle concentration at which the “solid” phase forms). From these calculations, the width and height of the ring deposit were approximately 123 *μ*m and 24.3 *μ*m, respectively, while the profilometer-measured values from the experiment were 100 *μ*m and 26 *μ*m. There is reasonably good agreement between the computed and experimental height and width of the ring deposit.

Figure 13 shows the mass of the deposited particles and suspended particles individually. The deposited particles are the ones in the “solid” phase regions where the particle concentration reach the maximum packing criterion (*C*_{P} = *C*_{P,max}), forming the close packing ring in Figure 12. At the end of the simulation, there are still a significant number of particles suspended throughout the drop phase. These suspended particles eventually deposit as the rest of the liquid evaporates forming the random clusters as seen in Figure 12.

The effect of the initial concentration of silica particles on the width and height of the ring deposit is shown in Figure 15. The height and width of the ring deposit increase as the initial concentration of the silica particles increases.

#### 4. Conclusion

The particle deposition pattern from an evaporating sessile droplet was calculated using the Galerkin finite element method. In addition to the convective-diffusive equation in the drop phase, a convective-diffusive equation along the drop/substrate interface was successfully added into the simulation for more accurate deposition profile calculations as well as calculating the profile when there is desorption. When there is no desorption, higher Pe and lower Da generated an edge deposition pattern, and the deposition rate decreased as Da decreased. With desorption, an edged deposition pattern is formed at higher Da_{−1} numbers. The rate of deposition decreases as Da_{−1} increases. These results gave parametric studies on the effects of convection, adsorption, and desorption of particles on the final deposition profiles.

In the case of very large Pe and very small Da, the maximum packing density criterion was applied in the simulation. This criterion accounts for large initial concentrations of particles that lead to the “solid” phase formation during evaporation of the drop. The final deposition profile from the simulation agrees with the experimental results for suspension of silica particles in a sessile droplet of water on a glass substrate (an edge pattern). A significant number of particles remaining suspended at the termination of the simulation provided a source for the discrete random clusters that were seen interior to the ring deposition pattern. However, this is only a start of a model with maximum packing. Other details, like fluid dynamics in a porous media, will be added to give a more valid model in our future work.

#### Data Availability

Data are available in logbooks and in computer files with the corresponding author.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

The authors thank the Department of Energy, Office of Basic Energy Sciences Division of Materials Science and Engineering, Biomolecular Materials Research Program (DEFG02-02-ER45975 and DEFG02-02-ER45976) and the Purdue University Shreve Fund for the financial support.

#### Supplementary Materials

The file We have changed “Method section” to “Section 2.2” here for correctness. Please confirm.is a supplementary of the methods for our research. It contains discussion on the boundary conditions of the vapor domain, the boundary conditions of the particle adsorption/desorption to the substrate, and a sketch of the “solid region” which is mentioned in Section 2.2.* (Supplementary Materials)*