## Abstract

During the drying of colloidal suspensions, the desiccation process causes the suspension near the air interface to consolidate into a connected porous matrix or crust. Fluid transport in the porous medium is governed by Darcy’s law and the equations of poroelasticity, while the equations of colloid physics govern processes in the suspension. We derive new equations describing this process, including unique boundary conditions coupling the two regions, yielding a moving-boundary model of the concentration and stress profiles during drying. A solution is found for the steady-state growth of a one-dimensional crust during constant evaporation rate from the surface. The solution is used to demonstrate the importance of the system boundary conditions on stress profiles and diffusivity in a drying crust.

## 1. Introduction

Crust formation during drying is an important phenomenon in many areas of research such as soil science, the drying of colloidal suspensions and the food sciences. Often, the crust characteristics can have controlling influence on the physics of a drying dispersion. For example in a soil the small permeability of the crust relative to the underlying wet soil means that the rate of drying of the soil is controlled by the crust thickness, as will be the rate of penetration of rainwater. Stress build ups in the crust can also lead to buckling or cracking of the crust (e.g. Brinker & Scherer 1990; Tsapis *et al.* 2005) and so crust formation is also of strong relevance to materials such as paints, concrete and ceramics, which are specifically designed to avoid cracking and the associated degradation of the material.

In this paper, we model a drying colloidal dispersion in which the solid particles are initially stably dispersed in the liquid phase as a colloidal sol. As liquid evaporates from the sol the concentration of particles increases until the particles gel into a saturated porous network at a critical concentration determined by inter-particle forces. In the sol phase, the concentration of solid particles can be described by a diffusion equation with non-constant diffusivity as previously shown by Russel *et al.* (1989). In the gel phase, we show that the concentration of solid particles is described by a similar diffusion equation, but with a modified diffusivity that depends upon the boundary conditions of the problem. When the gel is unstressed during drying the resulting diffusion equation is similar to that used by Brown *et al.* (2002). However when the gel is constrained during drying, the diffusivity is modified and tensile stresses arise which are the source of primary, ‘craquelure’ type cracks (Bohn *et al.* 2005). We use the model to examine the growth of a steadily thickening gel in the case of a constant evaporation rate from the drying surface, and this enables us to quantitatively examine the pressure, stress and solid phase profiles that arise during drying in both unstressed and constrained gels.

## 2. Governing equations in the colloidal system

We consider a sol consisting of hard-sphere particles with short-ranged attractive forces (‘sticky spheres’) that undergo dense coagulation upon gelation (see table 1 for a list of parameters relevant to the system). The sol is stable up to a solid phase fraction *ϕ*_{c} at which point it forms a gel of solid phase fraction *ϕ*_{g}. Gelation occurs at a much faster rate than the drying process so the sol is treated as immediately gelling once it has dried above the critical phase fraction *ϕ*_{c}. Owing to the short-ranged attractive interparticle forces and to the close-packed nature of the particles, we assume that the solid network behaves elastically under stress and does not compress irreversibly like a flocculated gel. An example of such a system is given by adhesive hard spheres proposed by Baxter (1968), and studied experimentally by Buzzaccaro *et al.* (2007).

We consider a system consisting initially of a homogeneous, stable sol in contact with air. If the vapour pressure of the air is less than the saturation vapour pressure at the sol/air interface then evaporation occurs which increases the concentration of colloid particles near the surface of the sol while particle diffusion increases the concentration of particles in the bulk of the sol. Provided the air is dry enough the concentration of particles at the interface eventually rises above the critical volume phase fraction *ϕ*_{c}, and a saturated gel ‘crust’ forms at the surface of the sol. This is shown schematically in figure 1. Further drying then causes thickening of the crust until finally air invades the pores of the gel. We consider the period before air invades the pores of the gel phase, so the gel is always saturated with liquid. This period is of interest as it comprises the time during drying in which tensile stresses arise (which can be responsible for cracking if sufficiently large). Stresses typically peak at the point at which air invades the pores of the gel and then drop off as the gel desaturates (e.g. Martinez & Lewis 2002).

