Capillary control of collapse in soft composite columns

Euler buckling is the elastic instability of a column subjected to longitudinal compressing forces at its ends. The buckling instability occurs when the compressing load reaches a critical value and an infinitesimal deflection leads to a large amplitude deflection. Since Euler’s original study, this process has been extensively studied in homogeneous, isotropic, linear-elastic solids. Here, we examine the nature of the buckling in inhomogeneous soft composite materials. In particular, we consider a soft host with liquid inclusions both large and small relative to the elastocapillarity length, which lead to softening and stiffening of a homogeneous composite respectively. However, by imposing a gradient of the volume fraction or varying the inclusion size we can deliberately manipulate the nature of Euler buckling.


I. INTRODUCTION
An elastic beam under a sufficiently large compressible axial load collapses, or buckles, when an infinitesimal deflection destroys the equilibrium. The critical load for the buckling of homogeneous, isotropic, linear-elastic rods with constant cross-section was derived by Euler in 1744 [1,2], and Lagrange analyzed the higher modes in 1770 [3]. From the mechanical failure of structural elements in civil engineering to the storage of information through controlled buckling of nanoscale beams for future nanomechanical computing [4], the buckling of slender structures has been a focus of studies in engineering, biology and physics for nearly 300 years [5].
The macroscopic response of a solid body to an external force lies at the heart of buckling, the details of which depend on the body shape, material composition, and internal structure. A vast range of distinct responses is displayed in materials with geometric inclusions of different elastic moduli [6], foams modeled by anisotropic Kelvin cells [7], porous and particle-reinforced hyperelastic solids with circular inclusions of variable stiffness [8], fiberreinforced elastomers with incompressible Neo-Hookean phases [9], long cylindrical shells with localized imperfections [10], finitely strained porous elastomers [11], and hyperelastic cylindrical shells [12], to mention but a few. Attempts to describe buckling have lead to, among other things, the celebrated theory of elasticity [5,13] and to finite-element simulation methods [14].
Kirchhoff [15,16] and Clebsch [17,18] described the basic theoretical analysis of elastic rods by replacing the stress acting inside a volume element by a resultant force and the moment vectors attached to a body defining * marc.sune.simon@su.se curve. These "Kirchhoff equations" relate the averaged forces and moments to the curve's strains [e.g., 19].
Recent work shows how the elastic response of soft materials with liquid inclusions is governed by interfacial stresses [20][21][22][23] and suggests the possibility that capillarity may play an important role in buckling instabilities. To that end, we reformulate the stability analysis of compressed rods viewed from the perspective of the theory of elasticity [5,13] to account for the surface tension effects of the inclusions. We incorporate the physics of capillarity into the Kirchhoff equations through the elastic moduli as given by a generalization of Eshelby's theory of inclusions [22][23][24]. Eshelby's theory describes how an inclusion of one elastic material deforms when it is embedded in an elastic host matrix [25]. However, it has recently been discovered that Eshelby's inclusion theory breaks down when the inclusion size R approaches the elastocapillary length L ≡ γ/E, where γ is the inclusion/host surface tension and E is the host Young's modulus [22][23][24]. Importantly, when R > L (R < L) the composite softens (stiffens). This basic physical process, wherein the inclusion size controls the Young's modulus of the the composite, E c , reveals the possibility of controlling the buckling process by controlling the properties and distribution of the inclusions.
A quantitative treatment of how the inclusion size, R, and volume fraction, φ, in soft composites influences their bulk mechanical properties underlies our understanding of their response under loads. In particular, by determining how the spatial variation of R and φ modify Euler buckling we provide a framework of either tailoring a material response or explaining observations in naturally occurring soft composites. Canonical examples of the latter include slender composite structures such as plant stems [26,27] and bones [28], bacterial biofilaments [29] or plant tendrils [30]. Indeed, these latter systems [29,30] very often grow into axisymmetric elongated structures by adding new material at a small active growing zone located near the tip, creating a varying composition along the growth axis. Although our analysis is confined to static elastic Kirchhoff rods, our results are clearly of use in interpreting the buckling instabilities prompted by tip growth [31,32] as well as the phenomenon of morphoelasticity induced by time-dependent compression [33,34].

