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 Reference Barnes1999; Derby Reference Derby2010; Thompson et al. Reference Thompson, Tipton, Juel, Hazel and Dowling2014; Mackay Reference Mackay2018). 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 Reference Oron, Davis and Bankoff1997; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Craster & Matar Reference Craster and Matar2009). 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. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Craster & Matar Reference Craster and Matar2009).
The spreading of Newtonian droplets over a thin layer has been addressed previously for a range of physical regimes (see Bergemann, Juel & Heil (Reference Bergemann, Juel and Heil2018), Jalaal, Seyfert & Snoeijer (Reference Jalaal, Seyfert and Snoeijer2019b), 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 (Reference Balmforth, Frigaard and Ovarlez2014), Coussot (Reference Coussot2014) and Bonn et al. (Reference Bonn, Denn, Berthier, Divoux and Manneville2017) 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 Reference Roussel and Coussot2005; Balmforth et al. Reference Balmforth, Craster, Perona, Rust and Sassi2007a; German & Bertola Reference German and Bertola2009; Jalaal, Balmforth & Stoeber Reference Jalaal, Balmforth and Stoeber2015; Liu et al. Reference Liu, Balmforth, Hormozi and Hewitt2016; Chen & Bertola Reference Chen and Bertola2017; Liu, Balmforth & Hormozi Reference Liu, Balmforth and Hormozi2018). 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 Reference Luu and Forterre2009; Saïdi, Martin & Magnin Reference Saïdi, Martin and Magnin2010; Luu & Forterre Reference Luu and Forterre2013; Blackwell et al. Reference Blackwell, Deetjen, Gaudio and Ewoldt2015; Sen, Morales & Ewoldt Reference Sen, Morales and Ewoldt2020) or fluid interfaces (Jalaal, Kemper & Lohse Reference Jalaal, Kemper and Lohse2019a).
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
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
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
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
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 Reference Roussel and Coussot2005; Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006).
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
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$.
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. Reference Oron, Davis and Bankoff1997; Craster & Matar Reference Craster and Matar2009), the key differences arising from the impact of the yield stress on the vertical profile of the radial velocity (Liu & Mei Reference Liu and Mei1989; Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006). 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 Reference Balmforth and Craster1999; Putz, Frigaard & Martinez Reference Putz, Frigaard and Martinez2009). 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
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
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. Reference Oron, Davis and Bankoff1997; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Craster & Matar Reference Craster and Matar2009).
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 Reference Balmforth, Ghadge and Myers2007b; Jalaal & Balmforth Reference Jalaal and Balmforth2016). We discuss them further below and in greater detail in appendix A.
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
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.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
(cf. Blake Reference Blake1990; Roussel & Coussot Reference Roussel and Coussot2005; Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006). 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
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 Reference Tanner1979; Tuck & Schwartz Reference Tuck and Schwartz1990; Jalaal et al. Reference Jalaal, Seyfert and Snoeijer2019b), 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.
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).
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
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.
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 Reference Popinet2003). 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
where
(cf. O'Donovan & Tanner Reference O'Donovan and Tanner1984). 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:
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. (Reference Liu, Balmforth, Hormozi and Hewitt2016, Reference Liu, Balmforth and Hormozi2018). 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.
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. (Reference Jalaal, Balmforth and Stoeber2015) 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).
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).
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 Reference Manglik, Wasekar and Zhang2001; Boujlel & Coussot Reference Boujlel and Coussot2013; Géraud et al. Reference Géraud, Jørgensen, Petit, Delanoë-Ayari, Jop and Barentin2014; Jørgensen et al. Reference Jørgensen, Le Merrer, Delanoë-Ayari and Barentin2015). 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. Reference Jalaal, Seyfert and Snoeijer2019b), 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.
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.
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 Reference King2001a,b; Rafaï, Bonn & Boudaoud Reference Rafaï, Bonn and Boudaoud2004). 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. Reference Coussot, Nguyen, Huynh and Bonn2002; Luu & Forterre Reference Luu and Forterre2009; Balmforth et al. Reference Balmforth, Frigaard and Ovarlez2014; Coussot Reference Coussot2014; Dinkgreve et al. Reference Dinkgreve, Paredes, Denn and Bonn2016; Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017), 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 Reference Geraud, Bocquet and Barentin2013).
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 Reference Saramito2007; Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016; Dimitriou & McKinley Reference Dimitriou and McKinley2019; Oishi, Thompson & Martins Reference Oishi, Thompson and Martins2019a,b), and the use of more advanced experimental tools to accurately measure droplet height and visualize the internal flow field.
Supplementary movie
A supplementary movie is available at https://doi.org/10.1017/jfm.2020.886.
Acknowledgements
M.J. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada through a Vanier Canada Graduate Scholarship.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Viscoplastic final shapes and contact lines
A.1. Computational details
We may streamline the path to the solution of (3.4) by introducing the new variables $\xi =r/R_*$ and $\eta (\xi )=h(r)/h(0)$. The task is then to solve
with eigenvalue and film thickness parameter
respectively, and the boundary conditions
From this solution, we may then find $\eta _*$, the first (scaled) radial position where $\eta (\xi _*)= {\hat {h}}_\infty$, corresponding to $r=R$. A simple rescaling then provides $h/h(0)$ as a function of $r/R$, from which we may compute the functions
We use MATLAB's boundary-value-problem solver BVP4C to solve (A1). To further ease the computation, we smooth out the switch of sign on the right-hand side using the function $\tanh [\varXi (\eta - {\hat {h}}_\infty )]$, with $\varXi =10^{6}$, which is sufficiently large to ensure the solution details are insensitive to the precise value. For the initial guess for the solver, we use either a final snapshot from the solution of the evolution equation, or continuation from a equilibrium profile with different parameter settings. For the solution shown in figure 3, the initial guess is taken from the final snapshot in figure 2(a), and contains a sufficient number of undulations at the edge such that the solver converges to an equilibrium profile with five switches of sign of $h-h_\infty$. The slight smoothing of those switches is visible in the plot of $h_{rr}$ at the right of figure 3(b). To map out the functions $\varLambda (R^{2} {\mathcal {B}}, {\hat {h}}_\infty )$ and $\mathcal {I}(R^{2} {\mathcal {B}}, {\hat {h}}_\infty )$, as required below, we accelerate the computations by using a shorter initial guess with only three sign switches.
The final step is to consider the volume constraint, which becomes
In conjunction with (A2), we may then write
The dependence of the function $\varLambda$ over ${\mathcal {I}}$ on the arguments $R^{2} {\mathcal {B}}$ and ${\hat {h}}_\infty \equiv h_\infty {\mathcal {L}}/ {\mathcal {H}}_f$ ensures that these relations constitute a pair of implicit algebraic equations for $R= {\mathcal {R}}_f/ {\mathcal {L}}$ and ${\mathcal {H}}_f/ {\mathcal {L}}$. Moreover, the relations in (A6) indicate that the functional dependence on the parameters of the problem is through the combinations ${\mathcal {B}} {\mathcal {J}}^{-2/7}$ and $h_\infty {\mathcal {J}}^{-2/7}$ (or $h_\infty / {\mathcal {B}}$), as in (3.7). Alternatively, one can map out the two functions on the $(R^{2} {\mathcal {B}}, {\hat {h}}_\infty )$ parameter plane, and then interpolate onto a grid of the physical parameters, ${\mathcal {B}} {\mathcal {J}}^{-2/7}$ and $h_\infty / {\mathcal {B}}$, as done in figure 5.
A.2. Edge structure
To analyse the structure at the edge, we consider the limit $h_\infty \ll 1$ and then resolve the narrow undulations there by introducing the new variables $\mathcal {G}(\zeta )=h/h_{\infty }$ and $\zeta = (r-R)( {\mathcal {J}}/h_{\infty }^{2})^{1/3}$. For finite ${\mathcal {B}}$ and to leading order, we then find
For $\mathcal {G} \rightarrow 1$, (A7) further simplifies to
The wavetrain therefore limits to a sequence of cubic polynomials, patched together to make the solution and its first two derivatives continuous, as illustrated in figure 10. The approach to the pre-wetted film $\mathcal {G}=1$ is thereby achieved by passing through an infinite sequence of switches in the sign of $\mathcal {G}-1$, with the second derivative $\mathcal {G}_{\zeta \zeta }$ taking a decaying sawtooth waveform.
We split up the wavetrain into the intervals $\{I_n\}$ and $\{I_n^{*}\}$ between the sign switches of $\mathcal {G}-1$. The $n$th interval over which $\mathcal {G}>1$ is $I_n=[\zeta _n,\zeta _n^{*}]$, and the following one, where $\mathcal {G}<1$, is $I_n^{*}=[\zeta _n^{*},\zeta _{n+1}]$. We then write
where primes denote the spatial derivatives with respect to $\zeta$ and the subscripts refer to the locations $\zeta _n$ where they are evaluated. Matching these solutions together at $\zeta _{n}^{*}$ leads to a system of algebraic equations:
For large $n$, we seek the solution in the form
where $\beta$ is an arbitrary amplitude that must be fixed by matching to the solution at the beginning of the wavetrain. The remaining constants satisfy
These equations have a single set of valid solutions with $\varrho < 1$ (so that the undulations do not diverge), with
Moreover, the distance along the wavetrain is given by
This geometric series has the finite limit
implying the wavetrain ends at finite radius.
The preceding analysis establishes that the solution for the wavetrain has a bounded solution in terms of the rescaled variables $\zeta$ and $\mathcal {G}$. In the limit $h_\infty \to 0$, the wavetrain therefore shrinks to a point. Moreover, given $(h,h_r,h_{rr})=(h_\infty \mathcal {G}, {\mathcal {J}}^{1/3}h_\infty ^{1/3}\mathcal {G}_\zeta , {\mathcal {J}}^{2/3}h_\infty ^{-1/3}\mathcal {G}_{\zeta \zeta })$, the solution for $r < R$ must approach the edge with $(h,h_r)\to 0$ for $r\to R$, but a diverging second derivative. Analysis of the singular point of (3.4) at $r=R$ then establishes the local solution used in § 3.2.2 to compute the solution for $h_\infty \to 0$.
Appendix B. Gerris simulations
The details of the Gerris simulations can be found in Jalaal (Reference Jalaal2016). We use a rectangular domain that is sufficiently large to eliminate the influence of the outer radial and top boundaries; see figure 11, which illustrates the adaptive gridding for an evolving solution. Symmetry and no-slip boundary conditions were applied along the centreline and bottom surface, respectively. ‘Outflow’ conditions were applied on the other two boundaries. We set the density and viscosity of the ambient medium above the viscoplastic fluid to be $10^{-2}\rho$ and $10^{-2}\mu$, respectively, which ensures that this fluid does not influence the dynamics (as confirmed through a number of other simulations, changing the density and viscosity ratios from $0.1$ to $0.002$). The computations are continued until the value of kinetic energy fell below $10^{-6} \sigma \mathcal {L}^{2}$, at which point a quasi-steady shape was largely established and before any residual spreading arose due to the regularized viscosity $\mu _{max}$.
To verify that the simulations were independent of the regularization parameter, we conducted numerical simulations for different $\mu _{max}$ and established that the changes to the final shapes were insignificant for $\mu _{max} / \mu > 10^{4}$. Figure 12 shows an example of final shapes for different regularization parameter.