Spreading of viscoplastic droplets

Abstract

The spreading under surface tension and gravity of a droplet of yield-stress fluid over a thin film of the same material is studied. The droplet converges to a final equilibrium shape once the driving stresses inside the droplet fall below the yield stress. Scaling laws are presented for the final radius and complemented with an asymptotic analysis for shallow droplets. Moreover, numerical simulations using the volume-of-fluid method and a regularized constitutive law, and experiments with an aqueous solution of Carbopol, are presented.

1. Introduction

The impact and spreading of droplets of complex fluids over surfaces occur in a wide variety of industrial applications. Examples include, but are not limited to, inkjet printing, spray coating and additive manufacturing (Barnes 1999; Derby 2010; Thompson et al. 2014; Mackay 2018). In many cases, such as three-dimensional (3-D) printers and paint sprays, liquid droplets or filaments are deposited on an existing layer of the same fluid. Hence, understanding the underlying fluid mechanics of droplet spreading on a thin film helps to improve the design of such systems. On the theoretical side, the removal of an advancing contact line has the extra advantage of simplifying the spreading problem substantially by removing the complicated physics associated with relieving the stress singularity that otherwise arises (e.g. Oron, Davis & Bankoff 1997; Bonn et al. 2009; Craster & Matar 2009). Precursor films are also expected to be drawn out ahead of spreading droplets by intermolecular forces, effectively emplacing a pre-wetted film even in situations in which one did not exist originally (e.g. Bonn et al. 2009; Craster & Matar 2009).

The spreading of Newtonian droplets over a thin layer has been addressed previously for a range of physical regimes (see Bergemann, Juel & Heil (2018), Jalaal, Seyfert & Snoeijer (2019b), and references therein). The present article aims to provide a discussion on the viscoplastic version of this problem. Viscoplastic or yield-stress fluids feature a mix of fluid and solid behaviour: if not sufficiently stressed, such materials behave more like an elastic material, but, above a critical yield stress, they flow like a viscous fluid. Yield stress is a common feature of many natural and industrial fluids, such as clay, cement, toothpaste, cosmetic creams, dairy products and waxy oil (see Balmforth, Frigaard & Ovarlez (2014), Coussot (2014) and Bonn et al. (2017) for reviews).

For a Newtonian droplet deposited on a thin film of the same fluid and spreading due to gravity and surface tension, flow continues until a completely flat film is formed. By contrast, when the driving stresses fall below the yield stress, flow ceases inside a viscoplastic droplet. Consequently, the final shape is not flat, as shown previously in different configurations (Roussel & Coussot 2005; Balmforth et al. 2007a; German & Bertola 2009; Jalaal, Balmforth & Stoeber 2015; Liu et al. 2016; Chen & Bertola 2017; Liu, Balmforth & Hormozi 2018). The final shape of the droplets is of particular importance in industrial processes such as 3-D printing, as it can determine the resolution of the printing process and the final quality of the product. Previous experiments have been reported for the problem of viscoplastic droplet impact on solid surfaces (Luu & Forterre 2009; Saïdi, Martin & Magnin 2010; Luu & Forterre 2013; Blackwell et al. 2015; Sen, Morales & Ewoldt 2020) or fluid interfaces (Jalaal, Kemper & Lohse 2019a).

In the present work, we explore the spreading of a viscoplastic droplet and its final shape, providing a theoretical framework for the problem complemented with experimental and computational results. The order of material in this paper is as follows. Section 2 presents simple scaling laws for the final shape of a viscoplastic droplet. Section 3 summarizes a viscoplastic lubrication theory suitable for shallow droplets. Section 4 presents the numerical simulations for a spreading viscoplastic droplet. Section 5 describes the experimental tests, and is followed by § 6, which summarizes the theoretical and experimental results. The appendices contain further technical details of the lubrication theory and numerical computations.

2. Scaling laws for the final shape

Consider a viscoplastic droplet deposited on a thin film of thickness $H_{\infty }$ at $t=0$. We imagine that the droplet yields entirely under capillary action or gravity, spreading and then braking to a halt due to the yield stress $\tau _0$. That is, we consider the situation in which the yield stress cannot localize flow and leave intact a substantial volume of the droplet to imprint a dependence of the final shape on the initial configuration. Global force balance over the entire droplet volume should then control the final radius $\mathcal {R}_f$ and height $\mathcal {H}_f$. If the rheology of the fluid only features through the yield stress, the physical parameters of the problem include $\tau _0$, the density of the droplet $\rho$, the surface tension coefficient $\sigma$, gravitational acceleration $g$ and the droplet volume $\mathcal {V}$.