In both the gel and the sol we define *ϕ* to be the local solid phase fraction and **v**_{1} and **v**_{2} to be the average velocities of the liquid and the particles relative to the laboratory frame. The volume-averaged velocity is given by2.1We take the solid and liquid phases to be incompressible, so conservation of volume of the liquid in the suspension means that the water content of a small volume can only vary if water is advected into or out of the volume. This implies that2.2while conservation of total volume implies that the total flux into or out of any volume must be zero, or equivalently2.3

### (a) Equations in the gel phase

In the gel phase, flow of the liquid through the solid matrix is given by Darcy’s law2.4where *k*_{g}(*ϕ*) is the permeability of the gel, *μ* is the dynamic viscosity of the liquid and *p* is the pore pressure of the liquid phase. Then by use of equations (2.1), (2.3) and (2.4) in the equation for conservation of liquid volume (2.2) we find that2.5Provided *p*=*p*(*ϕ*) this means that the equation governing the evolution of the solid phase fraction in the gel can be written as2.6where the gel diffusivity is given by2.7The form of the permeability *k*(*ϕ*) depends upon the qualities of the gel. For a regularly packed gel the permeability can be assumed to be given by the semi-empirical Carmen–Kozeny expression2.8where *R* is the radius of the particles, and *A* is the Kozeny constant.

### (b) Equations in the sol phase

For spherical particles in the sol phase, the motion of the particles relative to the suspending medium is given by2.9where *υ*_{p} is the volume of an individual particle and *Π* is the osmotic pressure of the particles in the sol: the pressure that particles in the sol would exert on a confining, water-permeable membrane (Russel *et al.* 1989; Peppin *et al.* 2005). This is the generalization of Einstein’s equation describing the sedimentation of dilute spheres relative to a fluid to more concentrated dispersions (Davis & Russel 1989) with a dimensionless friction factor *f* to account for particle–particle interactions. We define the effective permeability of the sol as2.10and then by inserting this into equation (2.9) and making use of equation (2.3), we find that equation (2.9) becomes equivalent to Darcy’s law (2.4), only modified by the pressure gradient of Darcy’s law (∇*p*) being replaced by the negative gradient of the osmotic pressure (−∇*Π*).

For the case of a suspension of hard spheres, the permeability can be described by the empirical expression2.11(Russel *et al.* 1989). Then using equations (2.1) and (2.3) and the modified Darcy’s law for the sol phase (2.9) in the equation for conservation of liquid volume (2.2), we find that the solid fraction in the sol satisfies the diffusion equation (2.6) with diffusivity2.12

### (c) Boundary conditions

As an example of the use of these equations we consider the one-dimensional system shown in figure 1. The co-ordinate system is chosen so that the gel–air interface corresponds to *z*=0, and the sol–gel interface is at *z*=−*h*(*t*). As the system is one dimensional, if we define the mass flux rate of liquid evaporating per unit area of the surface as *ρ*_{l}*E* (where *ρ*_{l} is the density of the liquid), then conservation of volume implies that **u**=*E***z** throughout the sol–gel system, where **z** is the unit vector in the *z* direction. Also at the gel–air interface, conservation of solid implies that we must set *v*_{2}=0 so2.13At the sol–gel interface there is a phase change from the dispersed sol to the aggregated gel. Thus2.14and2.15(subscripts g and s correspond to values in the gel phase and sol phase, respectively). Also at the sol–gel boundary conservation of solid gives that which becomes2.16Two further equations arising from matching chemical potentials of the solvent and particle phases at the sol–gel interface will determine *ϕ*_{c} and *ϕ*_{g} and these are discussed in the next section.

Finally, we consider a sol with a boundary far from the sol–gel interface so that the far-field boundary condition becomes2.17where *ϕ*_{0}<*ϕ*_{c} is the initial concentration of the sol. This final boundary condition can be simply modified to a zero solid-particle flux condition for the case of an impermeable boundary bounding a finite sol.

Note that we have chosen to consider the system as isothermal. At the evaporating interface cooling will occur owing to the latent heat of evaporation, while a small amount of heat will be released as the sol crystalizes into a gel. However, we shall see that the colloidal diffusivity *D* is typically much smaller than the thermal diffusivity of the system *κ*∼*O*(10^{−7} m^{2} s^{−1}), which means that we can neglect temperature variations across the film because heat diffuses away much faster than the gelation process occurs (e.g. Style & Wettlaufer 2007).

