2.1 Assumptions and mixture properties

For the binary liquid, we consider an ideal mixture. Non-ideal, partially miscible systems are not treated in this work, but it is possible to extend the theory if necessary. In addition, it is assumed that the surrounding atmosphere is an ideal gas. These assumptions enable the use of Raoult’s law, such that the vapour pressure of a liquid component is the product of that component’s vapour pressure when pure and its liquid-phase mole fraction. The total vapour pressure of the liquid phase, $p_{v}$ , is hence

where $x_{1}$ is the liquid-phase mole fraction of component $1$ on a solid-free basis, and ${p_{1}}^{\star }$ and ${p_{2}}^{\star }$ are the vapour pressures of each, pure, liquid.

Assuming that the concentrations of the volatile components in the surrounding atmosphere are zero and the system is well ventilated, the evaporation rates for the pure liquid of each component are given by

where ${k_{i}}^{\star }$ is the mass transfer coefficient for the transfer of pure liquid $i$ to the vapour phase.

By applying Raoult’s law, the total liquid evaporation rate, $E$ , is

For the viscosity of the mixture, a simple linear relation is adopted between the viscosities of the two liquids. It should be noted that the viscosity does increase with solids content; however, as we will discuss later, the flow within the droplet is dominated by surface tension effects. Later, we will show that the nature of this expression is unimportant, since surface tension effects dominate the flow profile.

To fit the surface tension of the binary liquid mixture as a function of liquid composition, one might crudely consider a linear relation. In order for the surface energy of the droplet to be minimised, however, there is a preferential partition of the liquid with lower surface tension to the droplet surface (Hildebrand & Scott Reference Hildebrand and Scott1964; Schmidt, Randall & Clever Reference Schmidt, Randall and Clever1966; Suri & Ramakrishna Reference Suri and Ramakrishna1968; Escobedo & Ali Mansoori Reference Escobedo and Ali Mansoori1998). For this reason and following Hildebrand & Scott (Reference Hildebrand and Scott1964) and Suri & Ramakrishna (Reference Suri and Ramakrishna1968), the surface tension is fitted using

where $A$ is the surface occupancy of the molecules ( $\text{cm}^{2}~\text{mol}^{-1}$ ) and is approximately the same for each component if the mixture is ideal, $R_{g}$ is the universal gas constant and $T$ is the temperature in kelvin. The $\pm$ term is used to ensure that the surface tension is always reduced compared to the linear relation. Note that the expression above does not include any impact of solids on surface tension. By performing experiments using inks containing different concentrations of light-emitting polymer, we did not observe a notable influence on surface tension.

2.2 Scaling terms and simplifications

Due to either surface inhomogeneities or adsorption of the solid material to the substrate, the three-phase contact line is assumed to be pinned. The model considers an axisymmetric droplet, in a cylindrical coordinate system (radius $r$ , height $h$ , azimuth ${\it\theta}$ ), such that velocity and gradients in the azimuthal direction are negligible ( $v_{{\it\theta}}=0,\partial /\partial {\it\theta}=0$ ). As can be seen in figure 1, the horizontal length scale is the droplet radius, $R$ , and the vertical length scale is the initial droplet height, $H$ . An extension to droplets with a non-circular cross-section, as in van Dam & Kuerten (Reference van Dam and Kuerten2008), could be an interesting topic for future work on P-OLED displays.

Figure 1.

Sketch of a droplet residing on a flat substrate.

Scaling is arbitrarily performed on the basis of the less volatile liquid, as outlined in table 1.

Table 1.

Scaling terms used to cast non-dimensional properties. A $^{\star }$ denotes saturated parameter values.

Typically, $R$ is about $125~{\rm\mu}\text{m}$ . The Bond number, $Bo$ , can be used to compare the relative magnitudes of gravity and surface tension on the droplet. It is found that typical values for the Bond number are very small, and consequently surface tension dominates at this length scale. Gravitational sagging of the droplet interface can hence be neglected.

where $g$ is the acceleration due to gravity. The liquid density ( ${{\it\rho}_{2}}^{\star }=1084~\text{kg}~\text{m}^{-3}$ ) (Garcia et al. Reference Garcia, Miranda, Leal, Ortega and Matos1991) and air–liquid surface tension ( ${{\it\gamma}_{2}}^{\star }=3.8\times 10^{-2}~\text{N}~\text{m}^{-1}$ ) (Wohlfarth Reference Wohlfarth and Lechner2008) are for methyl benzoate at 1 atm and $20\,^{\circ }\text{C}$ .

