The Plateau-Rayleigh instability in solids is a simple phase separation

TL;DR

This paper explains the Plateau-Rayleigh instability in elastic solids as a phase separation driven by surface tension, providing a theoretical model, finite-element verification, and analysis of interface properties.

Contribution

It introduces a simple energy-based model for the instability, predicts phase separation and hysteresis, and derives the interface length scale near the critical surface tension.

Findings

Instability occurs when surface tension exceeds a critical value.

The cylinder phase separates into segments with different stretches.

The interface length diverges near the critical surface tension.

Abstract

A long elastic cylinder, radius $a$ and shear-modulus $\mu$, becomes unstable given sufficient surface tension $\gamma$. We show this instability can be simply understood by considering the energy, $E(\lambda)$, of such a cylinder subject to a homogenous longitudinal stretch $\lambda$. Although $E(\lambda)$ has a unique minimum, if surface tension is sufficient ($\Gamma\equiv\gamma/(a\mu)>\sqrt{32}$) it looses convexity in a finite region. We use a Maxwell construction to show that, if stretched into this region, the cylinder will phase separate into two segments with different stretches $\lambda_1$ and $\lambda_2$. Our model thus explains why the instability has infinite wavelength, and allows us to calculate the instability's sub-critical hysteresis loop (as a function of imposed stretch), showing that instability proceeds with constant amplitude and at constant (positive) tension as the cylinder is stretched between $\lambda_1$ and $\lambda_2$. We use full nonlinear finite-element calculations to verify these predictions, and to characterize the interface between the two phases. Near $\Gamma=\sqrt{32}$ the length of such an interface diverges introducing a new length-scale and allowing us to construct a 1-D effective theory. This treatment yields an analytic expression for the interface itself, revealing its characteristic length grows as $l_{wall}\sim a/\sqrt{\Gamma-\sqrt{32}}$.

Equations

$$E=\gamma A+\int \frac{1}{2}\mu \left[\text{ Tr}\left(F\cdot F^T\right)\right]dV,$$

$$E(\lambda )$$

$$T=\frac{1}{L}\frac{dE}{d\lambda }=\pi \mu a^2\left(\frac{\Gamma }{\sqrt{\lambda }}+\lambda -\frac{1}{{\lambda }^2}\right),$$

$${\lambda }_{\pm }={\left[\frac{1}{4}\left(\Gamma \pm \sqrt{{\Gamma }^2-32}\right)\right]}^{2/3},$$

$$E^{\prime }({\lambda }_1)$$

$$E({\lambda }_1)+E^{\prime }({\lambda }_1)({\lambda }_2-{\lambda }_1)$$

$$A=\frac{a}{\sqrt{{\lambda }_1}}-\frac{a}{\sqrt{{\lambda }_2}}.$$

$${\lambda }_1=2^{1/3}-\frac{1}{\sqrt{3}}{2}^{7/12}\sqrt{ϵ},{\lambda }_2=2^{1/3}+\frac{1}{\sqrt{3}}{2}^{7/12}\sqrt{ϵ},$$

$$A=a\frac{1}{\sqrt{3}}\sqrt[12]{2}\sqrt{ϵ},T=\frac{\pi \mu a^2}{2^{2/3}}\left(\sqrt{2}ϵ+9\right).$$

$$E_{el}=\frac{1}{2}\mu \left(\frac{2}{{\lambda }^2}+{\lambda }^2+\frac{a^2{\lambda }^{\prime 2}}{16{\lambda }^3}\right)\pi a^2$$

$$E_s=\gamma \frac{2\pi a}{\sqrt{\lambda }}\sqrt{{\lambda }^2+\frac{a^2{\lambda }^{\prime 2}}{4{\lambda }^3}}.$$

$$E=\int E_s+E_{el}-T(\lambda -⟨\lambda ⟩)dz.$$

$$E=\int \mathrm{const}+\frac{17\pi a^4\mu ϵh^{\prime 2}}{242^{5/6}}+\frac{\pi a^2\mu {ϵ}^2}{12\sqrt[3]{2}}{h}^2\left(h^2-2\right)dz,$$

$$h(z)=\tanh \left(\frac{2^{3/4}\sqrt{ϵ}}{\sqrt{17}a}z\right),$$

Full text

Molecules at a condensed phase’s surface have fewer neighbors than those in the bulk, so all materials suffer an energy proportional to their exposed area. The resulting tendency to reduce area is familiar in fluids, explaining why droplets are round, taps drip and pond-skaters don’t drown. An isolated fluid body always form a sphere under surface tension, but the corresponding problem for a solid body is altogether more subtle, pitting surface tension against bulk elasticity. Here we give a full account of this competition in long solid cylinders.

