2.1. Diffusion

To make explicit how the theory outlined here compares with previous work, whilst (2.7) in § 2.2 was also used by Trueman et al. (Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b), (2.12)–(2.14) correct the hydrodynamics implementation. These equations are given here in their most general forms. Also, (2.8)–(2.11) add clarification regarding diffusion reference frames. Equations (2.18)–(2.23) are reproduced from the work of Trueman et al., but they are subsequently substituted into the general conservation equation derived in § 2.2. §§ 2.3, 3.1 and 3.2 follow the same approach as Trueman et al., although again following through with the general equations.

2.1.1. Explanation of hydrodynamics correction

Using the thermodynamic approach in this present work, it is found that the extension from the one-component case of Routh & Zimmerman (Reference Routh and Zimmerman2004) to the two-component case is non-trivial. The Stokes–Einstein diffusion coefficient, ${D_{SE}}$, was originally derived for a single particle in infinite solvent, so it needs to be determined how to use ${D_{SE}}$ in a system that is no longer infinitely dilute. According to Batchelor (Reference Batchelor1976), ${D_{SE}}$ applies to a force-free solvent. Hence instead of deriving that the diffusive flux of particles in a one-component film, ${\boldsymbol{j}_1}$, is

(2.1)\begin{equation}{\boldsymbol{j}_1} ={-} {D_{SE}}\boldsymbol{\nabla }{\mu _1},\end{equation}

where $- \boldsymbol{\nabla }{\mu _1}$ is the force acting on component one due to its chemical potential, ${\mu _1}$, Batchelor derives

(2.2)\begin{equation}{\boldsymbol{j}_1} ={-} {D_{SE}}\left( {\frac{1}{{1 - {\phi_1}}}} \right)\boldsymbol{\nabla }{\mu _1},\end{equation}

where ${\phi _1}$ is the volume fraction of component one. The factor of $1/(1 - {\phi _1})\; $ results from a force correction. This force correction could be expressed as

(2.3)\begin{equation}{(\boldsymbol{\nabla }{\mu _1})_{eff}} = \boldsymbol{\nabla }{\mu _1} - \frac{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}}{{{\nu _s}}}\boldsymbol{\nabla }{\mu _s},\end{equation}

where $(4/3){\rm \pi} R_1^3$ is the volume of a particle of component one, ${\nu _s}$ is the volume of a particle of the solvent and ${\mu _s}$ is the chemical potential of the solvent. The subscript $eff$ denotes the effective driving force for diffusion. Since the ratio of particle volumes is a constant, (2.3) could be equivalently written as

(2.4)\begin{equation}{\mu _{1,eff}} = {\mu _1} - \frac{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}}{{{\nu _s}}}{\mu _s}.\end{equation}

This correction appears in the equations of Routh & Zimmerman (Reference Routh and Zimmerman2004). Trueman et al. (Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b) incorrectly extended this by writing

(2.5)\begin{equation}{\mu _{1,eff}} = {\mu _1} - \frac{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}}{{{\nu _s}}}\left( {\frac{{1 - {\phi_1} - {\phi_2}}}{{1 - {\phi_1}}}} \right){\mu _s} - {\left( {\frac{{{R_1}}}{{{R_2}}}} \right)^3}\left( {\frac{{{\phi_2}}}{{1 - {\phi_1}}}} \right){\mu _2},\end{equation}

where ${\phi _2}$ is the volume fraction of component two, and ${\mu _2}$ is its chemical potential. Whilst this appears to be a logical extension of (2.4), (2.3) is still correct for a film with any number of components, following the reasoning of Batchelor (Reference Batchelor1976). Equation (2.3) should therefore be used to replace the erroneous (2.5) in the work of Trueman et al.

In a multicomponent solution, (2.3) is used to form an expression for ${\boldsymbol{j}_{\boldsymbol{i}}}$ by summing over the contributions from all components in the mixture

(2.6)\begin{equation}{\boldsymbol{j}_{\boldsymbol{i}}} = {\phi _i}\mathop \sum \limits_j \frac{{{B_{ij}}}}{{6{\rm \pi}\eta {R_i}}}{(\boldsymbol{\nabla }{\mu _j})_{eff}},\end{equation}

where ${B_{ij}}$ is the bulk mobility coefficient, which is generally a function of all ${\phi _j}$ and ratios of ${R_j}$. Substituting in expressions for ${B_{ij}}$ allows an equation of the form of (2.12) to be derived (Batchelor Reference Batchelor1983).

This present work will use the methodology of Trueman et al. (Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b), but with the corrected (2.12). The numerical results of Trueman et al. solved diffusion fluxes of the same leading order in $P{e_1}$ and $P{e_2}$ as the equations in this present work, and so the results in § 4.1 are expected to show the same qualitative behaviour.

2.2. Governing equations

The conservation equation for one type of particles can be written as

