The stressful way of droplets along single-fiber strands: A computational analysis

Droplets wetting and moving on fibers are omnipresent in both nature and industry—from spider web arrangement to engineering applications in the textile,^1,2 filters,^3,4 and the paper industry.

Plateau⁵ has been the first to study droplets spreading on fibers. He has observed that only due to surface tension, a liquid covering a fiber would turn into a string of beads, taking the shape of unduloids. Following up on these studies, many physicists have started investigating the flow of liquid films on threads, notably Lord Rayleigh,⁶ who has studied liquid film instabilities for falling cylindrical jets, which has been later expanded upon by Yarin et al.⁷ and conversely studied by Kliakhandler,⁸ Craster,⁹ Duprat.¹⁰ The analytical solution to the unduloidal shape studied by Plateau has been later derived by Carroll^11–13 who has analyzed droplets, which are molded by capillary forces. This analytical solution has been validated in the same paper for various systems and has been shown to be accurate for droplets that are not influenced by gravitational forces, are axisymmetrical, and present a so-called “barrel shape,” as opposed to a “clam-shell shape.”

The motion of single droplets on fibers has also been extensively studied, where the droplet is set in motion due to gravitational effects,^14,15 asymmetric slug in a tube,¹⁶ temperature gradients,¹⁷ and also due to radius gradients.¹⁸ However, to the best of the authors' knowledge, no studies exist which investigate the transient droplet motion on fibers by means of direct numerical simulations (DNS), that is, resolving all relevant spatiotemporal scales. Thus far, numerical studies have been concerned with the equilibrium shape of droplet at fibers only. Venkateshan et al.¹⁹ have studied the 3D shape of water droplets with different volumes on various fibrous structures, using the finite element code Surface Evolver. Also, using the same approach, Aziz et al.^20–22 have investigated wetting of a liquid droplet on a fiber and fibrous coatings. In particular, the authors have studied the shape of a droplet deposited on a single fiber, and the effect of the capillary force exerted on a fiber by a droplet and the volume residue after the droplet detachment. Modern manufacturing techniques allow now to produce fibers with controlled heterogeneous wettability,²³ which are of interest for various applications including oil–water separation²⁴ and catalysis.²⁵ The motion of a droplet on a fiber with controlled heterogeneous wettability due to a chemically inhomogeneous surface has not been studied numerically at all so far.

This contribution is concerned with a detailed numerical investigation of the transient spreading and motion of three-dimensional droplets on single fiber using direct numerical simulations. We aim to analyze the contribution of the local forces acting on a sliding droplet on a chemically homogeneous as well as heterogeneous (patterned) fiber strand. For this purpose, we deploy a phase-field method, which has been shown to be particularly suited for accurately describing transient wetting processes with contact line dynamics,^26–28 particularly since the underlying diffuse interface model inherently allows for a moving three-phase contact line in combination with the no-slip condition at solid walls.^29–31 

Our diffuse interface phase-field solver is used—phaseFieldFoam, implemented in OpenFOAM (FOAM-extend 4.0 and 4.1). The solver phaseFieldFoam has been extensively validated for several cases of static and dynamic wetting, on both simple and complex substrates, such as chemically patterned surfaces.^32–35 The objective of this work is twofold, viz.

• to detail on crucial methodological aspects to devise a high-fidelity simulation of wetting processes using the phase-field approach and finite-volume discretization with particular focus on non-planar substrates such as fibers,

• to provide local insight into the hydrodynamics of droplets spreading and moving on single fibers, in particular on the influence of wettability of chemically homogeneous substrates as well as of patterns of chemically heterogeneous substrates on local stresses and forces.

We assume two immiscible Newtonian fluid phases under isochoric and isothermal conditions. Then, the coupled Cahn–Hilliard Navier–Stokes equations read in semi-closed formulation (see Ref. 36)

where C denotes the phase-field order parameter, and $p̃$ is a modified pressure, since parts of the known Korteweg tensor term accounting for interfacial capillarity have been absorbed into the pressure gradient term. Assuming a liquid of Newtonian fluid rheology, $τ=2μdevD$⁠, where D is the rate-of-deformation (rate-of-strain) tensor, $D=12[∇u+(∇u)T]$⁠, and $devD$ denotes the deviatoric, trace-free part [$tr(devD)=0$]. Moreover, ρ and μ are the volumetric average density and dynamic viscosity of the fluids, given from the phasic densities and viscosities as $ρ=1+C2ρ1+1−C2ρ2$ and $μ=1+C2μ1+1−C2μ2$⁠, and $fb=ρg$ is the body force due to gravity. Note that herein we are using the contrast of volumetric phase fractions as order parameter C that is $C∈[−1,1]$⁠, opposed to concentrations in some other treatises.

With this, the set of governing equations (1) is closed—apart from J, which denotes the phase-field flux. Following Landau and Lifshitz,³⁷ it is governed by generalizing Fick's law as $J=−M ∇Φ$⁠, where $Φ$ is the chemical potential and M is the mobility. In other words, under the above assumptions, the dynamics of the phase-field system is fundamentally driven by the gradient in the local chemical potential $Φ$⁠. Note that the term $∇·(uJ)$ is required for thermodynamic consistency.^29,38

Within the diffuse-interface model framework, following Cahn and Hilliard,³⁹ the local chemical potential arises from the variational derivative of the total free energy F (Helmholtz free energy), which is considered as a functional of the order parameter C, viz.

where F_m denotes the bulk (mixing) contribution to the total free energy of the system

Herein, Ω is the fluid domain and $∂Ω$ is its boundary, a solid surface not permeable to fluids. The local mixing energy and the local wall energy densities are denoted f_m and f_w, respectively. The mixing energy density (Helmholtz free energy density)³⁹ can be written as