## 3. Crust formation with an unstressed gel

Comparing the diffusivities of the solid phase fraction *ϕ* in the gel and sol phases given by equations (2.7) and (2.12), respectively, it can be seen that the key difference is that the diffusivity in the sol phase is dependent upon the change in the *osmotic* pressure with *ϕ*, while the diffusivity in the gel phase is dependent upon the negative change in the pore pressure with *ϕ*. From the definition of the osmotic pressure when the system is in equilibrium with an external pressure *P*,3.1When the gel and sol are in equilibrium with the surrounding air, *P* in both sol and gel becomes the atmospheric pressure, *P*_{a}, ∂*p*/∂*ϕ*=−∂*Π*/∂*ϕ* and the expressions for the diffusivities (2.7) and (2.12) become equivalent. However, we shall see later that this equivalence only holds under certain conditions. In particular, when the gel is constrained by rigid boundaries so that anisotropic stresses on the gel are maintained, ∂*p*/∂*ϕ*≠−∂*Π*/∂*ϕ*. Therefore, the distinction between the form of the two diffusivities is an important one to make.

We initially consider the unstressed system where the diffusivity expressions (2.7) and (2.12) are identical. In order to solve the system of equations we need to know the relationship between *Π* and *ϕ*. For a system of sticky hard spheres the osmotic pressure is given by the expression3.2(Russel *et al.* 1989), where *k*_{B} is Boltzmann’s constant and *T* is the temperature, while *Z*(*ϕ*,*τ*) is the compressibility factor of the system which depends upon the solid phase fraction and a ‘stickiness parameter’ *τ* that represents the effect of changing the interparticle forces (Buzzaccaro *et al.* 2007). For a system of adhesive hard spheres the compressibility factor in the sol has been shown theoretically and experimentally to be given by3.3(Miller & Frenkel 2004; Buzzaccaro *et al.* 2007). The first term on the right-hand side is the Carnahan–Starling equation for the compressibility factor of a hard sphere suspension, while the second term is due to short-ranged attractive forces between the particles. For simplicity we choose *Z*_{adh}(*ϕ*,*τ*)=0 (weak adhesive forces between the particles), although non-zero *Z*_{adh} can readily be accounted for by insertion in the following equations.

In the gel it has been shown experimentally that the compressibility factor is well described by3.4where *λ*(*τ*) is a scaling factor between 0 and 1 that reduces with increasing ‘stickiness’ of the particles, and takes a value of 1 for a standard hard spheres colloidal crystal (Buzzaccaro *et al.* 2007). *ϕ*_{cp} is the close-packing fraction, and this can either refer to hexagonal close-packing fraction () or random close packing (*ϕ*_{cp}≈0.64) depending upon the conditions under which the gel is formed. Typically, slow aggregation leads to hexagonal close packing, while fast aggregation leads to random close packing. Here, we assume the drying rate is sufficiently slow that the gel is hexagonally close packed.

Inserting the above expressions for the osmotic pressures of the sol and the gel, along with the expressions for the permeabilities of the two phases (2.11) and (2.8), the diffusivities in each phase become
3.5and
3.6
while the advection–diffusion equation (2.6) in one dimension is3.7with *D*_{0}=*k*_{B}*T*/6*πμR*.