In general a surface energy, $\gamma$, can only compete with a bulk elastic shear-modulus, $\mu$, at scales below the elastocapiliary length, $l_{cap}=\gamma/\mu$, which is sub-Angtrom in crystalline materials but reaches up to microns or even millimeters in soft solids such as gels and biological tissues. Elastocapillary distortions Roman and Bico (2010) thus occur in the realm of the small and soft, including microfluidics Van Honschoten et al. (2009), the processing of polymer strands Zuo et al. (2005); Naraghi et al. (2007) and many aspects of biology Hazel and Heil (2005); Hannezo et al. (2011); Dervaux et al. (2011); Gonzalez-Rodriguez et al. (2015). In addition to the deformations of soft cylinders considered here Matsuo and Tanaka (1992); Barriere et al. (1996); Mora et al. (2010); Ciarletta and Ben Amar (2012); Taffetani and Ciarletta (2015a, b); Xuan and Biggins (2016), recent work has highlighted capillarity-driven bending of wet elastic rods and sheets Py et al. (2007); Bico et al. (2004); Kim and Mahadevan (2006); Mora et al. (2013), elastic modifications to wetting Style et al. (2013a, b, c), and the inhibitory role of surface tension in elastic creasing Amar and Ciarletta (2010); Yoon et al. (2010); Mora et al. (2011) and cavitation Gent (1990).

In fluid cylinders, surface tension famously causes disintegration into spherical droplets Strutt (Lord Rayleigh). The driver for this Plateau-Rayleigh instability is that an undulation along a cylinder reduces its surface area if its wavelength is longer than its circumference. How the instability proceeds depends in on the fluid Bhat et al. (2010); Eggers (1997), but this simple geometric origin means, inescapably, that all fluid cylinders ultimately disintegrate. A solid cylinder (radius $a$) is also subject to the same geometric reality, but in this case, the surface energy saved by undulation only outweighs the elastic cost of deforming the cylinder if $l_{cap}\gtrsim a$ Barriere et al. (1996), so only thin soft cylinders are unstable. Solid cylinders cannot break into spheres, so the instability produces stable undulating cylinders, fig. 1, with a beads-on-a-string morphology emerging at high surface tension. Previous groups have observed this instability Matsuo and Tanaka (1992); Mora et al. (2010) and studied its onset via linear stability theory Barriere et al. (1996); Mora et al. (2010), showing the first unstable mode has infinite wavelength. More recent work has focussed on stretched cylinders Ciarletta and Ben Amar (2012); Taffetani and Ciarletta (2015a, b), showing, with linear stability theory, that the instability occurs in a particular stretch-interval, $\lambda_{1}<\lambda<\lambda_{2}$, and, via finite element numerics, that, between $\lambda_{1}$ and $\lambda_{2}$, the bead amplitude is essentially constant while bead length falls Taffetani and Ciarletta (2015b). Here we demonstrate that this long-wavelength constant-amplitude beading arises because the instability is actually a longitudinal phase separation between a less and a more stretched phase. We are able to simply predict the instability’s full high-amplitude behavior, including its hysteresis loop, verify our predictions with finite elements, and identify and analyze the domain walls separating the two phases.

We consider a long elastic cylinder, radius $a$ and length $L$, that undergoes a displacement $\mathbf{u}$ and consequent deformation gradient $F=I+\nabla\mathbf{u}$. We model the cylinder as incompressible and neo-Hookean with shear-modulus $\mu$ and surface energy $\gamma$, so the cylinder’s total energy is

[TABLE]

where $A$ is the cylinder’s external area, and incompressibility requires $\textrm{ Det}(F)=1$. Previous studies have shown this energy first becomes linearly-unstable to an infinite wavelength perturbation when $\gamma=6\mu a$ Barriere et al. (1996); Mora et al. (2010). We first verify that long-wavelength behavior persists far beyond threshold by conducting full finite element calculations at $\gamma=8\mu a$ in long aspect ratio cylinders ($L=100a$, see numerical SI) and observe that, although instability arises with a finite wavelength, it coarsens (movie S1), saving energy, until, as seen in fig. 1, the longest possible mode is reached. To understand this, we imagine stretching the cylinder homogeneously to a length $\lambda L$, causing a contraction of its radius to $a/\sqrt{\lambda}$ so that the volume, $V=\pi a^{2}L$, is preserved. The cylinder’s surface area thus increases to $2\pi aL\sqrt{\lambda}$, and its deformation gradient is $F=\mathrm{diag}(1/\sqrt{\lambda},1/\sqrt{\lambda},\lambda)$, so its total energy is

[TABLE]