(2.7)\begin{equation}\frac{{\partial {\phi _1}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{N}_1} = 0,\end{equation}

where ${\boldsymbol{N}_1}$ is the total flux of particles of component one. Following the approach of Cussler (Reference Cussler2009), the total flux can be arbitrarily split up into a convective and diffusive flux

(2.8)\begin{equation}{\boldsymbol{N}_1} = \boldsymbol{j}_1^{\boldsymbol{a}} + {\phi _1}{\boldsymbol{U}^{\boldsymbol{a}}} = {\phi _1}{\boldsymbol{U}_1},\end{equation}

where the volume average velocity of type one particles is ${\boldsymbol{U}_1}$, and their diffusive flux $\boldsymbol{j}_1^{\boldsymbol{a}}$ relative to the reference velocity ${\boldsymbol{U}^{\boldsymbol{a}}}$ is

(2.9)\begin{equation}\boldsymbol{j}_1^{\boldsymbol{a}} = {\phi _1}({\boldsymbol{U}_1} - {\boldsymbol{U}^{\boldsymbol{a}}}) ={-} \frac{{D_{11}^a}}{{kT}}\boldsymbol{\nabla }{\mu _1} - \frac{{D_{12}^a}}{{kT}}\boldsymbol{\nabla }{\mu _2}.\end{equation}

The chemical potential of component i in the two-component film is denoted by ${\mu _i}$, and the diffusion coefficient of particles of type one in particles of type i, using a reference velocity ${\boldsymbol{U}^{\boldsymbol{a}}}$, is denoted by $D_{1i}^a$. It is convenient to choose the volume average velocity ${\boldsymbol{U}^{\boldsymbol{v}}}$as ${\boldsymbol{U}^{\boldsymbol{a}}}$, since ${\boldsymbol{U}^{\boldsymbol{v}}} = 0$ in a static film. Hence (2.7) becomes

(2.10)\begin{equation}\frac{{\partial {\phi _1}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\phi _1}{\boldsymbol{U}_1}) = \frac{{\partial {\phi _1}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_1^{\boldsymbol{v}} = 0,\end{equation}

and diffusion coefficients using a volume-fixed reference frame are required.

For a volume-fixed reference frame,

(2.11)\begin{equation}{\boldsymbol{N}_1} = \boldsymbol{j}_1^{\boldsymbol{v}} + {\phi _1}{\boldsymbol{U}^{\boldsymbol{v}}} ={-} \frac{{D_{11}^v}}{{kT}}\boldsymbol{\nabla }{\mu _1} - \frac{{D_{12}^v}}{{kT}}\boldsymbol{\nabla }{\mu _2} = {\phi _1}{\boldsymbol{U}_1},\end{equation}

so expressions for $D_{11}^v$ and $D_{12}^v$ are required. It is desired to relate these to the Stokes–Einstein diffusion coefficient of component one, ${D_{SE,1}} = kT/6{\rm \pi}\eta {R_1}$.

Hydrodynamic effects are taken into account via ${K_{11}}({\phi _1},{\phi _2})$ and ${K_{12}}({\phi _1},{\phi _2})$, which are referred to as sedimentation coefficients by analogy with sedimentation (Batchelor Reference Batchelor1976). It is expected that these will be functions of ${R_1}$ and ${R_2}$ as well as ${\phi _1}$ and ${\phi _2}$. This gives

(2.12)\begin{equation}{\boldsymbol{U}_1} ={-} \frac{1}{{6{\rm \pi}\eta {R_1}}}[{K_{11}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _1} + {\phi _2}{K_{12}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _2}].\end{equation}

Defining the sedimentation coefficients such that ${U_1}$ is a velocity relative to the volume average velocity, and so can be used in (2.10) and (2.11), the conservation equation for component one becomes

(2.13)\begin{equation}\frac{{\partial {\phi _1}}}{{\partial t}} = \frac{1}{{6{\rm \pi}\eta {R_1}}}\boldsymbol{\nabla }\boldsymbol{\cdot }[{\phi _1}{K_{11}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _1} + {\phi _1}{\phi _2}{K_{12}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _2}].\end{equation}

Hence, $D_{11}^v$ and $D_{12}^v$ can be found in relation to ${D_{SE,1}}$ by comparison with the sedimentation coefficients. Note that, in this work, the diffusive fluxes are written in terms of chemical potential gradients, as opposed to concentration gradients. This is why we divide ${D_{SE,1}}$ by $kT$ in the flux expressions, (2.12) and (2.13), and why a factor of ${\phi _1}$ appears in front of the $\boldsymbol{\nabla }{\mu _1}$ term in (2.13), and likewise for the ${\phi _2}$ term in front of the $\boldsymbol{\nabla }{\mu _2}$.

Expressions for ${K_{11}}({\phi _1},{\phi _2})$ and ${K_{12}}({\phi _1},{\phi _2})$ for rigid spheres with zero interaction potential in dilute solution can be found in Batchelor's work (Reference Batchelor1983). However, it is desired to use generalised forms of these expressions to run the model up to close packing. A factor of ${\phi _1}$ is included in front of the $\boldsymbol{\nabla }{\mu _2}$ term in (2.13) in order to match the leading-order coefficients in $({\phi _1},{\phi _2})$ for dilute solution. Batchelor's (Reference Batchelor1983) dilute expressions will need to be adapted such that the sedimentation coefficients fall to zero as the mixture approaches close packing, due to hydrodynamic hindrance (Russel et al. Reference Russel, Saville and Schowalter1989). Steric effects could be incorporated into the chemical potential terms if desired.

We consider the case of one-dimensional drying with the top surface decreasing normally to the substrate. By scaling with $\hat{z} = z/H$ and $\hat{t} = t\dot{E}/H$, (2.13), in one dimension, becomes

(2.14) \begin{equation}\frac{{\partial {\phi _1}}}{{\partial \hat{t}}} = \frac{\partial }{{\partial \hat{z}}}\left[ {\frac{1}{{6{\rm \pi}\eta {R_1}\dot{E}H}}\left[{\phi_1}{K_{11}}({\phi_1},{\phi_2})\frac{{\partial {\mu_1}}}{{\partial \hat{z}}} + {\phi_1}{\phi_2}{K_{12}}({\phi_1},{\phi_2})\frac{{\partial {\mu_2}}}{{\partial \hat{z}}}\right]}\right].\end{equation}

2.2.1. Evaporation rate

For the purpose of (2.14), $\dot{E}$ needs only to be a constant, characteristic evaporation velocity used for non-dimensionalising. When evaporation at the top surface is included in the boundary condition (§ 3.2), it is assumed that the evaporation rate is approximately constant. Evaporation is driven by the thermodynamic driving force between the vapour pressure of the solvent at the top of the film and the vapour pressure of the gas above the film. The vapour pressure of the gas above the film is affected by the time scale for diffusion from the layer of saturated vapour directly above the film to the bulk gas. This time scale is orders of magnitude smaller than the evaporation time (Popòv Reference Popòv2005), leading to a quasi-static problem (Routh Reference Routh2013).

Routh & Russel (Reference Routh and Russel1998) derived an expression for $\dot{E}$,

(2.15)\begin{equation}\dot{E} = {k_m}\frac{{{\nu _s}({p_{vap}} - p_{vap}^\infty )}}{{kT}},\end{equation}

where ${k_m}$ is the gas-side mass transfer coefficient, ${p_{vap}}$ is the vapour pressure of the solvent at the top of the film and $p_{vap}^\infty $ is the vapour pressure of the bulk gas.

The chemical potential of the solvent is related to the osmotic pressure, $\varPi $, as

(2.16)\begin{equation}\varPi ={-} \frac{{{\mu _s} - \mu _s^0}}{{{\nu _s}}} = \left( {\frac{{{\phi_1}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}} + \frac{{{\phi_2}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_2^3}}} \right)kTZ({\phi _1},{\phi _2}),\end{equation}

where $\mu _s^0$ is the solvent reference potential and $Z({\phi _1},{\phi _2})$ is the compressibility, which accounts for the non-ideality of the osmotic pressure at high volume fractions. Relating the chemical potential to the solvent pressure through ${\mu _s} - \mu _s^0 = kT\ln ({p_{vap}}/p_{vap}^0)$, where $p_{vap}^0$ is the solvent reference vapour pressure, and extending the approach of Routh & Russel (Reference Routh and Russel1998) to a two-component mixture, (2.16) is substituted into (2.15), giving

(2.17)\begin{equation}\dot{E} = {k_m}\frac{{{\nu _s}}}{{kT}}\left[ {p_{vap}^0exp \left[ { - {\nu_s}\left( {\frac{{{\phi_1}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}} + \frac{{{\phi_2}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_2^3}}} \right)Z({\phi_1},{\phi_2})} \right] - p_{vap}^\infty } \right].\end{equation}

Since ${\nu _s} \ll (4/3){\rm \pi} R_1^3$, and ${\nu _\textrm{s}} \ll (4/3){\rm \pi} R_2^3$, as long the film has not yet reached close packing, we obtain $\dot{E}\sim {k_m}{\nu _s}(p_{vap}^0 - p_{vap}^\infty )/kT$, which is independent of film composition. In other words, the driving force is dominated by the humidity of the gas. In this one-dimensional drying model, a large surface area of film is assumed, so geometric edge effects are not applicable (Routh Reference Routh2013). When the film is close packed, $Z({\phi _1},{\phi _2})$ diverges.

This work models drying up until the film is close packed throughout. There are stages of drying beyond this (deformation and aging) that lead to film formation (Sonzogni et al. Reference Sonzogni, Passeggi, Wedepohl, Calderón, Gugliotta, Gonzalez and Minari2018), but it is the first part of drying, where the particles can move throughout the film, that affects their arrangement in the dried film. If the particles are colloids, as they would be in many common examples of films (Routh Reference Routh2013), then they are assumed to be colloidally stable. However, the model would also be valid for particles smaller than the colloidal range until the composition at which they precipitate.

2.3. Derivation of the spatial derivative of the chemical potential

The chemical potential of component one can be written as a function of ${\phi _1}$ and ${\phi _2}$, such that $\partial {\mu _1}/\partial \hat{z} = f({\phi _1}(\hat{z}),{\phi _2}(\hat{z}))$, and likewise for the chemical potential of component two, ${\mu _2}$. These expressions could be directly inputted into (2.14). However, it is difficult to find expressions for the chemical potential of the solutes that are valid as the solution approaches close packing. One means of resolving this is to relate ${\mu _1}$ and ${\mu _2}$ to the chemical potential of the solvent, ${\mu _s}$, for which an expression that diverges at close packing can be found.

For a system which contains ${n_1} = {\phi _1}/(4/3){\rm \pi} R_1^3$ particles of component one, ${n_2} = {\phi _2}/(4/3){\rm \pi} R_2^3$ particles of component two and ${n_s} = (1 - {\phi _1} - {\phi _2})/{\nu _s}$ particles of solvent per unit volume, conservation of volume can be expressed as

(2.18)\begin{equation}{\textstyle{4 \over 3}}{\rm \pi} R_1^3{n_1} + {\textstyle{4 \over 3}}{\rm \pi} R_2^3{n_2} + {\nu _s}{n_s} = 1.\end{equation}

The Gibbs–Duhem equation relates the chemical potentials in the system at constant temperature and pressure as

(2.19)\begin{equation}{n_1}\boldsymbol{\nabla }{\mu _1} + {n_2}\boldsymbol{\nabla }{\mu _2} + {n_s}\boldsymbol{\nabla }{\mu _s} = 0.\end{equation}

The Gibbs–Duhem equation (2.19), in one dimension, can be rearranged, giving

(2.20)\begin{equation}\frac{{\partial {\mu _1}}}{{\partial \hat{z}}} ={-} \frac{1}{{{\nu _s}}}\frac{{\partial {\mu _s}}}{{\partial \hat{z}}}\frac{{(1 - {\phi _1} - {\phi _2})}}{{\left( {\frac{{{\phi_1}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}} + \frac{{{\phi_2}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_2^3}}\frac{{\partial {\mu_2}/\partial \hat{z}}}{{\partial {\mu_1}/\partial \hat{z}}}} \right)}},\end{equation}

and

(2.21)\begin{equation}\frac{{\partial {\mu _2}}}{{\partial \hat{z}}} ={-} \frac{1}{{{\nu _s}}}\frac{{\partial {\mu _s}}}{{\partial \hat{z}}}\frac{{(1 - {\phi _1} - {\phi _2})}}{{\left( {\frac{{{\phi_1}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_1^3}}\frac{{\partial {\mu_1}/\partial \hat{z}}}{{\partial {\mu_2}/\partial \hat{z}}} + \frac{{{\phi_2}}}{{{\textstyle{4 \over 3}}{\rm \pi} R_2^3}}} \right)}}.\end{equation}

Using the spatial derivative of (2.16), plus the relationship ${R_1}/{R_2} = P{e_1}/P{e_2}$, (2.20) and (2.21) become

(2.22)\begin{equation}\frac{{\partial {\mu _1}}}{{\partial \hat{z}}} = \frac{{(1 - {\phi _1} - {\phi _2})}}{{\left( {{\phi_1} + {{\left( {\frac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}\frac{{\partial {\mu_2}/\partial \hat{z}}}{{\partial {\mu_1}/\partial \hat{z}}}} \right)}}\frac{\partial }{{\partial \hat{z}}}\left[ {\left( {{\phi_1} + {{\left( {\frac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}} \right)kTZ({\phi_1},{\phi_2})} \right],\end{equation}

and

(2.23)\begin{equation}\frac{{\partial {\mu _2}}}{{\partial \hat{z}}} = \frac{{(1 - {\phi _1} - {\phi _2})}}{{\left( {{{\left( {\frac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1}\frac{{\partial {\mu_1}/\partial \hat{z}}}{{\partial {\mu_2}/\partial \hat{z}}} + {\phi_2}} \right)}}\frac{\partial }{{\partial \hat{z}}}\left[ {\left( {{{\left( {\frac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1} + {\phi_2}} \right)kTZ({\phi_1},{\phi_2})} \right].\end{equation}

The purpose of rewriting the chemical potential gradients as in (2.22) and (2.23) is to obtain the advantages, discussed below, of expressing them as functions of $\partial {\mu _s}/\partial \hat{z}$ and $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$. The chosen model formulation intends for the effects of concentrated solution to be addressed via the divergence in solvent chemical potential, ${\mu _s}$, and the particle interactions to be input via $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$. Hence, this extends approaches that are valid only for more dilute solution, for example, Zhou et al. (Reference Zhou, Jiang and Doi2017). Since the ratio of the chemical potential gradients is determined physically by the inter-particle interactions, the chemistry of the types of particles that are used to form the film, such as their surface charge, will be relevant (Atmuri et al. Reference Atmuri, Bhatia and Routh2012).

Note that the Gibbs–Duhem equation (2.19) relates ${\mu _1}$, ${\mu _2}$ and ${\mu _s}$. Equation (2.16) gives an expression for ${\mu _s}$ which diverges at close packing. Hence, only one of ${\mu _1}$ and ${\mu _2}$, or a relationship between them, can also be specified. Since (2.20) and (2.21) give $\partial {\mu _1}/\partial \hat{z}$ and $\partial {\mu _2}/\partial \hat{z}$ only as implicit functions to be inserted into the conservation equation, different rearrangements of the Gibbs–Duhem equation could have been chosen instead. The particular rearrangements in (2.20) and (2.21) were chosen since their right-hand sides depend on $\partial {\mu _s}/\partial \hat{z}$ and the ratio $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$. The solvent chemical potential is measurable up to close packing, where it diverges. This means that $\partial {\mu _1}/\partial \hat{z}$ and $\partial {\mu _2}/\partial \hat{z}$ also diverge at close packing, in an unknown fashion. An expression chosen for the ratio $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$ is more likely to be valid approaching close packing than expressions that could be put forward for each particle separately. Equations (2.20) and (2.21) are also of the same form as each other, as desired.

Equation (2.16) for the osmotic pressure is used to generate an expression for ${\mu _s}$ only, rather than being used as an equation of state to also generate the other chemical potentials. The thermodynamic consistency of this approach with respect to the Maxwell relations is further explained in § S6 of the supplementary material (SI) available at https://doi.org/10.1017/jfm.2021.800.

In summary, a thermodynamically consistent approach is achieved by satisfying the Gibbs–Duhem equation. In choosing to specify $Z({\phi _1},{\phi _2})$, as a proxy for $\partial {\mu _s}/\partial \hat{z}$, and $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$, the system chemical potentials are fully specified. As formulated, the right-hand sides of (2.22) and (2.23) require input of an expression for $(\partial {\mu _1}/\partial \hat{z})/(\partial {\mu _2}/\partial \hat{z})$. Although the chemical potential expressions given in (2.22) and (2.23) appear complicated, they are merely (2.16) combined with rearrangements of the Gibbs–Duhem equation, in the same way as Russel et al. (Reference Russel, Saville and Schowalter1989) and Trueman et al. (Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b).

2.4. Diffusiophoresis

The governing conservation equation (2.7) requires an expression for the total flux. In this section, the aim is to investigate the effect of diffusiophoresis. The total flux is written as the sum of two terms, one described as the diffusion term, and the other described as the diffusiophoresis term.

The deviation of osmotic pressure from ideality, and hence $Z({\phi _1},{\phi _2})$ from unity, accounts for the total excluded volume (Russel et al. Reference Russel, Saville and Schowalter1989). However, this model uses the same form of $Z({\phi _1},{\phi _2})$ for the diffusion term, for both the flux of component one and two. In not differentiating between the expressions used for the two fluxes, this form of equation will not give rise to diffusiophoresis. It will give a diffusional model in the sense of the fluxes of both particle types remaining inversely proportional to their particle radius. The flux contribution that is due to diffusiophoresis can be found from first principles for a dilute solution, and written as a separate term. The language of this paper refers to this as the diffusiophoresis term.

Excluded volume effects are more general than diffusiophoresis resulting from excluded volume, since the latter specifically refers to one component moving in response to the concentration gradient of another, whereas in general a component can have volume excluded by both itself and other components. An alternative approach could write the total flux as one term, with a form of $Z({\phi _1},{\phi _2})$ from an equation of state that would include the diffusiophoretic flux. However, the approach taken here allows the effect of diffusiophoresis to be identified separately.

With this approach, a model including diffusiophoresis is constructed, assuming ${R_2} \gg {R_1}$. The particles of type one are modelled using an extension of the Asakura–Oosawa model to concentrated solution. Under the Asakura–Oosawa model, the type one particles are assumed to be excluded from a layer of solvent of thickness ${R_{DP}}$ around each particle of type two. Using this model, in a dilute solution, the drift velocity of the large particles due to diffusiophoresis, ${\boldsymbol{U}_{\boldsymbol{P}}}$, in a static film is given by

(2.24)\begin{equation}{\boldsymbol{U}_{\boldsymbol{P}}} ={-} \frac{{3{\phi _1}kT}}{{8{\rm \pi}\eta }}\frac{{R_{DP}^2}}{{R_1^3}}\boldsymbol{\nabla }\ln {\phi _1} = {\varGamma _1}\boldsymbol{\nabla }\ln {\phi _1},\end{equation}

where ${\varGamma _1}$ is the diffusiophoretic drift coefficient (Anderson & Prieve Reference Anderson and Prieve1984; Sear & Warren Reference Sear and Warren2017). Following the approach of Marbach, Yoshida & Bocquet (Reference Marbach, Yoshida and Bocquet2017), this is extended to concentrated solution via invoking the osmotic pressure (SI, § S1.2), to obtain

(2.25)\begin{equation}{\boldsymbol{U}_{\boldsymbol{P}}} ={-} \frac{{3{\phi _1}}}{{8{\rm \pi}\eta }}\frac{{R_{DP}^2}}{{R_1^3}}\boldsymbol{\nabla }{\mu _1}.\end{equation}

This is taken to be the extra slip velocity, in addition to that due to diffusion, between the particles of component two and the solvent. For the case of hard spheres, which is considered here, ${R_{DP}} = {R_1}$. Since ${\boldsymbol{U}_{\boldsymbol{P}}}$ gives the relative diffusiophoretic velocity between type two particles and the solvent, this needs to be converted to a velocity relative to the volume average velocity, ${\boldsymbol{U}^{\boldsymbol{v}}}$. Note that the result in (2.24) is the correct diffusiophoretic drift velocity for dilute solution, expected to agree with explicit solvent models. This contrasts with the result of implicit solvent methods which neglect the solvent dynamics, and consequently may overpredict the diffusiophoretic velocity (Sear & Warren Reference Sear and Warren2017).

The velocity of each type of particle i is split up into a component due to diffusion (denoted by the superscript $D$), which has been found in § 2.2, and a component due to diffusiophoresis (denoted by the superscript $P$)

(2.26)\begin{equation}{\boldsymbol{U}_{\boldsymbol{i}}} = \boldsymbol{U}_{\boldsymbol{i}}^{\boldsymbol{D}} + \boldsymbol{U}_{\boldsymbol{i}}^{\boldsymbol{P}}.\end{equation}

The local average velocities of the type one, type two and solvent particles are denoted by ${\boldsymbol{U}_1}$, ${\boldsymbol{U}_2}$ and ${\boldsymbol{U}_{\boldsymbol{s}}}$, respectively.

It is necessary to obtain relationships between $\boldsymbol{U}_1^{\boldsymbol{P}}$, $\boldsymbol{U}_2^{\boldsymbol{P}}$ and $\boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{P}}$. First, note that for a static film, ${\boldsymbol{U}^{\boldsymbol{v}}} = 0$, so

(2.27)\begin{equation}{\boldsymbol{U}^{\boldsymbol{v}}} = {\phi _1}{\boldsymbol{U}_1} + {\phi _2}{\boldsymbol{U}_2} + (1 - {\phi _1} - {\phi _2}){\boldsymbol{U}_{\boldsymbol{s}}} = 0.\end{equation}

The components of the velocities due to diffusion already obey

(2.28)\begin{equation}{\boldsymbol{U}^{\boldsymbol{v},\boldsymbol{D}}} = {\phi _1}\boldsymbol{U}_1^{\boldsymbol{D}} + {\phi _2}\boldsymbol{U}_2^{\boldsymbol{D}} + (1 - {\phi _1} - {\phi _2})\boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{D}} = 0.\end{equation}

Hence, subtracting equation (2.28) from (2.27), provides

(2.29)\begin{equation}{\boldsymbol{U}^{\boldsymbol{v},\boldsymbol{P}}} = {\phi _1}\boldsymbol{U}_1^{\boldsymbol{P}} + {\phi _2}\boldsymbol{U}_2^{\boldsymbol{P}} + (1 - {\phi _1} - {\phi _2})\boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{P}} = 0.\end{equation}

The velocity of component two and the velocity of the solvent are related by

(2.30)\begin{equation}\boldsymbol{U}_2^{\boldsymbol{P}} - \boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{P}} = {\boldsymbol{U}_{\boldsymbol{P}}}.\end{equation}

Under the assumptions of this model, diffusiophoresis does not affect the velocity of component one relative to the solvent,

(2.31)\begin{equation}\boldsymbol{U}_1^{\boldsymbol{P}} = \boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{P}}.\end{equation}

This is the case in an idealised model and is therefore used as a starting point for a non-ideal model, with finite size ratios. Although deviation from this will be expected when the particle size ratio is finite (Howard & Nikoubashman Reference Howard and Nikoubashman2020), it would not be expected to be large. Solving (2.29)–(2.31) simultaneously, and including sedimentation coefficients, ${K_{P1}}({\phi _1},{\phi _2})$ and ${K_{P2}}({\phi _1},{\phi _2})$, obtains

(2.32)\begin{equation}\boldsymbol{U}_1^{\boldsymbol{P}} = \boldsymbol{U}_{\boldsymbol{s}}^{\boldsymbol{P}} ={-} {\phi _2}{\boldsymbol{U}_{\boldsymbol{P}}} = \frac{{3{\phi _1}{\phi _2}{K_{P1}}({\phi _1},{\phi _2})}}{{8{\rm \pi}\eta {R_1}}}\boldsymbol{\nabla }{\mu _1},\end{equation}

and

(2.33)\begin{equation}\boldsymbol{U}_2^{\boldsymbol{P}} = (1 - {\phi _2}){\boldsymbol{U}_{\boldsymbol{P}}} ={-} \frac{{3{\phi _1}(1 - {\phi _2}){K_{P2}}({\phi _1},{\phi _2})}}{{8{\rm \pi}\eta {R_1}}}\; \boldsymbol{\nabla }{\mu _1}.\end{equation}

To satisfy continuity, it is required that ${K_{P2}}({\phi _1},{\phi _2}) = {K_{P1}}({\phi _1},{\phi _2}) = {K_P}({\phi _1},{\phi _2})$. Note that this is not a general result; it is just a result of how this model chooses to write the total flux as two separate terms. Addressing the more general question of whether ${K_{12}}({\phi _1},{\phi _2}) = {K_{21}}({\phi _1},{\phi _2})$, these are not equal, as can be seen from analytical expressions for dilute solution (Batchelor Reference Batchelor1983). However, equating these expressions would be a reasonable approximation for generating example solutions. ${K_{ij}}({\phi _1},{\phi _2})$ and ${K_P}({\phi _1},{\phi _2})$ are both denoted with ‘$K$’ since they are both sedimentation coefficients, i.e. they represent the same physical phenomenon of hydrodynamic hindrance. The different subscripts identify the terms for which these are the coefficients.

Equation (2.32) can be interpreted as showing that component one and the solvent move backwards compared with the diffusiophoretic motion of component two, in order to conserve volume. Using (2.12), which gives $\boldsymbol{U}_1^{\boldsymbol{D}}$, and the analogous expression for $\boldsymbol{U}_2^{\boldsymbol{D}}$, the complete expressions for the component velocities including both diffusion and diffusiophoresis are

(2.34)\begin{equation}\begin{array}{c} {\boldsymbol{U}_1} ={-} \dfrac{1}{{6{\rm \pi}\eta {R_1}}}[{K_{11}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _1} + {\phi _2}{K_{12}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _2}]\\ \qquad\ \, + \dfrac{{3{\phi _1}{\phi _2}{K_P}({\phi _1},{\phi _2})}}{{8{\rm \pi}\eta {R_1}}}\boldsymbol{\nabla }{\mu _1}, \end{array}\end{equation}

and

(2.35)\begin{equation}\begin{array}{c} {\boldsymbol{U}_2}={-} \dfrac{1}{{6{\rm \pi}\eta {R_2}}}[{K_{22}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _2} + {\phi _1}{K_{21}}({\phi _1},{\phi _2})\boldsymbol{\nabla }{\mu _1}]\\ \qquad\ \, - \dfrac{{3{\phi _1}(1 - {\phi _2}){K_P}({\phi _1},{\phi _2})}}{{8{\rm \pi}\eta {R_1}}}\boldsymbol{\nabla }{\mu _1}. \end{array}\end{equation}

Including diffusiophoresis, the conservation equation for particles of type one is now

(2.36)\begin{equation}\begin{array}{c} \dfrac{{\partial {\phi _1}}}{{\partial t}} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left[ {\left[ {\dfrac{{{\phi_1}}}{{6{\rm \pi}\eta {R_1}}}{K_{11}}({\phi_1},{\phi_2}) - \dfrac{{3{{({\phi_1})}^2}{\phi_2}{K_P}({\phi_1},{\phi_2})}}{{8{\rm \pi}\eta {R_1}}}} \right]} \right.\boldsymbol{\nabla }{\mu _1}\\ \qquad\ \ \left. { + \dfrac{{{\phi_1}}}{{6{\rm \pi}\eta {R_1}}}{\phi_2}{K_{12}}({\phi_1},{\phi_2})\boldsymbol{\nabla }{\mu_2}} \right], \end{array}\end{equation}

and that for particles of type two becomes

(2.37)\begin{equation}\begin{array}{c} \dfrac{{\partial {\phi _2}}}{{\partial t}} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left[ {\left[ {\dfrac{{{\phi_2}}}{{6{\rm \pi}\eta {R_2}}}{\phi_1}{K_{21}}({\phi_1},{\phi_2}) + \dfrac{{3{\phi_1}{\phi_2}(1 - {\phi_2}){K_P}({\phi_1},{\phi_2})}}{{8{\rm \pi}\eta {R_1}}}} \right]\boldsymbol{\nabla }{\mu_1}} \right.\\ \qquad\ \ \left. { + \dfrac{{{\phi_2}}}{{6{\rm \pi}\eta {R_2}}}{K_{22}}({\phi_1},{\phi_2})\boldsymbol{\nabla }{\mu_2}} \right]. \end{array}\end{equation}

Note that the divergence of the total flux should reflect the total excluded volume, which is incorporated into $Z({\phi _1},{\phi _2})$, hence all the flux terms, including the diffusiophoresis terms, should diverge in the same form as $Z({\phi _1},{\phi _2})$. Since all of the flux terms contain a factor of $\boldsymbol{\nabla }{\mu _1}$ or $\boldsymbol{\nabla }{\mu _2}$, which diverge via $Z({\phi _1},{\phi _2})$ through the formulations for (2.22) and (2.23), (2.36) and (2.37) obtain suitable divergence at close packing.

3.1. Scaling

Equations (2.22) and (2.23), the spatial derivatives of the chemical potentials, are substituted into (2.36), and scaled in one dimension, to form the conservation equation for component one,

(3.1) \begin{align} \dfrac{{\partial {\phi _1}}}{{\partial \hat{t}}} & = \dfrac{1}{{P{e_1}}}\dfrac{\partial }{{\partial \hat{z}}}[{(1 - {\phi_1} - {\phi_2})} \notag\\ & \quad \left[ {\dfrac{{\left[ {{\phi_1}{K_{11}}({\phi_1},{\phi_2}) - \dfrac{9}{4}{{({\phi_1})}^2}{\phi_2}{K_P}({\phi_1},{\phi_2})} \right]}}{{\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}\dfrac{{\partial {\mu_2}/\partial \hat{z}}}{{\partial {\mu_1}/\partial \hat{z}}}} \right)}} + \dfrac{{{\phi_1}{\phi_2}{K_{12}}({\phi_1},{\phi_2})}}{{\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1}\dfrac{{\partial {\mu_1}/\partial \hat{z}}}{{\partial {\mu_2}/\partial \hat{z}}} + {\phi_2}} \right)}}{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}} \right]\notag\\ & \quad \left. {\dfrac{\partial }{{\partial \hat{z}}}\left[ {\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}} \right)Z({\phi_1},{\phi_2})} \right]\; } \right]. \end{align}

For component two, substituting equations (2.22) and (2.23) into (2.37) and scaling gives

(3.2) \begin{align} \dfrac{{\partial {\phi _2}}}{{\partial \hat{t}}} & = \dfrac{1}{{P{e_2}}}\dfrac{\partial }{{\partial \hat{z}}}\left[\vphantom{\left(\dfrac{{P{e_2}}}{{P{e_1}}}\right)}{(1 - {\phi_1} - {\phi_2})} \right.\notag\\ & \quad \left. \left[ {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}\left[ {\dfrac{{{\phi_1}{\phi_2}{K_{21}}({\phi_1},{\phi_2}) + \dfrac{9}{4}\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right){\phi_1}{\phi_2}(1 - {\phi_2}){K_P}({\phi_1},{\phi_2})}}{{\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}\dfrac{{\partial {\mu_2}/\partial \hat{z}}}{{\partial {\mu_1}/\partial \hat{z}}}} \right)}}} \right]\right.\right. \notag\\ &\quad + \left.\left. \dfrac{{{\phi_2}{K_{22}}({\phi_1},{\phi_2})}}{{\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1}\dfrac{{\partial {\mu_1}/\partial \hat{z}}}{{\partial {\mu_2}/\partial \hat{z}}} + {\phi_2}} \right)}} \right]\right. \notag\\ & \quad \left. {\; \dfrac{\partial }{{\partial \hat{z}}}\left[ {\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1} + {\phi_2}} \right)Z({\phi_1},{\phi_2})} \right]} \right]. \end{align}