Figure 2 shows the non-dimensional permabilities and diffusivities given by equations (2.8) and (2.11) and (3.5) and (3.6), respectively. The permeabilities, shown in figure 2*a*, have been non-dimensionalized with the particle radius squared so it can be seen that typical permeabilities for a system of particles of 10 nm particles are on the order of 10^{−18} *m*^{2}. The diffusivities, shown in figure 2*b*, are non-dimensionalised with *D*_{0}≈2.2×10^{−19}/*R*. Therefore typical diffusivities for 10 nm particles and *ϕ*<0.65 are on the order of 2.2×10^{−11} *m*^{2} *s*^{−1}.

The expressions for the osmotic pressure in the two different phases also allow us to determine *ϕ*_{c} and *ϕ*_{g} by equating solvent and particle chemical potentials at the sol–gel interface. Then3.8and3.9where Δ*s*(*ϕ*^{0}) is the entropy of fusion per particle upon crystallization at a concentration *ϕ*^{0} (e.g. Santore *et al.* 1990). These can be solved by choosing an arbitrary *ϕ*^{0} and requiring and as in the hard sphere limit in order to determine Δ*s*(*ϕ*^{0}) (these values have been observed experimentally by Russel *et al.* 1989). The results show that *ϕ*_{c} and *ϕ*_{g} vary relatively little for small amounts of ‘stickiness’ and so we take the values *ϕ*_{c}=0.494 and *ϕ*_{g}=0.545 in the following calculations.

### (a) Steady-state growth of a gel crust

We consider the initial period of evaporation in which the evaporation rate per unit area of the drying surface is a constant *E*, and which persists up to the ‘critical point’ at which the stiffness of the gel causes shrinkage to stop, and the liquid–air menisci penetrate into the interior of the gel causing the evaporation rate to fall (Brinker & Scherer 1990). We seek a steady-state solution where the interface is moving at a constant rate and so we move the co-ordinate system such that *z*=0 corresponds to the sol–gel interface and the gel–air interface is situated at *z*=*h*(*t*). The advection–diffusion equation then becomes3.10where a dot refers to a derivative with respect to time. We non-dimensionalize the system by setting *z*=(*D*_{0}/*E*)*ζ*, *t*′=(*D*_{0}/*E*^{2})*t*, and . Then by setting , the advection–diffusion equation becomes3.11with corresponding to the relevant non-dimensional diffusivity in the sol or gel phase.

In the gel phase as the osmotic pressure increases due to the divergent nature of the compressibility factor (3.4). So assuming the sol–air boundary is far from the sol–gel interface, we use the condition in equation (2.13) (which takes the same form in the new co-ordinate system) to find that3.12(we will prove that is an appropriate boundary condition *a posteriori*).

In the steady-state system, we let ∂*ϕ*/∂*t*′=0 and then integrate equation (3.11) once using equations (3.12) and (2.17) to find that3.13in the gel, while3.14in the sol. Then, inserting these results into the boundary condition for conservation of mass at the sol–gel interface (2.16) (after non-dimensionalisation and transformation into the new co-ordinate system), we find that3.15Physically, this equation expresses conservation of mass across the system during steady gel growth, as it is the result of equating the amount of solid advected towards the sol–gel boundary per unit time in the far-field sol, *E*(1+*β*)*ϕ*_{0}, to the amount of solid added to the gel per unit time, *Eβϕ*_{cp}. For a typical system with *ϕ*_{cp}≈0.74, and *ϕ*_{0}=0.2, the non-dimensional crust thickening rate is *β*≈0.37 so the gelation front will travel at 0.37 times the rate of evaporation, *E*.

We can solve these equations to give the solid phase fraction profile in the sol and the gel, using typical values for a colloidal system as given in table 2. The results are shown in figure 3. It can be seen that in the gel outside the boundary layer by the sol–gel interface, justifying our previous assumption that . In the sol, *ϕ* decays across a boundary layer to the initial value solid phase fraction *ϕ*_{0}.

### (b) Discussion