The Reynolds number, $Re$ , is very small. This enables simplification of the governing equations, by the creeping-flow assumption, in which inertial effects are neglected. Furthermore, we consider thin films, for which the lubrication approximation ( $Re\,H^{2}/R^{2}\ll 1$ ) holds.

The evaporative flux per unit area ( ${E_{2}}^{\star }=2.0\times 10^{-9}~\text{m}^{3}~\text{m}^{-2}~\text{s}^{-1}$ ) was measured by the mass loss from a petri dish containing pure methyl benzoate at $20\,^{\circ }\text{C}$ . The viscosity ( ${{\it\mu}_{2}}^{\star }=1.8\times 10^{-3}~\text{N}~\text{s}~\text{m}^{-2}$ ) (Garcia et al. Reference Garcia, Miranda, Leal, Ortega and Matos1991) is again for methyl benzoate at 1 atm and $20\,^{\circ }\text{C}$ .

The strength of diffusion in the horizontal direction will be a function not only of particle size but also of shape. We account for the relative strength of diffusion by using the Péclet number, $Pe_{p}$ , which is the ratio of the rates of convection and diffusion. Generally the diffusion is weak. As solids content increases, so the viscosity increases. This results in a smaller diffusion coefficient and a larger Péclet number. In our numerical solutions we will use the base case $Pe_{p}\rightarrow \infty$ , but we will handle a few cases of finite $Pe_{p}$ to investigate weak diffusion.

where $D_{p}$ is the Stokes–Einstein diffusion coefficient for the solid material in the liquid medium ( $1.2\times 10^{-12}~\text{m}^{2}~\text{s}^{-1}$ for a dilute suspension of rigid, 100 nm spheres) and a typical value for $H$ is ${\sim}10~{\rm\mu}\text{m}$ . We treat the solid as a suspension of non-interacting spheres. Note that for polymer chains the validity of the rigid-sphere assumption depends on the quality of the solvent. This could provide a topic for future research.

For the liquid–liquid diffusion, we use the liquid–liquid Péclet number, $Pe_{ll}$ .

The diffusion coefficient of the ‘solute’ liquid, $1$ , through the ‘solvent’ liquid, $2$ , is $D_{ll}$ . A typical order of magnitude for liquids is $10^{-9}~\text{m}^{2}~\text{s}^{-1}$ . For the numerical solutions we initially examine the case without diffusion ( $Pe_{ll}\rightarrow \infty$ ) in order to understand what would be expected to happen if there were no liquid–liquid diffusion. We then move on to consider finite $Pe_{ll}$ . Section 4.2 is devoted to the treatment of liquid–liquid diffusion and its significance.

The relative importance of viscous and surface tension forces is governed by the magnitude of the capillary number, $Ca$ . In the derivation of the governing equations, which follows in § 2.3, it is shown that the capillary number is

It is found under most conditions that surface tension effects dominate viscous ones. Even with an increase in viscosity of several orders of magnitude, due to increased solids content, the viscosity has no impact on film shape. Hence, all numerical solutions provided in this work will consider the limit of small $Ca$ . For the simulations, the limit of small capillary number is taken to be $Ca=10^{-3}$ ; this ensures that surface tension dominates whilst achieving an acceptable simulation time. Previous works (Fischer Reference Fischer2002; Eales et al. Reference Eales, Dartnell, Goddard and Routh2015a ) have shown that coffee-ring formation would be avoided if the capillary number could be made large.

As the volatile liquid mixture evaporates, the volume fraction of the solid component, ${\it\phi}$ , increases from an initial value of typically 1 % or less. Due to packing or gelling, a maximum value, ${\it\phi}_{max}$ , will eventually be reached. It is assumed that ${\it\phi}_{max}$ is the same for both liquids. The viscosity of the dispersion will increase with increasing solid volume fraction. The expression we use to account for this, following Krieger & Dougherty (Reference Krieger and Dougherty1959), is

Note that if polymer solutions are considered, a more specific, alternative relationship would be required. However, due to the very small capillary number, this viscosity correction term generally has no impact on the final film shape.

It is assumed that the presence of the solid material does not alter the surface tension significantly. Thus, the solutal Marangoni effect is neglected but it is possible to extend this work to account for it if required.