When the droplet spreads under capillary effects, the driving horizontal pressure force (given by the product of the pressure $p\sim \sigma \kappa$ and a typical vertical surface area $\mathcal {H}_f\mathcal {R}_f$) can be estimated as

(2.1)\begin{equation} p\mathcal{H}_f\mathcal{R}_f \sim \sigma \kappa \mathcal{H}_f\mathcal{R}_f \sim \sigma \frac{\mathcal{H}_f^{2}}{\mathcal{R}_f}, \end{equation}

where $\kappa \sim \mathcal {H}_f / \mathcal {R}_f^{2}$ is the curvature. On the other hand, if the droplet does not slip over the underlying surface and the yield stress acting over the base of the droplet provides the main resistance to flow, the opposing force is of the order of $\tau _0 \mathcal {R}^{2}_f$. By balancing the two, we arrive at

(2.2)\begin{equation} \sigma \mathcal{H}_f^{2} \sim \tau_0 \mathcal{R}^{3}_f. \end{equation}

Moreover, since $\mathcal {V}\sim \mathcal {H}_f\mathcal {R}^{2}_f$, if we define the length scale $\mathcal {L}=[3\mathcal {V}/(4{\rm \pi} )]^{1/3}$ (i.e. the radius of the corresponding spherical drop), then

(2.3)\begin{equation} \sigma \mathcal{V}^{2} \sim \tau_0 \mathcal{R}^{7}_f \quad \textrm{or} \quad \frac{\mathcal{R}_f }{\mathcal{L}} = \varOmega_c \mathcal{J}^{-1/7}, \quad \mathrm{where}\ \mathcal{J}=\frac{\tau_0 \mathcal{L}}{\sigma}\end{equation}

is a non-dimensional number that compares the yield stress and capillary pressure. Thus, as the strength of the plastic effect $\mathcal {J}$ increases, the final radius becomes correspondingly smaller. In (2.3), the prefactor $\varOmega _c$ encapsulates dependence on the remaining dimensionless groups in the problem. If gravity and any other effects are not important, the only remaining group is the scaled pre-wetted film thickness, $h_\infty \equiv H_\infty /\mathcal {L}$.

If we instead counter a driving hydrostatic pressure $p\sim \rho g \mathcal {H}_f$ by the resistance from the yield stress, the balance is

(2.4)\begin{align} \rho g \mathcal{H}_f^{2} \mathcal{R}_f \sim \tau_0 \mathcal{R}^{2}_f \quad \textrm{or} \quad \frac{\mathcal{R}_f }{ \mathcal{L}} = \varOmega_g \left( \frac{\rho g \mathcal{L}}{\tau_0} \right)^{1/5} \equiv \varOmega_g \mathcal{B}^{1/5} \mathcal{J}^{-1/5}, \quad \mathrm{where}\ \mathcal{B} = \frac{\rho g \mathcal{L}^{2}}{\sigma} \end{align}

is the Bond number, comparing the hydrostatic pressure and capillary pressure. Again, the prefactor $\varOmega _g$ contains the dependence on any other dimensionless groups. Note that the combination $\mathcal {B} / \mathcal {J}$ is independent of $\sigma$, eliminating surface tension from the right-hand side of (2.4) when capillary effects are not present and $\varOmega _g$ depends only on $h_\infty$. This limit is relevant to geophysical flows and rheometry with larger spatial scales (cf. Roussel & Coussot 2005; Balmforth et al. 2006).

More generally, we must take either $\varOmega _c$ or $\varOmega _g$ to depend on both $\mathcal {B}$ and $h_\infty =H_\infty /\mathcal {L}$, as well as any other dimensionless parameters stemming from further physical effects. In what follows, we assume that only $\mathcal {B}$ and $h_\infty$ are relevant and write the general relation

(2.5)\begin{equation} \frac{\mathcal{R}_f }{\mathcal{L}} = \varOmega (\mathcal{B},h_\infty) \mathcal{J}^{-1/7}, \end{equation}

with $\varOmega \to \varOmega _c$ in the capillary-dominated limit $\mathcal {B}\to 0$, and $\varOmega \to \varOmega _g \mathcal {B}^{1/5}\mathcal {J}^{-2/35}$ in the gravity-dominated limit $\mathcal {B}\gg 1$. In particular, we compare this scaling law with asymptotic analysis, numerical simulations and experiments, each of which provides more refined estimates for $\varOmega$.

3. Viscoplastic lubrication theory

