Fibrotaxis: gradient-free, spontaneous and controllable droplet motion on soft solids

Abstract

Most passive droplet transport strategies rely on spatial variations of material properties to drive droplet motion, leading to gradient-based mechanisms with intrinsic length scales that limit the droplet velocity or the transport distance. Here, we propose droplet fibrotaxis, a novel mechanism that leverages an anisotropic fiber-reinforced deformable solid to achieve spontaneous and gradient-free droplet transport. Using high-fidelity simulations, we identify the fluid wettability, fiber orientation, anisotropy strength and elastocapillary number as critical parameters that enable controllable droplet velocity and long-range droplet transport. Our results highlight the potential of fibrotaxis as a droplet transport mechanism that can have a strong impact on self-cleaning surfaces, water harvesting and medical diagnostics.

---

1 Introduction

The control and movement of small droplets is critical in many engineering and science applications, including microfluidics,^1–3 microfabrication,⁴ medical diagnostics,⁵ and drug delivery.⁶ While controllable droplet manipulation on solid surfaces has been achieved using, e.g., gradients of material properties or external force fields, a gradient-free mechanism that passively induces controllable droplet motion would open opportunities for transformative advances. Here, we propose droplet fibrotaxis—a novel, gradient-free transport mechanism for small droplets that enables droplet motion with controllable direction and speed without relying on external force fields.

Over the last few years, notable advances in droplet manipulation have been achieved using electric fields,⁷ light,⁸ vibrations,^9,10 as well as gradients of curvature,¹¹ topography,¹² confinement,¹³ wettability,^14,15 and surface charge density.¹⁶ A particularly interesting example is droplet durotaxis,¹⁷ which is a passive droplet transport mechanism that uses stiffness gradients of the solid substrate. Durotactic droplet motion relies on the ability of the droplet to deform the solid through elastocapillary interactions that emerge from forces at the fluid–fluid interface.¹⁸ Experimental¹⁷ and computational¹⁹ studies have shown that small droplets placed on a soft silicone gel with a stiffness gradient move spontaneously without application of external forces, travelling up or down the gradient depending on the droplet's wettability. More recent research has proposed new droplet motion mechanisms based on elastocapillarity, including systems based on gradients of strain (tensotaxis),²⁰ gradients of pressure (bendotaxis)²¹ or substrate stretching.²²

Fibrotaxis, the droplet transport mechanism proposed herein, emerges from elastocapillary interactions with soft, homogeneous, anisotropic solids. One way to understand soft anisotropic solids is by conceptualizing them as an isotropic matrix with an infinite number of fibers continuously distributed. Fibrotaxis draws inspiration from cell migration in an extracellular matrix, a highly anisotropic network of collagen,^23,24 but spontaneous droplet motion on soft anisotropic solids has not been explored. Our data show that droplet fibrotaxis features unique advantages, such as, independence of external forces, controllable droplet direction and regulable velocity. More importantly, while gradient-based mechanisms are limited by inherent length scales that restrict the maximum transport distance or velocity, fibrotaxis is a gradient-free transport process. To elucidate the mechanism underpinning fibrotaxis we develop a high-fidelity fluid-structure interaction model that accounts for the coupled dynamics of a nonlinear, anisotropic hyperelastic solid, and two immiscible fluids with surface tension at their interface. Our results for droplets moving on fiber-reinforced gels unveil a complex interaction between the droplet's wettability and the orientation of the fibers that permits to control the droplet direction and velocity. Together, our findings indicate that fibrotaxis is a robust mechanism for droplet transport with unique features that can open new opportunities in microfluidics, medical diagnostics and energy generation.

2 Model of droplet fibrotaxis

2.1 Elastocapillary wetting

A liquid droplet resting on a flat and rigid solid surrounded by air will form a contact angle at equilibrium that satisfies the Young–Dupré equation γ_SL + γ_LA cos θ = γ_SA. Here, θ is the static contact angle between the liquid-air interface and the solid, while γ_SL, γ_LA and γ_SA are the surface tensions at the solid–liquid, liquid-air, and solid-air interfaces, respectively; see Fig. 1A. When the solid is deformable, however, the process is controlled by elastocapillary wetting,²⁵ and the Young–Dupré equation breaks down. Capillary forces at the fluid–fluid interface pull up the solid at the contact line, creating a ridge. Simultaneously, the overpressure in the droplet's interior generates a dimple in the solid located under the droplet. The combined effect of these two actions produces a rotation of the interfaces near the contact line that gives rise to an apparent contact angle α that differs from θ; see Fig. 1B. We hypothesize that when a droplet is placed on a soft anisotropic solid, such as a fiber-reinforced polymer, the apparent contact angles on the left (α_L) and right (α_R) triple points differ, leading to an imbalance of the horizontal forces that produces droplet motion. To test this hypothesis, we propose a three-dimensional multi-component fluid–solid interaction model. We consider a system comprising of two immiscible fluids, i.e., liquid droplet and air, in contact with a solid. We use the phase-field method²⁶ to model the fluids, by replacing the liquid-air interface with a thin transition region. We model the fiber-reinforced anisotropic solid using a transversely isotropic hyperelastic material.

Fig. 1 Schematic of a droplet wetting on (A) rigid solid and a (B) fiber-reinforced deformable solid. In (B), gray colored arrow depicts the orientation of fibers. This depiction does not indicate that we model the fibers discretely; rather, the model assumes that there are infinite fibers continuously distributed.

2.2 Governing equations of fluid mechanics

We describe the dynamics of the droplet and air using the Navier–Stokes–Cahn–Hilliard (NSCH) equations, a phase-field model for flow of two immiscible fluids.²⁷ The NSCH equations assume that the two fluids are incompressible and share the same velocity field. The governing equations in Eulerian coordinates are

∇·v = 0, (1a)

ρ(∂_tv + v·∇v) = ∇·σ^f + ρg, (1b)