The first term on the right-hand side (RHS) of (4) is the gradient energy, which represents the interfacial energy density describing non-local interactions between the two components promoting complete mixing of the fluids, while the second term is the so-called bulk energy density, which models the counter tendency to separate.²⁸ 

The wall free energy density in (3) can be written as

where $Ψw(C)$ obeys the bulk limits $Ψw(C=−1)=0$ and $Ψw(C=1)=1$⁠, respectively.

A variational procedure, that is, variation in F_w and integration by parts, reveals that

where λ is the mixing energy coefficient and $fw′(C)=(σSG−σSL)Ψw′(C)$ denotes the derivative of the wall free energy density with respect to the order parameter. Effectively (5) postulates that the wall free energy density is a function only of the fluid composition next to the wall.²⁸ Equation (6) is commonly referred to as wetting condition. It describes the diffusive relaxation process of the dynamic contact angle locally constraining toward the equilibrium angle to leading order.^28,30

Now, choosing a particular double-well potential function $Ψ(C)$ allows closure. Using the phenomenological double-well function of polynomial form according Ginzburg and Landau, $Ψ(C)=14(C2−1)2$⁠, the equilibrium profile of the order parameter can be derived by minimization of the free energy: requiring $Φ:=δFδC=λε2Ψ′(C)−λΔC=!0$⁠, where $Φ$ denotes the chemical potential and, assuming a planar interface, leads to

where n denotes the coordinate normal to the interface, commonly referred to as interfacial capillary width. Within $3/2ε$⁠, the order parameter C varies from about −0.9 to 0.9. Under these assumptions, one can also show that the mixing energy coefficient λ within the diffuse interface context can be related to the surface tension coefficient employed in sharp interface models by

With this, we can exploit fundamental geometrical relations: if the interface makes a contact angle θ₀ with a wall (see Fig. 1), using a coordinate system, which holds the wall normal n_w as one coordinate and the wall tangent t_w

Differentiating the equilibrium profile C(n) along the coordinate normal to the wall as

allows us to recognize

from (7) and

from a simple geometrical observation. Thus, we finally arrive at the wetting boundary condition

Using the Young equation,⁴⁰ 

it can be shown that (5) can be written in closed form as³⁰ 

Here, we have used the wetting condition (6), since $Ψw(C)=C(C2−3)−24$ is a closure restoring (11), inserting λ according to (8) and $(σSG−σSL)$ according to (12) into (6). Note carefully how the particular choice of the double-well function $Ψ(C)$ according to Ginzburg and Landau has allowed us to consistently close for the phase-field transport equation and the wetting boundary condition. Both closures are indeed connected. Moreover, zero flux across the wall requires for the chemical potential that

The Navier–Stokes equations are supplemented by the no-slip wall boundary conditions on the solid substrate $∂Ωw$⁠, viz.

intuitively requiring that the velocity of the fluid at the wall $u|w$ is equal to the velocity of the wall $uw$⁠. Equivalently, one can demand for the normal velocity component that

that is, to fulfill a no-penetration condition due to impermeability of the wall, and for the tangential velocity component that

where $Pw:=I−nwnw$ is the projection operator onto the tangent plane of the boundary with outer unit normal $nw$⁠, and I denotes the identity or unit diagonal tensor. Moreover, it is

where K denotes the wall-tangential contribution of the viscous boundary force per unit area, $K:=μ(∇u+(∇u)T)nw$⁠, which amounts to the viscous wall stress $τw$ exerted to the wall.

In the following, we attempt to provide an overview of the numerical method of the phaseFieldFoam solver. One focus shall be on the discretization of the Cahn–Hilliard equation, since it is a non-linear fourth-order partial differential equation (PDE) and thus challenging to treat numerically in an accurate yet robust manner. Hence, special care must be taken in the solution method as set out in the reminder. Moreover, the Moving Reference Frame (MRF) technique and its implementation are described, since it is crucial for an efficient deployment of computational resources by significantly reducing the required size of the computation domain. Moreover, we detail on the requirement of a consistent no-slip boundary condition at fiber walls, which is different from the de facto standard Dirichlet boundary condition imposing a zero velocity field at walls in the linear momentum equation.

The Cahn–Hilliard equation is a non-linear fourth-order PDE and thus challenging to treat numerically in an accurate yet robust manner. The difficulty originates from its right-hand side. Illustratively, this can be seen by considering the mechanism of interfacial evolution, which within the diffuse-interface model theory is governed by dissipation of the free energy functional according to (3) as⁴¹ 

Note that (19) is void of the advection term, which is not needed here for making this point: in order to promote numerical stability, it is desirable to employ a temporal discretization, which satisfies the discrete energy stability condition, viz. the decay of energy with time as

where the superscripts n and o denote the new and old time levels, respectively. It is a common practice to devise a splitting technique to avoid implicit schemes, which come at the cost of solving a non-linear non-convex problem, that is, to incorporate the RHS of (19) at the new time level. While such schemes enjoy stability when solving the inherently stiff phase-field equations, there are alternative schemes offering a better performance. Depending on the approximation of the potential terms in $fm(Cn, Co)$ and $fw(Cn, Co)$⁠, see (3), one arrives at different schemes with distinct stability properties. One distinguishes linear and non-linear approximations of the potential terms. We shall continue to focus on linear schemes, since they do not suffer from the significant computational overhead of non-linear schemes, which require iterative methods. More particular, in the remainder we consider the well-established linear Eyre approximation, dating back to the work of Elliott and Stuart⁴² but named after Eyre,⁴³ and the optimal dissipation approximation.⁴⁴ For the sake of brevity, we shall refer to these schemes as stable and optimal, according to their properties, namely, unconditional energy stability due the existence of numerical dissipation, and conditional energy stability but optimal dissipation, respectively. For an overview and detailed discussion, the reader is referred to the work of Tierra and Guillén-González.⁴¹ Since Ref. 41 details on the mixing energy potential f_m term, we shall further restrict ourselves in the remainder to the wall energy potential f_w and its numerical treatment in the boundary condition for the order parameter. It is important to note that both terms have to be approximated consistently, that is, using the same scheme, so as to simulate wetting processes at high physical fidelity by phase-field methods.