As sketched in figure 1, and assuming that the droplet remains axisymmetric, we employ cylindrical polar coordinates $(r,z)$ to describe the geometry of a shallow viscoplastic droplet. The top surface of the fluid lies at $z=h(r,t)$, and the droplet is emplaced upon an existing fluid layer of thickness $H_\infty$. There is a rigid plane at $z=0$ over which the fluid cannot slip. To simplify our discussion, we used the Bingham constitutive law, which combines the yield stress $\tau _0$ with a plastic viscosity $\mu$.

Figure 1. Sketches of the geometry and anatomy of a spreading viscoplastic droplet. The top surface is $z=h(r,t)$; below $z=Y(r,t)$ the fluid is fully yielded and flows in a plug-like manner over $Y < z < h$, where fluid is held near the yield stress (the ‘pseudo-plug’).

Lubrication theory applies when droplets are relatively shallow and inertia is negligible. In this instance, the hydrostatic and capillary pressures are largely independent of $z$ and drive spreading, which is primarily countered by the vertical shear stress. The analysis for viscoplastic fluid follows along similar lines to that for Newtonian fluid (e.g. Oron et al. 1997; Craster & Matar 2009), the key differences arising from the impact of the yield stress on the vertical profile of the radial velocity (Liu & Mei 1989; Balmforth et al. 2006). In particular, as illustrated in figure 1, the velocity field adopts a distinctive anatomy in which a lower slice of the fluid in $0 < z < Y(r,t)$ is fully yielded and the radial velocity has a parabolic profile; over the region $Y(r,t) < z < h(r,t)$, the radial velocity becomes plug-like and independent of $z$ to leading order. The upper region possesses a shear stress that lies below the yield stress, an observation that has previously led to the incorrect conclusion that the fluid here is unyielded, contradicting the radial expansion and appearing to imply a paradox in lubrication theory. In fact, over the plug-like region, or ‘pseudo-plug’, the extensional stresses are of the same order as the shear stress (a feature demanded by the proper 3-D form of the Bingham constitutive law when deformation rates are relatively small) and their conspiracy holds the stress slightly above $\tau _0$ to permit the radial expansion (Balmforth & Craster 1999; Putz, Frigaard & Martinez 2009). In the limit $\tau _0\to 0$, the pseudo-plug disappears ($Y\to h$) and we recover the fully parabolic Newtonian flow profile. Conversely, when $Y\to 0$, the pseudo-plug reaches the base, bringing the entire fluid layer to a halt.

Given the shallow-layer velocity field of the lubrication analysis, the expression of depth-integrated mass conservation provides an evolution equation for the local fluid depth. We express this equation in the dimensionless form

(3.1a,b)\begin{equation} h_t = \frac{1}{6 r} [r p_r Y^{2} ( 3h-Y ) ]_r, \quad p = {\mathcal{B}} h - \frac{1}{r} (r h_r)_r, \end{equation}

after scaling lengths by $\mathcal {L}$, pressure $p(r,t)$ by $\sigma / \mathcal {L}$, velocity by $U = \sigma /\mu$ and time by $\mathcal {L}/U$; the surface bordering the pseudo-plug is given by

(3.2)\begin{equation} Y = \max\left( 0, h- \frac{ {\mathcal{J}}}{ \vert p_{r} \vert} \right).\end{equation}

For ${\mathcal {J}}=0$, one has $Y=h$ and (3.1a,b) reduces to the standard evolution equation for a Newtonian film and therefore recovers known spreading laws for viscous droplets under the action of gravity, capillarity or a combination of both (Oron et al. 1997; Bonn et al. 2009; Craster & Matar 2009).

3.1. Sample spreading solutions

To provide sample solutions for spreading viscoplastic droplets, we solve (3.1a,b) numerically using centred finite differences to approximate radial derivatives and a stiff integrator to step the surface profile forwards in time. The length of the domain is chosen sufficiently large that the droplet never reaches the border. We begin from the initial condition $h(r,0) = \max (0, 1- 3r^{2} /16 )^{3} + h_{\infty }$, which smoothly interpolates between the pre-wetted film of scaled thickness $h_\infty$ and an initial ‘bump’ with a dimensionless radius $R(0)=4/\sqrt {3}$ (chosen so that the dimensional volume is ${\mathcal {V}}=4{\rm \pi} {\mathcal {L}}^{3}/3$, given the scaling of lengths by $\mathcal {L}$); the results are insensitive to the precise initial shape of the bump provided the droplet becomes fully yielded during spreading. We also replace (3.2) by the regularization $Y=\max (\varepsilon ,h- {\mathcal {J}}/|p_r|)$ to ease computations, with $\varepsilon =10^{-6}$ (having verified that the precise value makes no significant difference to the results).