It has been argued in previous works (Ozawa et al. Reference Ozawa, Nishitani and Doi2005; Tarasevich et al. Reference Tarasevich, Vodolazskaya and Bondarenko2013) that evaporation will be completely suppressed when the volume fraction reaches ${\it\phi}_{max}$ . If this were the case, we would expect the final deposit to be wet. In other words, a fraction ( $1-{\it\phi}_{max}$ ) of the final film volume would still contain liquid. In practice this is not observed, meaning that evaporation, whilst significantly hindered, cannot be completely suppressed. This is accounted for by using a hindrance fraction, $f_{h}$ , such that the evaporation in close-packed/gelled regions is hindered in the following manner. In regions where ${\it\phi}={\it\phi}_{max}$ :

where $0\leqslant f_{h}<1$ , $E_{h}$ is the evaporation rate in the close-packed/gelled region and $E$ is the evaporation rate if there were no hindrance to evaporation. In the base case, the hindrance fraction, $f_{h}$ , is chosen to be constant at 0.95.

Since the purpose of this work is to predict trends and provide guidance in selection of binary liquid inks for inkjet printing applications, we will not concern ourselves greatly with the subtleties of edge-enhanced evaporation profiles. ${E_{1}}^{\star }$ and ${E_{2}}^{\star }$ are assumed to be spatially uniform. The evaporation rate, $E$ , is taken to be a function of position and time, solely because $x_{1}$ is a function of position.

2.3 Derivation of the governing equations

For pure liquid droplets, the governing equations for droplet height and volume fraction of solid material have been derived by Routh & Russel (Reference Routh and Russel1998), Salamanca et al. (Reference Salamanca, Ciampi, Faux, Glover, Mcdonald, Routh, Peters, Satguru and Keddie2001), Ozawa et al. (Reference Ozawa, Nishitani and Doi2005), van Dam & Kuerten (Reference van Dam and Kuerten2008), Tarasevich et al. (Reference Tarasevich, Vodolazskaya and Bondarenko2013), Wray et al. (Reference Wray, Papageorgiou, Craster, Sefiane and Matar2014), Eales (Reference Eales2015) and Eales et al. (Reference Eales, Dartnell, Goddard and Routh2015a ,Reference Eales, Dartnell, Goddard and Routh b ). Here, we extend the work to binary liquid mixtures. In addition, we derive the governing equation for the liquid composition and allow for evaporation hindrance upon gelation.

Starting with the incompressible, axisymmetric, steady-state Navier–Stokes equations, the non-dimensionalisation can be performed using the scalings outlined in table 1. Applying the lubrication assumption, we find that for the thin film there is negligible depthwise gradient in pressure. The horizontal velocity and the surface pressure gradient are related by

where $\bar{{\it\mu}}$ is the non-dimensional viscosity, neglecting the influence of polymer, given by

An expression for the horizontal velocity can be determined by integrating (2.15) twice, with respect to the vertical coordinate. The no-slip boundary condition is used at the substrate ( $\bar{z}=0,\bar{v_{r}}=0$ ) and the tangential stress at the droplet surface must be equal to the surface tension gradient, following equation (2.17). Note that the nonlinear term in (2.6) is neglected in this balance. This is because an order-of-magnitude analysis conducted using a typical value for the surface occupancy of molecules, $A$ , from Suri & Ramakrishna (Reference Suri and Ramakrishna1968) shows that the nonlinear contribution is of negligible magnitude:

The pressure gradient is found by differentiating the Young–Laplace equation ( $p=-{\it\gamma}{\it\kappa}$ ), where ${\it\kappa}$ is the mean curvature of the droplet surface. A standard expression for ${\it\kappa}$ can be derived using differential geometry. Within the lubrication approximation, the pressure is

where $\bar{{\it\gamma}}$ is the non-dimensional surface tension at the air–liquid interface, given by

Subsequent differentiation yields the pressure gradient as

Hence, using (2.18), the horizontal velocity is

In the lubrication regime, we are interested in the horizontal flow within the droplet and ignore depthwise gradients. It is convenient to consider the vertically averaged horizontal velocity:

where

By examining the boundary condition at the droplet surface, defined in (2.17), two limits for the vertically averaged horizontal velocity can be considered, as follows.

Velocity limit 1: $\bar{{\it\mu}}$ and $\partial \bar{v_{r}}/\partial \bar{z}|_{\bar{z}=\bar{h}}$ are of order unity.

In this limit, the air–liquid interface has approximately zero shear stress. This is because a strong surface flow rapidly alters the surface composition to ensure there is no gradient in surface tension. This is a localised effect and the liquid composition can still vary in the bulk liquid below the surface.

Velocity limit 2: $\partial x_{1}/\partial \bar{r}|_{\bar{z}=\bar{h}}$ is of order unity.

In this limit, the vertical gradient in horizontal velocity is very large at the air–liquid interface. The curvature-induced capillary flow is negligible compared to the Marangoni flow.