For the discretization of the temporal term (ddtSchemes), a first-order accurate Euler scheme is employed. Regarding the advection term discretization (divSchemes), for the divergence term of the scalar transport of the order parameter a high-resolution scheme (Gamma) is used, while the convection term of the momentum equation is discretized using limitedLinearV scheme, which is second-order accurate. The gradient term (gradSchemes) is discretized with the second-order accurate (linear) scheme and the Laplacian term (laplacianSchemes) with (linear corrected).

The wetting boundary condition (11) for the phase-field evolution equation can be also stated in general notation as so-called convective boundary condition (or boundary condition of third kind), viz.

For implementation, we cast this as Robin or mixed-type boundary condition, viz.

where $ϕb$ refers to the boundary patch value (i.e., the value on the boundary face b) and $ϕP$ is the internal cell value adjacent to the boundary face, cp. Fig. 2. Moreover, $ωb∈[0,1]$ are weights, $dn$ is face-normal component of the distance vector between the boundary face center and the adjacent cell center, and $ϕref$ and $gref$ denote a reference value and reference gradient, respectively. To deploy a Robin (mixed) boundary condition, appropriate values for $ϕref, gref$⁠, and ω_b must be found. For this, (21) is linearized as

With this, we obtain

From comparison with (22), we find $ωb=hk/|dn|+h$ and $gref=0$⁠. It remains to determine k, h, and $ϕref$ as well as $ϕP$ from comparing with (6). Immediately, one identifies $k=λ$⁠. To take advantage from a semi-implicit discretization in a linear scheme, we note that $ϕb=Cbn$⁠, where n denotes the new time level. Moreover, with the wall free energy density, cf. (5), which is repeated here for convenience,

where the non-linear wall potential is $Ψw(C)=C3−3C4$ for the chosen phenomenological Ginzburg–Landau potential, and where $σSG−σSL=σ cos θ0$ [cf. (12)], the first derivative of $fw(C)$ with respect to the order parameter C can be written as

Adding dissipation to this yields the stable scheme, which is of first order in time,

and choosing the parameter $β=Co$ leads to the optimal scheme, a second order in time scheme,

Comparison between (6) and (21) then allows to eventually determine h and $ϕref$⁠. The parameters are summarized in Table I. It is noteworthy that these parameters are valid only for the chosen phenomenological Ginzburg–Landau potential and must be changed upon choice of other potentials.

Tierra and Guillén-González⁴¹ have noted that the energy stability condition is only satisfied for appropriate values of β, viz. $β≥(3(Cn)2−1)/2$⁠. With such a condition, a value of β can be determined large enough to assure that the scheme is energy-stable. However, this has limits: it has been shown in Ref. 44 that wrong equilibrium solutions can be the consequence of introducing too much numerical dissipation.

From Sec. II C, it is evident that special attention has to be devoted when discretizing the viscous stress term in the momentum equation. In particular, (18) must be enforced on discrete level.

For Newtonian fluids, the viscous stress tensor can be closed as $τ=μ(∇u+(∇u)T)$⁠. Then, deploying the finite volume method

where i and b denote inner and boundary faces, respectively. Here, we have used that $∇·(μ(∇u)T)=(∇u)·∇μ+μ∇(∇·u)$⁠. For a single boundary face b, finite volume discretization yields

where the last identity shows the split into explicit contribution to the source vector (first term) and the implicit contribution to the system matrix (second term), respectively.

However, while this is accurate for orthogonal Cartesian boundary cells, for non-planar boundary patches (skewed meshes) one must ensure that the applied boundary condition guarantees that the wall shear stress $τw$ is tangent to the wall along with $uw$⁠,⁴⁵ cp. (18),

In order to reduce computational costs, a Moving Reference Frame (MRF) technique has been employed. Following the droplet's center of mass during the simulation, such a technique allows to significantly reduce the computational domain size.

Here, we are dealing with a non-inertial non-rotating reference frame. This method has been successfully used in rising bubble simulations^46,47 and is convenient to implement due to its simplicity.

Figure 3 depicts the situation for a droplet sliding on the surface of a vertically aligned fiber: the droplet slides with respect to the inertial reference frame (x, y, z) but remains centered within the computational domain in the non-inertial frame of reference $(x̂,ŷ,ẑ)$⁠.

For this, we add the frame acceleration $aF$ to the momentum equation (1) as

The acceleration $aF$ and velocity $uF$ of the moving reference frame relative to the inertial reference frame are given by

Herein, $xF$ is the droplet center of mass relative to the moving reference frame.

The velocity of the moving reference frame is adjusted at the beginning of each time step so as to keep the droplet barycenter at a fixed location within the computational domain. This is done by applying a velocity correction $ΔuF$ for each time step. This correction is computed using a PID controller where the control variable is the droplet's center position vector. This vector is compared to a given constant target position. The applied velocity correction is

where

Following the proposal in Ref. 46, a PID controller suffers significantly less from oscillations compared to a PD controller, which has been used in Ref. 47.

At every time step, the frame velocity field is then updated

and used according to (32) for the calculation of the reference frame acceleration.

Since there also exists a solid wall in contact with the droplet, we impose a wall velocity, which is equal in value to the frame velocity but opposite in sign, viz.

Effectively, this ensures that the droplet barycenter is kept in place.