Here, v is the fluid velocity; ρ is the fluid density; ∂_t denotes partial differentiation with respect to time; σ^f is the fluid Cauchy stress tensor; g is the acceleration of gravity; c ∈ [0,1] is the phase-field denoting the volume fraction of liquid; ε is the diffuse interface length scale; and M is the mobility coefficient. The fluid Cauchy stress tensor is defined as , where p is the pressure and ∇^s is the symmetrization of the spatial gradient operator, η is the dynamic viscosity of the mixture given by η = η₁c + η₂(1 − c), where η₁ and η₂ are the dynamic viscosities of the liquid and air, respectively. Because all of our simulations are performed at length scales [small script l] < 1 mm, gravity forces are negligible and we will ignore them in our calculations.

2.3 Governing equations of solid mechanics

We describe the dynamics of the solid using the linear momentum balance equation in Lagrangian coordinates,

ρ^s₀∂_t²u|_X = ∇_X · P, (2)

where ρ^s₀ is the mass density of the solid in the undeformed configuration, u is the solid displacement and P is the first Piola–Kirchhoff stress tensor. In eqn (2), we use the spatial derivative with respect to the material coordinates X in ∇_X. The subscript X in ∂_t²u|_X indicates that the second-order time derivative is taken by holding X fixed.

We model the solid as an anisotropic material. Our model is based on the theory of fiber-reinforced composites^28,29 and has been widely used for anisotropic gels. We consider a material composed of an isotropic matrix, in which one or several families of fibers are homogeneously distributed along specific directions. We do not model the fibers discretely; rather, the model assumes that there are infinite fibers continuously distributed. Here, we use a single family of fibers (i.e., all fibers have the same orientation) to represent a transversely isotropic material,³⁰ the simplest form of an anisotropic material. The fibers in the undeformed configuration are oriented along the direction defined by the unit vector d. We define the orientation tensor A = d ⊗ d. The strain energy density of the solid is given by ref. 30 as,

W = W_vol + W_iso + W_aniso, (3)

where W_vol and W_iso are the volumetric and isochoric isotropic contributions, respectively. Additionally, W_aniso is the anisotropic strain energy density contribution from the fibers. The isotropic contributions W_vol and W_iso are given by a compressible Neo–Hookean model³¹where κ and μ are the bulk and shear moduli of the solid, J is the Jacobian determinant of the deformation gradient F = I + ∇_Xu, I is the identity tensor, tr(·) is the trace operator, C = F^TF is the right Cauchy-Green deformation tensor and d is the number of spatial dimensions. The Young's modulus and the Poisson's ratio of the matrix are defined from μ and κ as and . The anisotropic strain energy density contribution is given by ref. 30 as,where H = tr(CA) is an invariant associated with the fibers, k₁ is a material parameter with dimensions of stress that defines the strength of anisotropy and k₂ is a dimensionless parameter. We assume that the fibers do not contribute to the material's mechanical response under compression. In eqn (2), P is computed from W as .

2.4 Dimensionless governing equations of the fluid–solid interaction problem

We solve the multi-component fluid–solid interaction problem by rescaling the units of length, time, and mass with respect to the droplet radius R, visco-capillary time-scale and ρR³, respectively. With this non-dimensionalization, the coupled system described by eqn (1) and (3) can be characterized by eight dimensionless numbers: (a) Ohnesorge number of the liquid which is the ratio of inertio-capillary to inertio-viscous time scale; (b) viscosity ratio of liquid to air ; (c) Cahn number which is the ratio of the diffuse interface length scale to the droplet radius; (d) Péclet number which is the ratio of advection to diffusion; (e) elastocapillary number , which quantifies the strength of elastocapillary effects; (f) dimensionless anisotropy strength ; (g) Poisson's ratio ν; and (h) k₂.

We use a superscript * to denote dimensionless quantities. For consistency, we also apply the superscript * to quantities that were already dimensionless before non-dimensionalization, such as c, F and C. We neglect the gravity forces and write the dimensionless governing equations for the fluid–solid interaction as,

∇*·v* = 0, (6a)

whereand . The dimensionless linear momentum balance for the solid iswhere
We solve the governing equations using a body-fitted fluid-structure algorithm. First, we rewrite the fluid equations in Arbitrary Lagrangian–Eulerian form.³² We then recast the governing equations in weak form and discretize them using Isogeometric Analysis.³³ We use a space of splines that is conforming at the fluid–solid interface and perform time integration using the generalized-α method.³⁴ We provide more details of our computational method in the Appendix.

3 Results

3.1 Fibrotaxis enables gradient-free, spontaneous droplet motion

We study fibrotaxis with droplet and solid properties similar to those used for durotaxis¹⁷ and elastocapillary wetting experiments.³⁵ A droplet surrounded by air is placed on a soft solid. Before deformation, the solid geometry is that of a rectangular prism of size 1000 μm × 50 μm × 500 μm; see Fig. 2A. The selected properties of the droplet and solid yield an elastocapillary length γ_LA/E equal to 9.2 μm, which is within the elastocapillary wetting range.³⁵ The elastic properties of the solid are similar to that of the spin-coated silicone gel used in ref. 17 while the solid's material properties defining the anisotropy strength are similar to those observed in soft biological tissues.³⁶ In our simulation, we use a transversely isotropic material with one family of fibers oriented in the direction d = cos βe₁ + sin βe₂, where e₁ and e₂ are the unit vectors in the x− and y− directions and β is the angle of the fibers with respect to e₁. The droplet, shown in a semi-transparent color is initially at rest (Fig. 2A), but the compliance and anisotropy of the solid substrate induce droplet motion immediately after the simulation starts. The droplet, shown in blue color is at time t = 1.38 s (Fig. 2A), and has moved parallel to the x-axis towards the right-hand side of the solid surface. As the droplet moves, it displaces the surrounding air, giving rise to a vortical structure, illustrated in Fig. 2A by streamlines. Fig. 2B shows a scaled view of the solid's deformed configuration at t = 1.38 s with the image colored by vertical solid displacements. We observe a depression of the solid under the droplet, a ridge at the triple line, and negligible displacements everywhere else. Because the solid fibers are parallel to the xy-plane and the droplet is a spherical cap, the solid deformation is symmetric with respect to the z-axis. However, the solid's anisotropy breaks the symmetry with respect to the x-axis resulting in a much larger solid deformation on the left-hand side—this symmetry breaking drives droplet motion.