For any scenario, the vertically averaged, horizontal velocity will be between the above limits.

Now that the horizontal velocity has been determined, the governing equation for the droplet height can be derived using a material balance, yielding a partial differential equation for the droplet height as a function of horizontal position and time:

where $\bar{E}(\bar{r})=x_{1}({E_{1}}^{\star }/{E_{2}}^{\star })+(1-x_{1})$ . There are two contributions to the change in height. The first is a decrease due to evaporation of the volatile liquid. The second term relates to flow driven by surface tension.

The solid volume fraction, ${\it\phi}$ , can be similarly derived to follow from

where the Péclet number, $Pe_{p}$ , is defined as $Pe_{p}=(R^{2}{E_{2}}^{\star })/(HD_{p})$ .

There are three contributions to the change in solid volume fraction: an increase due to evaporation of the volatile liquid, a decrease due to convection from regions of smaller volume fraction and a spatially variable diffusion.

A means to track the evolving liquid composition is also necessary. On a solid-free basis, one obtains

where

There are three contributions to the liquid compositional change. The first is a decrease in mole fraction of component 1, due to evaporation of component 1 and the increase in solid volume fraction. The second is an increase due to convection from areas of larger concentration. The third is a spatially variable, liquid–liquid diffusion contribution.

2.4 Boundary and initial conditions

The governing equation for height is a fourth-order partial differential equation. Four boundary conditions are hence required. The edge of the droplet has zero height and the horizontal velocity at this point is zero, to ensure that the droplet remains pinned. The system is symmetrical at the central boundary $(\partial \bar{h}/\partial \bar{r}|_{\bar{r}=0},\partial {\it\phi}/\partial \bar{r}|_{\bar{r}=0},\partial x_{1}/\partial \bar{r}|_{\bar{r}=0}=0)$ and there is zero flux across the centreline.

As the volatile liquid evaporates, the volume fraction of solid material increases. This is most marked towards the edge of the droplet where the height is smallest. Soon, a close-packed or gelled region forms at the periphery, when ${\it\phi}={\it\phi}_{max}$ . The innermost position of this region is termed the liquid/gel front, since its position changes with time. In the gelled region, the mole fraction of the liquid components may still change, provided evaporation occurs, but the height and solid volume fraction are fixed. The governing equations for height and solid volume fraction derived in § 2.3 are only applicable for the liquid domain, between the centre and the liquid/gel front. When this situation arises, the outer boundary conditions need to be changed. In this scenario the height is fixed at the location of the liquid/gel front, $\bar{r}_{f}$ , and the horizontal velocity is found by equating the flux through the liquid/gel front with the liquid evaporated from within the gelled region as shown in figure 2. For close-packed solids this assumption is very good. We note that shrinkage will occur for polymer gels, and the inflow of liquid from the liquid domain might be insufficient to ensure that the height is fixed in the gelled region. However, our present work focuses more on the final, rather than intermediary, film profile. Therefore, for our purposes, whether a polymer gel would start to shrink immediately, or after the entire droplet has gelled, would not be important.

Figure 2.

Sketch of the balance used to find the horizontal velocity at the position of the liquid/gel front. The flux into the gelled region is matched to the total evaporation rate from the gelled region.

This is a volumetric average and hence there is no need to include a term for the rate of front progression, $\text{d}\bar{r}_{f}/\text{d}\bar{t}$ . The integral in the following equation can be numerically evaluated at each time step:

The partial differential equation for the liquid composition (2.31) does not hold in the gelled region. It was previously assumed that volume fraction and height were functions of time and horizontal position. Instead the equation for the liquid composition in the gelled region becomes

The notable difference, compared to (2.31), is that there is no contribution from the change in polymer volume fraction.

Initially, the droplet is assumed to have a spherical cap shape, $\bar{h}=1-\bar{r}^{2}$ , and homogeneous distribution of solid material, at volume fraction ${\it\phi}_{0}$ . The liquid is initially homogeneously distributed, with mole fraction ratio ${x_{1}}_{0}$ : $(1-{x_{1}}_{0})$ .

2.5 Numerical implementation

A mixed scheme is chosen to solve the governing equation for height. Provided the time step is sufficiently small, the height change will be slower than changes to the spatial derivatives of height. This approximation enables the partial differential equation, (2.29), to be linearised. This is then solved using Newton’s method. An implicit formulation is used for the solid volume fraction and the mole fraction of liquid component 1.

2.6 Groups to examine