Figure 2(a) shows a solution with ${\mathcal {J}}=5/4$ and ${\mathcal {B}}=0$, using an initial fluid film of thickness $h_\infty =0.01$. The dimensionless yield stress is sufficiently small that the entire droplet yields under capillary action at $t=0$ and then spreads much like a Newtonian droplet. Subsequently, however, the yield stress comes into play as driving stresses decline, and eventually the droplet brakes to rest. Distinctive spatial oscillations appear near the edge of the droplet, becoming frozen into the final shape as flow ceases. These undulations also appear in the Newtonian problem and have been reported previously for viscoplastic films (Balmforth, Ghadge & Myers 2007b; Jalaal & Balmforth 2016). We discuss them further below and in greater detail in appendix A.

Figure 2. Numerical solutions of the evolution equation (3.1a,b). (a) Snapshots of $h(r,t)$ (grey lines) for ${\mathcal {J}}=5/4$ and ${\mathcal {B}}=0$, with the inset showing a magnification of the edge of the droplet; the final snapshot is plotted in black. The surface $z=Y(r,t)$ at an intermediate time (corresponding to the curve of $h(r,t)$ in blue) is included as the red dashed line. (b,c) Time series of the edge (defined as the first radial position $R(t)$ where $h(R,t)=h_\infty$) for (b) $( {\mathcal {J}},h_\infty )=(1.25,0.01)$ with ${\mathcal {B}}=0$, 3/4 and 2, and (c) $( {\mathcal {J}}, {\mathcal {B}})=(1.25,0)$ with $h_\infty =10^{-3}$, $3\times 10^{-3}$ and $0.01$. The red dotted lines show corresponding Newtonian solutions ($\mathcal {J}=0$), and the stars indicate the times of the snapshots in panel (a). The relatively sharp features in $R(t)$ for $t\approx 1$ in panels (b) and (c) correspond to the creation of the undulating wavetrain at the edge of the droplet.

The solutions are much the same, though wider and flatter, with gravity (${\mathcal {B}}>0$). Figures 2(b) and 2(c) show time series of the edge $R(t)$ for further solutions with different parameter settings. Here, the edge is measured for numerical convenience as the first location $R(t)$ where $h(R,t)=h_\infty$. With ${\mathcal {J}}=0$, the edge continues to expand; the yield stress, however, inevitably brings fluids to rest.

3.2. Final shapes

A viscoplastic droplet comes to rest when $Y \rightarrow 0$, implying $h |p_r| = {\mathcal {J}}$, which comprises a third-order differential equation for the final profile in view of (3.1a,b). A complication in solving this equation is the presence of the factor $|p_r|$, which leaves open the sign of the pressure gradient. To determine this sign, we appeal to the evolution equation (3.1a,b) and the limit of its solution as the fluid comes to rest. In particular, over the bulk of the droplet, the sign of $-p_r$ corresponds to the sense of the flux, which must be positive. Near the edge, however, the spatial oscillations complicate matters. There, as is clear from the magnification in figure 2(a), the undulations correspond to a travelling wavetrain such that

(3.3)\begin{equation} h_t \to - R_t h_r \sim \tfrac{1}{6} [ p_r Y^{2} (3h-Y)]_r . \end{equation}

Integrating this equation and observing that $h\to h_\infty$ for $p_r Y^{2} (3h-Y)\to 0$, we find that the sign of $-p_r$ must be given by the sign of $h-h_\infty$. Thus,

(3.4)\begin{equation} h_{rrr} + \frac{1}{r} h_{rr} - \frac{1}{r^{2}} h_r - {\mathcal{B}} h_r = \frac{ {\mathcal{J}}}{h}\textrm{sgn}(h-h_\infty). \end{equation}

3.2.1. Gravity-dominated limit

In the gravity-dominated limit, the higher derivatives disappear from the left-hand side of (3.4), leaving $- hh_r \sim {\mathcal {J}}/ {\mathcal {B}}$ (since $h\geq h_\infty$). Thence

(3.5)\begin{equation} h = \left[ h_\infty^{2} + \frac{2 {\mathcal{J}}}{ {\mathcal{B}}} (R-r) \right]^{1/2}\end{equation}

(cf. Blake 1990; Roussel & Coussot 2005; Balmforth et al. 2006). The corresponding droplet volume must be $\mathcal {V}=2 {\rm \pi}\mathcal {L}^{3} \int _0^{R} [h(r)-h_\infty ] r\,\textrm {d}r$, which demands that the final radius satisfy the algebraic problem