Fig. 2 Three-dimensional simulation of droplet fibrotaxis. (A) Shows the positions of the droplet at different times. The initial position of the droplet is indicated by a semi-transparent spherical cap, while the droplet in blue is at time t = 1.38 s. The droplets are represented by isosurfaces of the phase field. We plot the streamlines of the fluid velocity along the mid-plane at t = 1.38 s, which are colored with the velocity magnitude. The black solid arrow depicts the fibers' orientation. Note that we do not model the fibers discretely; rather, the model assumes that there are infinite fibers continuously distributed. (B) Shows the solid deformation (magnified 3×) colored by the dimensionless vertical solid displacement u_y. We use a non-wetting droplet with θ = 120° and a radius of 160 μm. The fluids properties are γ_LA = 46 mN m⁻¹, ρ = 1260 kg m⁻³, η₁ = 1.412 Pa s and η₂ = 1.85 × 10⁻⁵ Pa s. The solid properties are ρ^s₀ = 1000 kg m⁻³, E = 5 kPa, ν = 0.25, k₁ = 40 kPa, k₂ = 7 and β = 60°. The initial solid geometry is a rectangular prism of size 1000 μm × 500 μm × 50 μm. The numerical parameters are M = 2 × 10⁻¹¹ m³ s kg⁻¹, ε = 20 μm and Δt = 100 μs. The values we choose correspond to the following dimensionless parameters: Oh = 14.67, [small eta, Greek, circumflex] = 7.6 × 10⁴, Cn = 0.125, Pe = 906.52, ζ = 0.058 and ϒ = 3.77 × 10⁴. The initial solid geometry is a rectangular prism of size 1000 μm × 50 μm × 500 μm. We perform computations in a box of size 1000 μm × 400 μm × 500 μm that includes the solid and fluid domains.

To better understand the mechanism of fibrotaxis, we perform an analogous two-dimensional simulation that permits a more detailed analysis; see Fig. 3. The droplet, initially at rest on the left-hand side of the solid surface, spontaneously moves to the right. Fig. 3A shows the system at a very early time—at this point the droplet's motion is negligible, but approximately 6 s later, it has reached the right end of the solid surface (Fig. 3B). The insets in Fig. 3A show that the right triple point has a much smaller displacement than the left one. This is a result of capillary forces pulling up the solid more parallel to the fibers on the right triple point, which leads to a stiffer response of the solid. This trend remains throughout the entire simulation as illustrated in Fig. 3C, which shows the vertical displacement of the solid along the fluid–solid interface. Thus, although our solid is homogeneous, the solid response is stiffer on one side of the droplet and at this point motion occurs due to the same mechanism as in durotaxis; see ref. 17. However, a critical advantage of fibrotaxis is that it is a gradient-free mechanism, unlike other transport processes that require gradients of stiffness (durotaxis),¹⁷ temperature, or wettability.¹⁴ The gradient-free mechanism is advantageous because it is devoid of spatial length scales that would set hard limits on the maximum transport length or maximum velocity. The fibrotaxis transport length is, in principle, limited only by dissipative mechanisms such as droplet evaporation or contact angle hysteresis, which we have neglected here. Additionally, the gradient-free nature of fibrotaxis allows simple control of the droplet velocity. As shown in Fig. 3D, the droplet velocity (v_d) remains steady for the length-scale considered in our simulations. Fig. 3D also shows that the variations in the left and right contact angles (α_L and α_R) are smaller than ∼1%. We have confirmed from our simulations that the total free energy (which is the analogue of entropy for an isothermal system) of the coupled fluid–solid system is consistently non-increasing over time.

Fig. 3 Mechanism of droplet fibrotaxis. (A) and (B) show droplet positions at two times. The insets in (A) show the solid deformation (10× magnified) and the orientation of the capillary forces acting at the wetting ridges. The black solid arrows in (A) and (B) depict the fibers' orientation. (C) Shows the vertical displacement of the solid u_y along the fluid–solid interface. (D) Shows the time evolution of droplet's velocity (v_d) and the apparent contact angles at the left (α_L) and right (α_R) contact lines. We plot the droplet velocity using the data of droplet position over time. The variation of apparent contact angles here is a moving average that filters out noise in the data resulting from inaccurate measurements of the contact line due to the diffuse interface. We use a droplet of radius 160 μm surrounded by air. The droplet is non-wetting with θ = 105°. The fluids properties are identical to those used in Fig. 2. For the solid, we choose ρ^s₀ = 1000 kg m⁻³, E = 5 kPa, ν = 0.25, k₁ = 50 kPa, k₂ = 7, β = 60° and a thickness of 50 μm. We also use M = 2 × 10⁻¹¹ m³ s kg⁻¹, ε = 20 μm and Δt = 25 μs. The values we choose correspond to Oh = 14.67, [small eta, Greek, circumflex] = 7.6 × 10⁴, Cn = 0.125, Pe = 906.52, ζ = 0.058 and ϒ = 4.71 × 10⁴. We perform our computations in a box of size 1000 μm × 500 μm.

To better illustrate the advantages of a gradient-free mechanism, we perform a durotaxis simulation of a glycerol droplet; see Fig. 4. The substrate has a stiffness gradient and is softer on the left-hand side. The droplet moves towards the right end of the solid. The absence of left-to-right symmetry can be seen in Fig. 4B, which shows the vertical solid displacements at the fluid–solid interface at two time instants. As the droplet moves towards the stiffer part, the displacements on both triple points become smaller, leading to a lower driving force that results in a decreasing droplet velocity (Fig. 4C). When the droplet reaches parts of the solid that are very stiff, the problem reverts to that of a droplet on a rigid solid and the droplet motion stops. Therefore, because the solid's stiffness on the left-hand side cannot be made arbitrarily small, the maximum transport distance or the maximum velocity are severely limited by constraints emanating from a gradient-based transport mechanism.