The results illustrated in figure 3 show that in the gel the solid phase fraction of the gel becomes close packed so that *ϕ*=*ϕ*_{cp} outside of a boundary layer of thickness *O*(*d*), where *d*=*D*_{0}/*E*. Similarly in the sol the solid phase fraction drops to the initial phase fraction *ϕ*_{0} outside of a similar boundary layer of thickness *O*(*d*). Colloidal particles are usually defined as being between 1 nm and 1 μm in radius, so using the values in table 2 we find that the boundary layer thickness *d* lies between 0.2 μm for large micron-sized colloids and 220 μm for small, nanometre-sized colloids. For many typical drying experiments (e.g. Dufresne *et al.* 2003) the thickness of the crust is on a scale of millimetres rather than micrometres and so this means that in a densely gelling colloidal system we can treat the bulk of the crust and the sol as being at constant phase fraction *ϕ*_{cp} or *ϕ*_{0} respectively. In particular, we note that for the larger, micron-sized particles the boundary layer thickness *d* is less than the radius of the particles so that we can effectively treat the sol–gel interface as a sharp interface between a close packed gel and a sol with phase fraction *ϕ*_{0}.

In this system, the boundary layer thicknesses are both *O*(*D*_{0}/*E*) because both diffusivities are of *O*(*D*_{0}). Smaller diffusivities, which might occur due to an increased particle stickiness, will result in thinner boundary layers. On the other hand, larger diffusivities, which can arise from including charges on the particles, will result in thicker booundary layers. Such effects can be analysed by substitution of appropriate expressions for *Z*_{s}(*ϕ*) and *Z*_{g}(*ϕ*) in the above equations.

In addition to calculating the solid-fraction profile, we can obtain the pore pressure profiles in the sol and the gel from the relations *p*_{g}=−*Π*_{g}(*ϕ*)+*P*_{a} and *p*_{s}=−*Π*_{s}(*ϕ*)+*P*_{a}. These profiles are shown in figure 4*a*,*b* which show the non-dimensional pore pressure in the gel and the sol, respectively. Here *p*_{i} is the pressure at the sol–gel interface. It can be seen in figure 4*a* that outside of the boundary layer the pressure gradient in the gel asymptotes to a constant value. This can be seen from equation (2.5) by letting . Then3.16and the gradient of *p* asymptotes to a constant value. Allowing in equation (3.13), and by use of equations (2.7) and (2.8) we find that3.17outside the boundary layer in the gel. This asymptotic limit is shown as the dashed line in figure 4*a* and is in good agreement with the pore pressure for *z*>*d*. For typical parameter values from table 2 this implies that the pressure gradient is given by ∂*p*/∂*z*=1.04×10^{8} Pa m^{−1}.

Figure 4*b* shows the pore pressure profile in the sol phase. In the sol outside of the boundary layer the pressure asymptotes to a fixed value which corresponds to the pore pressure in the bulk sol . The dashed line in the figure shows this limit which is in good agreement with the pore pressure in the sol for *z*<−3*d*.

From comparing figure 4*a*,*b* it can be seen that the bulk of the pressure change in the system occurs in the gel phase in which the pressure can reach very large negative values. The gel remains saturated until the pore pressure reaches the capillary pressure *p*≈2*γ*/0.15*R* at which air invades the pores of the gel (e.g. Haines 1925). This capillary pressure represents the minimum pressure that is attainable in the gel, which non-dimensionally is3.18

Using the parameter values from table 2 we find that , and so it can be seen from figure 4*a* that for most colloids (with 1 *nm*≤*R*≤1 μm), at the point at which the pressure reaches its minimum the boundary layer thickness is small relative to the thickness of the crust and so the pressure gradient is approximately constant throughout the gel crust.

In the model we have assumed the evaporation rate is constant. This is justified by experiments on hard particles that show that the constant rate period (during which *E* is approximately constant) persists up to the point of air invasion of the pores (Wedin *et al.* 2004). However, if *E* does change, these results suggest that we can approximate the pressure distribution in the gel as a linear function provided the evaporation rate changes slowly.

Dufresne *et al.* (2006) conducted an experiment where a gel was allowed to form from evaporation of a colloidal silica sol. At a certain point evaporation was halted by covering the evaporating surface of the gel, and the sol–gel interface was examined. It was noted that after evaporation ceased particles in a region of the first several 100 μm of the compact region of the sol–gel system diffused back into the sol, leaving a sharp interface. This is consistent with the boundary layer in the sol found in our model and with the two phase idea in our model of a gel phase bonded together by adhesive forces and a stable sol. However we note that the addition of chemical bonding and concentration-based electrostatic double layer forces in colloidal silica will presumably yield more complicated expressions for the osmotic pressures than those considered here (e.g. Bergna 2006).

### (c) Flocculated gels