A key observation is that the only dimensionless groups appearing in (3.1) and (3.2) are $P{e_1}$ and $P{e_2}$ i.e. no further dimensionless groups appear due to diffusiophoresis, as would be expected from dimensional analysis.

3.2. Coordinate transform

In order to create a static top boundary, the conservation equations are transformed using $\xi = \hat{z}/(1 - \hat{t})$ and $\tau = \hat{t}$. This results in

(3.3) \begin{align} & \dfrac{{\partial {\phi _1}}}{{\partial \tau }} + \dfrac{\xi }{{1 - \tau }}\dfrac{{\partial {\phi _1}}}{{\partial \xi }} = \dfrac{1}{{P{e_1}{{(1 - \tau )}^2}}}\dfrac{\partial }{{\partial \xi }}[{(1 - {\phi_1} - {\phi_2})} \notag\\ & \quad \left[ {\dfrac{{\left[ {{\phi_1}{K_{11}}({\phi_1},{\phi_2}) - \dfrac{9}{4}{{({\phi_1})}^2}{\phi_2}{K_P}({\phi_1},{\phi_2})} \right]}}{{\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}\dfrac{{\partial {\mu_2}/\partial \xi }}{{\partial {\mu_1}/\partial \xi }}} \right)}} + \dfrac{{{\phi_1}{\phi_2}{K_{12}}({\phi_1},{\phi_2})}}{{\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1}\dfrac{{\partial {\mu_1}/\partial \xi }}{{\partial {\mu_2}/\partial \xi }} + {\phi_2}} \right)}}{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}} \right]\notag\\ & \quad \left. {\dfrac{\partial }{{\partial \xi }}\left[ {\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}} \right)Z({\phi_1},{\phi_2})} \right]\; } \right], \end{align}