A. Stability analysis
Euler's laws describe the conservation of linear and angular momenta and thereby underlie the dynamics of elastic rods. Consider two adjacent cross-sections of a rod defined by an infinitesimal element of length dl. In equilibrium, Euler's laws yield the balance of the total forces and the total couple exerted on the reference segment, and hence in the continuous limit are where F is the resultant internal stress on a cross-section, K is the external bending force on the rod per unit length,t = dl/dl is the unit vector tangential to the rod, and M is the resultant moment of the internal stresses on the cross-section (Fig. 1). The constitutive relationship between the resultant stresses σ ij and strains ǫ ij in cartesian coordinates is

Sfrag replacements
where again E c is Young's modulus of the composite and the z-axis coincides with the axis of the rod. The former relation holds for a slightly deflected thin rod in absence of an external twisting moment, so that the bending occurs in a single plane. Therefore, the strain-energy is a quadratic function of the transverse coordinates, and the constitutive relations for the bending moments are where X, Y give the displacement of points on the elastic line from their corresponding undeformed positions, and I x and I y are the second moments of area about the x and y axes respectively, viz., Assuming a constant cross-section with area S along the rod, I x and I y are constant throughout despite any inhomogeneous mass distribution along z that we may consider.
In the absence of external transverse bending forces (K = 0), when external tensile (or compressive) forces are applied to the rod ends, F z ≡ T , the equation of equilibrium for a rod with circular cross-section (I x = I y ≡ I) is obtained by combining Eqs. (1), (2) and (4) and is where primes denote differentiation with respect to z.
In classical Euler buckling one solves the following boundary value problem for Eq. (6). A rod of length l compressed with axial force T such that its ends cannot be displaced X(z = 0) = X(z = l) = 0 but can freely rotate X ′′ (0) = X ′′ (l) = 0. This problem has the trivial solution X(z) = 0, which corresponds to the rod remaining straight. When T reaches a critical value T cr the trivial solution becomes unstable and the rod buckles. As the compressive force increases, successive higher order buckling modes appear. We now ask how these classical solutions are modified for soft composites with structured effective Young's modulus, E c .

B. Composite mechanics
The theory of effective elastic moduli of solid composites is generally ascribed to Eshelby [25]. Initially conceived to treat composites of host materials such as glass or steel, with E = O(GPa), containing dilute inclusions, Eshelby's theory has been extended to non-dilute composites [35][36][37]. However, Eshelby's approach does not account for the energy between the inclusion and the host. Although this is quantitatively valid when the inclusion size is much larger than the elastocapillary length R ≫ L, as defined above, such is not the case otherwise. Recently, the other limit where R L and surfacetension effects in soft solids are important has been a focus of research [38][39][40][41], and a generalized version of Eshelby's theory that accounts for the effects of interfacial surface tension has been derived [22,23]. This theory provides an expression for the effective elastic modulus for soft composite solids, E c , in terms of the elastic moduli of the host material (Young's modulus E and Poisson's ratio ν), the dimensionless number γ′ ≡ L/R, and the inclusion volume fraction φ as follows where we have assumed an incompressible solid; ν = 1/2. In the appropriate limit for the original Eshelby theory [25] the composite softens as φ increases, which is the behavior clearly recovered when γ′ < 2/3 in Eq. (7) [22,23]. However, surface energy effects become important when γ′ > 2/3, a situation of particular relevance when the host is soft and the inclusions softer, and then the composite becomes stiffer, for example when liquid droplets of size O(100µm) are embedded in compliant materials, such as gels, with E = O(kPa) [22,23]. Stiffer materials with E = O(MPa), like elastomers, may also exhibit enhanced stiffness due to capillarity when liquid inclusion size is O(0.1µm). These results hold in the dilute and non-dilute limits [42,43].