Although only densely coagulating gels have been considered so far, the equations presented here can also be extended to the consideration of loosely flocculating gels in the case of an unstressed gel. A commonly used model for the compression of a loosely flocculated gel is that of Buscall & White (1987). In this model, for a gel initially of solid fraction *ϕ*_{1}, if the stress *Π*=*P*_{a}−*p* on the solid part of a gel becomes greater than the compressive yield stress *P*_{y}(*ϕ*_{1}), where *P*_{y} is an increasing function of *ϕ*, the solid phase fraction increases irreversibly at a specified rate until the compressive yield stress is such that *Π*=*P*_{y}(*ϕ*) where *ϕ*>*ϕ*_{1}. In practice, this means that for slowly compressing gels the stress on the solid structure of the gel is always given by *Π*=*P*_{y}(*ϕ*) and this assumption has previously been shown to be appropriate and used by Brown *et al.* (2002) in their model of a compacting gel. Therefore, in order to consider the evolution of a flocculating gel, we need only replace the expression *Π*_{g}(*ϕ*) from the case of a densely coagulating gel with *Π*_{g}(*ϕ*)=*P*_{y}(*ϕ*) for an appropriate expression for the yield stress (e.g. Usher *et al.* 2006).

## 4. Stress in the gel phase

We now consider the case when a gel dries under lateral constraint so that the gel is under uniaxial strain. For an ideal poroelastic material, in order to analyse the stress on the gel during drying we would treat the gel as a linear poroelastic system (e.g. Brinker & Scherer 1990). However in reality gel parameters such as the bulk and shear moduli are highly nonlinear and so we must adapt the linear theory to accurately reflect the gel properties. Elastic stresses can build up in the gel; however, we assume that the sol remains unstressed so that the equations describing the evolution of the sol remain the same as described in the previous section. Equivalently we assume that the sol does not undergo a glass transition so its shear modulus is always zero. We also assume that the gel is not viscoelastic so it does not flow in response to applied stresses.

For many gels it is observed that the elastic strain response to an imposed stress depends only upon the effective stress which represents the stress on the solid phase of the gel. Here *σ*_{ij} is total stress, *δ*_{ij} is the Kronecker delta and *α* is the Biot–Willis constant and the solid is incompressible so *α*=1 and becomes the Terzaghi effective stress *σ*_{ij}+*pδ*_{ij} (Muir Wood 1990). In line with this observation we choose the simplest possible extension of the linear poroelastic, isotropic stress–strain relationship, making the assumptions that the strain response of the gel depends only upon the Terzaghi effective stress so that4.1where *G* and *ν* are the shear modulus and Poisson’s ratio of the gel, respectively, which are functions of the effective stress. *ϵ*_{ij} is the natural or Hencky strain, and a *δ* before a variable implies a small change in the variable. By summing over the indices of this constitutive relationship, we find that4.2where the bulk modulus4.3