The numerical solutions in this paper will consider the relative volatility of the two liquid components, their initial mixing ratio, viscosity ratio and surface tension ratio, as well as the magnitude of the solid and liquid–liquid Péclet numbers.

4.1 Relative volatility

To examine the influence of the relative volatilities of the liquid components, the viscosity and surface tension ratios are set to unity. In addition, the convection-dominated regime is initially considered, so that both the polymer and liquid–liquid Péclet numbers tend to infinity. The final film shape is shown in figure 3 as a function of relative volatility. It should be noted that the case of ${E_{1}}^{\star }/{E_{2}}^{\star }=1$ is equivalent to a single-liquid system.

Figure 3.

Numerical results for the final film shape as a function of the relative volatility. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

The coffee-ring effect is hindered when the relative volatility is increased. The peak ring height, at the edge, is reduced and more material is deposited centrally. Soon after evaporation commences, the composition within the droplet depends on the horizontal position. This can be understood by examining the denominator in (2.31). At the edge of the droplet, the height of the droplet is very small, therefore the change in composition, due to the different evaporation rates of the two liquid components, is largest at this point. The edge becomes rich in the less volatile liquid and depleted in the more volatile liquid. Hence, the evaporation rate is smaller at the edge and larger at the centre when compared to a single-liquid system. The effect, therefore, is a reduction in the strength of the outward capillary flow required to maintain an equilibrium, spherical cap shape. Thus, less polymer is transported towards the contact line and less coffee-ring formation is observed. With increasing relative volatility, the less volatile-rich region becomes less confined to the periphery, as can be seen in figure 4. This reduces the coffee-ring effect further.

Figure 4.

The liquid composition profile after a dimensionless time $\bar{t}=0.01$ . The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

As a measure of the extent of the coffee-ring effect, we consider the ratio of the peak height of the deposited ring to the final height at the centre of the dried film. Figure 5 shows this measure as a function of the relative volatility. It shows that the coffee-ring effect is reduced when the relative volatility is increased. There are, however, diminishing returns at higher values of ${E_{1}}^{\star }/{E_{2}}^{\star }$ .

Figure 5.

The ratio of the peak height of the deposited ring to the final central height, as a function of the relative volatility. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

Depending on the desired final film shape, the relative volatility of the two liquids can be optimised. For example, if one wishes to achieve as flat a film as possible, there will be an optimal relative volatility. Figure 3 suggests that this would be at a relative volatility of greater than $10$ . If it was desired to obtain a deposit in the shape of a spherical cap, a much larger relative volatility would be appropriate.

4.2 Liquid–liquid Péclet number

In figure 6, the film shape resulting from droplets with a variety of liquid–liquid Péclet numbers can be seen. At infinite Péclet number, diffusion is negligible. The mechanism outlined in § 4.1 occurs and the coffee-ring effect is reduced when compared to a single-liquid system or a system in which the liquid volatilities are matched. In the presence of weak liquid–liquid diffusion, i.e. $Pe_{ll}=10^{1}$ , liquid material is redistributed within the droplet. The result is that the less volatile-rich region at the edge diminishes, due to the increased presence of the more volatile liquid. Subsequently, there will be a larger evaporation rate near the edge than in the $Pe_{ll}\rightarrow \infty$ case. It will take less time for the polymer to reach the gelation volume fraction and so there is less time for the outward capillary flow to convect polymer before the height is fixed by polymer gelation. The coffee-ring effect is further reduced, albeit only slightly.

Figure 6.

Final film shape as a function of the liquid–liquid Péclet number. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ and $Pe_{p}\rightarrow \infty$ .

The optimal improvement in film shape is achieved at approximately $Pe_{ll}=1$ . Below this, diffusion dominates and causes equilibration in liquid composition throughout the droplet. This removes the mechanism for the improvement in film shape by ensuring a spatially uniform evaporation profile across the droplet surface. At $Pe_{ll}=0$ , the system essentially behaves like a single-liquid system but with an evaporation rate that is initially the midpoint of ${E_{1}}^{\star }$ and ${E_{2}}^{\star }$ , and gradually decreases with time, as the mole fraction of more volatile liquid decreases.

4.3 Initial liquid composition

The influence of initial liquid composition on final film shape is considered, for a system in which the viscosities and surface tensions of the two liquid components are matched. The trend with surface tension ratio is discussed later in § 4.5. It should be noted that the initial liquid composition has an influence on this trend also. In figure 7, the final shape is shown for a variety of initial mixtures, with a relative volatility of 10 and for a polymer content ${\it\phi}_{0}=0.032$ and ${\it\phi}_{max}=0.080$ .

Figure 7.

