Hostname: page-component-76dd75c94c-sgvz2 Total loading time: 0 Render date: 2024-04-30T09:42:21.274Z Has data issue: false hasContentIssue false

Spreading of viscoplastic droplets

Published online by Cambridge University Press:  05 March 2021

Maziyar Jalaal*
Affiliation:
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, CambridgeCB3 0WA, UK Department of Mechanical Engineering, University of British Columbia, Vancouver, BCV6T 1Z4, Canada
Boris Stoeber
Affiliation:
Department of Mechanical Engineering, University of British Columbia, Vancouver, BCV6T 1Z4, Canada Department of Electrical and Computer Engineering, University of British Columbia, Vancouver, BCV6T 1Z4, Canada
Neil J. Balmforth
Affiliation:
Department of Mathematics, University of British Columbia, Vancouver, BCV6T 1Z2, Canada
*
Email addresses for correspondence: mazi@alumni.ubc.ca, mj547@cam.ac.uk

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.

Type
JFM Papers
Copyright
© The Author(s), 2021. Published by Cambridge University Press

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

(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 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

(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. 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

(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. 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.

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 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

(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 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.

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 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

(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 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:

(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. (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.

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 (af, 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. (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).

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 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.

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 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

(A1)\begin{equation} \eta''' + \xi^{-1} \eta'' - \xi^{-2} \eta' - R_*^{2} {\mathcal{B}} \eta' = \frac{\varLambda_*}{\eta}\,\textrm{sgn} ( \eta - {\hat{h}}_\infty ) \end{equation}

with eigenvalue and film thickness parameter

(A2a,b)\begin{equation} \varLambda_* (R_*^{2} {\mathcal{B}}, {\hat{h}}_\infty) = \frac{R_*^{3} {\mathcal{J}}}{[h(0)]^{2}} \quad \textrm{and} \quad {\hat{h}}_\infty = \frac{h_\infty}{h(0)} \end{equation}

respectively, and the boundary conditions

(A3ad)\begin{equation} \eta(0)=1, \quad \eta'(0)=0, \quad \eta(1) = {\hat{h}}_\infty \quad \textrm{and} \quad \eta_{rr}(1)=0 . \end{equation}

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

(A4a,b)\begin{equation} \varLambda ( R^{2} {\mathcal{B}}, {\hat{h}}_\infty ) = \varLambda_* ( R_*^{2} {\mathcal{B}}, {\hat{h}}_\infty ) \quad \textrm{and} \quad \mathcal{I} ( R^{2} {\mathcal{B}}, {\hat{h}}_\infty ) = \frac{\int_0^{R} [h(r)-h_\infty] r\,\textrm{d}r }{R^{2} h(0)}. \end{equation}

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

(A5)\begin{equation} \mathcal{V}= 2 {\rm \pi}{\mathcal{R}}_f^{2}\mathcal{H}_f \mathcal{I} (R^{2} {\mathcal{B}}, {\hat{h}}_\infty). \end{equation}

In conjunction with (A2), we may then write

(A6a,b)\begin{equation} \frac{ {\mathcal{R}}_f}{ {\mathcal{L}}} = \left( \frac{4\varLambda}{9 {\mathcal{I}}^{2}}\right)^{1/7} {\mathcal{J}}^{-1/7} \quad \textrm{and} \quad \frac{ {\mathcal{H}}_f}{ {\mathcal{L}}} = \left( \frac{8}{27\varLambda^{2} {\mathcal{I}}^{3}}\right)^{1/7} {\mathcal{J}}^{2/7}. \end{equation}

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

(A7)\begin{equation} \mathcal{G}_{\zeta \zeta \zeta} = \frac{1}{\mathcal{G}}\,\mathrm{sgn}(\mathcal{G}-1). \end{equation}

For $\mathcal {G} \rightarrow 1$, (A7) further simplifies to

(A8)\begin{equation} \mathcal{G}_{\zeta \zeta \zeta} = \mathrm{sgn}(\mathcal{G}-1). \end{equation}

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.

Figure 10. A sketch of the decaying undulations at the edge.

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

(A9)\begin{equation} \left.\begin{gathered} I_{n}{:}\quad \mathcal{G} = 1+\mathcal{G}^{\prime}_{n}(\zeta-\zeta_{n}) + \tfrac{1}{2} \mathcal{G}^{\prime \prime}_{n}(\zeta-\zeta_{n})^{2} +\tfrac{1}{6} (\zeta-\zeta_{n})^{3}, \\ I_{n}{}^{*}{:}\quad \mathcal{G} = 1+\mathcal{G}^{\prime}_{n+1}(\zeta-\zeta_{n+1}) + \tfrac{1}{2} \mathcal{G}^{\prime \prime}_{n+1}(\zeta-\zeta_{n+1})^{2} -\tfrac{1}{6} (\zeta-\zeta_{n+1})^{3}, \end{gathered}\right\} \end{equation}

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:

(A10)\begin{equation} \left.\begin{gathered} \mathcal{G}^{\prime}_{n}+\tfrac{1}{2} \mathcal{G}^{\prime \prime}_{n} z_{A} + \tfrac{1}{6} z_{A}^{2} = 0, \\ \mathcal{G}^{\prime}_{n+1}+\tfrac{1}{2} \mathcal{G}^{\prime \prime}_{n+1} z_{B} - \tfrac{1}{6} z_{B}^{2} = 0, \\ \mathcal{G}^{\prime}_{n}+ \mathcal{G}^{\prime \prime}_{n} z_{A} + \tfrac{1}{2} z_{A}^{2} - \mathcal{G}^{\prime}_{n+1}- \mathcal{G}^{\prime \prime}_{n+1} z_{B} + \tfrac{1}{2} z_{B}^{2}=0, \\ \mathcal{G}^{\prime \prime}_{n} + z_{A} + z_{B} = 0,\\ z_{A} = \zeta^{*}_{n} - \zeta_{n}, \quad z_{B} = \zeta_{n+1} - \zeta^{*}_{n}. \end{gathered}\right\} \end{equation}

For large $n$, we seek the solution in the form

(A11)\begin{equation} \left.\begin{gathered} \mathcal{G}_{n}^{\prime} \sim \varrho^{2n} a \beta^{2},\quad \mathcal{G}_{n}^{\prime \prime} \sim \varrho^{n} \beta, \quad z_{A} \sim \varrho^{n} \beta b,\\ \mathcal{G}_{n+1}^{\prime} \sim \varrho^{2n+2} a \beta^{2}, \quad \mathcal{G}_{n+1}^{\prime \prime} \sim \varrho^{n+1} \beta ,\quad z_{B} \sim \varrho^{n+1} \beta c, \end{gathered}\right\} \end{equation}

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

(A12)\begin{equation} \left.\begin{gathered} a+\tfrac{1}{2}b+\tfrac{1}{6}b^{2}=0,\quad a+\tfrac{1}{2}c-\tfrac{1}{6}c^{2}=0,\\ a+b+\tfrac{1}{2}b^{2}-a\varrho^{2}-cr^{2}+\tfrac{1}{2}\varrho^{2}c^{2}=0,\quad 1+b-\varrho+\varrho c=0. \end{gathered}\right\} \end{equation}

These equations have a single set of valid solutions with $\varrho < 1$ (so that the undulations do not diverge), with

(A13ad)\begin{equation} \varrho = \tfrac{7}{2} - \tfrac{3}{2}\sqrt{5} \approx 0.14,\quad a = 0.312, \quad b =-1.38, \quad c =3.618. \end{equation}

Moreover, the distance along the wavetrain is given by

(A14)\begin{equation} \zeta_{n} \sim \beta ( b + c \varrho) \sum_{m=0}^{n} \varrho^{m}. \end{equation}

This geometric series has the finite limit

(A15)\begin{equation} \mathop{\textrm{Lim}}\limits_{{n \rightarrow \infty}} \zeta_{n} = \frac{\beta ( b + c \varrho)}{1-\varrho} = -\beta, \end{equation}

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}$.

Figure 11. An example of numerical grids, where the maximum level in the grid generation was eight over the interface of the droplet. Inside the droplet, the level of grids was always larger than six. The interface corresponds to the case of $\mathcal {J}=0.144$ and $\mathcal {B}=0$ at $t=1.$ The magnified view of the grids around the edge of the droplet is shown on the right.

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.

Figure 12. Effect of regularization parameter on the Gerris simulations. Results correspond to $\mathcal {J}=0.144$, $\mathcal {B}=0$ and $h_\infty =0.0175$.

References

Balmforth, N.J. & Craster, R.V. 1999 A consistent thin-layer theory for bingham plastics. J. Non-Newtonian Fluid Mech. 84 (1), 6581.CrossRefGoogle Scholar
Balmforth, N.J., Craster, R.V., Perona, P., Rust, A.C. & Sassi, R. 2007 a Viscoplastic dam breaks and the bostwick consistometer. J. Non-Newtonian Fluid Mech. 142 (1–3), 6378.CrossRefGoogle Scholar
Balmforth, N.J., Craster, R.V., Rust, A.C. & Sassi, R. 2006 Viscoplastic flow over an inclined surface. J. Non-Newtonian Fluid Mech. 139 (1), 103127.CrossRefGoogle Scholar
Balmforth, N.J., Frigaard, I.A. & Ovarlez, G. 2014 Yielding to stress: recent developments in viscoplastic fluid mechanics. Annu. Rev. Fluid Mech. 46, 121146.CrossRefGoogle Scholar
Balmforth, N.J., Ghadge, S. & Myers, T. 2007 b Surface tension driven fingering of a viscoplastic film. J. Non-Newtonian Fluid Mech. 142 (1), 143149.CrossRefGoogle Scholar
Barnes, H.A. 1999 The yield stress—a review or ‘$\pi \alpha \nu \tau \alpha \ \rho \varepsilon \iota$’—everything flows? J. Non-Newtonian Fluid Mech. 81 (1), 133178.CrossRefGoogle Scholar
Bergemann, N., Juel, A. & Heil, M. 2018 Viscous drops on a layer of the same fluid: from sinking, wedging and spreading to their long-time evolution. J. Fluid Mech. 843, 128.CrossRefGoogle Scholar
Blackwell, B.C., Deetjen, M.E., Gaudio, J.E. & Ewoldt, R.H. 2015 Sticking and splashing in yield-stress fluid drop impacts on coated surfaces. Phys. Fluids 27 (4), 043101.CrossRefGoogle Scholar
Blake, S. 1990 Viscoplastic models of lava domes. In Lava Flows and Domes (ed. J.H. Fink), pp. 88–126. Springer.CrossRefGoogle Scholar
Bonn, D., Denn, M.M., Berthier, L., Divoux, T. & Manneville, S. 2017 Yield stress materials in soft condensed matter. Rev. Mod. Phys. 89 (3), 035005.CrossRefGoogle Scholar
Bonn, D., Eggers, J., Indekeu, J., Meunier, J. & Rolley, E. 2009 Wetting and spreading. Rev. Mod. Phys. 81 (2), 739.CrossRefGoogle Scholar
Boujlel, J. & Coussot, P. 2013 Measuring the surface tension of yield stress fluids. Soft Matter 9 (25), 58985908.CrossRefGoogle Scholar
Chen, S. & Bertola, V. 2017 Morphology of viscoplastic drop impact on viscoplastic surfaces. Soft Matter 13 (4), 711719.CrossRefGoogle ScholarPubMed
Coussot, P. 2014 Yield stress fluid flows: a review of experimental data. J. Non-Newtonian Fluid Mech. 211, 3149.CrossRefGoogle Scholar
Coussot, P., Nguyen, Q.D., Huynh, H.T. & Bonn, D. 2002 Avalanche behavior in yield stress fluids. Phys. Rev. Lett. 88 (17), 175501.CrossRefGoogle ScholarPubMed
Craster, R.V. & Matar, O.K. 2009 Dynamics and stability of thin liquid films. Rev. Mod. Phys. 81 (3), 1131.CrossRefGoogle Scholar
Derby, B. 2010 Inkjet printing of functional and structural materials: fluid property requirements, feature stability, and resolution. Annu. Rev. Mater. Res. 40, 395414.CrossRefGoogle Scholar
Dimitriou, C.J. & McKinley, G.H. 2019 A canonical framework for modeling elasto-viscoplasticity in complex fluids. J. Non-Newtonian Fluid Mech. 265, 116132.CrossRefGoogle Scholar
Dinkgreve, M., Paredes, J., Denn, M.M. & Bonn, D. 2016 On different ways of measuring ‘the’ yield stress. J. Non-Newtonian Fluid Mech. 238, 233241.CrossRefGoogle Scholar
Fraggedakis, D., Dimakopoulos, Y. & Tsamopoulos, J. 2016 Yielding the yield stress analysis: a thorough comparison of recently proposed elasto-visco-plastic (evp) fluid models. J. Non-Newtonian Fluid Mech. 236, 104122.CrossRefGoogle Scholar
Geraud, B., Bocquet, L. & Barentin, C. 2013 Confined flows of a polymer microgel. Eur. Phys. J. E 36 (3), 30.CrossRefGoogle ScholarPubMed
Géraud, B., Jørgensen, L., Petit, L., Delanoë-Ayari, H., Jop, P. & Barentin, C. 2014 Capillary rise of yield-stress fluids. Europhys. Lett. 107 (5), 58002.CrossRefGoogle Scholar
German, G. & Bertola, V. 2009 Impact of shear-thinning and yield-stress drops on solid substrates. J. Phys.: Condens. Matter 21 (37), 375111.Google ScholarPubMed
Jalaal, M. 2016 Controlled spreading of complex droplets. PhD thesis, University of British Columbia.Google Scholar
Jalaal, M. & Balmforth, N.J. 2016 Long bubbles in tubes filled with viscoplastic fluid. J. Non-Newtonian Fluid Mech. 238, 100106.CrossRefGoogle Scholar
Jalaal, M., Balmforth, N.J. & Stoeber, B. 2015 Slip of spreading viscoplastic droplets. Langmuir 31 (44), 1207112075.CrossRefGoogle ScholarPubMed
Jalaal, M., Kemper, D. & Lohse, D. 2019 a Viscoplastic water entry. J. Fluid Mech. 864, 596613.CrossRefGoogle Scholar
Jalaal, M., Seyfert, C. & Snoeijer, J.H. 2019 b Capillary ripples in thin viscous films. J. Fluid Mech. 880, 430440.CrossRefGoogle Scholar
Jørgensen, L., Le Merrer, M., Delanoë-Ayari, H. & Barentin, C. 2015 Yield stress and elasticity influence on surface tension measurements. Soft Matter 11 (25), 51115121.CrossRefGoogle ScholarPubMed
King, J.R. 2001 a The spreading of power-law fluids. In IUTAM Symposium on Free Surface Flows (ed. A.C. King & Y.D. Shikhmurzaev), pp. 153–160. Springer.CrossRefGoogle Scholar
King, J.R. 2001 b Thin-film flows and high-order degenerate parabolic equations. In IUTAM Symposium on Free Surface Flows (ed. A.C. King & Y.D. Shikhmurzaev), pp. 7–18. Springer.CrossRefGoogle Scholar
Liu, K.F. & Mei, C.C. 1989 Slow spreading of a sheet of bingham fluid on an inclined plane. J. Fluid Mech. 207, 505529.CrossRefGoogle Scholar
Liu, Y., Balmforth, N.J. & Hormozi, S. 2018 Axisymmetric viscoplastic dambreaks and the slump test. J. Non-Newtonian Fluid Mech. 258, 4557.CrossRefGoogle Scholar
Liu, Y., Balmforth, N.J., Hormozi, S. & Hewitt, D.R. 2016 Two-dimensional viscoplastic dambreaks. J. Non-Newtonian Fluid Mech. 238, 6579.CrossRefGoogle Scholar
Luu, L.-H. & Forterre, Y. 2009 Drop impact of yield-stress fluids. J. Fluid Mech. 632, 301327.CrossRefGoogle Scholar
Luu, L.-H. & Forterre, Y. 2013 Giant drag reduction in complex fluid drops on rough hydrophobic surfaces. Phys. Rev. Lett. 110 (18), 184501.CrossRefGoogle ScholarPubMed
Mackay, M.E. 2018 The importance of rheological behavior in the additive manufacturing technique material extrusion. J. Rheol. 62 (6), 15491561.CrossRefGoogle Scholar
Manglik, R.M., Wasekar, V.M. & Zhang, J. 2001 Dynamic and equilibrium surface tension of aqueous surfactant and polymeric solutions. Exp. Therm. Fluid Sci. 25 (1), 5564.CrossRefGoogle Scholar
O'Donovan, E.J. & Tanner, R.I. 1984 Numerical study of the bingham squeeze film problem. J. Non-Newtonian Fluid Mech. 15 (1), 7583.CrossRefGoogle Scholar
Oishi, C.M., Thompson, R.L. & Martins, F.P. 2019 a Impact of capillary drops of complex fluids on a solid surface. Phys. Fluids 31 (12), 123109.CrossRefGoogle Scholar
Oishi, C.M., Thompson, R.L. & Martins, F.P. 2019 b Normal and oblique drop impact of yield stress fluids with thixotropic effects. J. Fluid Mech. 876, 642679.CrossRefGoogle Scholar
Oron, A., Davis, S.H. & Bankoff, S.G. 1997 Long-scale evolution of thin liquid films. Rev. Mod. Phys. 69 (3), 931.CrossRefGoogle Scholar
Popinet, S. 2003 Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 190 (2), 572600.CrossRefGoogle Scholar
Putz, A., Frigaard, I.A. & Martinez, D.M. 2009 On the lubrication paradox and the use of regularisation methods for lubrication flows. J. Non-Newtonian Fluid Mech. 163 (1), 6277.CrossRefGoogle Scholar
Rafaï, S., Bonn, D. & Boudaoud, A. 2004 Spreading of non-Newtonian fluids on hydrophilic surfaces. J. Fluid Mech. 513, 77.CrossRefGoogle Scholar
Roussel, N. & Coussot, P. 2005 “Fifty-cent rheometer” for yield stress measurements: from slump to spreading flow. J. Rheol. 49 (3), 705718.CrossRefGoogle Scholar
Saïdi, A., Martin, C. & Magnin, A. 2010 Influence of yield stress on the fluid droplet impact control. J. Non-Newtonian Fluid Mech. 165 (11), 596606.CrossRefGoogle Scholar
Saramito, P. 2007 A new constitutive equation for elastoviscoplastic fluid flows. J. Non-Newtonian Fluid Mech. 145 (1), 114.CrossRefGoogle Scholar
Sen, S., Morales, A.G. & Ewoldt, R.H. 2020 Viscoplastic drop impact on thin films. J. Fluid Mech. 891, A27.CrossRefGoogle Scholar
Tanner, L.H. 1979 The spreading of silicone oil drops on horizontal surfaces. J. Phys. D: Appl. Phys. 12 (9), 1473.CrossRefGoogle Scholar
Thompson, A.B., Tipton, C.R., Juel, A., Hazel, A.L. & Dowling, M. 2014 Sequential deposition of overlapping droplets to form a liquid line. J. Fluid Mech. 761, 261281.CrossRefGoogle Scholar
Tuck, E.O. & Schwartz, L.W. 1990 A numerical and asymptotic study of some third-order ordinary differential equations relevant to draining and coating flows. SIAM Rev. 32 (3), 453469.CrossRefGoogle Scholar
Figure 0

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’).

Figure 1

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.

Figure 2

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).

Figure 3

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$.

Figure 4

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).

Figure 5

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 (af, respectively).

Figure 6

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.

Figure 7

Table 1. Properties of the experimental fluids.

Figure 8

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).

Figure 9

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.

Figure 10

Figure 10. A sketch of the decaying undulations at the edge.

Figure 11

Figure 11. An example of numerical grids, where the maximum level in the grid generation was eight over the interface of the droplet. Inside the droplet, the level of grids was always larger than six. The interface corresponds to the case of $\mathcal {J}=0.144$ and $\mathcal {B}=0$ at $t=1.$ The magnified view of the grids around the edge of the droplet is shown on the right.

Figure 12

Figure 12. Effect of regularization parameter on the Gerris simulations. Results correspond to $\mathcal {J}=0.144$, $\mathcal {B}=0$ and $h_\infty =0.0175$.

Jalaal et al. supplementary movie

Extrusion of a viscoplastic droplet (see panel a in figure 8.)

Download Jalaal et al. supplementary movie(Video)
Video 61.3 KB