For a gel constrained between rigid side walls (such as in the experiments of Dufresne *et al.* (2003), or for a gel of infinite lateral extent, it is reasonable to treat the system as undergoing uniaxial strain so that *δϵ*_{xx}=*δϵ*_{yy}=0 as contraction or expansion in the lateral directions is prevented. Also *σ*_{zz}=−*P*_{a} at the sol–gel and gel–air interfaces and we can assume that there are no changes in vertical stress in the gel and *δσ*_{zz}=0. Under these assumptions, and using the definition of the effective stress in equation (4.1), we find that4.4where4.5(cf. the discussion of uniaxial strain by Wang 2000).

Now from the definition of the solid phase fraction *ϕ*=*V*_{s}/*V* with *V*_{s} being the volume of solid particles in the network, we find that4.6and in the case of an incompressible solid phase, the first term on the right-hand side disappears so that using equation (4.2) we find that4.7Combining equations (4.7) and (4.4) we finally obtain that4.8and this equation can then be inserted into equation (2.7) in order to determine the diffusivity of the gel.

All that remains in order to obtain a complete set of equations is to determine the parameters *K* and *ν* which we assume to depend only upon the effective stress through its three invariants during elastic compression. In the absence of detailed knowledge of their effective stress dependence, we assume that bulk modulus depends only upon the first invariant of the effective stress (the mean effective stress) and that Poisson’s ratio is a constant. This form of bulk modulus can be justified analytically (Coussy 2004), and has been experimentally observed to be a good approximation in many materials (e.g. Buscall & White 1987; Muir Wood 1990; Coussy 2004). Although the assumption of a constant Poisson’s ratio with a dependent bulk modulus implies that the gel is non-conservative (Zytynski *et al.* 1978), we note that this approximation is often used, for instance in critical state soil mechanics (Muir Wood 1990), for its good agreement with experiment.

Letting , we see from equation (4.7) that *K* is only a function of *ϕ* so that4.9In order to calculate *K*(*ϕ*) we can consider the special case where the stress is given by *σ*_{ij}=−*Pδ*_{ij} so that the mean effective stress is equivalent to the negative of the gel ‘osmotic pressure’ −*Π* and *K*(*ϕ*)=*ϕ*∂*Π*/∂*ϕ* so that we can use experimental observations of *Π*(*ϕ*) in the gel to calculate the bulk modulus dependence on *ϕ*. For a dense gel of hard spheres, such as that observed by Buzzaccaro *et al.* (2007), the osmotic pressure in the gel is given by equations (3.4) and (3.2), and so using these equations, and the relationship between the mean effective stress and the osmotic pressure in the above expression for *K* we find4.10

Figure 5 shows the bulk modulus as a function of the solid phase fraction *ϕ*_{cp}. We can also calculate the shear modulus from the definition (4.3) and the Poisson’s ratio *ν*. For a Poisson’s ratio of 0.2, we see that *G*(*ϕ*)=3*K*(*ϕ*)/4. This shows that the shear rigidity takes a small, finite value when the particles first form a gel. However, as with the bulk rigidity, the shear rigidity of the gel will also increase substantially as the particles approach close packing.

### (a) Discussion

By the use of equation (4.4), we see immediately the impact that changing the boundary conditions during crust formation has upon the gel diffusivity. Using the equation for gel diffusivity (2.7) in the unstressed case4.11while in the stressed case4.12Interestingly, a simple alteration in the boundary conditions also changes the governing equations in the bulk. A typical Poisson’s ratio for a close-packed colloidal crystal is approximately 0.2 (Frenkel & Ladd 1987) so that the diffusivity in the stressed case *D*_{g}^{s}(*ϕ*) is approximately twice as large as the diffusivity in the unstressed case. In the case of steady-state crust growth this implies that the boundary layer thickness in the gel at the sol–gel interface will be twice the thickness in the stressed case compared with the unstressed case. Outside the boundary layer, the solid phase fraction will still be an asymptote to *ϕ*_{cp} and from equation (3.17) it can be seen that the pressure gradient will still be an asymptote to the same value that it takes in the unstressed case. This must be true so that the flux of the liquid to the evaporating surface is the same in both stressed and unstressed cases.

As the boundary layer thickness is small in both stressed and unstressed cases, in terms of pressure and solid phase fraction there will be little macroscopically to distinguish between the two cases. However, the key difference occurs in the presence of a tensile stress across the crust. From equation (4.4) and the definition of the effective stress, we find that *δσ*_{kk}=−4*ηδp*, and so4.13where we have used lateral symmetry, *σ*_{xx}=*σ*_{yy} and the fact that *σ*_{zz}=−*P*_{a} to obtain this relation between lateral tensile stress and pressure. The dependence of these relations on Poisson’s ratio is shown in figure 6 for physically allowable positive values of Poisson’s ratio. Importantly, ∂*σ*_{xx}/∂*p* is always negative while is always positive. During drying the pore pressure becomes large and negative so that the lateral stress *σ*_{xx} will be tensile, while the lateral effective stress is compressive. Therefore, during drying because the effective stress is compressive even for weakly adhesive gels, the gel particles will not pull apart even though large tensile stresses can be formed across the gel. The lateral tensile stresses then cause the growth of cracks in the gel (Atkinson & Craster 1991; Scherer 1992).

If we assume a constant Poisson’s ratio we can integrate equation (4.13) using continuity of pressure and stress at the sol–gel interface to find that the lateral stress in the gel is given by *σ*_{xx}=−*P*_{a}−2*η*(*p*−*p*_{i}). Because the pressure gradient is a constant in the bulk of the gel the tensile stress is approximately a linear function of the distance from the sol–gel interface; from use of equation (3.17) we find that *σ*_{xx}≈−*P*_{a}+2*ηEμz*/*k*_{g}(*ϕ*_{cp}) in the gel. Also noting from the calculation of the velocity of the gelation front that the position of the evaporating interface is given approximately by *z*=*βEt*, we find that the stress at the evaporating interface is well approximated by4.14which increases linearly with time as has been seen in experimental observations (e.g. Martinez & Lewis 2002; Wedin *et al.* 2004). Typical values of parameters from table 2 give a surface stress of *σ*_{xx}=−*P*_{a}+8.6×10^{3}*t*.

The greatest possible tensile stress occurs when the pore pressure reaches its minimum value *p*=*P*_{a}−2*γ*/0.15*R* where4.15The approximation follows because the capillary pressure term is much larger than the interface pressure and atmospheric pressure. This is an important result because it is commonly assumed that the tensile stress in a gel that drives cracking is the same as the maximum capillary pressure in the liquid 2*γ*/0.15*R* (e.g. Dufresne *et al.* 2003; Lee & Routh 2004). However, equation (4.15) shows that when the gel crust is constrained during drying, the maximum tensile stress is altered from the maximum capillary pressure by a factor of 2*η*. For a Poisson’s ratio of 0.2, this means that the maximum tensile stress on the gel is only three quarters of the maximum capillary pressure. This difference will have implications in the context of interpretation of experimental measurements of tensile stress arising during drying.

## 5. Conclusions

In this paper, we studied the growth of a gel crust at the surface of a drying sol. We derive equations and boundary conditions which couple poroelastic effects in the gel with colloidal effects in the sol. In the sol and the gel phases, the solid phase fraction *ϕ* satisfies similar advection–diffusion equations with nonlinear diffusivitiesThe dependence of *D*_{g} on the gradient of pore pressure with *ϕ* gives rise to a dependence of the gel diffusivity on the constraints placed on the system during drying. For example, we have shown that in the case of uniaxial strain the gel diffusivity approximately doubles over its value in the unstressed case.

During steady-state growth of the gel due to evaporation, *ϕ* has been shown to be constant outside of thin boundary layers on either side of the gel–sol interface, while the pore pressure at the gel–air interface scales with the thickness of the crust. This means that macroscopically the system can be thought of as two adjacent bulk phases of constant solid phase fractions *ϕ*_{0} and *ϕ*_{cp}, with the pore pressure in the gel increasing linearly from the sol–gel interface to the evaporating surface.

By considering the non-hydrostatic stress state in the gel we find that the stress at the surface of the gel increases linearly with time in agreement with experimental observations. We also find that the maximum tensile stress responsible for cracking is not equal to the maximum capillary pressure attainable at the surface before air invasion of the pores occurs, as has been previously assumed. Instead, integration of equation (4.13) in the case of constant Poisson’s ratio shows that the maximum tensile stress at the surface of the gel is given bywhere *P*_{cap}^{max} is the maximum capillary pressure. Thus, for a realistic Poisson’s ratio of 0.2, the maximum tensile stress is 0.75*P*_{cap}^{max}. This result is important in interpretation of the results of stress measurements in drying films and in the analysis of fracture toughness of drying gels.

This work is applicable to the modelling of any system involving crust formation in a drying colloidal suspension. In addition, the equations for consolidation in a drying gel can also be used separately to analyse the evolution of a drying gel. Such systems are commonly found and so this model will be of use in understanding the behaviour of gel crusts and stress evolution in the many fields of colloid science, including soil desiccation, xerogel formation and crack formation in drying films such as ceramics and paints. Potential extensions of this work include more complex interparticle potentials, the effect of additional stress regimes, non-steady interface growth and modelling of unsaturated crusts.

## Acknowledgements

This publication is based on work supported by Award No. KUK-C1-013-04, made by King Abdullah University of Science and Technology (KAUST).

## Footnotes

- Received January 26, 2010.
- Accepted June 1, 2010.

- © 2010 The Royal Society