Fig. 4 Droplet motion driven by durotaxis. (A) Shows the droplet positions at two times. The initial position of the droplet is indicated by a black spherical cap, while the droplet in blue is at time t = 3.1 s. (B) Shows the vertical solid displacement u_y at the fluid–solid interface. (C) Shows the time evolution of the droplet velocity. We use a non-wetting droplet with θ = 120° and radius of 180 μm, and place it in contact with a heterogeneous solid of thickness 50 μm. The fluids properties are identical to those used in Fig. 2. We assume the solid is isotropic and its stiffness varies linearly from E = 5 kPa to E = 50 kPa. We also use ν = 0.25, M = 2 × 10⁻¹¹ m³ s kg⁻¹, ε = 25 μm and Δt = 25 μs. The values we choose correspond to Oh = 13.82, [small eta, Greek, circumflex] = 7.6 × 10⁴, Cn = 0.14, Pe = 1147.4, ζ ∈ [0.05,0.005] and ϒ = 0. We perform our computations in a box of size 1000 μm × 500 μm.

3.2 The interplay of wettability and fiber orientation regulates droplet fibrotaxis

The effect of wettability on droplet durotaxis was elucidated in ref. 19. When wetting droplets are placed on a substrate with a stiffness gradient, they migrate to the soft part. However, non-wetting droplets undergoing durotaxis move to the stiffer part of the substrate. This implies that if we keep the solid's mechanical properties unchanged, motion can be inverted going from a wetting to a non-wetting droplet. In fibrotaxis, wetting (respectively, non-wetting) droplets also move to the part of the substrate with a softer (respectively, stiffer) response. However, this implies that if we keep the solid surface unchanged, the droplet motion cannot be inverted by changing the fluid's wettability. This is illustrated in Fig. 5. The two fluids and the mechanical properties of the solid are identical to those in Fig. 3, but now the solid has more affinity to the fluid, which results in a wetting droplet. The simulation results show motion towards the right-hand side just like in the non-wetting droplet; see Fig. 5B. However, although the droplet moves towards the right in Fig. 3 and 5, the two processes are different. Indeed, two events are inverted in Fig. 5 with respect to Fig. 3 that have opposing effects and result in invariant motion direction: first, the displacements in the wetting case are larger on the right triple point while in the non-wetting case are larger on the left triple point; second, in the wetting case the droplet moves toward the side with a softer response while in the non-wetting case the droplet moves toward the side with stiffer response. In both the wetting and the non-wetting cases, the motion is in the direction of a lower apparent contact angle.

Fig. 5 (A) and (B) show droplet positions at two different times. The insets in (A) show the solid deformation (10× magnified) at the wetting ridges. The black solid arrows in (A) and (B) depict the fibers' orientation. (C) Shows the vertical displacement of the solid u_y along the fluid–solid interface. (D) Shows the time evolution of the apparent contact angles at the left (α_L) and right (α_R) contact lines. We use a droplet of radius 160 μm. The droplet is wetting with θ = 75°. The properties of fluid, solid, numerical and geometrical parameters are identical to those used in Fig. 3. The values we choose correspond to Oh = 14.67, [small eta, Greek, circumflex] = 7.6 × 10⁴, Cn = 0.125, Pe = 906.52, ζ = 0.058 and ϒ = 4.71 × 10⁴. We perform our computations in a box of size 1000 μm × 500 μm.

3.3 Control of droplet fibrotaxis

The proposed mechanism of droplet fibrotaxis and the roles of fiber orientation and wettability that we identified suggest possible control strategies for droplet motion that include directionality and speed. We illustrate this possibility by constructing a droplet velocity diagram on the θ − β phase space. Fig. 6 shows the results of 77 high-fidelity simulations (represented by circles in the plot), each corresponding to a θ − β pair. All simulations were performed using the same value for the droplet radius and the same fluid properties. We positioned the droplet at the center of the domain and assumed that its shape corresponds to that of a circular sector at contact angle θ with the solid—this would represent a quasi-equilibrium solution for a rigid solid. Our simulations reveal that the magnitude of the droplet velocity is maximum when θ ≈ β or θ ≈ β + 90°, which represent cases for which the solid's response is stiffest and softest, respectively, on one of the triple points. The maximum velocity magnitudes are not exactly achieved at θ = β and θ = β + 90° because the orientation of the fibers changes with deformation and can only be defined unambiguously in the undeformed configuration. When β = 0°, 90° or 180°, the droplet velocity is negligible because the solid's response is the same on both triple points for any contact angle θ. For the same reason, a neutrally-wetting droplet (θ = 90°) shows negligible velocity for any fiber orientation. Fig. 6 also indicates that the velocity magnitude is symmetric with respect to β = 90°, but the motion direction is flipped. Thus, for a fixed θ, the fiber orientations β and 180 − 2β lead to the same velocity magnitude, but opposite motion directions. Another important conclusion that emerges from Fig. 6 is that, for a given fiber orientation, the direction of droplet motion cannot be reversed changing wettability only. Overall, the data show that fibrotaxis permits to achieve a wide range of droplet velocities with controllable direction in microfluidic applications.

Fig. 6 Control of droplet fibrotaxis. (A) Shows the dependence of droplet velocity on fiber angle β and wettability θ. Each circle in the plot indicates a numerical simulation for a θ − β pair. We use a droplet with a radius of 180 μm surrounded by air. The material parameters of the solid are E = 5 kPa, ν = 0.25, k₁ = 50 kPa and k₂ = 7. (B) Shows the variation of droplet velocity with the strength of anisotropy measured by the parameter k₁ (see eqn (5)). We use a non-wetting droplet with θ = 120° and radius of 160 μm, and place it in contact with air. The material parameters of the solid are E = 10 kPa, ν = 0.25, β = 40° and k₂ = 7. (C) Shows the variation of droplet velocity with the elastocapillary number ζ. We use θ = 120° and β = 60°. The material parameters of the solid are E ∈ [2.5,100] kPa, ν = 0.25, k₁ = 50 kPa and k₂ = 7. (D) Shows the solid deformation (10× magnified) at the left and right wetting ridges for two different values of ζ from (C): 0.068 and 0.0026. The black solid arrows in (D) depict the fibers' orientation. For (A)–(C), we use fluids properties identical to those used in Fig. 2. Additionally, we use a solid thickness of 50 μm, M = 2 × 10⁻¹¹ m³ s kg⁻¹, ε = 25 μm and Δt = 25 μs. We perform our computations in a box of size 1000 μm × 500 μm.