(3.6a,b)\begin{align} R\equiv\frac{\mathcal{R}_f}{\mathcal{L}} = \left( \frac{25 {\mathcal{B}}}{8 {\mathcal{J}}}\right)^{1/5} \left[ (1+A)^{5/2} - \frac{1}{2} A^{3/2} (2A+5)- {\frac{15}{8}}A^{1/2}\right]^{-2/5},\quad A = \frac{ {\mathcal{B}} h_\infty^{2}}{2 {\mathcal{J}} R} . \end{align}

The solution can be written formally as ${\mathcal {R}}_f/\mathcal {L} = \varOmega _g(h_\infty \sqrt { {\mathcal {B}}/ {\mathcal {J}}}\,) {\mathcal {B}}^{1/5} {\mathcal {J}}^{-1/5}$ in the manner of (2.4). Notably, when $h_\infty \to 0$, one obtains $\varOmega _g\to (25/8)^{1/5}\approx 1.26$.

3.2.2. Finite ${\mathcal {B}}$

Away from the gravity-dominated limit, we must attack (3.4) numerically. For a pre-wetted film with finite thickness, the boundary conditions are the symmetry condition $h_r(0)=0$ and the far-field condition, $h\to h_\infty$ for large radii. The latter requires the imposition of two conditions in order that the droplet profile meets the pre-wetted film continuously. Given the form of the solution at the edge, we choose $h(R_*)=h_\infty$ and $h_{rr}(R_*)=0$, where $R_*$ denotes a radius well into the decaying undulations. Applying this boundary condition corresponds to pinning the solution at a point further along the wavetrain. In addition, we must also arrive at the correct droplet volume. Thus, the third-order equation (3.4) must be solved subject to three boundary conditions and the volume constraint, demanding that the edge position $R_*$ be found as part of the solution (i.e. an eigenvalue). Appendix A describes further details of the numerical construction of the final profile, as well as a more detailed consideration of the undulations at the edge.

Figure 3 shows a sample numerical solution with ${\mathcal {B}}=0$. This particular example corresponds to the solution of the evolution equation in figure 2(a), and is compared with the final snapshot of that computation in figure 3. The decaying undulations converge to a sawtooth wave in $h_{rr}$, with the corners corresponding to the sign switches of $h-h_\infty$. Unlike the decaying capillary waves of moving Newtonian contact lines, which have fixed wavelength (Tanner 1979; Tuck & Schwartz 1990; Jalaal et al. 2019b), the viscoplastic undulations shorten with distance along the wavetrain. In appendix A, we outline how the waveform converges to piecewise-cubic polynomials, with an accumulation point at a finite outer radius.

Figure 3. (a) A final equilibrium profile for $( {\mathcal {J}}, {\mathcal {B}},h_\infty )=(1.25,0,0.01)$. (b) A magnification near the edge is shown, along with $h_{rr}$. The dashed lines show the final snapshot of the numerical solution of the evolution equation from figure 2(a).

The analysis of the edge behaviour in appendix A also indicates that the wavetrain shrinks to a point in the limit $h_\infty \to 0$. Moreover, the limiting solution corresponds to solving (3.4) with outer boundary conditions based on the local solution $h \sim \mathcal {C}(R-r)^{3/2}$ for $r\to R$ and an unknown constant $\mathcal {C}$. Figure 4 illustrates such limiting profiles for various values of the gravity parameter. Evidently, since $h_r(R)=0$, the limiting contact angle is zero here, implying spreading over a perfectly wetting surface (although one can also solve (3.4) with a prescribed contact angle).

Figure 4. (a) Equilibrium profiles for $h_\infty \to 0$ with $R^{2} {\mathcal {B}}=0$, 10, 30, 100 and 300, together with the limit ${\mathcal {B}}\to \infty$ given by (3.5). The dots show the approximation $(1-r^{2})^{3/2}$. (b) The coefficient $\varOmega (x= {\mathcal {B}} {\mathcal {J}}^{-2/7},0)$ in (3.7), along with the limits for $x=0$ and $x\gg 1$.

From the computed solutions for final equilibrium profiles, we may evaluate the final radius and depth, which are given by the (fairly complicated) algebraic problem outlined in appendix A. For the final radius, the formal solution may be written as

(3.7)\begin{equation} \frac{\mathcal{R}_f}{\mathcal{L}} = \varOmega\left(\frac{ {\mathcal{B}}}{ {\mathcal{J}}^{2/7}}, \frac{h_\infty}{ {\mathcal{B}}}\right) {\mathcal{J}}^{-1/7}, \end{equation}