Note that our phase-field solver has been previously subject to significant validation efforts regarding wetting and de-wetting processes⁴⁸ and applied to different complex wetting physics.^32,49 Here, we focus on three-dimensional cases relevant to droplets wetting on fibers such as the evolution of the spreading diameter and equilibrium shape of a droplet on flat and particularly spherical surfaces,⁵⁰ as well as the spreading of a droplet on single fiber strand, where the dimensions of the equilibrium shape of an axisymmetric droplet on a cylinder are compared with literature-known reference data.^15,51,52

For model calibration, we have initially pursued a parameter study based on the first case (see Sec. IV A) in order to determine free model parameters of the phase-field method, which have been then fixed for all further simulations in this study. Because all the cases examined are capillary-dominated transient flows with low Reynolds numbers, we use an adaptive time step criteria, where the time steps are adjusted based on the maximum Brackbill–Kothe–Zemach (BKZ) criterion $ΔtmaxBKZ=1+52ρσh3$⁠, as proposed in Ref. 53. In this work, we set $CFLBKZ=1$⁠. For time accuracy, we further restrict the maximum time step size computed according this BKZ criterion and require $Δtmax=2×10−6s$⁠.

The equilibrium shape of an oil droplet in water on a flat surface with static contact angle, θ₀, is investigated. A schematic representation of the initial and equilibrium shapes of the droplet on a flat surface with static contact angle is shown in Fig. 4. Table II lists the physical properties of the system.⁵⁰ 

For initialization, a droplet with radius $R0=1×10−3m$ is placed on a smooth flat surface such that its center is at the distance of the initial droplet radius R₀ from the surface. The computational domain is of size $0.004×0.004×0.004 m3$⁠. Since the initial contact angle is different from the equilibrium contact angle, the droplet will start spreading until it reaches its equilibrium shape. Due to absence of gravity in this case, the Eötvös number, which describes the ratio between gravitational and capillary forces, is zero, $Eo:=Δρ gR02/σ=0$⁠. The no-slip boundary condition is applied at the bottom boundary with free-slip boundary conditions being applied on every other boundary.

In order to calibrate the underlying diffuse-interface model for all subsequent simulations, we determine the minimum number of cells $NC:=4εdx$ (where d_x is the cell size at the interface) required to sufficiently resolve the diffuse interface and the highest Cahn number $Cn:=εD0$⁠, which relates the capillary width to reference length scale, that correctly models interface dynamics. Here, the spreading diameter is investigated for various values of the above-mentioned parameters.

Too low numbers of cells N_C to resolve the diffuse interface result with nonphysical interface evolution, and consequently wrong spreading diameters as can be seen in Fig. 5. We show the evolution of the dimensionless spreading diameter $S/D0$ as a function of the dimensionless time $t/τ$ with $τ:=ρR03/σ$⁠, where σ is the surface tension coefficient, for a contact angle $θ0=60°$ and $Cn=0.02$⁠, with varying N_C.

From this, we conclude that at least six cells are required to resolve the interface sufficiently, which is in accordance our previous study⁴⁸ where simulations of a rising bubble were performed.

Fixing N_C = 6, the effect of the Cahn number can be studied. We have varied the Cahn number as $Cn∈{0.01,0.02,0.04}$⁠. Results are shown in Fig. 6. As the Cahn number decreases, the simulations start to converge. More specifically, Yue³⁰ showed that for $Cn≤0.02$⁠, the simulated results are independent of $Cn$⁠. Thus, we choose $Cn=0.02$ to calculate the value of the capillary width ε that is used in simulations.

Fixing $Cn=0.02$⁠, the final parameter to be determined is the scalar mobility M, cp. (1). Yue et al.³⁰ have shown that the mobility should be adjusted by a scaling formula. The mobility used in our simulations can be calculated based on $Cn=4S=4Mμe/D0$⁠, where $μe=μ1μ2$ is the effective viscosity, which is used due to the highly dissimilar viscosities. Here, S reflects the diffusion length scale at the contact line.³⁰ 

The spreading of a droplet has two stages—an initial stage where inertial–capillary forces are controlling droplet motion and a second stage where spreading is dominated by viscous forces. It has been observed that the spreading radius follows a power-law scaling in time⁵⁴ 

where S is the spreading diameter, c is a scaling factor, which depends on equilibrium contact angle, and n is a scaling exponent, which has been observed by Winkels et al.⁵⁵ to not be influenced by wetting properties. Other authors^54,56 have seen that the contact angle does influence the scaling exponent, and report that it decreases with increasing contact angle. Das et al.⁵⁷ have noted that this discrepancy may be due to the fact that system's properties are different, and as such, the contribution of the viscous forces is changed with each study.

We investigate the variation of the normalized spreading diameter for five values of equilibrium contact angles, $θ0∈{30°,60°,90°,120°,150°}$⁠. The equivalent spreading diameter S/D is displayed in Fig. 7(a). Figure 7(b) shows the variation of n and c when changing the equilibrium contact angle.

By fitting Eq. (37) to the monotonic regime of data points obtained used in Fig. 7(a), the values of wetting pre-factor c and the scaling exponent n have been computed for the various contact angles, both of which are influenced by the changes in equilibrium contact angle θ₀, shown in Fig. 7(b).

We have found that n decreases linearly with increasing θ₀, except for contact angles lower than $90°$⁠, where the exponent's value stabilizes at $n≈1/2$⁠. This agrees with the observations by Biance et al.,⁵¹ who conducted experiments of water droplets spreading on treated glass plates. Bird et al.⁵⁴ observed a similar behavior, where they have seen the scaling exponent monotonically decreases with increasing equilibrium contact angle.