Final film shape as a function of the initial mole fraction of liquid 1. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

At a fixed relative volatility and polymer content, there is an optimal initial mixing ratio that gives rise to the most improvement in film shape. In figure 8, the ratio of peak ring height to final central height is plotted as a function of the initial mole fraction of the more volatile liquid. When the initial mole fraction is $0$ or $1$ the system is essentially a single-liquid problem and no improvement in film shape due to the relative volatility mechanism occurs. In this example, the optimal initial composition is ${x_{1}}_{opt}=0.48$ for a relative volatility of 10 and ${x_{1}}_{opt}=0.54$ for a relative volatility of 20.

Figure 8.

The ratio of the peak height of the deposited ring to the final central height, as a function of the initial mole fraction of liquid 1 and the relative volatility. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

The initial and maximum volume fractions of the polymer are also important. In figure 9, the final shape is shown as a function of the initial liquid composition, at a much smaller polymer volume fraction ratio, ${\it\phi}_{0}/{\it\phi}_{max}$ , with ${\it\phi}_{0}=0.0032$ and ${\it\phi}_{max}=0.1600$ . The evaporation suppression factor for these simulations, $f_{h}$ , is set to unity, in order to avoid a complication concerning polymer volume fraction to be discussed later, in § 4.7. As the polymer volume fraction ratio is decreased, the optimal initial mole fraction of the more volatile component increases. For smaller polymer contents, a larger ${x_{1}}_{0}$ causes polymer gelation to occur much sooner. There is hence less time for the outward capillary flow to drive coffee-ring formation before the height becomes fixed.

Figure 9.

Final film shape as a function of the initial mole fraction of liquid 1. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.0032$ , ${\it\phi}_{max}=0.160$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

4.4 Viscosity ratio

As discussed in § 2.2, the capillary number is small: surface tension dominates over viscous effects. For this reason, the viscosity ratio is unimportant. It does alter the horizontal velocity slightly, but it does not alter the magnitude of the horizontal velocity and so the resultant effect on film shape is marginal. Figure 10 shows the film shape for various values of the viscosity ratio, which results in an appreciable change to the final film only when it is very large, i.e. the viscosity of the more volatile component is three orders of magnitude larger than the viscosity of the less volatile liquid. This is because the more volatile-rich central region has a reduced mobility, so the outward capillary flow is hindered. Situations with such a large differential in liquid viscosity are probably physically unreasonable.

Figure 10.

Final film shape as a function of the viscosity ratio. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

4.5 Surface tension ratio

By using typical values for the surface occupancy of molecules in (2.6), quoted by Suri & Ramakrishna (Reference Suri and Ramakrishna1968), we observe no appreciable difference in predicted film shape between using a linear fit and the nonlinear model of equation (2.6) for the surface tension of the binary liquid mixture.

There are two different situations. Either the more volatile liquid has a higher surface tension than the less volatile liquid, ${\it\gamma}_{1}/{\it\gamma}_{2}>1$ , or it has a lower surface tension, ${\it\gamma}_{1}/{\it\gamma}_{2}<1$ . As can be seen in figure 11, there is little difference in film shape in the latter case. When the more volatile liquid has a higher surface tension, however, we observe some interesting phenomena. There is a rapid decrease in height towards the centre of the droplet and a large increase in height towards the periphery. This leads to a crash in the code, which cannot track the height of the droplet accurately. It is believed that the reason for this is an instability in the height caused by the compositional and resulting surface tension gradients within the droplet driving a large internal flow. This will be investigated further in § 6.

Figure 11.

Final film shape as a function of the surface tension ratio. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

4.6 Solid Péclet number

The results shown up to this point consider the case of infinite polymer Péclet number. We extend this to consider weak diffusion. When the Péclet number decreases, diffusion of the polymer becomes more appreciable, and there is inward diffusion of polymer, down the gradient in polymer volume fraction. This results in a decrease in the coffee-ring effect. Less material is deposited towards the edge and more is deposited centrally. Figure 12 shows the final film shape as a function of the polymer Péclet number.

Figure 12.

Final film shape as a function of the polymer Péclet number. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=0.95$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1$ , ${x_{1}}_{0}=0.5$ and $Pe_{ll}\rightarrow \infty$ .

4.7 Comparison of the two horizontal velocity limits

The results in §§ 4.1–4.6 consider only the first horizontal velocity limit, stated in (2.26). It is important to consider what happens if the horizontal velocity is closer to the second limit, given in (2.28). In any situation, the actual horizontal velocity will be between these two extremes.

The final film shape is shown in figure 13 for a typical scenario with surface tension ratio ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }<1$ .