where we have introduced $\Gamma=\gamma/(\mu a)=l_{cap}/a$, the elastocapillary number for the cylinder. We plot $E(\lambda)$ for several values of $\Gamma$ in fig. 2a. The elastic part of the energy is a convex function minimized at $\lambda=1$. The surface energy never dominates at small or large $\lambda$, and the the total energy always has a single minimum, though it moves to $\lambda<1$ as $\Gamma$ increases Mora et al. (2013). However, the surface tension term, proportional to $\sqrt{\lambda}$, is a concave function and, when $\Gamma$ is sufficient, it introduces a concave region into the total energy. This concave region ($\frac{d^{2}E}{d\lambda^{2}}<0$) is mechanically unstable since therein the tension,

[TABLE]

falls with increasing stretch, as seen in fig.  2b. The limits of this mechanically unstable region are

[TABLE]

where $\frac{d^{2}E}{d\lambda^{2}}=0$. The concave region arises when $\Gamma>\sqrt{32}$, first appearing at $\lambda\pm=2^{1/3}$ and spreading out reaching $\lambda_{-}=1$ at $\Gamma=6$, both in agreement with previous linear stability thresholds Mora et al. (2010); Taffetani and Ciarletta (2015a). In general, given a pair of stretch ratios, $\lambda_{1}$ and $\lambda_{2}$, the cylinder can achieve on average any intermediate stretch $\lambda_{1}<\lambda<\lambda_{2}$ by phase separating into a $\lambda_{1}$ region and a $\lambda_{2}$ region with appropiate length fractions, at an energy cost that corresponds to the chord connecting $E(\lambda_{1})$ and $E(\lambda_{2})$. In the concave region such chords lie below the original energy, so phase-separation saves energy. As understood by Maxwell Clerk-Maxwell (1875), the optimal phase-separation arises when the chord is tangent to the energy at both $\lambda_{1}$ and $\lambda_{2}$:

[TABLE]

The solutions for $\lambda_{1}$ and $\lambda_{2}$ as functions of $\Gamma$ are shown in fig. 2c. As the cylinder is stretched from $\lambda_{1}$ to $\lambda_{2}$, it lowers its energy by phase separating onto the common tangent line, shown as a dashed line in fig. 2a, stretching from $\lambda_{1}$ to $\lambda_{2}$ by adjusting the length-fraction of the two phases. Since the resulting effective energy is a straight line, the tension is constant during the transition as shown by the dotted line in fig. 2b. Furthermore, an obvious conclusion of this phase separation approach is that no amount of surface tension will destabilize a solid cylinder unless it is also stretched to the transition tension $T=(1/L)E^{\prime}(\lambda_{1})$, plotted in fig. 3a. Following the common-tangent yields a convex effective energy, so there is no additional instability.

When it is phase-separated, the difference in final radius between the $\lambda_{1}$ segment (the “bead”) and the $\lambda_{2}$ segment (the “string”) is simply

[TABLE]

This amplitude does not change as the cylinder is stretched from $\lambda_{1}$ to $\lambda_{2}$. Such constant amplitude beading between two stretch ratios, (movie S2), has been observed without explanation in rat sciatic nerve Markin et al. (1999) and elasto-capiliary finite elements Taffetani and Ciarletta (2015b). We compare experimental amplitudes with our finite-element and phase-separation amplitudes in fig. 3b. The three agree well, with finite elements exactly reproducing the phase-separation amplitudes, verifying our approach.

If instability were produced by increasing $\Gamma$ it would be continuous (supercritical) but if it is controlled at fixed $\Gamma$ by stretching from $\lambda_{1}$ to $\lambda_{2}$ it will be observed as a discontinuous (subcritical) with constant amplitude $A(\Gamma)$ throughout. In a stretching cycle the cylinder will be stretched from the energy minimizing contraction $\lambda_{min}$ to $\lambda_{1}$ homogeneously. Upon passing $\lambda_{1}$ it would save energy by phase separating, but the cylinder nevertheless remains mechanically stable until it reaches the end of convexity, $\lambda_{-}$, only then phase separating into $\lambda_{1}$ and $\lambda_{2}$. As stretching continues further, the cylinder shifts material from the $\lambda_{1}$ phase to the $\lambda_{2}$ phase, causing bead width to fall at constant amplitude, until all the material is at $\lambda_{2}$, concluding the instability. Similarly in unloading, the instability should occur at $\lambda_{+}$, but then persist down to $\lambda_{1}$. We thus predict the sub-critical hysteresis loop for a stretch cycle shown in fig. 3c. Hysteresis also slightly complicates the tension-extension behavior: we expect the tension to remain on the homogeneous curve in loading until $\lambda_{-}$ then drop to the phase separation value, and remain there until $\lambda_{2}$.