III. COMPRESSING INHOMOGENEOUS COMPOSITES
The dimensionless form of Eq. (6) that we will treat is where we scale length with I/S and hence the dimensionless axial force is Γ ≡ T /(E S), the dimensionless composite Young's modulus is E rel ≡ E c /E, and we have not changed the notation for the independent and dependent variables. Given the expression for the effective Young's modulus of a soft composite, Eq. (7), an axially inhomogeneous elastic modulus can be achieved by varying either the inclusion volume fraction, φ(z), or the ratio of the elastocapillary length to the inclusions radius, γ′(z). To this end, we analyze a circular cross section column, or rod, with the linear inclusion volume fraction profile, where φ(z = 0) ≡ φ 0 and l is the dimensionless column length.
We now solve the Euler buckling boundary value problem for Eq. (8), where X(0) = X(l) = X ′′ (0) = X ′′ (l) = 0, but we use Eq. (9) in E rel . Eq. (8) becomes where We use standard variation of parameters (e.g., [44]) to solve Eq. (10), up to a constant C, in terms of Airy functions Ai, Bi as whereã = Γ cr a andb = Γ cr b. When the compression exerted on the inhomogeneous column ends exceeds a critical value, Γ cr , the column deflections are given by Eq. (11) as a function of the strength of the gradient in inclusion volume fraction, |φ 0 /l|. The values of Γ cr are given by the non-trivial solutions of Eq. (10) and hence are the roots of the transcendental equation which we solve numerically to determine the failure modes. Next we study the first two modes of the heterogeneous composite rod in both the stiffening and the softening regimes.

A. Stiffened composites
As noted in §II B, when γ′ > 2/3 liquid inclusions in a soft host will stiffen the composite and hence E rel > 1. Therefore, when the inclusion volume fraction is governed by Eq. (9), the stiffness of the column will decrease along the axis: the composite Young's modulus decreases from z = 0 towards the bulk host elasticity at the opposite end where E rel (z = l) = 1 (upper panel in Fig. 2a).
To examine the behavior of the soft composite scenarios we compare them to a column with a homogeneous constant stiffness, or the Kirchhoff case, the first two sinusoidal buckling modes of which are For n = 1 the maximum deflection of the Kirchhoff case occurs at z = l/2, which we compare in Fig. 2a to the heterogeneous stiffened composite rod, obtained by substituting the numerical result for the first critical compression Γ cr into Eq. (11) [45]. There is clearly an asymmetry that distinguishes the responses of the homogeneous and the heterogeneous elastic columns. In particular, the apex of the deflection for the latter is shifted towards the compliant end, z = l, which is an anticipated consequence of Eq. (9).
As shown in Fig. 2b, the asymmetry becomes more evident for the one-cycle sinusoidal second buckling mode; n = 2 in Eq. (13). Physically, unlike the symmetric bending of a rod with constant elastic modulus, the rod with  varying Young's modulus responds unevenly, with the stiffened region driving the compliant region towards the compliant boundary. This soft composite behavior can be brought out by examining the mode-deviations, ǫ ≡ X n (z) − X(z), as shown in Fig. 3. In the first buckling mode, n = 1, ǫ changes from positive in the stiffened end to negative towards the compliant side (Fig. 3a). Accordingly, the heterogeneous rod buckles less (more) than the homogeneous column where it is stiffer (softer). These differ- ences are reflected in terms of the behavior of ǫ(z) near the midpoint of the column. For example, in the inset of Fig. 3a we see that the sum of ǫ(z) and its specular reflection with respect to z = 0.5, ǫ(1 − z), is non-zero. The symmetry breaking of the deviations arises from the spatial variation of E rel and shown in the upper panel of Fig. 2a, a phenomenon enhanced at higher compressions as seen for the second buckling mode shown in Fig. 3b. Finally, Fig. 3 demonstrates the two consequences of the ratio of elastocapillary length to inclusion radius, γ′, and the gradient of inclusion volume fraction, |φ 0 /l| in this stiffening regime. First, the stiffer the rod (large γ′ and/or large φ 0 ), the greater the deviation relative to the Kirchhoff constant elastic modulus case, and the more prominent the spatial asymmetry. Second, the system is more sensitive to the inclusion volume fraction gradient than to the inclusion size variation.