and

(3.4) \begin{align}& \dfrac{{\partial {\phi _2}}}{{\partial \tau }} + \dfrac{\xi }{{1 - \tau }}\dfrac{{\partial {\phi _2}}}{{\partial \xi }} = \dfrac{1}{{P{e_2}{{(1 - \tau )}^2}}}\dfrac{\partial }{{\partial \xi }}\left[\vphantom{\left(\dfrac{{P{e_2}}}{{P{e_1}}}\right)}{(1 - {\phi_1} -{\phi_2})} \right. \notag\\ &\quad\left. \left[ {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}\left[ {\dfrac{{{\phi_1}{\phi_2}{K_{21}}({\phi_1},{\phi_2}) + \dfrac{9}{4}\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right){\phi_1}{\phi_2}(1 - {\phi_2}){K_P}({\phi_1},{\phi_2})}}{{\left( {{\phi_1} + {{\left( {\dfrac{{P{e_1}}}{{P{e_2}}}} \right)}^3}{\phi_2}\dfrac{{\partial {\mu_2}/\partial \xi }}{{\partial {\mu_1}/\partial \xi }}} \right)}}} \right]\right.\right. \notag\\ &\quad + \left. \left. \dfrac{{{\phi_2}{K_{22}}({\phi_1},{\phi_2})}}{{\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1}\dfrac{{\partial {\mu_1}/\partial \xi }}{{\partial {\mu_2}/\partial \xi }} + {\phi_2}} \right)}} \right]\right. \notag\\ & \quad \left. {\; \dfrac{\partial }{{\partial \xi }}\left[ {\left( {{{\left( {\dfrac{{P{e_2}}}{{P{e_1}}}} \right)}^3}{\phi_1} + {\phi_2}} \right)Z({\phi_1},{\phi_2})} \right]} \right]. \end{align}