Despite this, the results obtained by us and the previously mentioned authors do not agree with the observation made by Legendre et al.⁵⁶ where numerical simulations of spreading droplets were performed and Das et al.⁵⁷ who investigated the spreading behavior of droplets over a bundle of fibers: they found that n is independent of θ₀.

Das et al.⁵⁷ have used the Laplace number $La:=σρR0/μ2$⁠, where μ is the dynamic viscosity of the oil droplet, which relates the inertial and surface tension forces to viscous forces, to explain these discrepancies. They have adopted the hypothesis that for $La∼101$ and lower, the effect of the contact angle on the spreading dynamics was negligible. Legendre et al.⁵⁶ showed that n varies from 2/3 to 1/2 when La varies from 10⁰ to 10⁴, for a contact angle of $65°$⁠.

For our simulation setup, $La∼102$⁠, which could indicate that for these Laplace numbers and for non-wetting liquids (⁠$θ0<90°$⁠), the scaling exponent value is close to 1/2, as it is seen in Fig. 7(b).

Moreover, in the studies of Bird et al.,⁵⁴ it was observed the scaling exponent gives values close to 1/2 for $La∼102$⁠, much like the results we obtain.

Regarding the pre-factor c, it demonstrated that it also decreases with increasing contact angle, which was also noted in Refs. 54–57.

One should note that the droplet's spreading diameter is not zero upon initialization. This is due to grid resolution in our numerical simulation, causing an overlap between the droplet and the surface. This can be fixed by simply increasing the grid resolution such that $S/D(t/τ=0)$ tends to zero. Nevertheless, we have shown that our simulated results agree with experimental observations, by capturing the scaling $S/D∼τn$ where n varies from $∼1/2$ to $∼1/4$⁠, across the various equilibrium contact angles chosen.

The equilibrium shape of the droplet spreading in the absence of gravity can be derived analytically, where the equilibrium contact radius r_f, curvature radius R_f, and droplet height h_f are calculated from⁵⁰ 

Figure 8 displays the characteristic radii and heights of the droplets from simulations against the ones obtained analytically, for the same equilibrium contact angles. The results show a very good agreement between simulated and analytical results, also substantiating that the parameter values for N_C, $Cn$⁠, and M are appropriate.

For further validation, we consider a case of particular relevance for spreading on non-planar substrates such as fibers: we perform simulations to study the equilibrium shapes of droplets spreading on a spherical surface.

As shown in Fig. 9, a solid sphere with a radius equal to the droplet radius is introduced in the center of the previous computational domain. The droplet is initialized in such a way that its center is at a distance equal to R₀ from the spherical surface, just touching it. Since the initial contact angle is different from the static contact angle, the droplet will spread until it reaches equilibrium. Physical properties of the system are the same as in Sec. IV.

The equilibrium shapes of droplets on spherical surfaces can be obtained by solving the following non-linear system:⁵⁰ 

For the five contact angles considered, the simulated results of the characteristic dimensions of the spreading droplet are compared with the ones obtained analytically, as shown in Fig. 10. Simulation results are in very good agreement with analytical results.

In a final validation study, we investigate the equilibrium shape of droplets spreading on a fiber. The unduloidal shape of a wetting droplet on a fiber cannot be described with simple analytical equations. Despite this, looking at the asymptotic regimes of very large droplets (⁠$Ω≫dv3$⁠) and very small ones (⁠$Ω≪dv3$⁠), where Ω is the droplet volume and d_v is the fiber diameter, one can derive analytical results of the full transition region between these two cases.^15,52

Large droplets tend to keep their spherical shape when they are placed on a fiber, thus $W≈R0$⁠, with W being the system width and R₀ the deposited droplet diameter. For small droplet sizes, the droplet spreads on the fiber in such a way that the curvature of its interface is only slightly lower than the fiber curvature $2/dv$⁠, for which a shape approximates to that of a cylinder. Gilet et al.¹⁵ have derived asymptotics for these two scenarios, which are succinctly shown in Table III. Herein, X is the total length of the droplet.

With the asymptotics from Table III, the dimensions of both large and small droplets on a vertical fiber can be assessed quantitatively. Gilet et al.¹⁴ have demonstrated an excellent agreement between the theoretical and experimental results for various $Ω/dv3$⁠.

We perform simulations of droplets spreading on fibers for various $Ω/dv$⁠, where both the width W and extension X are compared to the asymptotic solutions—see Fig. 11, for a contact angle of $15°$⁠. It can be seen that the simulations results are in very good agreement with the theoretical predictions, substantiating the capability of our solver to predict the equilibrium shapes of large and small droplets spreading on a fiber.

Varying the droplet's volume and viscosity, we study the terminal velocities of droplets moving on a vertical fiber strand, their accelerations, and also the viscous energy dissipation rates during the early stages of motion. Table IV lists the physical properties of the oil–air system.

The computational setup for the present problem is shown in Fig. 12, where a droplet is placed onto a fiber. The same boundary conditions have been used for the sides of the domain, and for the fiber, a no-slip boundary condition is employed, as shown in Fig. 12. We consider droplets placed on smooth, chemically homogeneous and heterogeneous fibers of diameter $dv=200μm$⁠. Local insights are gained by evaluating the volumetric energy densities in particular, regarding viscous dissipation. To reduce the computational effort, the solver has also been enhanced to use the moving reference-frame technique as described in Sec. III C. Furthermore, an adaptive mesh refinement technique has been employed, where we use refinement criteria based on the order parameter so as to dynamically refine the mesh at the interface and inside the liquid phase.