which identifies the dependence of the coefficient on gravity and the prewetted film thickness. The prefactor $\varOmega ({ {\mathcal {B}}}{ {\mathcal {J}}^{-2/7}},{h_\infty }/{ {\mathcal {B}}})$ is shown in figure 5. In the limit $h_\infty \to 0$, the function $\varOmega (x,0)$ is shown in more detail in figure 4(b); one can see that $\varOmega _c\equiv \varOmega (0,0)\approx 1.74$ and $\varOmega (x,0)\to (25x/8)^{1/5}$ for $x\gg 1$, which aligns with the result in § 3.2.1.

Figure 5. The function $\varOmega ({ {\mathcal {B}}}{ {\mathcal {J}}^{-2/7}},{h_\infty }/{ {\mathcal {B}}})$ in (3.7) plotted as a surface over the $({ {\mathcal {B}}}{ {\mathcal {J}}^{-2/7}},{h_\infty }/{ {\mathcal {B}}})$ plane. The dashed line shows the result from figure 4(b).

4. Numerical simulations

To complement the lubrication analysis, we solve the spreading problem numerically away from the shallow limit using the open-source code Gerris (Popinet 2003). The code employs a volume-of-fluid scheme to deal with the interface and an adaptive grid to achieve high resolution inside the droplet and along its interface (see appendix B for more details). For the rheology, we use the regularized Bingham model with

(4.1a,b)\begin{equation} \tau_{ij} = \mu_1 \dot\gamma_{ij},\quad \mu_1 = \min \left(\frac{\tau_0}{\dot\gamma} + \mu,\mu_{max} \right) \dot\gamma_{ij}, \end{equation}

where

(4.2a,b)\begin{equation} \dot\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i},\quad \dot\gamma\equiv \left( \frac{1}{2}\sum_{i,j}\dot\gamma_{ij}\dot\gamma_{ji} \right)^{1/2} \end{equation}