Fig. 6B shows that the anisotropy strength k₁ (see eqn (5)) can be used as a third independent control parameter for the droplet velocity. The plot indicates that the droplet remains immobile for isotropic solids (k₁ = 0) and moves faster as we increase the strength of anisotropy. At small values of k₁, the droplet velocity increases sharply, for example, by about ∼200% when k₁ is increased from ∼1 to ∼6 kPa. For large values of k₁, the droplet's velocity plateaus at about ∼280 μm s⁻¹. Although we use fixed values of θ and β in Fig. 6B, we have verified from our numerical simulations that a similar trend in the variation of droplet velocity with k₁ remains across a wide range of θ and β values considered in Fig. 6A.

Fig. 6C shows that the elastocapillary number ζ serves as a fourth independent control parameter for the droplet velocity. To plot Fig. 6C, we vary E and keep θ, β and k₁ fixed. The plot indicates that the droplet velocity increases as we increase the elastocapillary number. At high values of elastocapillary numbers, the droplet velocity decreases sharply, for example by ∼46% when the elastocapillary number is lowered from 0.102 to 0.068. At very low elastocapillary numbers, the droplet remains nearly stationary. The trend in Fig. 6C is because the droplet velocity is directly proportional to the magnitude of the relative solid displacement at the left and right triple points, a behavior that has been previously observed in durotaxis^17,19—this is illustrated in Fig. 6D. In Fig. 6D, at an elastocapillary number of 0.068, the solid displacement at the left triple point is about ∼1300% more than that on the right triple point, resulting in a high droplet velocity; but, at an elastocapillary number of 0.0026, the solid displacement at the left triple point is only ∼40% more than that on the right triple point, leading to a low droplet velocity. We can also conclude from Fig. 6C that when the elastocapillary number is theoretically zero, the droplet will remain stationary because the problem corresponds to that of a droplet on a rigid solid.

4 Conclusions

The ability to move, mix and process small droplets without external energy supply can open many new opportunities in science, engineering and industry, including, but not limited to, (i) new ways to perform multiple or sequential processes in microfluidic and lab-on-a-chip devices, (ii) innovative mechanisms for enhanced heat transfer that ensure that fresh coolant is continually in contact with the hot surface, (iii) better designs for self-cleaning surfaces that use droplet motion to remove dirt and contaminants, (iv) efficient methods to guide condensation and harvest water from air in arid regions, and (v) pump-free and smaller devices for biofluid manipulation in medical diagnostics.

This paper proposed fibrotaxis—a simple and robust mechanism for transporting droplets on fiber-reinforced deformable solids. Notably, fibrotaxis is a gradient-free transport mechanism that sets it apart from most droplet motion techniques. Using high-fidelity numerical simulations, we demonstrated that fibrotaxis exhibits three key features essential for efficient droplet transport: controllable droplet velocity and droplet direction, long-distance droplet transport, and spontaneity of droplet motion. We show that adjusting four control parameters—fiber angle, wettability, anisotropy strength, and elastocapillary number allows us to control droplet velocity. By constructing a velocity diagram, we additionally show that fiber angle and wettability can control droplet direction. Our results indicate that fibrotaxis is capable of moving droplets of ∼200 μm radius at speeds of ∼200 μm s⁻¹ using anisotropic solids that can be fabricated through reinforcement of standard silicone gels with silica or polymeric³⁷ fibers. This velocity range is appropriate for most microfluidics, diagnostics, water harvesting and heat exchange applications, but higher velocities could be achieved, for example, by coating the solid surface with a thin layer of lubricant or by creating surface textures in the solid.³⁸ These methods reduce viscous dissipation, which facilitates faster droplet motion. Anisotropy is ubiquitous in soft gels. Even those gels which are isotropic prior to deformation can exhibit anisotropy under sufficient deformation due to their nonlinear response or structural changes in the material.^39,40

We investigated fibrotaxis in its simplest form, which emerges from the use of a single family of rectilinear fibers embedded in a homogeneous solid matrix. Our theory, extends to other forms of anisotropy, including multiple fiber families with different orientations (see the ESI† for more details) and curvilinear fibers. These forms of anisotropy could facilitate the design of predefined droplet trajectories and velocities. Our computational insights into fibrotaxis will offer valuable guidance for designing future experiments to control droplet motion on anisotropic solids. Fibrotaxis can be combined with other forms of droplet control, such as wettability patterns^41,42 or transport mechanisms based on the use of magnetic fields.⁴³ Future efforts could also focus on enhancing the present model by incorporating a nonlinear poro-elastic solid, offering a more realistic depiction of anisotropic gels.

Author contributions

Sthavishtha R. Bhopalam: conceptualization (supporting); data curation (lead); formal analysis (lead); investigation (equal); methodology (equal); software (lead); validation (lead); visualization (lead); writing – original draft (equal); writing – review & editing (equal). Jesus Bueno: software (supporting); validation (supporting). Hector Gomez: conceptualization (lead); formal analysis (supporting); funding acquisition (lead); investigation (equal); methodology (equal); project administration (lead); resources (lead); supervision (lead); writing – original draft (equal); writing – review & editing (equal). All authors reviewed and approved the final version of the manuscript.

Data availability

Data for this article are available at https://doi.org/10.4231/7AXT-KP80.

Conflicts of interest

There are no conflicts to declare.

A Appendix

A.1 Computational method

Here, we briefly present the key ingredients of our computational method, which has been successfully applied and validated through numerical experiments on various elastocapillary wetting problems.^20,44,45 In what follows, we use the dimensional form of the governing equations from eqn (1) and (2) to derive the variational formulation for the FSI problem.

Arbitrary Lagrangian–Eulerian formulation of fluid mechanics equations. We solve the governing equations using a body-fitted fluid-structure algorithm similar to that used by the authors in ref. 44 and 45. The fluid equations (see eqn (1)) are re-written in Arbitrary Lagrangian–Eulerian form,³² which are given by,

∇·v = 0, (9a)

where v̂ is the fluid mesh velocity, obtained as a mapping between the Eulerian and referential coordinates x̂. The subscript x̂ in eqn (9b) and (9c) indicates that the time derivative is taken by holding x̂ fixed. For the body-fitted fluid-structure algorithm, we update the fluid mesh by solving successive fictitious linear elasticity problems.⁴⁶ Because density variations do not make a significant difference at this scale, we assumed the fluid density to be constant.