The motion of droplets is governed by viscous and capillary forces (neglecting charge effects here), where the viscous force $Fv$ arises due to viscous fluid flow inside the droplet and the capillary forces $Fc$ originates due to the capillarity of the interface and the difference of its advancing and receding contact angles at the contact line. With this, one can analyze the equation of motion of a sliding droplet. On a vertical fiber strand, the motion of a droplet with volume Ω and density ρ_L is driven by the gravitational force $Fg=Fgêg$ of constant magnitude $Fg=ρLΩg$⁠, where $êg$ is the unit vector in direction of gravity and $g=‖g‖=9.81ms2$ is the magnitude of the gravitational acceleration. The component of the full equation of motion in direction of gravity is

Here, $Fv=−Fv·êg$ and $Fc=Fc·êg$ are components of the viscous and capillary forces, respectively. Note that the sign of the viscous force in Eq. (40) is taken negative by convention as this force is expected to act opposite to the direction of droplet motion and gravity. The sign of the capillary force F_c can be positive or negative, as the capillary force may act along or against the direction of gravity.

To gain further insight, we investigate the rate at which work is done on a fluid element in the droplet while moving on the fiber strand changing its shape and volume. In general,

where $σ=−pI+τ$ with $p=−13trσ$ denotes the total stress tensor, with p being the mechanical pressure. Since we assume isochoric flow, $∇·u≡0$⁠, we are left with the so-called dissipation function

Assuming a liquid of Newtonian fluid rheology, we can compute the total dissipation rate within the moving droplet, viz.

With this, the viscous force in (40) can be approximated as

where U_p is the droplet barycenter velocity.

The capillary force F_c in (40) is computed from

that is, the force the diffuse interface exerts on the substrate surface $∂Ω$ due to its capillarity.⁵⁸ 

Due to the small size of the fiber (and thus small contact line perimeter), the effect of capillary force is assumed to be negligible and the droplet motion is controlled by the effect of viscous forces. Based on this, Gilet et al.¹⁵ have proposed a correlation for the droplet's velocity as a function of the normalized volume, for a highly wetting liquid, viz.

where $Cν$ is a proportionality factor that accounts for the effect of surface tension. For a liquid spreading on a dry surface, Hoffman⁵⁹ showed that α = 15, which is used here to predict the droplets' terminal velocity. One should note that despite the complex dependence of W/X with $Ω/dv3$⁠, for sufficiently large ratios of $Ω/dv3$⁠, one finds that $W/X→0.5$⁠.⁵² 

Simulation results of the terminal droplet velocities are shown to be in very good agreement with the correlation (46), for all considered droplet volumes and viscosities, as seen in Fig. 13 (error bars indicate 10% deviation). For our simulations, the equilibrium contact angle has been set to $θ0=15°$⁠.

Figure 14 displays the temporal evolution of the droplets' barycenter velocities. Overshoots in velocities are observed at $t∼1ms$⁠, after which the velocities decrease and reach quasi-steady-state values.

The times it takes for droplets to reach their respective quasi-steady-state velocities have also been evaluated. Gilet et al. have proposed that the terminal velocity can be estimated using (46), where $ts−s∼Up/g$⁠. Our results show that the actual times to achieve steady-state velocities are higher: we observe that $ts−s∼(Up/g)n$⁠, with $n≈0.7$ resulting in a much better agreement with simulated results. This discrepancy may be caused by the way the droplet is initially placed onto the fiber in experiments compared to how it is initialized for simulations. Moreover, the equilibrium contact angle used in simulations may not the same as the one observed experimentally, which would affect the time to reach quasi-steady-state motion.

Gilet's assumption proposes that the droplet motion is governed by viscous forces, making it vital to study the viscous dissipation processes during the droplets' motion at the fiber strands in more detail. Several correlations to calculate the viscous force have been proposed in Refs. 15, 52, 60, and 61. We decide to use the one by Ref. 15 due to its simplicity

where U_p is the droplets' terminal velocity.

Figure 15 shows the viscous force as processed from simulation results over the droplets' velocities for kinematic viscosities of $10cSt$ and $20cSt$⁠, comparing it with the values obtained from (47), with the error bars indicating relative deviations of 10%.

A good agreement is observed. However, for the higher droplet velocities and for the lower viscosities, the simulation results start to deviate slightly from the expected values.

To identify the processes which govern dynamic wetting of droplets on fibers, the viscous dissipation rate has been compared between simulations. Here, the dissipation arising from contact line relaxation is neglected.⁶² Figure 16 shows the rate of dissipation for a droplet moving on a fiber strand with an equilibrium contact angle $θ0=15°$⁠, for a dynamic viscosity of $10cSt$ and $20cSt$⁠.

Initially, the droplets experience very rapid spreading over the fiber, driven by the difference between the equilibrium and its initial contact angle. In this initial stage, the viscous dissipation increases very rapidly due to the formation of a small cusp at the contact line, which has been also observed for the case of droplet spreading on flat substrates.⁶³ Here, we confirm the same phenomena extend to cylinder-shaped surfaces. One can also see that this behavior is quite similar for the various droplet volumes studied, only changing when the viscosity is also varied, which may indicate that the cusp formation is not significantly affected by the droplet volume as it happens over a very short span of time and it may be only influenced by the physical properties of the fluid and equilibrium contact angle.

After the initial spreading stage, the dissipation rate will increase up to a maximum and then decreases to a lower constant value, when the droplet reaches quasi-steady-state. We observe that the dissipation rate is larger the larger the droplets, which is expected due to their larger terminal velocities. Interestingly, simulations with the higher viscous oil (⁠$20cSt$⁠) exhibit lower dissipation rates. This is because viscosity influences the fluid velocity [cp. (46)], and a more viscous fluid will have a lower velocity in average, leading to lower dissipation rates at the bulk and wedge. Furthermore, this lower velocity promotes lower X/W ratios, which also contributes to the lower dissipation rate.