This is the coordinate transform that is classically used in solving problems with this geometry, as in Routh & Zimmerman (Reference Routh and Zimmerman2004), Trueman et al. (Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b) and Howard et al. (Reference Howard, Nikoubashman and Panagiotopoulos2017b). The boundary conditions are no flux of particles at both the top and bottom boundaries (see Section S1.1 of the SI). For the top boundary condition, the evaporation rate is assumed to be constant, as explained in § 2.2.

Setting ${K_P}({\phi _1},{\phi _2}) = 0$ in (3.3)–(3.4) gives a system of PDEs to describe a diffusion-only model that neglects diffusiophoresis. This can be compared with the diffusion–diffusiophoresis model, with non-zero ${K_P}({\phi _1},{\phi _2})$, to observe the effect of diffusiophoresis. Discussion regarding satisfying the Onsager reciprocal relations with each model is provided in the SI, § S3.3.

3.3. Numerical method

The system of (3.3)–(3.4) can be solved if ${\phi _{1,t = 0}}$, ${\phi _{2,t = 0}}$, $P{e_1}$, $P{e_2}$ and the maximum volume fraction, ${\phi _m}$, which could be a function of ${\phi _1}$ and $\; {\phi _2}$, are specified. Appropriate forms of $K({\phi _1},{\phi _2})$ and $Z({\phi _1},{\phi _2})$ can be taken from the literature. Also needing to be specified is the ratio of the chemical potential gradients, as explained in § 2.3.