Our hysteresis predictions are compared with a finite element stretch-cycle in fig. 3c, showing good agreement with one discrepancy: the finite-element loading cycle amplitude drops off slightly early. The reason for this is made evident by plotting (fig. 4a) the longitudinal stretch ratio $\lambda(z)$ (on the axis of the cylinder) as a function of length along the cylinder, running from the middle of the bead ($z=0$) to the end of the cylinder ($z=L/2$), for various values of the total imposed stretch $<\lambda>$ that traverse the instability. In these plots we see clearly that the cylinder has separated into $\lambda_{1}$ and $\lambda_{2}$ regions, but also that they are connected by a transition region—a domain wall—with finite length. The energy cost and length-fraction of a domain wall is negligible in an infinite cylinder, but matters in a finite one: upon stretching a finite cylinder towards $\lambda_{2}$ the amplitude drops early because the vanishing length of the $\lambda_{1}$ segment becomes comparable to the domain wall, so the domain-wall raises the phase-separation energy above the homogeneous energy. As expected, in longer cylinders, the instability persists closer to $\lambda_{2}$. In contrast, at the onset of instability ($\lambda_{-}$) the energy minimizing configuration already requires an extensive length fraction of both phases, so instability proceeds even in modest length cylinders.

Around the critical point, $\Gamma=\sqrt{32}$ and $\lambda=2^{1/3}$, we can solve the phase separation model analytically, and model the domain wall. Expanding (5-6) in $\epsilon=\Gamma-\sqrt{32}$, we see the leading order behavior of $\lambda_{1}$ and $\lambda_{2}$ is

[TABLE]

and hence the amplitude and tension are

[TABLE]

To model the domain wall we require an energy that penalizes stretch gradients along the cylinder. Such a term could arise from additional elastic shear or additional surface area. We capture both effects by allowing the longitudinal stretch, $\lambda(z)$, to vary slowly along the cylinder’s length so a point initially at $(r,z)$ is moved to $(r/\sqrt{\lambda(z)},\int\lambda(z)dz)$. This introduces a new shear term into the deformation gradient, $F_{rz}=-r\lambda^{\prime}/(2\lambda^{3/2})$, increasing the elastic energy per unit-length to

[TABLE]

and the surface energy per-unit length to

[TABLE]

Finally, we use $T$ as a Lagrange multiplier constraining the average stretch to $\left<\lambda\right>$, leading us to the total energy

[TABLE]

Near the critical point we can take the above form for $T$ and write $\lambda=2^{1/3}+\frac{1}{\sqrt{3}}2^{7/12}\sqrt{\epsilon}h(z)$, where $h(z)$ varies between $\pm 1$ through a domain wall. Expanding in $\epsilon$ and keeping only the lowest order terms $h^{\prime}$ and $h$ terms yields

[TABLE]

where the $\epsilon^{2}$ term is minimized by $h=\pm 1$, while the $\epsilon$ term penalizes gradients. Minimizing $E$ w.r.t variations in $h$ requires $17a^{2}h^{\prime\prime}-4\sqrt{2}\epsilon h\left(h^{2}-1\right)=0$, and hence

[TABLE]

an explicit form for the domain-wall. The walls are exponentially localized with a characteristic length $l_{wall}=2\times\sqrt{17}a/(2^{3/4}\sqrt{\epsilon})$ that diverges near the critical point. In fig. 4b we plot $\lambda(z)$ through finite-element domain walls at various $\Gamma$, showing very good agreement with eqn. (14) near $\Gamma=\sqrt{32}$. In fig. 4c we show the length of the finite element domain walls (defined as the $z$ interval required for $\tanh(1)\approx 0.76$ of the $\lambda$ variation) diverges near the critical point, in accord with our prediction. Although this analysis captures the essential physics, it assumes $z$ displacements are independent of $r$; relaxing this assumption leads to a small decrease in $l_{wall}$ by a factor of $\sqrt{34/33}$ (see theoretical supplement).

In summary, the solid Plateau-Rayleigh instability is a phase separation between regions of different but homogeneous longitudinal stretches connected by exponentially localized domain walls. This structure, reminiscent of elastic necking Antman (1973); Ericksen (1975), balloon instabilities Gent (1999); Meng et al. (2014) and volume transitions in swelling gels Tanaka (1978); Matsuo and Tanaka (1992); Boudaoud and Chaïeb (2003), provides a picture of the instability that is both simple and complete. Since phase-separation originates solely in the non-convexity of surface energy, this picture will certainly generalize to more sophisticated elastic constitutive laws, Voigt viscoelasticity and more elaborate extruded prismatic shapes, and, since it naturally produces a high amplitude beads-on-a-string morphology, may underpin the more famous bead/string pattern formed during the breakup of viscoelastic threads Bhat et al. (2010). Mastering these instabilities is essential for the better sculpting, spinning and processing of polymer fibers Zuo et al. (2005); Naraghi et al. (2007).