To verify or falsify the assumption, the capillary force was negligible when compared to the viscous force, we compute the capillary force according to (45). Figure 17 shows the capillary force that has been computed from our simulation results using (45), for all cases previously studied.

Looking at the capillary force, it is several orders of magnitude smaller than the viscous force, and thus, its contribution to the overall droplet motion is indeed negligible, as assumed by Gilet et al. Certainly, this is because the contact line perimeter is small, which is then reflected on the capillary force.

Since the perimeter of the contact line is kept the same, because the fiber diameter is also the same throughout simulations, the observed changes in the capillary force are only due to the difference in the advancing and receding contact angles, which increases with increasing droplet volume (and thus velocities), leading to higher capillary forces.⁶⁴ We also observe that lower viscosities in general lead to higher capillary forces.

The focus of this study is to investigate the influence that locally varying wettability and also stripe periodicity have on the droplets' velocity profile, dissipation rate, viscous force and capillary force (which is assumed to be negligible compared to the viscous force), and quasi-steady-state motion.

To investigate the motion of droplets on chemically heterogeneous fibers, the same computational domain is used as depicted in Fig. 12. However, instead of applying a homogeneous wetting boundary condition with one equilibrium contact angle, the fiber's substrate is considered patterned with stripes of alternating wettability, where the equilibrium contact angle values have been locally varied between $θ0=55°$ and $θ0=105°$⁠, respectively, as shown in Fig. 18. For the simulations, a droplet with volume $Ω=2.1μl$ and $νL=10cSt$ is chosen.

Figure 19 shows the velocity and dissipation rate profiles over time, as well as the gravitational, viscous, and capillary forces. Figure 20 shows the evolution of the dissipation rate density and velocity vectors in the moving reference frame for the case with $θ0∈{55°,105°}$⁠.

Comparing the dissipation rates depicted in Fig. 19 with the ones from Fig. 16, they are similar in the early stages of droplet motion—one observes a rapid increase in the dissipation rate due to the formation of a cusp during the initial spreading stage. The dissipation rate is biased when the transition from philic to phobic surfaces (and vice versa) happens. Here, the capillary force may contribute to both droplet acceleration and deceleration, depending on the local contact angle and droplet shape at the contact line.

Initially, the droplet is moving exclusively on a phobic stripe, where an overshoot of the velocity and dissipation rate can be seen, at around $t/τ=1$⁠. The change in the viscous force, being a function of the dissipation rate and velocity, will be similar. The capillary force at this stage is orders of magnitude smaller than the viscous force. Following this overshoot, the droplet velocity will decrease until $t/τ=4$⁠, where the front of the droplet passes through a philic stripe, forcing the front edge of the droplet to spread further. This can be seen as a rapid increase in both the velocity and dissipation rate as well as the viscous and capillary force. The magnitude of the capillary force is now much larger than before due to the deformation of the contact line. Since the front and rear of the droplet now have very different contact angles, contributions from the capillary force are larger. This increase in capillary force acts in the direction of the gravitational force, accelerating the droplet. At around $t/τ=6.5$⁠, the rear edge of the droplet finally approaches the philic region. Since the spreading occurs in the direction that is opposite to the droplet motion, the droplet decelerates, which also leads to a decrease in the dissipation rate and viscous force. The capillary force will also decrease since the “hysteresis” effect is now much smaller than before, reaching values that are once again much smaller than the viscous force. At this stage, the capillary force contributes to the deceleration of the droplet, as the trailing edge is pulled in the direction opposite of gravity. At around $t/τ=18.5$⁠, the front of the droplet will pass through a phobic region forcing the contact line region to bend, further decreasing the droplet's velocity and dissipation rate and increasing the capillary force due the increasingly different contact angles at the rear and front edges. The capillary force has contributed to the deceleration of the droplet.

This process repeats as the droplet is moving over the stripped patterns. Thus, its velocity will oscillate around its mean quasi-steady-state value. Thus, when compared to the homogeneous cases, the influence of the capillary force is indeed not negligible. Local deformations at the contact line as well as the difference in rear and front contact angles lead to larger values of the capillary force, which significantly affects the overall motion of the droplet.

To further investigate under which conditions the capillary force acts with or against the gravitational force, we decompose the contribution of the capillary force as

where $Fca=Fca·êg$ and $Fcr=Fcr·êg$ are the capillary force contributions at the advancing and receding edge of the droplet. Figure 21 shows the time evolution of both contributions to the capillary force. As the front edge passes to a philic region, the capillary force at the advancing and receding edges increases, which leads to a positive contribution to the capillary force in the direction of the droplet motion—accelerating the droplet. As the trailing edge of the droplet is also passing the philic region, starting at around $t/τ=6$⁠, the receding contact line will move in the direction opposite of the droplet motion, leading to the decrease in the receding capillary force. This leads to a decrease in the capillary force and thus the droplet decelerates, until it reaches a quasi-steady state. At $t/τ=18.5$⁠, the front edge passes to the phobic region. The decrease in the advancing capillary force decelerates the droplet once more.

We can infer that the value of the capillary force will be positive when the contact line accelerates in the direction of motion, making the droplet accelerate, and negative for a contact line that accelerates opposite to the direction of motion, decelerating the droplet. Moreover, one can see that since the droplet is initialized in the phobic region, the contribution of the capillary force at the front of the droplet is negative, and the capillary force at the rear is positive—this is the opposite in the philic region.