Figure 13.

Final film shape for the two limits of horizontal velocity. The results are shown for $Ca=10^{-3}$ , ${\it\phi}_{0}=0.032$ , ${\it\phi}_{max}=0.080$ , $f_{h}=1.0$ , ${E_{1}}^{\star }/{E_{2}}^{\star }=10$ , ${{\it\mu}_{1}}^{\star }/{{\it\mu}_{2}}^{\star }=1$ , ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=0.9$ , ${x_{1}}_{0}=0.5$ , $Pe_{p}\rightarrow \infty$ and $Pe_{ll}\rightarrow \infty$ .

In the second horizontal velocity limit, a strong outward Marangoni flow convects solid to the edge and results in an enhanced coffee-ring formation.

It is important to note that a numerical solution is also not possible for the second horizontal velocity limit, when the surface tension ratio ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }>1$ . It is thought that this is due to the onset of an instability and will be examined further in the linear stability analysis of § 6. We note that prior to the crash in the code, the solids fraction falls dramatically near the edge. This could be consistent with a strong inward Marangoni flow and associated depinning of the contact line, resulting in a deposit only at the droplet centre, as widely reported in the literature, one example being Babatunde et al. (Reference Babatunde, Hong, Makaso and Fukai2013).

5.1 Initial liquid composition

Figure 14 shows the film shape resulting from inks A–E. The initial composition appears to have a subtle effect. This is because liquid–liquid diffusion is significant and the tendency is to equilibrate the liquid composition within the droplet, thus removing or at least reducing the extent of the relative volatility mechanism proposed in § 4.1. Previous works, for example Fischer (Reference Fischer2002) and Eales et al. (Reference Eales, Dartnell, Goddard and Routh2015a ), have shown that coffee-ring formation is accentuated by edge-enhanced evaporation. Despite the capillary number being small in both cases, the evaporation rate profiles across droplets of methyl benzoate and anisole could be different. This could explain the repeatable difference in the final shapes for methyl benzoate and anisole, shown in figure 14. It should be noted that the error bars are not shown in figure 14, for clarity, but the experimental error is $\pm 10\,\%$ by height and cannot explain the trend.

Figure 14.

Experimental results for the final film profile from anisole/methyl benzoate inks A–E, as a function of the initial mole fraction of anisole. The properties of inks A–E can be found in table 2.

5.2 Relative volatility

The plots in figure 15 represent the final film profiles for inks A, C and F. There appears to be an improvement in film shape as the relative volatility is increased. Less material is deposited at the edge and more is deposited at the centre. The extent of coffee-ring formation is reduced. As the relative volatility is increased, we expect there to be a larger differential in evaporation rate between the centre and edge.

Figure 15.

Experimental results for the final film profile as a function of the relative volatility. The inks shown are A, C and F. The properties of these inks can be found in table 2.

This agrees with the work of Teichler et al. (Reference Teichler, Perelaer and Schubert2013), who found that the surface roughness of films decreases when the difference in liquid boiling points is increased.

It was previously argued that the capillary number is very small, such that surface tension effects dominate. This was true for the majority of droplets under consideration. In the case of ink F, however, the droplet wetted the ITO much better and the capillary number was of order $10^{-1}$ . Whilst surface tension effects are still important, this could make a small contribution to the improvement in shape observed in figure 15, as predicted by our previous work (Eales et al. Reference Eales, Dartnell, Goddard and Routh2015a ).

Figure 16.

Experimental results for the final film profile as a function of the surface tension ratio. The inset shows the white-light interferometer data for the anisole/pentyl benzene droplet (ink H). A primary coffee ring was observed, with a satellite at a seemingly random orientation. The inks shown are G and H. The properties of these inks can be found in table 2.

5.3 Surface tension ratio

In figure 16, the final film shape is shown for the anisole/methyl benzoate mixture of ink G and for the anisole/pentyl benzene mixture of ink H. For the former, the surface tension ratio ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=0.93$ at $80\,^{\circ }\text{C}$ and the resultant film shape is similar to those previously observed. For the latter, however, the surface tension ratio, ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }$ , is $1.19$ and an extreme coffee-ring formation results. There is minimal material deposited centrally within the droplet and a region of intense deposition at the periphery. These results were repeatable. This would appear to be consistent with the dynamics of the numerical solution under these conditions (see figure 11). Interestingly, the cross-section for ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }=1.19$ in figure 16 shows a larger coffee-ring height on one side and a smaller one on the other, with what appears to be a satellite coffee ring (see the cross-section of the final film profile shown in the inset at the top right of figure 16). Looking at the 3D deposit, this occurred repeatably but at a seemingly random orientation. The secondary crescent-shaped ring of much smaller height suggests that there is an overspill of the droplet, beyond the original confines of the edge. The instability and the associated strong outward flow are a possible cause of this phenomenon. Henceforth we will say that an instability has occurred if we observe a satellite coffee ring.