Boundary conditions and fluid–solid interface conditions for the FSI problem. In our numerical simulations, we assume that the initial shape of the droplet corresponds to that of a circular sector at contact angle θ with the solid—this would represent a quasi-equilibrium solution for a rigid solid. We perform our simulations on prismatic domains. We position the outer boundary of our computational domain far away from the droplet. We impose the following boundary conditions at the lateral fluid boundaries: (a) zero normal velocity, (b) zero tangential traction and (c) zero phase-field flux. However, at the top fluid boundary, we set all fluid velocity components to zero. On the lateral solid boundaries, we impose zero normal displacement and zero tangential traction. But, at the bottom solid boundary, we set all solid displacement components to zero.

For the coupled fluid-structure interaction problem, we impose the following coupling conditions at the fluid–solid interface: (a) wettability condition ∇c·n^f = |∇c|cos θ, (b) kinematic condition and (c) traction balance σ^fn^f = σ^sn^f where, n^f is the unit normal vector at the fluid–solid interface pointing in the direction from fluid to solid, and σ^s = J⁻¹PF^T is the solid Cauchy stress tensor. For weak imposition of the wetting boundary condition, we split eqn (9c) into two second-order equations by defining the chemical potential .

Variational formulation of the FSI problem. We multiply the governing equations with weight functions and integrate them over the computational domains in the spatial (Ω^f_t) and referential (Ω^s₀) configurations to obtain the weak form. In our FSI problem, we consider the material and referential configurations to be the same. We denote [scr V, script letter V]^f, [scr V, script letter V]^s and [scr V, script letter V]^m to be the trial function spaces for the fluid, solid and mesh motion problems. We also denote [scr W, script letter W]^f, [scr W, script letter W]^s and [scr W, script letter W]^m to be the weight function spaces in a similar fashion. We state the weak form as follows: find {p,v,c,μ} ∈ [scr V, script letter V]^f, u ∈ [scr V, script letter V]^s and u^m ∈ [scr V, script letter V]^m such that ∀ {w^p,w^v,w^c,w^μ} ∈ [scr W, script letter W]^f, ∀w^s ∈ [scr W, script letter W]^s and ∀w^m ∈ [scr W, script letter W]^m, B^f({w^p,w^v,w^c,w^μ},{p,v,c,μ};v̂) + B^s(w^s,u) + B^m(w^m,u^m) = 0, whereΓ_t^sf is the fluid–solid interface in the spatial configuration,B^m(w^m,u^m) is the variational formulation of the fluid mesh motion and u^m is the mesh displacement. For numerical implementation of eqn (10), we adopt the variational multiscale method (VMS).⁴⁷

Spatial discretization of the FSI problem. We spatially discretize the weak form using Isogeometric Analysis (IGA).³³ We use IGA because of its proven success in computational phase-field modeling and its ability to use basis functions with controllable inter-element continuity. To derive the semi-discrete variational formulation, we substitute the trial and weight function spaces defined earlier with finite-dimensional trial and weight function spaces such that [scr W, script letter W]_h^f ⊂ [scr W, script letter W]^f, [scr W, script letter W]_h^s ⊂ [scr W, script letter W]^s, [scr W, script letter W]_h^m ⊂ [scr W, script letter W]^m, [scr V, script letter V]_h^f ⊂ [scr V, script letter V]^f, [scr V, script letter V]_h^s ⊂ [scr V, script letter V]^s and [scr V, script letter V]_h^m ⊂ [scr V, script letter V]^m. The solution variables and corresponding weight functions are defined aswhere p_A, w^p_A, u_A, w^s_A are the control variables, A is a control variable index, I_f and I_s are the index sets of the fluid and solid control variables, respectively. Additionally, N̂_A (and N_A) represent spline basis functions defined on referential (and spatial) configuration of the combined fluid and solid. In our simulations, we use quadratic spline basis functions with [capital script C]¹ – inter-element continuity everywhere except along four parametric lines that enclose the solid domain and include the fluid–solid interface.

Time integration and solution strategy. We perform time-integration of our equations using the generalized-α method³⁴ and solve the coupled FSI problem using a quasi-direct solution strategy.⁴⁸ We use the Newton-Raphson method with a backtracking linesearch for solving the nonlinear system of equations while we solve the linear solver in each Newton iteration using Generalized Minimal Residual method (GMRES) with a Jacobi preconditioner. We define a residual vector for each solution variable by substituting the discrete solutions into the variational formulation (see eqn (10) and (11)). We monitor the convergence of the nonlinear solver by ensuring that at least one of the following conditions is met: (a) the relative tolerance of each residual vector is smaller than 10⁻³, (b) the absolute tolerance of each residual vector is smaller than 5 × 10⁻⁶ and (c) the norm of the change in the solution between each Newton step is smaller than 10⁻³. Also, we set the linear iterative solver to converge until the relative reduction in the preconditioned residual norm reaches 10⁻⁵. We develop our code using the open-source frameworks, PetIGA⁴⁹ and PETSc.⁵⁰

A.2 Validation of computational model

Here, we present two validation examples performed using our FSI model and computational method. The first example demonstrates the capability of our model to quantitatively capture the elastocapillary FSI phenomenon. The second example shows that our model can quantitatively reproduce the contact line velocity observed in capillary tube experiments.

Validation of FSI model. We verify that our FSI code and formulation produce quantitatively accurate results by simulating the static wetting experiment of a glycerol droplet on a thin sheet of silicone gel.³⁵ By setting k₁ = 0 in our model (see eqn (5)), we perform numerical simulations in a two-dimensional box such that the length of the box is twice the width. We position the outer boundary of our computational domain far from the droplet.

In our simulations, we initially position the droplet at the center of the domain such that the droplet-air interface intersects the undeformed solid. We assume that the initial shape of the droplet is a circular sector that makes a contact angle of 96.24° with the solid at equilibrium. We impose zero tangential traction and zero phase-field flux along the boundaries of our domain. At the lateral boundaries of the domain, we impose zero velocity in the normal direction. At the bottom boundary, we impose zero normal displacements for the solid. At the top boundary, we impose zero velocity along both the normal and tangential directions.