(cf. O'Donovan & Tanner 1984). Here, the divergence of the viscosity at low shear rates is controlled by the regularization parameter, $\mu _{max}$, chosen sufficiently large to render the simulations insensitive to the precise value (see appendix B). We use a parabolic initial shape for the droplet, merged to a flat pre-wetted film, motivated by the experimental profiles observed in § 5:

(4.3)\begin{equation} h(r,0) = h_\infty +(2/3)^{1/3} \max (0,1-(r/R_0)^{2}), \quad \mathrm{with}\ R_0=(2/3)^{1/3}. \end{equation}

The simulation includes inertia, allowing us to take the velocity field to be zero at $t=0$. However, for the physical parameters studied here, the effect of inertia is always small.

Figure 6 shows an example for $\mathcal {J}=0.144$, ${\mathcal {B}}=0$ and $h_\infty =0.0175$, plotting snapshots of the interface with superposed density maps of $\dot {\gamma }$. Initially, the surface tension arising due to the high curvature at the edge of the droplet drives flow, flattening the entire profile over later times. Throughout, a plug remains at the core of the droplet (marked as zone I in figure 6b), much like in the gravity-driven problems explored by Liu et al. (2016, 2018). Once the droplet becomes shallower, a further region of low strain rates forms close to the interface (denoted as zone II in figure 6c), resembling the pseudo-plug region of the lubrication analysis. As the droplet brakes to rest, the regions of small strain rate broaden to span the droplet, with a train of undulations appearing at the edge as in the lubrication analysis (see appendix B and figure 12). We calculate final shapes for a range of $\mathcal {J}$ and $\mathcal {B}$, and later, in § 6, we compare them with those obtained from the lubrication analysis and the experiments.

Figure 6. Distribution of $\log _{10}(\mathcal {\dot {\gamma }})$ in a spreading viscoplastic droplet for $\mathcal {J}=0.144$, $\mathcal {B}=0$ and $h_\infty =0.0175$, at the times $t=0.29$, 87, 290, 580, 1450, 2030 (a–f, respectively).

5. Experiments

We experimentally study the spreading of viscoplastic droplets by extruding them from a syringe onto a pre-wetted surface. The experimental apparatus consisted of a hydrophobic nozzle (with inner diameter of 0.3 mm) connected to a syringe pump (KD Scientific-Legato 111), where the vertical location of the nozzle could be adjusted using a translation stage. Chemically treated glass slides were used to suppress any slip over the underlying surface (see Jalaal et al. (2015) for details). To form the pre-wetted film, spacers of a given height (adhesive tape of ${\sim }60\ \mathrm {\mu }$m thickness) were placed on either side of the surface, and a mound of the fluid was spread out evenly with a flat blade. The experimental set-up is sketched in figure 7(a).

Figure 7. (a) Schematic picture of the set-up for the extrusion tests, where the nozzle tip is placed just above the substrate. (b) Flow curves of the fluids used in the tests on a log–log scale. The black curves show the Herschel–Bulkley fits.

To minimize inertial effects, we slowly extruded the droplets on the surfaces. The nozzle tip was placed $200\ \mathrm {\mu }$m above the film, and droplets of different volumes ($\mathcal {V}=0.0042$ to 1.4367 ml) were deposited. In all experiments, the extrusion flow rate was 2 ml min$^{-1}$. The shapes of the droplets during spreading were recorded using side imaging, a cold light-emitting diode illuminating the test section, and images were recorded with a high-speed camera attached to a microscope. The shape and volume were obtained through image processing.

The working fluids were aqueous suspensions of Carbopol Ultrez 21 (by Lubrizol), neutralized with triethanolamine. The rheology of the fluids was characterized using controlled shear-rate tests with an Anton Paar (Physica MCR-302) rheometer fitted with sand-blasted (PP25-S) parallel plates (roughness of ${\sim }4\ \mathrm {\mu }$m). Figure 7(b) shows the flow curves of the five concentrations of Carbopol used, plus Herschel–Bulkley fits of the form $\tau = \tau _0 + K\dot {\gamma }^{n}$ to the measured shear stress $\tau$ and shear rate $\dot \gamma$. The fitting parameters (yield stress, consistency index $K$ and flow index $n$) are provided in table 1, supplemented by estimates of the elastic modulus $G$ found from constant, low-shear-rate tests (monitoring the stress growth with strain).

Table 1. Properties of the experimental fluids.

The density of the solutions is very close to that of water. Measuring the surface tension of yield-stress fluids is challenging. For Carbopol solutions, reported values vary over a range of ${\sim }51\text {--}69$ mN m$^{-1}$ (Manglik, Wasekar & Zhang 2001; Boujlel & Coussot 2013; Géraud et al. 2014; Jørgensen et al. 2015). Here, we did not measure the surface tension of our samples and used the fiducial value of 63 mN m$^{-1}$ as the rough average of all the reported values.

Figure 8(a) displays an example extrusion test. A small droplet forms at the beginning of the injection and grows in time as the pumping continues. When the pump is switched off, the droplet keeps spreading freely until reaching the final state. Figures 8(b) and 8(c) show results from tests for a given $\mathcal {B}$ at different $\mathcal {J}$ (same volume with different yield stress). As expected, an increase in the magnitude of $\mathcal {J}$ results in a smaller final radius of the droplet. Note that the photographs of the final shapes were taken 30 s after the extrusion commenced, to ensure the equilibrium shape was reached. Although some success has been enjoyed with Newtonian fluid (Jalaal et al. 2019b), we were also unable to clearly visualize any finer structure at the droplet edge such as the undulations that emerge in the lubrication theory and numerical simulations.

Figure 8. (a) An example of an extrusion test for $\mathcal {B}=1.286$ and $\mathcal {J}=0.463$. Droplets start to grow with the pumping and reach a final shape quickly after the extrusion is over. Snapshots are at $t=0$, $t=1.5$, $t=3$ and $t=4.5$ s, respectively (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2020.886). The surface is pre-wetted with $h_{\infty } \sim 300\ \mathrm {\mu }$m. (b) The variations of the radius of droplets over time. Pumping is terminated at $\sim$3.5 s. (c) The final state of the droplets, corresponding to coloured dots in panel (b).

Experiments were conducted for different $\mathcal {J}$ (by changing the yield stress; cf. figure 7b) and $\mathcal {B}$ (by changing the droplet size), achieving 18 different parameter settings in total, as restricted by the limitations outlined below (cf. inset in figure 9). The experiments were repeated (at least five times) for each data point, and the average values were calculated. The standard deviations (resulting from the accuracy of the image processing and the repeatability of the tests) are small (${\sim }$5–7 %) for the majority of the extrusions. The standard deviation is, however, larger, $\sim$15 % at the lowest values of $\mathcal {B}$ since the drops are smaller and controlled deposition is harder. Note that the dimensional groups have limited ranges: to attain small values of $\mathcal {B}$, we had to extrude a small volume, promoting the effect of the nozzle. For larger values of $\mathcal {B}$, with our lowest yield stress, the droplets acquired a very large footprint that exceeded the boundary of the biggest chemically treated substrate available to us.

Figure 9. Comparison of the experimental and theoretical results of the final radius for different $\mathcal {J}$ and $\mathcal {B}$. Grey lines are asymptotic for $h_{\infty } \rightarrow 0$. The inset shows the dimensional final radius for different yield stress and volume. The asymptotic and simulation results for $\mathcal {B}=0$ and $\mathcal {B}=0.14$ were so close that we just show the former for clarity.

6. Discussion and conclusion

Figure 9 summarizes our results for the dimensionless final radius $\mathcal {R}_f/\mathcal {L}$ obtained from the lubrication theory, numerical simulations and experiments. In the simulations, the value of the pre-wetted film is set at $h_{\infty }\approx 0.0175$. In experiments, we measure this value to be $h_{\infty } \approx 0.07 \pm 0.03$. The results of the lubrication theory are shown for $h_\infty \to 0$. The three sets of results show broad agreement with the predictions of scaling theory, and, in particular, the trends ${\mathcal {J}}^{-1/7}$ and ${\mathcal {J}}^{-1/5}$ predicted in the capillary- and gravity-dominated limits.

Although the theoretical and experimental results mostly overlap in figure 9, there are discrepancies. First, for the range of $\mathcal {B}$ explored here, from figure 5, we see that the prefactor increases by approximately 10 % if one increases the pre-wetted film thickness from the limit $h_\infty \to 0$ up to $h_\infty =0.0175$. Consequently, the results from the lubrication theory should be approximately 10 % higher in figure 9 to match up properly against the simulations. Evidently, the asymptotics overpredict the final radii, a discrepancy that must originate in the shallow approximation of the lubrication analysis.

Besides this, the experimental data also fall somewhat below the results from the simulations. This is surprising given that the pre-wetted film thickness in the experiments is larger than that used in the simulations, and a larger film thickness is expected to furnish larger final radii. In other words, the experimental drops definitely do not spread as far as suggested by the simulations. This second discrepancy might have several origins. For instance, in our simulations, we modelled a freely spreading droplet. In the experiments, however, the droplets are extruded on the surface from a syringe. The deposition method might therefore be responsible. The difference is, in fact, more pronounced when the yield stress is stronger (larger values of $\mathcal {J}$), as also notable in figure 9. Indeed, in the simulations with the largest yield stresses, not all of the droplet yields during spreading, leaving intact significant plugged regions. By contrast, the action of extruding the Carbopol through the syringe forces the fluid to yield everywhere.

The fact that the flow history may impact the final state through the evolution of the plugged regions raises a second potential source for the discrepancy: the simulations exploit the Bingham model whereas the Carbopol suspensions clearly have a nonlinear plastic viscosity (see the Herschel–Bulkley fits in table 1). The shear-thinning viscosity may impact the final state by affecting the evolution of the plugs, even if that rate-dependent effect is not expected to contribute to the final force balance, or by affecting the dynamics at a nearly singular contact line (cf. King 2001a,b; Rafaï, Bonn & Boudaoud 2004). Worse, Carbopol is not an ideal yield-stress fluid. In particular, this material has been reported previously to be viscoelastic and sometimes thixotropic (Coussot et al. 2002; Luu & Forterre 2009; Balmforth et al. 2014; Coussot 2014; Dinkgreve et al. 2016; Fraggedakis, Dimakopoulos & Tsamopoulos 2016; Bonn et al. 2017), all of which might contribute to the discrepancy in final radii.

In addition to the points above, there are some other technical difficulties in the experiments that might be responsible. At the lowest $\mathcal {J}$, the nozzle may have affected the spreading and final shape, especially when the droplets were small. Moreover, there are uncertainties in the experimental values of the surface tension coefficient and yield stress that may affect the dimensionless parameters. Additionally, these values are measured for a liquid bulk and it is not clear if they hold for small and confined geometries (e.g. Geraud, Bocquet & Barentin 2013).

To conclude, in this paper, we have studied the spreading of a viscoplastic droplet on a thin film. In contrast to Newtonian droplets, a viscoplastic droplet spreading on a pre-wetted surface reaches a final shape, suggesting a practical means to control the final radius by tuning the yield stress. To gauge that final radius as a function of the physical parameters of the problem, we have provided simple scalings laws, a lubrication theory for shallow droplets, numerical simulations for deeper droplets and experiments with a Carbopol gel. Our study has direct applications in many industries, such as 3-D printing and coating processes, in which the spreading of droplets of yield-stress fluid plays a key role. Possible extensions of our work include the development of frameworks for more complicated constitutive models such as elastoviscoplastic models (e.g. Saramito 2007; Fraggedakis et al. 2016; Dimitriou & McKinley 2019; Oishi, Thompson & Martins 2019a,b), and the use of more advanced experimental tools to accurately measure droplet height and visualize the internal flow field.