The compressibility is taken to be $Z({\phi _1},{\phi _2}) = {\phi _m}{({\phi _m} - {\phi _1} - {\phi _2})^{ - 1}}$ (Trueman et al. Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b). The canonical concentrated extension of the dilute form of $K(\phi )$ for the one-component case, $1 - 6.55\phi $, is ${(1 - \phi )^{6.55}}$ (Russel et al. Reference Russel, Saville and Schowalter1989). For the two-component case, this is extended to ${K_{ij|i,j = 1\;\textrm{or}\;2}}({\phi _1},{\phi _2}) = {(1 - {\phi _1} - {\phi _2})^{6.55}}$ (Trueman et al. Reference Trueman, Lago Domingues, Emmett, Murray and Routh2012b). This is adopted for the example results in this work. The particle concentration up to which this is applicable is discussed in § S3.1 of the SI. It is shown in § S3.2 of the SI that the stratification results are minimally affected by how the form of K is extended.

A code was written in MATLAB to solve the PDEs numerically. The governing PDEs are of a similar form to diffusion–advection equations, where $\xi /(1 - \tau )$ is the effective advection velocity, but the effective diffusion coefficient is a complicated function of ${\phi _1}$ and ${\phi _2}$. The PDEs are solved numerically using a finite volume method for the transformed diffusion term. This conserves mass by construction and is an appropriate choice for results containing sharp fronts. Due to the $\xi $-dependence of its coefficient, the quasi-convective term cannot be formulated using the finite volume method, so finite differences are employed instead. Since the $\partial \phi /\partial \xi $ term represents advection due to the stretching of the coordinate frame, this term is coded using backward finite differences. This ensures that information travels from upstream to downstream, avoiding potential stability issues that might result from using forward or central finite differences instead. The Euler method is used for time-stepping. The numerical results are verified by checking the final mass balance for each component, and by comparison with asymptotic solutions for high $Pe$ (§ 5).