In order to also assess the accuracy of our force computation from DNS data, the right-hand side (RHS) and left-hand side (LHS) of (40) are calculated and compared—see Fig. 22. A good agreement between the terms in the RHS and LHS of (40) is observed, substantiating that the computation of the forces from DNS is reliable. A maximum deviation of about 38% occurs at the very early stages of motion, at around $t/τ=0.15$—this can be attributed to the initialization setup, where an un-physical pressure field is initially used (assumed to be uniform and with a value of zero) and so the diffuse interface must first relax, which leads to a larger error at the onset of motion. Despite this, one more notable deviation occurs at $t/τ=4$⁠, when the front edge of the droplet passes to the philic stripe. This deviation between the acceleration term and the sum of the forces can be attributed to the deviation from the tangent hyperbolic profile of the phase field: as the front edge of the droplet passes to the philic stripe, the profile of the diffuse interface is deteriorated. Therefore, the interfacial profile can no longer be assumed to be of equilibrium tangent–hyperbolic shape. However, our analysis using Eq. (45) is based on this assumption. This explains the deviation in Fig. 22 to be largest at the point where the contact line passes the wettability step on the fiber surface. The faster this transition happens the more the deterioration of the equilibrium profile, and the greater the deviation. This is substantiated by the later droplet transition, at time $t/τ=18$⁠, where the droplet velocity is lower and so is the deviation between LHS and RHS. Certainly, one could use the new relaxation model accounting for highly dynamic diffuse interface deformations,³⁶ but this is not within the scope of this work and shall be subject to future research.

Since the focus here is to study the influence that chemical heterogeneities have on the droplet motion and the forces acting on the droplet, another case with $θ0∈{65°,95°}$ has been investigated, for the same stripe length in order to disclose the influence of the wettability range. Moreover, the periodicity of the stripes has been doubled, to investigate its effect on motion. Figure 23 shows the sliding velocity and dissipation rate profiles during the simulation. Both the magnitude of the capillary force and viscous dissipation force profiles are shown.

From Fig. 23, it is evident that changing the wettability range influences sliding velocity, dissipation rate, and the forces acting on the droplet significantly. For the case with $θ0∈{65°,95°}$⁠, it can be seen that the profile of the plots resembles the ones for the case of $θ0∈{55°,105°}$⁠. Because the wettability range is now smaller when compared to the case with $θ0∈{55°,105°}$⁠, the spreading motion of the droplet as it passes through the philic stripe is not as intense and the increase in sliding velocity is not as large. Conversely, this also leads to lower contact angle hysteresis and so the capillary force is slightly smaller as well. Since the contact angle of the philic region is larger as well, the droplet velocity becomes larger for the case with $θ0∈{65°,95°}$⁠, which becomes noticeable at around $t/τ=8$⁠. Thus, the front edge of the droplet will also reach the phobic region faster for the case of $θ0∈{65°,95°}$⁠, at around $t/τ=17$⁠. After this, the capillary force will once again increase to the values seen before (since the difference in contact angles is almost the same), although slower than previously because the contact line is also moving slower.

Another case where the same wetting range is used, but the periodicity of the stripes is doubled (by halving their length) is investigated. Initially, the rear and front edges of the droplet are in different regions, and the capillary force is initially large and has the same value as for the case initially described. Hence, the initial overshoot of the sliding velocity is smaller. The front of the droplet quickly encounters a philic stripe, forcing the advancing edge of the droplet to spread in the same way as for the case with the regular periodicity, just at a shorter time. Around the time the droplet velocity starts decreasing, the front edge passes through a phobic region, a decrease in velocity, dissipation energy, and viscous force is observed. The capillary force steadily increases as the front edge slowly passes from the philic to the phobic region. Thus, the terminal droplet velocity is lower when compared to the other cases. Since the droplet velocity is so low, the viscous force will be smaller than the capillary force, which opposite to the assumption by Gilet et al.¹⁴ This showcases that the choice of surface chemistry can greatly affect the droplet's motion and shape and correlations must be applied with care.

We discuss the motion of droplets moving along vertically arranged, chemically homogeneous and heterogeneous fibers based on local insights gained by direct numerical simulations. For this, we have utilized a diffuse interface phase-field method implemented in our in-house solver, phaseFieldFoam, in OpenFOAM (FOAM-extend 4.0/4.1). For validation, we compare (inter alia) the equilibrium shapes of droplets spreading on fibers with results obtained from analytical expressions, for various volume-to-fiber-diameter ratios.

For the homogeneous fibers, we compare our simulation results regarding the terminal velocity and the viscous force with existing correlations. We observe that simulation results are in very good agreement with these correlations. Moreover, we study the motion until the droplets have reached terminal velocities, where we observe an overshoot in the droplet velocity in the beginning for all situations. From our computational analysis, we could improve a correlation by Gilet et al.¹⁵ for the time until the droplets reach their terminal velocities. We further investigate the effect that the capillary force has on the overall droplet motion, and confirm the assumption made by Gilet et al. for cases of chemically homogeneous surfaces—that its influence is negligible when compared to the viscous dissipation. Eventually, the dissipation rates during the droplets' motion have been studied in detail. We observe a rapid increase in the dissipation rates, due to the formation of a cusp at the contact line. This observation holds for all volumes. For the less viscous liquid, the dissipation rate is found to be larger than for the more viscous liquid.

For the heterogeneous fibers, we investigate the evolution of the dissipation rates, droplet velocities and viscous and capillary forces as the droplet slides through the phobic and philic regions. Here, surface chemistry has a significant influence on both the droplet's velocity and shape, as well as the forces that act on the droplet. Notably, we observed that the capillary force, which we showed to be negligible when compared to the viscous force for a droplet moving on a chemically homogeneous fiber, can reach values that are at the same order of magnitude, and in some cases larger than, as the viscous force. Thus, correlations must be applied with care in cases of chemically heterogeneous fiber substrates.

This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project No. 265191195-SFB 1194.

The authors have no conflicts to disclose.

Francisco Bodziony: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Martin Wörner: Conceptualization (equal); Formal analysis (equal); Visualization (equal); Writing – review & editing (equal). Holger Marschall: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.