Fig. 7 shows our simulation results at t = 300 μs (assumed to be the steady state). Fig. 7A shows the formation of two wetting ridges at the fluid–solid interface. In Fig. 7B, we plot the vertical displacement of the fluid–solid interface along the horizontal direction measured from the center of the droplet x_c. From Fig. 7B, it is evident that the data from our simulations are in excellent agreement with experiments.³⁵

Fig. 7 Static wetting of glycerol droplet on silicone gel. (A) Shows how the droplet deforms the solid at steady state. (B) Shows the vertical displacement of the solid along the horizontal direction at t = 300 μs. The fluids properties are γ_LA = 46 mN m⁻¹, ρ = 1260 kg m⁻³, η₁ = 1.4 Pa s, η₂ = 1.4 Pa s. The material properties of the solid are E = 3.0 kPa, ν = 0.49 and ρ^s₀ = 1260 kg m⁻³. The numerical parameters are M = 5 × 10⁻⁸ m³ s kg⁻¹, and ε = 3 μm. We perform computations in a two-dimensional box of size 1000 μm × 500 μm and spatially discretize it with a mesh of 400 × 200 [capital script C]¹ – quadratic elements.

Validation of contact line dynamics. To validate the contact line dynamics captured by our computational method, we simulate the steady-state displacement of glycerol in a rigid capillary tube when air is injected from one side of the tube.⁵¹ Since this example serves as a validation case for the fluid mechanics formulation alone, we do not solve the equations of solid dynamics and mesh motion.

Fig. 8A shows the computational setup of our capillary tube. We solve the governing equations given by eqn (1) (or eqn (6) for the dimensionless governing equations) in axi-symmetric cylindrical coordinates. We perform numerical simulations in a two-dimensional box (see Fig. 8A) such that the length of the box is five times larger than the width of the box. In our simulations, we initially assume that the shape of the air–glycerol interface is flat. We inject air into the inlet boundary of the capillary tube such that the axial velocity is given by , where R_c is the radius of the capillary tube and V̄ is the average velocity of air at the inlet boundary. At the inlet boundary of our computational domain, we additionally impose a zero radial velocity. At the outlet boundary of our domain, we specify zero radial velocity, zero tangential traction, and zero phase-field flux. At the symmetric boundary, we impose the following symmetry boundary conditions: zero radial velocity, zero tangential traction, and zero phase-field flux. At the solid wall, we impose a no-penetration condition for the velocity, i.e., v = 0 and a dynamic wettability condition , where D_w is the dynamic wall mobility that accounts for solid wall relaxation and bulk diffusion near the solid wall while ϑ is the contact angle made by the glycerol–air interface with the solid.

Fig. 8 Displacement of glycerol by air in a capillary tube. (A) Shows the computational domain. (B) Shows the glycerol-air interface profile for three different Ca obtained using experiments⁵¹ and present simulations. The radius of capillary tube R_c is 375 μm. The fluids properties are γ_LA = 65 mN m⁻¹, ρ = 1260 kg m⁻³, η₁ = 1.4 Pa s, η₂ = 1.85 × 10⁻⁵ Pa s and ϑ = 68°. The numerical parameters are M = 5 × 10⁻⁹ m³ s kg⁻¹, D_w = 3.327 kg m⁻¹ s⁻¹ and ε = 45 μm. We perform computations in a two-dimensional box of size 1875 μm × 375 μm and spatially discretize it with a mesh of 96 × 480 [capital script C]¹ – quadratic elements.

Fig. 8B shows the present simulation results of the air–glycerol interface profile in the capillary tube for different , which are in excellent agreement with the experiments. It is evident from Fig. 8B that as Ca increases, the air–glycerol interface deforms significantly and bends closer to the symmetry axis. For the range of Ca captured in Fig. 8B, experiments⁵¹ indicate that Ca_cl = Ca and Ca_tip = Ca, where Ca_cl and Ca_tip are computed from the temporal data of contact line positions at the solid wall and symmetry axis respectively (see Fig. 8A). Table 1 shows the data for Ca_cl and Ca_tip from our simulations. The deviation of the values of Ca_cl and Ca_tip between simulations and experiments is less than 2%.

Table 1 Data of Ca_cl and Ca_tip predicted from simulations of air–glycerol displacement in a capillary tube

Ca	Cacl	Catip
3.00 × 10^−3	3.01 × 10^−3	3.07 × 10^−3
6.00 × 10^−3	6.02 × 10^−3	6.04 × 10^−3
1.2 × 10^−2	1.19 × 10^−2	1.21 × 10^−2

---

Acknowledgements

S. R. B. thanks Mario de Lucio for valuable discussions on the solid model. This research was supported by the National Science Foundation (Award no. CBET 2012242). This work uses the Bridges-2 system at the Pittsburgh Supercomputing Center (PSC) through allocation MCH220014 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296. The opinions, findings, and conclusions, or recommendations expressed are those of the authors and do not necessarily reflect the views of the National Science Foundation.

References