B. Softened composites
In the limit of small elastocapillary scale or large inclusion size we have γ′ → 0, in which case the effective Young's modulus, Eq. (7), reduces to Eshelby's result for liquid droplets in an elastic solid; E c = E/(1 − φ). This differs from the stiffening case in how the dilute theory (DT) [22,23] in the softening case deviates from nondilute approaches as φ increases [24,42,43].
Here we employ a particular non-dilute theory that generalizes the multiphase scheme introduced by Mori and Tanaka [46] by treating the effects of interfacial stress between the inclusion and the bulk host due to Mancar-ella et al. [24] (MSW). The effective Young's modulus of the composite is We take this approach because the expression is simpler than the quantitatively equivalent one obtained using the three-phase theory of soft composite solids [42,43]. We again consider the linear model for the distribution of inclusions given by Eq. (9), within the MSW approach the equilibrium equation (8) for a column with hinged ends becomes where a = 2 + (4 The solution of the boundary value problem for Eq. (15) is where U (n(Γ), 0; g(z, Γ)) is the confluent hypergeometric function, L (−1) n(Γ) (g(z, Γ)) are the associated Laguerre polynomials; g(z, Γ) ≡ −2 (c+d·z) Substituting the roots of Eq. (17) into Eq. (16), we find the corresponding buckling modes of the softened composite column. In Fig. 4 we compare the first two modes with the homogeneous reference cases of Eq. (13), and with the DT result Eq. (11). The asymmetric buckling of the heterogeneously softened rod in Fig. 4 is analogous to the stiffened case in Fig. 2. Namely, the most acute bending occurs at the compliant end of the rod, in the softened case this is towards the origin (see the upper plot in Fig. 4a), with the (larger) bulk modulus being reached at the end of the rod where the boundary conditions for all three cases must be satisfied and thus both the DT and MSW deflections must be smaller than the homogeneous case. Although the buckling modes plotted in Fig. 4 are qualitatively similar, there is some deviation between those results obtained by DT and MSW. Fig. 4 shows that the asymmetry about the midplane of the column is amplified for the MSW case as expected from the larger deviation of E rel relative to the DT and homogeneous cases.
Finally, we compare and contrast the stiffening and softening regimes for both the DT and MSW cases, by plotting in Fig. 5 the first two critical loads (Γ cr ) associated with the numerical roots of Eqs. (12), and (17) respectively. Here the stiffening and softening regimes are in clear evidence. Firstly, relative to the homogeneous case (φ 0 = 0), the critical load for φ 0 = 0 is larger (smaller) in the stiffening (softening) regime. These differences are substantially enhanced as the magnitude of the inclusion gradient, φ 0 , increases.  Also, Fig. 5 shows the expected convergence of both approaches in the dilute limit for the softening regime as well as the capillary force dependence of their deviation with increasing volume fraction. Moreover, and perhaps of use for other modeling purposes, in the softening regime Fig. 5a shows that the critical stress depends linearly upon the magnitude of the inclusion gradient in the MSW framework whereas in the stiffening regime the DT   In both §III A and III B we either linearly stiffened or softened a column using Eq. (9) for a given ratio of elastocapillary length to inclusion radius, γ′. In consequence, the bulk modulus decreased or increased with distance along the column. This suggests the possibility of a "polar" configuration in which one side of a column will be stiffer than the bulk host and the other will be softer. The most straightforward treatment of this is a linear variation of γ′(z), akin to Eq. (9), viz., where γ 0 ′ and γ l ′ are the values of γ′ ≡ L/R at z = 0 and z = l, respectively. The resulting effective Young's modulus, for both the DT and the MSW approaches, can be expressed as a ratio of two linear polynomials. The equilibrium equation corresponding to Eq. (8) is thus mathematically analogous to Eq. (15), that is the MSW theory with φ(z) obeying Eq. (9) for the MSW approach. Parameters b and d for both models are negative when γ 0 ′ > γ l ′ (0 ≤ φ ≤ 1), and hence evaluating Eqs. (16) and (17) yields complex numbers. However, we note that in both of these equations y 1 (z) ≡ exp(−g(z, Γ)/2) U (n(Γ), 0; g(z, Γ)) and y 2 (z) ≡ exp(−g(z, Γ)/2) L (−1) −n(Γ) (g(z, Γ)) are eigenfunctions of the linear differential operator acting on the function X(z) and have eigenvalue 0, corresponding to the equilibrium equation (15). It can be shown that the associated linear differential operator is Hermitian. Therefore, Eqs. (16) and (17) can be rewritten in terms of linear combinations of the eigenfunctions y 1 (z) and y 2 (z), such that they will take on real values. Nonetheless, since the deflection of a polar compressed rod given by Eq. (16) is defined up to an undetermined constant, without loss of generality we consider the real part of Eq. (16), which is our approach for the remainder of this section.
We note that in the softening regime, as expected, the DT deviates significantly from the non-dilute theories as φ increases as seen in Figs. 4 and 5. In the stiffening regime both non-dilute theories, the MSW and threephase generalized self-consistent theory [42], differ as φ increases beyond the dilute limit. Therefore, Eqs. (16) and (17), for φ constant and γ′ linear, capture the behavior in the dilute limit. Moreover, the two non-dilute theories accounting for interfacial tension [42,43] apply both for stiffening and softening regardless of the volume fraction. Unfortunately, however, the expression for the effective Young's modulus in the strongly nondilute regime is so complex that Eq. (8) cannot be easily analyzed mathematically.
As shown in Fig. 6, the effective Young's modulus as defined by Eqs. (7) and (14)  dence on γ′. Therefore, in the dilute limit with a linear model such as that in Eq. (18), the range of Young's modulus is limited. Such a constraint could be overcome by considering a higher order functional variation of γ′(z). However, the complications introduced by such a variation introduce more structure than is warranted by present observational constraints. The buckling shapes of a compressed polar elastic compressed rod using Eq. 18 are shown in Fig. 7. We chose an intermediate volume fraction to tailor the transition of the Young's modulus between the rod's ends (see Fig. 7aupper plot), yet without going too deeply into the nondilute regime. We observe the same qualitative features as for the previous cases, when only stiffening or softening was considered. Although strictly speaking dilute conditions do not apply for φ = 0.4, we point out that there is near perfect agreement between the DT and MSW approaches.