5.1. Diffusion

For the case where ${K_{12}}({\phi _1},{\phi _2}) = {K_{21}}({\phi _1},{\phi _2}) = 0$ and the chemical potentials are just determined by entropy, $\textrm{d}{\mu _1}/\textrm{d}{\mu _2}\; = \textrm{d}\ln {\phi _1}/\textrm{d}\ln {\phi _2}$, and employing the limit where ${(P{e_2}/P{e_1})^3} \gg 1$, it is shown in the SI that there is a unique relation linking ${\phi _1}$ and ${\phi _2}$

(5.2)\begin{equation}\frac{{{\phi _1} - {\phi _{1,\tau = 0}}}}{{\frac{{{\phi _m}{\phi _{1,\tau = 0}}}}{{({\phi _{1,\tau = 0}} + {\phi _{2,\tau = 0}})}} - {\phi _{1,\tau = 0}}}} = {\left[ {\frac{{{\phi_2} - {\phi_{2,\tau = 0}}}}{{\frac{{{\phi_m}{\phi_{2,\tau = 0}}}}{{({\phi_{1,\tau = 0}} + {\phi_{2,\tau = 0}})}} - {\phi_{2,\tau = 0}}}}} \right]^{P{e_1}/P{e_2}}}.\end{equation}

This can be compared with numerical results. Note that (5.2) is independent of ${K_{11}}({\phi _1},{\phi _2})$, ${K_{22}}({\phi _1},{\phi _2})$, $Z({\phi _1},{\phi _2})$ and $\tau $.

The result in (5.2) can be substituted into the expression for $\textrm{d}{\phi _1}/\textrm{d}X$ and $\textrm{d}{\phi _2}/\textrm{d}X$, giving $X({\phi _1})$ or $X({\phi _2})$. For the boundary condition, we arbitrarily set $X = 0$ at $({\phi _{1,X = 0}} + {\phi _{2,X = 0}}) = [({\phi _{1,\tau = 0}} + {\phi _{2,\tau = 0}}) + ({\phi _{1,X \to \infty }} + {\phi _{2,X \to \infty }})]/2$. This sets the asymptotic coordinate system to be centred at the transition.

The resulting asymptotic solution for an example set of Péclet numbers is shown in figure 8. In the transition region, there is excellent agreement between the asymptotic and numerical solutions. This validates the numerical code used to generate figures 3–5. When ${\phi _{1,\tau = 0}} = {\phi _{2,\tau = 0}}$, the asymptotic solution predicts that ${\phi _{1,X \to \infty }} = {\phi _{2,X \to \infty }}$ (SI). Whilst there is therefore disagreement in the values at the top surface, these do tend to the values given by the asymptotic solution over time. The small excess of large particles at the top surface, a diffusive effect, counterbalances the excess of small particles in the rest of the film (see (5.2)), to satisfy the mass balance.

Figure 8. Volume fraction of each component as a function of height for $P{e_1} = 10$ and $P{e_2} = 20$, with ${\phi _{1,\hat{t} = 0}} = {\phi _{2,\hat{t} = 0}} = 0.10$ and ${\phi _m} = 0.64$. Both the numerical and asymptotic solutions are shown for comparison.

5.2. Diffusiophoresis

As for the diffusion-only results, the diffusiophoresis results at high $P{e_1}$ and $P{e_2}$ appear to be curves of constant shape moving through the film. The numerical code is verified using the equations with enhanced diffusiophoresis. This is because the main purpose of this asymptotic solution is to specifically verify the coding of the diffusiophoresis term.

The approach to deriving the asymptotic solution for the non-enhanced diffusiophoresis case is very similar to that described in § 5.1. However, there are two key transitions in the curves for the enhanced diffusiophoresis case: The first is the position of the particle front i.e. the front of total volume fraction. This is denoted by $P(\tau )$, as in the diffusion-only case. The second is between the region of constant high ${\phi _1}$ and low ${\phi _2}$ values and the region of constant high ${\phi _2}$ and low ${\phi _1}$. The position of this transition is denoted by $Q(\tau )$. The transitions are annotated on figure 9. This section seeks to construct an asymptotic solution for all parts of the film:

• • the position of both transitions, $P(\tau )$, as already described, and $Q(\tau )$;

• • the concentration profiles around $P(\tau )$; and

• • the values of ${\phi _1}$ and ${\phi _2}$ either side of $Q(\tau )$.

The full derivation is shown in the SI from § S2.2. Briefly, the change of variable $X = [\xi - P(\tau )]P{e_2}$ forms differential equations for $X({\phi _1},{\phi _2})$, which can be solved to give the profile in the film region around $P(\tau )$, from the bottom of the film up to just beneath $Q(\tau )$.

Figure 9. Volume fraction of each component as a function of height for $P{e_1} = 10$ and $P{e_2} = 20$, with ${\phi _{1,\hat{t} = 0}} = {\phi _{2,\hat{t} = 0}} = 0.10$ and ${\phi _m} = 0.64$. Both the numerical and asymptotic solutions are shown for comparison. The model includes enhanced diffusiophoresis. The positions of $P(\tau )$ and $Q(\tau )$ at $\tau = 0.34$ are indicated as examples.