6.1 Description and methodology

We wish to determine any situations in which the ink droplets, printed onto our substrates, will have an instability in the height, due to the growth of perturbations. In the eventuality that they become unstable, it is desired to understand the mechanism driving the growth of perturbations. The presence of polymer is ignored, as are the polymer and liquid–liquid diffusion, because these are likely to have a stabilising effect. A detailed analysis of how to avoid instabilities may include these also, but this falls beyond the scope of the present work. The linear stability analysis is performed on (2.29) and (2.31), using base states $\bar{h}_{b}$ and $\bar{x}_{1,b}$ and small perturbation terms $\hat{\bar{h}}$ and $\hat{\bar{x}}_{1}$ , such that

In (6.1) and (6.2), $k$ is the wavenumber and ${\it\omega}$ is the growth rate. To determine the stability, it is necessary to determine the sign of ${\it\omega}$ . Expressions can be found for each of the terms in the partial differential equations (2.29) and (2.31), ignoring all nonlinear perturbation terms. The expressions resulting from consideration of the real part of the perturbation terms alone can be summarised as

where the matrix coefficients $\unicode[STIX]{x1D608}_{ij}$ represent the coefficient of growth for perturbation $i$ with respect to the introduced perturbation $j$ . For conciseness, the full form of these matrix coefficients is reported in supporting information available at http://dx.doi.org/10.1017/jfm.2016.163. Generally, the coefficients $\unicode[STIX]{x1D608}_{ij}$ are functions of $Ca$ , $k$ , $\bar{r}$ , $({E_{1}}^{\star }/{E_{2}}^{\star }-1)$ , $({{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }-1)$ , $\bar{x}_{1,0}$ , $\bar{h}(\bar{r})$ and $\bar{x}_{1}(\bar{r})$ .

Equation (6.3) is solved for ${\it\omega}$ . Two solutions are obtained and the system will be unstable if either solution is positive, at any physically realisable wavenumber, $k$ .

It is known that the height at the edge of the droplet is zero. Hence, only discrete values of $k=(m{\rm\pi})/2$ are physically meaningful for the system, where $m$ is a positive integer.

6.2 Solution methodology

In practice, the capillary number is small. The base case is therefore taken to be $Ca=10^{-3}$ throughout this analysis. We consider the situation soon after evaporation has commenced. The droplet has a spherical cap shape $\bar{h}=1-\bar{r}^{2}$ . For any ${x_{1}}_{0}$ , the liquid composition profile $\bar{x}_{1}(\bar{r})$ is taken from the numerical solution after one time step (with ${\rm\Delta}\bar{t}=10^{-3}$ ).

As the wavenumber $k$ gets larger, the solutions for ${\it\omega}$ tend to a steady-state value. For this reason, we solve for only the first $50$ discrete values of $k$ .

For any given parameter values the resulting ${\it\omega}$ as a function of $k$ is plotted at various radial positions, $r$ . Any positive value for ${\it\omega}$ indicates an instability, and hence a global stability locus can be generated as a function of ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }$ , ${E_{1}}^{\star }/{E_{2}}^{\star }$ and ${x_{1}}_{0}$ .

6.3 Results and interpretation

For conciseness, the results shown in §§ 6.3.1–6.3.5 and discussed below consider the first horizontal velocity limit, outlined in § 2.3 and (2.26). This corresponds to the condition of zero shear stress on the free surface. Comparative results are shown in § 6.3.6 for the second limit (2.28), a dominant Marangoni flow.

One solution for ${\it\omega}$ is always negative. Hence, only the one solution for ${\it\omega}$ , termed ${\it\omega}_{2}$ , is shown in subsequent figures. The stability results for the binary liquid droplets depend on the position under consideration, $\bar{r}$ , the surface tension ratio, ${{\it\gamma}_{1}}^{\star }/{{\it\gamma}_{2}}^{\star }$ , the initial mole fraction of the more volatile liquid, ${x_{1}}_{0}$ , and the relative volatility, ${E_{1}}^{\star }/{E_{2}}^{\star }$ . When considering the global stability, the situation under consideration is only stable if it is stable at every location within the droplet. For this reason, the linear stability calculations must be repeated many times.