1. M. Joanicot and A. Ajdari, Science, 2005, 309, 887–888.
2. H. A. Stone, A. D. Stroock and A. Ajdari, Annu. Rev. Fluid Mech., 2004, 36, 381–411.
3. R. Seemann, M. Brinkmann, T. Pfohl and S. Herminghaus, Rep. Prog. Phys., 2011, 75, 016601.
4. M. Srinivasarao, D. Collings, A. Philips and S. Patel, Science, 2001, 292, 79–83.
5. M. Serra, D. Ferraro, I. Pereiro, J.-L. Viovy and S. Descroix, Lab Chip, 2017, 17, 3979–3999.
6. M. Li, M. Brinkmann, I. Pagonabarraga, R. Seemann and J.-B. Fleury, Commun. Phys., 2018, 1, 23.
7. J. Li, N. S. Ha, T. L. Liu, R. M. van Dam and C.-J. C. J. Kim, Nature, 2019, 572, 507–510.
8. W. Li, X. Tang and L. Wang, Sci. Adv., 2020, 6, eabc1693.
9. K. John and U. Thiele, Phys. Rev. Lett., 2010, 104, 107801.
10. P. Brunet, J. Eggers and R. Deegan, Phys. Rev. Lett., 2007, 99, 144501.
11. C. Lv, C. Chen, Y.-C. Chuang, F.-G. Tseng, Y. Yin, F. Grey and Q. Zheng, Phys. Rev. Lett., 2014, 113, 026101.
12. J. Li, Y. Hou, Y. Liu, C. Hao, M. Li, M. K. Chaudhury, S. Yao and Z. Wang, Nat. Phys., 2016, 12, 606–612.
13. R. Dangla, S. C. Kayi and C. N. Baroud, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 853–858.
14. S. Daniel, M. K. Chaudhury and J. C. Chen, Science, 2001, 291, 633–636.
15. U. Thiele, K. John and M. Bär, Phys. Rev. Lett., 2004, 93, 027802.
16. Q. Sun, D. Wang, Y. Li, J. Zhang, S. Ye, J. Cui, L. Chen, Z. Wang, H.-J. Butt and D. Vollmer, et al. , Nat. Mater., 2019, 18, 936–941.
17. R. W. Style, Y. Che, S. J. Park, B. M. Weon, J. H. Je, C. Hyland, G. K. German, M. P. Power, L. A. Wilen, J. S. Wettlaufer and E. R. Dufresne, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 12541–12544.
18. J. Bico, É. Reyssat and B. Roman, Annu. Rev. Fluid Mech., 2018, 50, 629–659.
19. J. Bueno, Y. Bazilevs, R. Juanes and H. Gomez, Soft Matter, 2018, 14, 1417–1426.
20. J. Bueno, Y. Bazilevs, R. Juanes and H. Gomez, Extreme Mech. Lett., 2017, 13, 10–16.
21. A. T. Bradley, F. Box, I. J. Hewitt and D. Vella, Phys. Rev. Lett., 2019, 122, 074503.
22. K. Smith-Mannschott, Q. Xu, S. Heyden, N. Bain, J. H. Snoeijer, E. R. Dufresne and R. W. Style, Phys. Rev. Lett., 2021, 126, 158004.
23. A. Saez, M. Ghibaudo, A. Buguin, P. Silberzan and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 8281–8286.
24. R. S. Sopher, H. Tokash, S. Natan, M. Sharabi, O. Shelah, O. Tchaicheeyan and A. Lesman, Biophys. J., 2018, 115, 1357–1370.
25. R. W. Style, A. Jagota, C.-Y. Hui and E. R. Dufresne, Annu. Rev. Condens. Matter Phys., 2017, 8, 99–118.
26. D. M. Anderson, G. B. McFadden and A. A. Wheeler, Annu. Rev. Fluid Mech., 1998, 30, 139–165.
27. D. Jacqmin, J. Comput. Phys., 1999, 155, 96–127.
28. G. A. Holzapfel, T. C. Gasser and R. W. Ogden, J. Elasticity Phys. Sci. Solids, 2000, 61, 1–48.
29. T. C. Gasser, R. W. Ogden and G. A. Holzapfel, J. R. Soc., Interface, 2006, 3, 15–35.
30. D. R. Nolan, A. L. Gower, M. Destrade, R. W. Ogden and J. McGarry, J. Mech. Behav. Biomed. Mater., 2014, 39, 48–60.
31. J. C. Simo and T. J. R. Hughes, Computational Inelasticity, Springer, New York, 1998, vol. 7.
32. J. Donea, A. Huerta, J.-P. Ponthot and A. Rodrguez-Ferran, Arbitrary Lagrangian–Eulerian Methods, in Encyclopedia of Computational Mechanics, Fluids, 2004, vol. 3, ch. 14.
33. T. J. R. Hughes, J. A. Cottrell and Y. Bazilevs, Comput. Methods Appl. Mech. Eng., 2005, 194, 4135–4195.
34. K. E. Jansen, C. H. Whiting and G. M. Hulbert, Comput. Methods Appl. Mech. Eng., 2000, 190, 305–319.
35. R. W. Style, R. Boltyanskiy, Y. Che, J. S. Wettlaufer, L. A. Wilen and E. R. Dufresne, Phys. Rev. Lett., 2013, 110, 066103.
36. G. A. Holzapfel, T. C. Gasser and M. Stadler, Eur. J. Mech. A/Solids, 2002, 21, 441–463.
37. Z. Wang, C. Xiang, X. Yao, P. Le Floch, J. Mendez and Z. Suo, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 5967–5972.
38. M. Coux and J. M. Kolinski, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 32285–32292.
39. D. Vader, A. Kabla, D. Weitz and L. Mahadevan, PLoS One, 2009, 4, e5902.
40. R. Baker and C. S. Desai, Int. J. Numer. Anal. Methods Geomech., 1984, 8, 167–185.
41. L. Sun, F. Bian, Y. Wang, Y. Wang, X. Zhang and Y. Zhao, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 4527–4532.
42. A. F. Demirörs, S. Aykut, S. Ganzeboom, Y. A. Meier and E. Poloni, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2111291118.
43. A. Li, H. Li, Z. Li, Z. Zhao, K. Li, M. Li and Y. Song, Sci. Adv., 2020, 6, eaay5808.
44. S. R. Bhopalam, J. Bueno and H. Gomez, Comput. Methods Appl. Mech. Eng., 2022, 400, 115507.
45. J. Bueno, H. Casquero, Y. Bazilevs and H. Gomez, Meccanica, 2018, 53, 1221–1237.
46. T. Wick, Comput. Struct., 2011, 89, 1456–1467.
47. T. J. R. Hughes, G. Scovazzi and L. P. Franca, Multiscale and stabilized methods, in Encyclopedia of Computational Mechanics Second Edition, Wiley Online Library, 2017, pp. 1–64.
48. Y. Bazilevs, V. M. Calo, T. J. R. Hughes and Y. Zhang, Comput. Mech., 2008, 43, 3–37.
49. L. Dalcin, N. Collier, P. Vignal, A. Côrtes and V. Calo, Comput. Methods Appl. Mech. Eng., 2016, 308, 151–181.
50. S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang and H. Zhang, PETSc Web page, https://www.mcs.anl.gov/petsc, 2019.
51. B. Zhao, A. Alizadeh Pahlavan, L. Cueto-Felgueroso and R. Juanes, Phys. Rev. Lett., 2018, 120, 084501.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm01022a

---