Unlike the diffusion case, in the limit ${(P{e_2}/P{e_1})^3} \gg 1$, the ODEs cannot be approximated further. Since, from the numerical results, $\textrm{d}{\phi _2}/\textrm{d}{\phi _1}$ appears to be large, the product ${(P{e_1}/P{e_2})^3}{\phi _2}(\textrm{d}{\mu _2}/\textrm{d}X)/(\textrm{d}{\mu _1}/\textrm{d}X)$ cannot be neglected. For the same reason, it cannot be assumed that $|{(P{e_2}/P{e_1})^3}{\phi _1}(\textrm{d}{\mu _1}/\textrm{d}X)/(\textrm{d}{\mu _2}/\textrm{d}X)|\gg {\phi _2}$. For the case ${K_{12}}({\phi _1},{\phi _2}) = {K_{21}}({\phi _1},{\phi _2}) = 0$ and ${K_{11}}({\phi _1},{\phi _2}) = {K_{22}}({\phi _1},{\phi _2}) = {K_P}({\phi _1},{\phi _2})$, and no enthalpic contribution to the chemical potentials, $\textrm{d}{\mu _1}/\textrm{d}{\mu _2}\; = \textrm{d}\ln {\phi _1}/\textrm{d}\ln {\phi _2}$, equating the asymptotic expressions for $\textrm{d}{\phi _1}/\textrm{d}X$ and $\textrm{d}{\phi _2}/\textrm{d}X$ yields a quadratic in $\textrm{d}{\phi _1}/\textrm{d}{\phi _2}$. This equation is independent of ${K_{ij}}({\phi _1},{\phi _2})$, ${K_P}({\phi _1},{\phi _2})$, $Z({\phi _1},{\phi _2})$ and $\tau $. The quadratic equation can be solved for $\textrm{d}{\phi _1}/\textrm{d}{\phi _2}$. The resulting solution for $\textrm{d}{\phi _1}/\textrm{d}{\phi _2}$ is integrated numerically to give ${\phi _1}({\phi _2})$ around $P(\tau )$.

Finding the values of ${\phi _1}$ and ${\phi _2}$ either side of the transition $Q(\tau )$ will provide a boundary condition for integrating $\textrm{d}{\phi _1}/\textrm{d}{\phi _2}$ around $P(\tau )$, as well as allowing the asymptotic solution for the upper part of the film to be constructed. The change of variable $Y = [\xi - Q(\tau )]P{e_2}$ is used to expand around this transition. It is assumed that $Q(\tau )$ takes the form

(5.3)\begin{equation}Q(\tau ) = \frac{{1 - v\tau }}{{1 - \tau }},\end{equation}

which is a front travelling at constant speed v in $(\hat{z},\; \hat{t})$ coordinates. The speed v at which $Q(\tau )$ travels can be found by integrating the equation for $\textrm{d}{\phi _1}/\textrm{d}Y$ or $\textrm{d}{\phi _2}/\textrm{d}Y$ across the transition (SI, S2.2.6). This is equivalent to a mass balance on ${\phi _1}$ or ${\phi _2}$ across $Q(\tau )$. Physically, both the diffusion and diffusiophoretic fluxes work to propel this front, changing $({\phi _1},{\phi _2})$ from the values on one side of the front to those on the other. Note that it is expected that $v < {\phi _m}/[{\phi _m} - ({\phi _{1,\tau = 0}} + {\phi _{2,\tau = 0}})]$, since $Q(\tau )$ is located higher in the film than $P(\tau )$.

The ODEs for $\textrm{d}{\phi _1}/\textrm{d}Y$ and $\textrm{d}{\phi _2}/\textrm{d}Y$ can be integrated, and are combined for the case where ${K_{12}}({\phi _1},{\phi _2}) = {K_{21}}({\phi _1},{\phi _2}) = 0$, ${K_{11}}({\phi _1},{\phi _2}) = {K_{22}}({\phi _1},{\phi _2}) = {K_P}({\phi _1},{\phi _2})$ and the particle chemical potentials are purely entropic, $\textrm{d}{\mu _1}/\textrm{d}{\mu _2}\; = \textrm{d}\ln {\phi _1}/\textrm{d}\ln {\phi _2}$. In the upper region of the film, ${\phi _1} + {\phi _2} \approx {\phi _m}$. Substituting in ${\phi _2} = {\phi _m} - {\phi _1}$ gives a quadratic equation for the values of ${\phi _1}$ either side of $Q(\tau )$.

Using ${K_{11}}({\phi _1},{\phi _2}) = {(1 - {\phi _1} - {\phi _2})^{6.55}}$ and $Z({\phi _1},{\phi _2}) = {\phi _m}{({\phi _m} - {\phi _1} - {\phi _2})^{ - 1}}$, figure 9 is obtained, which compares the resulting asymptotic solution with the numerical results. Figure 9 uses $P{e_1} = 10$ and $P{e_2} = 20$ i.e. the same parameters as figure 8, but now with the diffusiophoresis term included.

There is reasonably good agreement between the numerical and asymptotic solutions, although the asymptotic solution slightly overpredicts the sharpness of the transition around $P(\tau )$. Using the numerical results for ${\phi _1}$ to reconstruct the ${\phi _2}$ curves with the asymptotic solution, and vice versa, gives excellent agreement. Therefore, the small discrepancy that does arise is mostly from the step of finding $X({\phi _1})$.

At fixed values of ${\phi _{1,\tau = 0}}$, ${\phi _{2,\tau = 0}}$ and ${\phi _m}$, the asymptotic solution for v shows minimal variation with $P{e_1}$ and $P{e_2}$. This allows a schematic diagram to be drawn for the regions in the film, as exemplified in figure 10 for the set of ${\phi _{1,\hat{t} = 0}}$, ${\phi _{2,\hat{t} = 0}}$ and ${\phi _m}$ values that were used in figure 9.

Aside from verifying the numerical code, the asymptotic solutions shed insight into the large $Pe$ regime. Remembering the definition of the Péclet number, $Pe = 6{\rm \pi}\eta R\dot{E}H/kT$, this corresponds to evaporation being fast compared with diffusion. Highly stratified films are predicted in this regime, with the asymptotic solutions allowing prediction of the composition of each region of the film. Figure 10 shows the depths of both (i) the region between the top surface and $Q(\tau )$, and (ii) the region between $Q(\tau )$ and $P(\tau )$, increasing linearly with $\hat{t}$. The depth of the region between the top surface and $Q(\tau )$ is the predicted thickness of the layer of small particles at the top surface. It can be seen in figure 10 that this small particle-rich layer is thinner than the large particle-rich layer beneath it, as a result of the relative speeds at which the top surface, $Q(\tau )$ and $P(\tau )$ recede. These fronts recede at speeds of 1, v and ${\phi _m}/[{\phi _m} - ({\phi _{1,\tau = 0}} + {\phi _{2,\tau = 0}})]$, respectively, in $(\hat{z},\hat{t})$ coordinates. The speeds are presented as the magnitudes of the gradients of the diagonal lines separating the regions in figure 10. Hence the schematic diagram may provide an explanation for why a thin layer of small particles can form at the top surface of a drying film.

Figure 10. Schematic of the regions of the asymptotic solution for large $Pe$, when ${\phi _{1,\hat{t} = 0}} = {\phi _{2,\hat{t} = 0}} = 0.10$ and ${\phi _m} = 0.64$. The unreachable region is shaded in grey.