Modelling the evaporation of nanoparticle suspensions from heterogeneous surfaces

TL;DR

This paper introduces a Monte Carlo grid-based model for simulating nanoparticle suspension evaporation on heterogeneous surfaces, accounting for surface roughness and contact-line pinning, with applications to inkjet-printed electrical connections.

Contribution

The model accurately reproduces dewetting phenomena and relates material interactions to drying behavior, aiding in optimizing surface designs for better electrical connections.

Findings

The model predicts dewetting behavior on heterogeneous surfaces.

Surface roughness influences contact-line pinning and drying dynamics.

Positioning hydrophobic regions below metal height reduces dewetting risk.

Abstract

We present a Monte Carlo (MC) grid-based model for the drying of drops of a nanoparticle suspension upon a heterogeneous surface. The model consists of a generalised lattice-gas in which the interaction parameters in the Hamiltonian can be varied to model different properties of the materials involved. We show how to choose correctly the interactions, to minimise the effects of the underlying grid so that hemispherical droplets form. We also include the effects of surface roughness to examine the effects of contact-line pinning on the dynamics. When there is a `lid' above the system, which prevents evaporation, equilibrium drops form on the surface, which we use to determine the contact angle and how it varies as the parameters of the model are changed. This enables us to relate the interaction parameters to the materials used in applications. The model has also been applied to drying on heterogeneous surfaces, in particular to the case where the suspension is deposited on a surface consisting of a pair of hydrophilic conducting metal surfaces that are either side of a band of hydrophobic insulating polymer. This situation occurs when using inkjet printing to manufacture electrical connections between the metallic parts of the surface. The process is not always without problems, since the liquid can dewet from the hydrophobic part of the surface, breaking the bridge before the drying process is complete. The MC model reproduces the observed dewetting, allowing the parameters to be varied so that the conditions for the best connection can be established. We show that if the hydrophobic portion of the surface is located at a step below the height of the neighbouring metal, the chance of dewetting of the liquid during the drying process is significantly reduced.

Equations

$$s={\gamma }_{\text{sg}}-({\gamma }_{\text{sl}}+{\gamma }_{\text{lg}}).$$

$${\gamma }_{\text{lg}}\cos \theta ={\gamma }_{\text{sg}}-{\gamma }_{\text{sl}}.$$

$$\mathbf{i}=(i,j,k)$$

$$E=$$

$$-\mu \sum _{\mathbf{i}}{l}_{\mathbf{i}}+\sum _{\mathbf{i}}{V}_{\mathbf{i}}^l{l}_{\mathbf{i}}+\sum _{\mathbf{i}}{V}_{\mathbf{i}}^n{n}_{\mathbf{i}}.$$

$$c_{\mathbf{ij}}=\{\begin{array}{ll}1& \mathbf{j}\in \{\text{NN }\mathbf{i}\}\\ \frac{3}{10}& \mathbf{j}\in \{\text{NNN }\mathbf{i}\}\\ \frac{1}{20}& \mathbf{j}\in \{\text{NNNN }\mathbf{i}\}\\ 0& \text{otherwise}\end{array}$$

$$V_{\mathbf{i}}^l=\{\begin{array}{ll}\infty & k<1\\ -12{\epsilon }_{wl}/5& k=1\\ 0& \text{otherwise,}\end{array}$$

$$P(S_{\alpha },t+1)$$

$$+\sum _{\beta \ne \alpha }{w}_{\beta \to \alpha }P(S_{\beta },t)$$

$$\frac{w_{\alpha \to \beta }}{w_{\beta \to \alpha }}=e^{-\Delta E/k_{B}T}$$

$$M=\xi \frac{V-N{\sigma }^3}{N{\sigma }^3}$$

$$⟨r^2⟩=6Dt$$

Full text

I Introduction

How ink or paint dries, i.e. how liquid droplets containing nanoparticles deposited on a surface evolve in time as the liquid evaporates has significant relevance in modern manufacturing. Inkjet deposition is increasingly used during the manufacture of functional nano-structured materials. An innovative recent example is the application described in Crozier et al. (2014) which uses inkjet printing as part of a new method for constructing solar panels. This includes making the electrical interconnections after the various different layers that form a solar cell have been laid down on the glass substrate and then scribed using laser ablation Crozier et al. (2014). The benefits of using inkjet printing include reduced costs, wastage and potentially improved performance. Several parts of the structure are inkjet printed. An insulating polymer layer can be formed by depositing an ink that consists of a polymer solution. Another ink consists of a suspension of conducting metal nanoparticles. As the solvent evaporates, the goal is for the remaining nanoparticles to form an electrically conducting connection over the surface.

The surface onto which the nanoparticle ink is printed consists of two materials: (i) a metal conducting surface that is either side of (ii) a strip of a polymer insulating material. The metallic part of the surface is hydrophilic and the polymeric part is hydrophobic. When the liquid is deposited onto these two materials side-by-side, there is a tendency for the liquid to dewet from the surface of the insulator and move onto the metal, since this reduces the energy of the system. In the solar cell manufacturing process Crozier et al. (2014), this insulating polymer strip is created by inkjet printing into a trench created on the surface by laser ablation, at a previous stage — see Fig. 1.

If the nanoparticle suspension deposited perpendicular to the polymer strip is to dry to form an electrical connection, it is crucial that the ink does not dewet from the hydrophobic surface. The aim of the present work is to understand when this dewetting occurs and also to determine if there are processes that can be done during manufacturing to prevent dewetting.

The specific example considered is a particular case of a more general class of problem: that of modelling the evaporation of a nanoparticle suspension from a heterogeneous surface. The deposition and drying of the ink involves processes that occur over a huge range of time and length scales. The procedure can be roughly split into two parts: (i) the process of the ink being sprayed from the print head and arriving at the surface and (ii) the behaviour of the ink as it dries, once it is on the surface. In our work, we focus solely on stage (ii), in which there are still processes that occur over a great range of length and time scales.

The nanoparticles move throughout the liquid with a diffusive dynamics, where changes occur on a time scale much larger than the time scale for rearrangements of the solvent molecules. Drop shape changes occur on a time scale that is very much larger than the molecular time scale and also the nanoparticle diffusive time scale. There are also several disparate length scales, ranging from the solvent molecular diameter scale, to the size of the nanoparticles, the scale of any surface structures and then largest of all, the liquid drop size. Because of this, modelling such a multi-scale system has many challenges. Mesoscopic thin-film partial differential equation based models can be used Oron et al. (1997); Kalliadasis and Thiele (2007); Thiele et al. (2009); Moosavi et al. (2009); Frastia et al. (2011); Fraštia et al. (2012) but relating properties of the microscopic interactions between the particles and the structures they form in the liquid is not straightforward, because this type of model describes the distribution of the nanoparticles over the surface via a height-averaged concentration profile. This does not allow a description of the variations in the nanoparticles density distribution in the direction perpendicular to the surface. A fully microscopic approach, such as molecular dynamics (MD) does include every aspect of the motion of the particles and can be used to describe small liquid drops on a surface Ingebrigtsen and Toxvaerd (2007); Lu et al. (2013); Tretyakov et al. (2013); Becker et al. (2014). Generally MD is computationally infeasible even for moderate system sizes due to the long time scales over which the evaporative drying occurs. Similarly, classical density functional theory (DFT) Evans (1979, 1992); Hansen and McDonald (2013) and dynamical DFT Marini et al. (1999); Archer and Evans (2004); Archer (2006, 2009) can describe in great detail the density profile of the liquid at the interface and the structure down to the scale of individual particles Evans (1992); Hansen and McDonald (2013); Nold et al. (2014); Hughes et al. (2015, 2017) but the level of detail makes this also computationally very expensive.

We require a coarse-grained model to describe the fluid dynamical processes of interest here but not to the degree of coarse-graining as is present in the thin-film equation based models. Thus, we develop a lattice model for the system using Monte Carlo (MC) to capture the non-equilibrium dynamics and model the system time evolution as a series of discrete events. We model the nanoparticles individually, incorporating in the model their diffusion through the liquid over time thus enabling a description of the structures they may form on the surface. However, instead of modelling every solvent molecule individually, we group them together and statistically model them by a single, larger effective ‘particle’ of the same size as the nanoparticles, also residing on a lattice. MC models of this type have been used before, initially by treating the system effectively in two dimensions Rabani et al. (2003); Pauliac-Vaujour et al. (2008); Vancea et al. (2008); Stannard (2011). However, more recently models that are fully three-dimensional have been used Sztrum et al. (2005); Kim et al. (2011); Jung et al. (2014); Crivoi and Duan (2014); Tewes et al. (2016). Our model here is of this kind but differs from previous studies in the manner in which we describe the particle interactions, allowing for correct modelling of the (hemispherical) liquid drop shape. Additionally, the effect of surface roughness is incorporated.

How a liquid wets a surface is characterised by the spreading parameter $s$ De Gennes et al. (2004). It is defined as the difference in the surface tensions between the liquid, gas and the substrate:

[TABLE]

The first term, $\gamma_{\text{sg}}$, is the excess free energy per unit area of the substrate when dry (i.e. in contact with the gas phase), referred to as the solid-gas interfacial tension. The second term is the excess free energy per unit area of the substrate when it is wet by a thick film of the liquid and is the sum of the solid-liquid interfacial tension $\gamma_{\text{sl}}$ and the liquid-gas interfacial tension $\gamma_{\text{lg}}$. When $s>0$ the liquid seeks to spread over the surface. In contrast, when $s<0$ the liquid only partially wets the substrate, forming a drop with contact angle $\theta$. Young’s equation De Gennes et al. (2004) relates the contact angle $\theta$ to the interfacial tensions

[TABLE]

Therefore, the contact angle and spreading parameter are related by $s=\gamma_{\text{lg}}(\cos\theta-1)$. From our simulation results we can calculate the contact angle and also determine how this depends on the parameters in our model. Thus, to model a particular experiment, we have to find the contact angle of the solvent on the particular material(s) in the substrate (many are available in the literature) and then we select our model parameters to match the experiments.

The remainder of this paper is laid out as follows: In Sec. II we describe our model and the MC algorithm for the dynamics. This section also presents results for the model when no nanoparticles are present, to illustrate the wetting behaviour of the pure solvent liquid on a uniform planar surface. We determine the dependance of the contact angle on the model parameters, to enable selecting values to match experiments. In Sec. III we briefly present the bulk solvent fluid phase diagram. In Sec. IV we present results for droplets containing nanoparticles evaporating from a smooth planar surface and also show how to include the effect of surface roughness by changing the fluid dynamics in the vicinity of the surface. Sec. V presents results for the drying of the nanoparticle suspension from a heterogeneous surface, with emphasis on the drying of liquid bridges spanning a hydrophobic patch. Finally, in Sec. VI we make a few concluding remarks.

II Lattice model for the system

The system is discretised onto a regular three dimensional grid with lattice spacing $\sigma$ and with periodic boundary conditions in the $x$ and $y$ coordinates. The surface of the substrate onto which the nanoparticle suspension (ink) is deposited is perpendicular to the $z$ direction. Ink cannot penetrate the surface. At the top of the simulation box, at $z=L$, we apply various different boundary conditions, discussed below. A typical starting configuration is displayed in Fig. 2. Each lattice site above the surface can be in one of three states: (i) empty, (ii) containing a nanoparticle or (iii) containing liquid. We refer to a lattice site containing liquid as containing a liquid ‘particle’, but it should be borne in mind that this does not mean an individual solvent molecule but rather many of them grouped together in a volume $\sigma^{3}$. The lattice constant $\sigma$ is most easily envisaged as being the diameter of the nanoparticles but this does not have to be so: one can also consider $\sigma$ to be a larger coarse-graining length scale, in which case when a lattice site is said to be ‘containing’ a nanoparticle, we mean ‘contains mostly nanoparticles’.

We define $n_{\mathbf{i}}$ and $l_{\mathbf{i}}$ to be the occupation numbers of lattice site $\mathbf{i}$ for nanoparticles and liquid respectively, where

[TABLE]

is the discrete position vector (c.f. Fig. 2). If site $\mathbf{i}$ is occupied by liquid, then $l_{\mathbf{i}}=1$, otherwise $l_{\mathbf{i}}=0$. Similarly, if site $\mathbf{i}$ is occupied by a nanoparticle then $n_{\mathbf{i}}=1$ and $n_{\mathbf{i}}=0$ if there is no nanoparticle. Liquid and a nanoparticle cannot occupy the same site.

We model the total energy of the system $E$ by the following sum:

[TABLE]

The first term, a sum over pairs of lattice sites, is the contribution from particle interactions. The overall strength of the interactions between pairs of nanoparticles is determined by the parameter $\varepsilon_{nn}$, between liquid and nanoparticles by $\varepsilon_{nl}$ and between pairs of liquid particles by $\varepsilon_{ll}$. The precise value of the interaction energy between pairs of particles at sites $\mathbf{i}$ and $\mathbf{j}$ is determined by the dimensionless coefficient $c_{\mathbf{i}\mathbf{j}}$, which decreases in value as the distance between the pair of particles increases. We use the following values

[TABLE]

where NN $\mathbf{i}$, NNN $\mathbf{i}$ and NNNN $\mathbf{i}$ stand for nearest neighbours, next nearest neighbours and next-next nearest neighbours, respectively Thus, we truncate all interactions for $|\mathbf{i}-\mathbf{j}|>\sqrt{3}\sigma$. The influence on the wetting behaviour of truncating the range of the interactions is discussed in Refs. Hughes et al. (2015, 2017).

The choice of particular values in Eq. (5) is important, as this leads to liquid droplets on the surface having a hemispherical shape. For example, if instead we set $c_{\mathbf{i}\mathbf{j}}=0$ for $\mathbf{j}\in$ NNN $\mathbf{i}$ and $\mathbf{j}\in$ NNNN $\mathbf{i}$, (i.e. just nearest neighbour interactions) then the system forms unrealistic rectangular shaped droplets, particularly at low temperatures. Thus, with the values in Eq. (5) the dependence of the gas-liquid surface tension on the orientation of the interface with respect to the grid is minimised. That one should select the particular values in Eq. (5) comes from noting that the sum over neighbours in Eq. (4) has the same form as a finite difference approximation for the Laplacian Robbins et al. (2011). It can be shown that the values for $c_{\mathbf{i}\mathbf{j}}$ given in Eq. (5) minimise the errors from discretising the Laplacian on the grid Kumar (2004), dictating the choice in Eq. (5) 111 In addition to leading to hemispherical droplets, we expect the particular values for $c_{\mathbf{i}\mathbf{j}}$ in Eq. (5) to result in the roughening transition Weeks (1980); Abraham (1980) temperature to be suppressed. Indeed, at all the temperatures for which we have performed simulations, we have seen no evidence of a roughening transition..

The second term in Eq. (4) is the contribution from treating the liquid as being coupled to a reservoir, which is the vapour above the surface. $\mu$ is the chemical potential of the vapour. The value of $\mu$ determines the rate at which the liquid evaporates from the surface. The last two terms of Eq. (4) are the contribution from the interaction with the surface, where $V_{\mathbf{i}}^{l}$ and $V_{\mathbf{i}}^{n}$ are the external potentials due to the surface exerted on the liquid and the nanoparticles, respectively. Assuming that the surface is composed of particles interacting with the fluid with interaction strength $\varepsilon_{wl}$ and a pair potential of the same form as the pair potentials in Eq. (4), then for a flat structureless surface the potential takes the form

[TABLE]

where $k$ is the perpendicular distance from the surface. Similarly, if the interaction strength with the nanoparticles is $\varepsilon_{wn}$, then the external potential for the nanoparticles takes the same form as (6) but the suffix $l$ is replaced by $n$.

The external potential in Eq. (6) is modified when the surface varies in height or if the material changes. For example, to model the situation illustrated in Fig. 1, since the polymer hydrophobic section of the surface is inkjet printed at an earlier stage, its surface height can be controlled. We denote the step in height from the polymer to the metal part of the surface, to be $h$.

In all that follows below, we non-dimensionalise and set $\varepsilon_{ll}$ to be the unit of energy and the lattice spacing $\sigma$ to be our unit of length. All other parameters are given in terms of these.

The Monte Carlo Algorithm

We denote a particular state of the system as $S_{\alpha}\equiv\{n_{\mathbf{1}},n_{\mathbf{2}},\ldots,l_{\mathbf{1}},l_{\mathbf{2}},\ldots\}$, i.e. a particular set of values of the occupation numbers, which we index with the label $\alpha$. We also denote the probability of the system being in this state at time $t$ as $P(S_{\alpha},t)$. The time evolution of this probability is given by the master equation

[TABLE]

where $w_{\alpha\to\beta}$ is the transition rate from state $S_{\alpha}$ to state $S_{\beta}$.

In equilibrium, where $P(S_{\alpha},t+1)=P(S_{\alpha},t)$, we have Landau and Binder (2014):

[TABLE]

where $T$ is the temperature, $k_{B}$ is Boltzmann’s constant and $\Delta E=E(S_{\beta})-E(S_{\alpha})$, with the energy $E$ given in Eq. (4). The following Metropolis Monte Carlo algorithm satisfies this Landau and Binder (2014):

Pick a random particle. 2. 2.

Pick a random neighbouring site. 3. 3.

Calculate the change in energy, $\Delta E$, from swapping these particles using Eq. (4). 4. 4.

If $\Delta E<0$, perform the swap. Otherwise, swap the particles with probability $e^{-\Delta E/k_{B}T}$.

The assumption here is that even when the system is out of equilibrium, the transition rates given by Eq. (8) still hold, allowing us to use the algorithm to determine the non-equilibrium dynamics of the liquid.

This algorithm is refined for a system with both nanoparticles and liquid. To evolve the liquid, a random site on the lattice is picked. The MC algorithm described above is used, however any move involving a nanoparticle is forbidden during a liquid step.

The nanoparticles are modelled differently. Instead of selecting any random particle from the system, we keep an explicit list of nanoparticle positions and select a nanoparticle from this list to attempt the following dynamics: Firstly, to prevent nanoparticles from being left floating when the surrounding liquid moves away, after selecting a nanoparticle, we first check if there is a vacancy in the lattice site below the chosen nanoparticle. If there is, the nanoparticle is moved down to that empty site, finishing the move. If the site below is non-empty or is part of the surface, we then perform a weighted sum over the neighbouring lattice sites to determine how much liquid there is in the vicinity of the nanoparticle. We calculate the quantity $\bar{l}_{\mathbf{i}}=\sum_{\mathbf{j}}c_{\mathbf{i}\mathbf{j}}l_{\mathbf{j}}$, where the coefficients $c_{\mathbf{i}\mathbf{j}}$ are the same as those used to calculate the energy, given in Eq. (5). If $\bar{l}_{\mathbf{i}}<5/3$ then the nanoparticle move is rejected. Only if $\bar{l}_{\mathbf{i}}\geq 5/3$ we do allow the nanoparticle to move, swapping with one of the neighbouring liquid particles, as per steps 2–4 above. This is done to prevent excessive nanoparticle movement once most of the liquid has evaporated, since the physical origin of the nanoparticle dynamics is the Brownian motion due to being suspended in the liquid. If there are not enough liquid particles neighbouring the nanoparticle, then it remains stationary. The threshold value $5/3$ was determined empirically; a lower value makes the nanoparticles too mobile on the dry surface, but higher values leads to the formation of immobile nanoparticle clusters. Note that the algorithm described above for evolving the nanoparticles violates detailed balance. This is in keeping with previous MC models for systems of this kind Rabani et al. (2003); Pauliac-Vaujour et al. (2008); Vancea et al. (2008) and is a consequence of the facilitated dynamics of the nanoparticles. Of course, for the liquid there is detailed balance.

Liquid and nanoparticles evolve at different rates. We perform $M$ liquid steps for every nanoparticle step. This ratio determines the diffusion coefficient of the nanoparticles in the liquid Rabani et al. (2003); Pauliac-Vaujour et al. (2008); Vancea et al. (2008). We set the value of $M$ to depend on the ratio of nanoparticles to non-nanoparticles in the system as:

[TABLE]

where $V$ is the volume of the system and $N$ is the total number of nanoparticles. For all simulations in this paper, we use a value of $\xi=0.2$. For typical systems this corresponds to a value of $M\approx 30$. Eq. (9) is required to prevent the nanoparticles “speeding up” as the liquid evaporates from the system, which decreases the ratio of liquid to nanoparticles.

Diffusion coefficient

In what follows the system is referred to as having evolved for a time of $x$ Monte Carlo steps (MC steps), which means that there has been an attempted move on average $x$ times per lattice site. To relate MC steps to the physical time scales, the diffusion coefficient of a single nanoparticle moving though the bulk liquid is determined.

This is calculated by running multiple simulations with a single nanoparticle in a system full of liquid. The distance $r$ that the nanoparticle travels is recorded at certain time intervals. A plot of $\langle r^{2}\rangle$ against the number of MC steps is then made. Note that $\langle{\cal P}\rangle$ denotes the statistical average of any quantity ${\cal P}$. Using the relation Einstein (1905)

[TABLE]

where $D$ is the diffusion coefficient and $t$ is time, the value of $D$ can be determined from the gradient of the plot.

For a system with $\mu/\epsilon_{ll}=6$ (a system filled with liquid), $k_{B}T/\epsilon_{ll}=0.6$ and averaging over $10,000$ simulations, a value of $D=2.6\times 10^{-4}\,\sigma^{2}\,{\text{MC step}}^{-1}$ was found. Thus, the Brownian timescale $\tau_{B}\equiv\sigma^{2}/6D=650$ MC steps. $\tau_{B}$ is the time it takes on average for a nanoparticle to diffuse a distance of order its own diameter. We obtain a similar value at the higher temperature $k_{B}T/\epsilon_{ll}=1.0$, since the value of $D$ only starts to change when the temperature is high enough or the chemical potential is low enough that the density of the vacancies in the liquid becomes sizeable. Although we specify times below in units of MC steps, knowing the value of $D$ allows to easily relate to the true timescales in a given system.

Determination of contact angles

Once the system, such as that illustrated in Fig. 2, has reached equilibrium we can measure the contact angle. This is done by taking an average along the length of the liquid ridge in the $y$-direction (c.f. Fig. 2). We average over the configurations of a liquid ridge instead of a hemispherical drop because this is easier to measure and gives us more samples to average over. This average calculates a density profile $\rho_{\mathbf{i}}=\langle l_{i}\rangle$. From this density profile, we define the liquid drop to be where $\rho_{\mathbf{i}}\sigma^{3}>0.5$. We then fit a circle to the top portion of the boundary of the drop using the Taubin circle fitting method Taubin (1991), illustrated in Fig. 3. From this circle, it is then straightforward to determine the contact angle, which is the angle made with the surface. The density profile in Fig. 3 is for a system with temperature $k_{B}T/\varepsilon_{ll}=1.0$ and wall attraction strength $\varepsilon_{wl}/\varepsilon_{ll}=0.7$. The 0.7 value corresponds to a weakly hydrophilic interface and so the liquid does not spread and forms a drop with a contact angle $\approx 75^{\circ}$. Increasing $\varepsilon_{wl}$ decreases the contact angle, corresponding to the surface becoming more hydrophilic. On the other hand, decreasing $\varepsilon_{wl}$ makes the surface hydrophobic.

In Fig. 4 we display a plot of the contact angle as the surface-liquid attraction strength $\varepsilon_{wl}$ is varied, for the temperature $k_{B}T/\varepsilon_{ll}=1.0$. In addition to results for this ‘smooth’ surface, we also include in Fig. 4 results for a ‘rough’ surface, discussed further below. We see that as the attraction strength increases, the contact angle decreases, until eventually at $\varepsilon_{wl}\approx\varepsilon_{ll}$ there is a wetting transition to a state where the liquid wets the surface, with contact angle $0^{\circ}$. For small values of $\varepsilon_{wl}$ the surface only weakly attracts the liquid, corresponding to a strongly hydrophobic surface on which the drop takes a shape that is close to a full circle, with a large contact angle. Owing to the way we define the wall potential, the contact angle plot in Fig. 4 varies only weakly with the temperature in the range $0.6<k_{B}T/\varepsilon_{ll}<1.2$, the range in which most of our results are calculated. At higher temperatures, one should expect the wall attraction strength for wetting to be lower. However, at higher temperatures the interfacial fluctuations become significant and the system is no longer in the regime relevant to modelling the drying of inkjet printed drops. At lower temperatures (results not displayed), the simulations become slow and the system becomes hard to equilibrate.

We also display in Fig. 4 the contact angle obtained for the liquid on a rough surface. This surface is physically rough on the scale of the lattice, modelled by randomly raising and lowering respectively one third of the blocks on the surface by one lattice spacing $\sigma$. When the wall is sufficiently attractive, for $\varepsilon_{wl}/\varepsilon_{ll}>0.5$, this generates a surface that contains many pits, into which liquid condenses (from the vapour) and becomes trapped. This makes the surface effectively more attractive and so the contact angle in this regime is decreased, compared to the smooth surface. However, for $\varepsilon_{wl}/\varepsilon_{ll}<0.5$ the surface roughness makes the surface more hydrophobic and with a larger contact angle than the smooth surface with the corresponding value of $\varepsilon_{wl}$. This is the well-known lotus effect used to create superhydrophobic surfaces via surface roughness Wenzel (1936); Cassie and Baxter (1944); Bico et al. (1999); De Gennes et al. (2004).

III Bulk solvent phase behaviour

Understanding the behaviour of the liquid in equilibrium gives us insight into how the liquid behaves out of equilibrium. Calculating the binodal allows us to pick parameters that correspond to a suitably high density liquid phase coexisting with low density vapour phase.

The binodal gives the coexisting density values for a system in equilibrium. Two coexisting phases have the same chemical potential, temperature and pressure in each phase.

Since we do not need to know the binodal densities with great accuracy we calculate the binodal by performing simulations in a long, narrow box of size $10\sigma\times 10\sigma\times 80\sigma$, with periodic boundary conditions, treated in the canonical ensemble. Initially one end of the box of filled with liquid particles, with the other half being empty. The simulation then equilibrates in a state with half the system in the liquid phase, coexisting with the other half containing the vapour.

To estimate the density of the two coexisting phases we first calculate the mean density $\rho_{\mathbf{i}}$ of each $10\sigma\times 10\sigma$ layer of the box. The layer densities are then split into two groups: those with $\rho_{\mathbf{i}}>0.5$ and those with $\rho_{\mathbf{i}}<0.5$. In each of these groups the statistical outliers are eliminated, since these are layers that correspond to the interface between the gas and the liquid. Then the mean of the remainder in each group is used as the density on the binodal. The result of this approach, over a range of temperatures, yields the binodal displayed in Fig 5. For example, when $k_{B}T/\varepsilon_{ll}=0.9$ the density of the coexisting liquid and vapour is $\rho_{l}=0.99$ and $\rho_{g}=0.01$.

In the vicinity of the critical temperature $T_{c}$, the binodal can not be calculated with any great accuracy using the approach described above, due to the fact that the average densities of the two coexisting phases become rather similar and also because the system is strongly fluctuating with a diverging correlation length at $T_{c}$ Hansen and McDonald (2013); Landau and Binder (2014). The top of the binodal curve displayed in Fig. 5 was estimated by inspecting the histogram of densities $\rho_{\mathbf{i}}$ in each layer. Below the critical temperature, this histogram has two distinct maxima, corresponding to the two phases. We use the density value at each of these maxima as our estimate for the densities of the two coexisting phases. For $T>T_{c}$, the density histogram has only one maximum, at $\rho_{\mathbf{i}}=0.5$. Based on this method we find that the critical temperature is $k_{B}T_{c}/\varepsilon_{ll}=2.08\pm 0.02$.

IV Evaporating droplets

For evaporation to occur, the statistical mechanics of the system must be done in the semi-grand canonical ensemble: the liquid is treated grand canonically, whilst the nanoparticles are dealt with canonically. The liquid is treated grand canonically because the vapour above the drop acts as a reservoir with chemical potential $\mu$, with which the system can exchange particles, allowing the number of liquid particles in the system to vary over time. This is achieved by periodically setting the density of the top layer of the system to the low density result $\rho_{\mathbf{i}}\approx e^{\beta\mu}/(1+e^{\beta\mu})$ Hughes et al. (2014, 2015). This enables the removal of particles from the system as the liquid drop evaporates. In contrast, the nanoparticles are treated canonically, since the number of nanoparticles in the system is fixed over time. In contrast, as discussed above, in order to determine the contact angle of a drop of liquid on the surface, we must treat it canonically, with a fixed number of liquid particles in the system.

We initiate the system with a fraction $\phi$ of the liquid particles replaced by nanoparticles. Fig. 6 shows the evaporation of a droplet containing nanoparticles with initial concentration $\phi=0.15$. The substrate area is $120\sigma\times 120\sigma$ and the height of the top of the simulation box above the substrate is $80\sigma$. The initial droplet consists of a hemisphere with a radius of $40\sigma$ with the vertical part linearly scaled to have a height of $24\sigma$. The chemical potential is $\mu/\varepsilon_{ll}=-9$ and temperature $k_{B}T/\varepsilon_{ll}=0.8$ which corresponds to an equilibrium vapour with a density $\rho_{g}=0.001$. The interaction parameters are $\varepsilon_{lw}/\varepsilon_{ll}=\varepsilon_{nw}/\varepsilon_{ll}=0.8$, $\varepsilon_{ln}/\varepsilon_{ll}=1.25$ and $\varepsilon_{nn}/\varepsilon_{ll}=1.5$.

The drop in Fig. 6 initially spreads to cover a greater area of the surface, since the starting configuration does not have the equilibrium contact angle. However, over time, liquid evaporates and the drop reduces in volume and so subsequently the area of the surface covered by the drop decreases – i.e. the contact line initially advances and then later recedes. Owing to the smooth surface, the drop retains a dynamic contact angle that is close in value to the equilibrium contact angle throughout most of the time evolution. Since the nanoparticles are attracted to the liquid they generally follow the liquid.

After most of the liquid has evaporated there is then a further spreading of nanoparticles over the surface. Because of the smoothness of the surface and the non-zero vapour density, the residual liquid facilitates a diffusive dynamics that allows the nanoparticles to spread out over the surface to a state where the average distance of the nanoparticles from the centre of the system is larger than when the liquid is present.

Surface roughness

The roughness of surfaces is known to play an important role in how liquids spread. Surface roughness can hinder contact line motion over the surface and can lead to significant differences between the advancing, receding and equilibrium contact angles De Gennes et al. (2004). We consider two different methods of modelling the effect of surface roughness. The first is to introduce a dynamic rule that forbids moves parallel to the surface for all particles in contact with the surface. This is equivalent to a no-slip boundary condition. Thus, for a contact line to advance, particles in the second layer of lattice sites or higher above the surface must advance and then drop down to wet the dry surface ahead of the spreading droplet.

Fig. 7 shows snapshots from a simulations with the same parameter values as the evaporation simulation in Fig. 6 but with the no-slip dynamical rule forbidding moves across the surface. The droplet still spreads to a cover an area similar to that in the case with the smooth surface – i.e. to a value similar to that dictated by the equilibrium contact angle for this particular value of $\varepsilon_{wl}$. We then find that once most of the liquid has evaporated, the nanoparticles are left in an almost uniform circle which has a slightly larger radius than the original drop. There is also no further spreading out over the surface, even though the vapour density is still non-zero.

We have also investigated droplet evaporation from the rough surface considered at the end of Sec. II that is physically rough on the scale of the lattice, made by randomly setting the hight of the surface to be 0 or $\pm\sigma$, each with equal probability. Results for this surface are displayed in Fig. 8. Recall that for $\varepsilon_{wl}/\varepsilon_{ll}>0.5$ the contact angle is less than on the corresponding flat surface (see Fig. 4). This second approach to modelling surface roughness generates a wall that contains many pits, within which liquid becomes trapped. This leads to a much higher amount of liquid remaining adsorbed on the surface than in the cases in Figs. 6 and 7. The adsorbed liquid facilitates the spreading of the nanoparticles over the surface out to distances well beyond where the liquid droplet was located. Whilst this facilitated dynamics is interesting, it is not what is observed on the experimental surfaces of interest here.

In Fig. 9 we display a plot of the mean distance $\langle r\rangle$ of the nanoparticles from the centre of the system (the centre of where the droplet was initiated) as a function of time for the three different surface roughness cases. Initially the average radius increases due to the droplet spreading in order to try and reach the equilibrium contact angle. However, except in the case where we model the surface roughness via the no-slip dynamic rule, once the liquid has evaporated, $\langle r\rangle$ further increases because the nanoparticles continue to spread out over the surface, facilitated by the vapour of liquid particles. The plateau value of $\langle r\rangle$ for the physically rough surface (Fig. 8) is even greater than the smooth surface case (Fig. 6) due to the higher amount of liquid adsorbed on the surface, in the surface pits. Eventually, $\langle r\rangle$ tends to a constant value as the spreading nanoparticles become trapped in the pits in the surface. For the case with the no-slip dynamic rule modelling surface roughness, as the droplet spreads and evaporates, $\langle r\rangle$ reaches a maximal value at around $t=4\times 10^{11}$. It then decreases slightly as the droplet contact line starts to recede, due to the droplet volume being decreased by the evaporation.

The results in Figs. 7 and 9 show that incorporating the effects of surface roughness via the no-slip dynamical rule seems to model the required physics. It also has the additional advantage of not introducing an additional length scale to be considered, namely the length scale of the surface roughness. Thus, this is the model we adopt henceforth to model the effects of surface roughness.

V Modelling the ink drying process

Evaporating liquid bridge over a hydrophobic strip

In Fig. 10 we display snapshots as the liquid evaporates from a surface containing a hydrophobic strip. The interaction parameters are given in the figure caption. The surface is smooth – i.e. we do not implement the no-slip dynamical rule. The chemical potential in the vapour is set to be $\mu/\varepsilon_{ll}=-6$, which corresponds to the vapour phase being the thermodynamic equilibrium state, so the liquid seeks to evaporate from the surface. The substrate is made of a central hydrophobic strip of width 20$\sigma$ (coloured yellow) which we denote region A, with $\varepsilon^{\rm A}_{wl}/\varepsilon_{ll}=\varepsilon^{\rm A}_{wn}/\varepsilon_{ll}=0.4$, i.e. only weakly attracting the liquid and the nanoparticles. From Fig. 4 we see that on this part of the surface the liquid has contact angle $\theta\approx 110^{\circ}$. Either side of this (coloured grey) is region B, where the surface is hydrophilic, having attraction strength parameters $\varepsilon^{\rm B}_{wl}/\varepsilon_{ll}=\varepsilon^{\rm B}_{wn}/\varepsilon_{ll}=0.8$, corresponding to $\theta\approx 60^{\circ}$. When there is no step in height from region B to region A ($h=0$), then Fig. 10 shows that during the drying, the liquid dewets from the hydrophobic part of the surface, breaking the nanoparticle bridge at time $t\approx 1.2\times 10^{12}$ MC steps. When there is a small step of height $h=\sigma$ (results not displayed), then the behaviour is similar, although the breaking of the bridge is slightly delayed. In contrast, a step of height $h\geq 2\sigma$ enables the liquid bridge to remain intact as it dries, so that a nanoparticle bridge is formed, spanning the hydrophobic part of the surface. The nanoparticle density is even slightly increased on the hydrophobic part of the surface when $h\geq 2\sigma$ (see Fig. 10).

The reason a step enables the liquid bridge to remain is that a corner is created into which the liquid is strongly attracted. The ability of corners and wedges to promote wetting by a liquid is well known Concus and Finn (1969); Pomeau (1986); Hauge (1992); Rejmer et al. (1999); Malijevský and Parry (2013a, b). Since surface roughness can also modify the wetability of surfaces, a combination of steps and roughness can be used to control dewetting.

Fig. 11 shows results from a case when the nanoparticles are less strongly attracted to one another, which enhances the spreading over the hydrophilic part of this (smooth) surface, compared to the case in Fig. 10. With no step present ($h=0$), the bridge of liquid breaks at the time $t\approx 2.4\times 10^{11}$ MC steps and the nanoparticles temporally group together with the remaining liquid, but eventually spread out over the hydrophilic region. With a step of height $h=\sigma$, the bridge still breaks at $t\approx 2.4\times 10^{11}$ MC steps. More nanoparticles remain at the corner formed from the step but the end result is similar to the case with no step. Although not shown here, when $h=2\sigma$, the connection breaks at $t\approx 2.1\times 10^{11}$ MC steps but the break occurs in the hydrophilic region and the nanoparticles collect in the hydrophobic region, scattering randomly as the rest of the liquid evaporates, due to the smoothness of the surface.

Fig. 12 shows results from a simulation where the parameters are the same as in Fig. 11, except here we assume the surface is rough, i.e. we implement the no-slip dynamical rule. When there is no step ($h=0$), the bridge breaks at $t\approx 4.4\times 10^{11}$ MC steps. When the step height $h=\sigma$, the bridge connection almost holds, but eventually breaks at $t\approx 5.2\times 10^{11}$ MC steps. Interestingly, however, due to the attractive step from the hydrophobic to the hydrophilic region, most of the nanoparticles are stabilised in a cluster on the hydrophobic region.

Evaporating film over an hydrophobic strip

Figs. 13 and 14 illustrate the situation when a film of nanoparticle suspension that initially has uniform thickness evaporatively dewets from the same surface considered already, i.e. with both hydrophobic and hydrophilic parts. All the parameters for the simulations in Figs. 13 and 14 are the same as in Fig. 11 except for the temperature which is increased from $k_{B}T/\varepsilon_{ll}=0.6$ to $k_{B}T/\varepsilon_{ll}=0.76$, which slightly speeds up the simulations.

Fig. 13 for the smooth surface shows there are differences between $h=0$ and $h=\sigma$. As the liquid evaporates, in both cases holes appear in the film during the drying, at around $t\approx 2.0\times 10^{11}$ MC steps. We are not able to determine conclusively whether these holes are nucleated or are formed via spinodal dewetting, which is expected to occur when the film thickness decrease below a critical value Kalliadasis and Thiele (2007); Thiele et al. (2009); Robbins et al. (2011). In the $h=0$ case, the holes appear first in the hydrophobic region. This leads to a dewetting of the liquid from off the hydrophobic region, moving many of the nanoparticles onto the hydrophilic region. In contrast, for the $h=\sigma$ case, since the film is thicker over the hydrophobic region, the holes instead appear first over the hydrophilic region. Thus, in the $h=\sigma$ case, initially the dewetting from the hydrophilic part of the surface leads to a clear increase in the density of the nanoparticles on the hydrophobic region. However, they then subsequently move back onto the hydrophilic part of the surface as the evaporation continues. In both cases, after most of the liquid has evaporated, the nanoparticles are distributed inhomogeneously over the surface, having a greater density on the hydrophilic part of the surface. However, for the $h=\sigma$ case, because the nanoparticles congregate at the corner of the steps, there is therefore slightly more bare patches on the hydrophilic part of the surface, compared to the $h=0$ case.

Fig. 14 shows snapshots from two simulations with the no-slip dynamical rule, which prevents horizontal movement of particles that are in contact with the surface. Holes in the film appear in a manner similar to that observed in the early stages of the dynamics when the surface is smooth (Fig. 13). However, once the holes are formed, the dynamics is changed significantly. The surface roughness results in the nanoparticles becoming congregated in clumps and they spread far less than in the case with the smooth surface. For the $h=0$ case in Fig. 14(a), the final state consists of the nanoparticles being clustered into two mounds with fewer lone nanoparticles than observed on the smooth surface. Surprisingly, one of the nanoparticle clusters spans the hydrophobic region of the surface. We believe this stems from the interplay of the no-slip dynamics and the fact that the attraction of the nanoparticles to each other is stronger than their attraction to the surface.

In the case with a step of height $h=\sigma$ displayed in Fig. 14(b), the dewetting initiates in the thinner film on the hydrophilic part of the surface, similar to the smooth surface case in Fig. 13(b). This leads to the nanoparticles becoming deposited on the hydrophobic region, similar to in the smooth surface case, except that fewer nanoparticles remain on the hydrophilic region. However, in contrast to the smooth surface case, ultimately the relative lack of mobility leads to the nanoparticles remaining on the hydrophobic region, forming a large cluster that is stabilised at the edges by the step.

VI Conclusion

In this paper we have presented a MC model of the drying of a nanoparticle suspension on heterogeneous surfaces. This mixture of liquid and nanoparticles is a simple model for the ink that is used in the ink jet printing manufacturing process described in Crozier et al. (2014). The model contains parameters which can be determined from experiments. Measuring the equilibrium contact angle of drops of the liquid on the relevant surfaces, in conjunction with the present work, allows the determination of the required values of the liquid-liquid and liquid-wall attraction parameters. Similarly, knowledge of the diffusion coefficient allows to relate the MC time step to the experimental time scales. The model can include the effects of surface roughness via a simple no-slip dynamical rule that forbids the motion of all particles that are in contact with the surface.

A key finding of the present work is the observation that when printing a bridge over a hydrophobic region to connect hydrophilic strips either side, adhesion is improved when the hydrophobic strip is at a lower level than the surrounding hydrophilic regions. We find that when the bridge does not properly form, generally the break occurs over the hydrophobic strip. However, for some parameter values occasionally the counter-intuitive result occurs, whereby the film breaks, but with nanoparticles congregating in the hydrophobic strip. This effect generally occurs when considering the evaporation of a film of liquid, rather than a bridge. That said, evaporating films can still result in clumps of nanoparticles distributed over the two regions.

The results have shown the necessity to choose ink and surface parameters carefully to obtain the best connections when ink jet printing. For example, it may be possible to enhance particle bridge formation by adjusting the surface chemistry of the nanoparticles to make them favour the hydrophobic portion of the surface. This aspect has not been explored here. However, such enhancement might also instead lead to results such as that in Fig. 12, where the bulk of the nanoparticles are deposited on the hydrophobic part of the surface and the bridge is broken. Further work will directly relate the parameters in the model Hamiltonian to the properties of the specialist materials used in the printing process. Our work here shows that to fully understand the observed phenomena requires knowledge of both the fluid dynamics and the thermodynamics.

Acknowledgements

The authors would like to thank Adam Brunton of M-SOLV for useful discussions concerning ink jet printing and also Dmitri Tseluiko for insightful comments on our work.