IV. CONCLUSIONS
We have studied the buckling of inhomogeneous soft composite columns, or rods, with axially varying elasticity. Their spatial structure is tailored either by changing the volume fraction of inclusions (φ) or the ratio of the elastocapillarity length to the inclusion size (γ′) along the column axis. We have extended the classical theoretical description of compressed elastic rods by incorporating these inclusion/host surface tension effects on the effective elastic modulus of the mixture. The resulting equilibrium equation that accounts for the rod's response to a compressive force has analytical solutions for a linear  model for the variation of either φ or γ′. This provides a framework of broad relevance to soft composite materials and could be tested by considering a distribution of sizes of the liquid droplets embedded in a soft solid host, which might be possible by suitable variations of the experiment described by Style et al. [22].
We studied three different problems of heterogeneous rods: stiffening (softening) by increasing the volume fraction of small (large) inclusions (i.e., of radius R < (3/2)L (R > (3/2)L)) as a function of distance along the axis as described by Eq. (9), and polar elasticity by introducing a gradient of the inclusion size by Eq. (18) at a constant volume fraction. Their principal common feature is the intuitive result that a compressed column of variable elasticity bends most easily where it is softer. Accordingly, the symmetric buckling characteristic of a rod with homogeneous properties is broken whenever a rod with inhomogeneous stiffness is considered. In both the dilute and non-dilute regimes, these three general cases, mutatis mutandis, include the overall behavior of the of the buckling instabilities in linearly heterogeneous soft composite rods.
We have studied these rather simple models with the hope that considering an elementary architecture might motivate experimental testing of these ideas as well as their potential for biological relevance [31,32]. Indeed, recent advances in understanding multi-scale morphomechanical effects in biological systems [47] and experimental realizations of soft-solid cavitation [48] and of stiffness-modulated elastic solids [49], provide useful techniques and settings to implement and explore our models.