## Abstract

This work investigates the morphological stability of a soft body composed of two heavy elastic layers attached to a rigid surface and subjected only to the bulk gravity force. Using theoretical and computational tools, we characterize the selection of different patterns as well as their nonlinear evolution, unveiling the interplay between elastic and geometric effects for their formation. Unlike similar gravity-induced shape transitions in fluids, such as the Rayleigh–Taylor instability, we prove that the nonlinear elastic effects saturate the dynamic instability of the bifurcated solutions, displaying a rich morphological diagram where both digitations and stable wrinkling can emerge. The results of this work provide important guidelines for the design of novel soft systems with tunable shapes, with several applications in engineering sciences.

This article is part of the themed issue ‘Patterning through instabilities in complex media: theory and applications.’

## 1. Introduction

Shape transitions in soft solids result from a bifurcation of the elastic solutions driven by either geometrical or constitutive nonlinearities. The characterization of the emerging morphologies is the subject of morpho-elasticity, a recent branch of continuum mechanics at the interface between finite elasticity and perturbation theory. This vibrant research field has rapidly developed in the last decade, pushed by the technological availability of experimental devices controlling the extreme deformations of soft incompressible materials, such as hydrogels [1–3] and elastomers [4,5].

Although their boundary value problems are intrinsically different, the study of pattern formation in soft solids has highlighted some similarities, as well as several relevant differences, to the instability characteristics of hydrodynamic systems. For example, if the surface tension in a thin fluid filament triggers the formation of droplets, which spontaneously break down [6], such a dynamics can be stabilized by elastic effects in soft solid cylinders [7], thus driving the emergence of stable beads-on-a-string patterns [8]. Similarly, while fingering at the interface of two immiscible viscous fluids is an unstable process [9], stable digitations may occur after a subcritical bifurcation for a fluid pushing against an elastic surface [10] and at the interface between a thin elastic layer and a glass plate to which it is adhered [11].

Another interesting example is the gravity-induced instability in an elastic layer attached to a rigid substrate with a traction-free surface facing downwards. In contrast to gravity waves in a fluid layer, the free surface experiences fluctuations that eventually saturate when the large deformations store an elastic free energy of the same order as the corresponding variation in the potential energy. The linear stability analysis of this problem has been recently performed [12], then refined to consider the effect of an applied strain on the elastic layer [13]. Nonetheless, this problem had been previously solved using numerical techniques [14], and is often used as a test case to study the stability of discrete solutions obtained by means of mixed finite-element techniques [15,16].

Since Rayleigh [17] and Taylor [18], it has been well known that the horizontal interface between one fluid layer and a lighter one below it is unstable to sinusoidal perturbations, forming protrusions that grow with a characteristic time. However, if one takes surface tension into account, the growth of small wavelength protrusions is inhibited by capillary effects, thus larger wavelength drops grow and eventually drip [19]. In this work, we aim at studying this kind of gravity instability in a soft system made of two heavy elastic layers attached at one end to a rigid surface. In particular, we are interested in characterizing both the pattern formation and its nonlinear evolution, determining the interplay between elastic and geometric effects for the emergence of a given pattern.

The article is organized as follows. In §2, we define the nonlinear elastic problem and identify its basic solution. In §3, we perform the linear stability analysis of the problem, deriving the marginal stability curves as a function of the elastic and geometric dimensionless parameters. In §4, we perform numerical simulations using finite elements to study pattern formation in the fully nonlinear regime. In §5, we finally discuss the theoretical and numerical results, adding a few concluding remarks.

## 2. The nonlinear elastic problem and its basic solution

In a Cartesian coordinate system with unit base vectors *E*_{i}, with *i*=(*X*,*Y*,*Z*), we consider a soft body made of two hyperelastic layers, as sketched in figure 1.

Let be the three-dimensional Euclidean space, the body occupies a domain , having a thickness *H* along the *Y*-axis and a length *L* along the *X*-axis, with *L*≫*H*. We also consider that the body is infinitely long along the *Z* direction, so that a plane strain assumption can be made, hence

The body is clamped to a rigid substrate at *Y* =0, so that its volume *Ω* can be split in the two subdomains *Ω*_{a} and *Ω*_{b} occupied by the constituting layers, such that
where ** X** is the material position vector and

*H*

_{a}and

*H*

_{b}are the thicknesses of the layers.

Indicating that ** x**=

**(**

*x**X*,

*Y*) is the spatial position vector, the kinematics is described by the geometrical deformation tensor

*F*=Grad

**. We also assume that the layers behave as incompressible neo-Hookean materials and the strain energy density of each layer is given by 2.1where**

*x**I*

_{1}is the trace of the right Cauchy–Green tensor

*C*=

*F*

^{T}

*F*and

*p*is the Lagrangian multiplier enforcing the internal constraint of incompressibility.

Using the constitutive assumption given by equation (2.1), the first Piola–Kirchhoff stress tensor *S* reads

Assuming quasi-static conditions, the balance of linear momentum for the elastic body subjected to its own weight reads
2.2where Div is the material divergence, *ρ*_{a} and *ρ*_{b} are the densities of the layers and ** g**=

*g*

*E*_{Y}is the gravity acceleration vector.

In the following, we aim to provide a unified analysis of the two configurations depicted in figure 1. For the sake of notational compactness, we consider that a positive *g* represents the body hanging down a rigid wall (figure 1*a*), and a negative *g* the body placed on top of a rigid substrate (figure 1*b*).

The two boundary conditions at the fixed substrate and at the free surface read
2.3where ** u**=(

**−**

*x***) is the displacement vector field. The elastic boundary value problem is finally complemented by the following displacement and stress continuity conditions at the interface between the two layers, respectively: 2.4**

*X*The boundary value problem given by equations (2.2)–(2.4) admits a basic solution given by
2.5so that no basic deformation is allowed by the incompressibility constraint, and the body is subjected to a hydrostatic pressure linearly dependent on *Y*. We also highlight that the pressure field given by equation (2.5) is discontinuous across the material interface if *μ*_{a}≠*μ*_{b} or *ρ*_{a}≠*ρ*_{b}.

## 3. Linear stability analysis of the basic solution

### (a) Incremental equations

We now aim at investigating the stability of the basic elastic solution given by equation (2.5) using the method of incremental deformations superposed on a finite strain [20].

Let us perturb the basic configuration by applying an incremental displacement *δ*** u**. If we set

*δ*

*F*=Grad

*δ*

**, the linearized incremental Piola–Kirchhoff stress tensor is where is the tensor of instantaneous elastic moduli,**

*u**I*is the identity tensor,

*δp*is the increment of the Lagrangian multiplier

*p*and the two dots operator (:) denotes the double contraction of the indices, namely

Recalling that the basic solution is undeformed, the incremental incompressibility and equilibrium equations read, respectively, 3.1and 3.2

The incremental counterparts of two boundary conditions at the fixed substrate and at the free surface may be rewritten as, respectively, 3.3 3.4 3.5 3.6

Similarly, the incremental versions of the displacement and stress continuity conditions at the interface read 3.7and 3.8

In the following, we derive the solution of the incremental boundary value problem given by equations (3.1)–(3.8).

### (b) Solution of the incremental boundary value problem

Let us now assume an ansatz by variable separation in the expression of the incremental displacement, namely
3.9where *k* is the horizontal spatial wavenumber. We recall that such a functional dependence along the *X* direction suitably describes both the infinite geometry, for which *k* is a continuous variable, and a finite length *L*, so that *k*=2*πn*/*L* with integer mode *n*.

From equation (3.2), we get that 3.10

From the first component of equation (3.1), we obtain the expression for *δp* as
3.11

By substituting equations (3.10) and (3.11) in the second component of equation (3.1), we obtain the following ordinary differential equation: 3.12which is valid for both layers and whose solution is given by 3.13

Hence, setting
we impose the conditions given in equations (3.3)–(3.8). We find eight linear algebraic equations in the unknowns *v*_{j}, , so that we can write such a system in the compact form *M*** v**=

**0**, where

*M*is the 8×8 coefficient matrix. Hence, we find that a non-null solution of such a linear system exists if and only if 3.14The full form of

*M*is reported in appendix A.

### (c) Results of the linear stability analysis

Let us now discuss the results of the linear stability analysis by making use of the following dimensionless parameters:

A great simplification arises if we set both *α*_{ρ}=1 and *α*_{μ}=1 or if we impose *α*_{H}=0 in equation (3.14), so that the body is made of a single homogeneous slab. In particular, we recover the same expression reported in [12]
highlighting that an elastic bifurcation occurs for the critical value *ρgH*/*μ*_{a}≃6.22 with critical wavenumber *kH*≃2.11.

Let us now analyse the resulting solutions when *α*_{ρ}=1, namely assuming that the body force is the same for both layers. In the case in which *α*_{μ}≠1 or *α*_{H}≠1, we find only one root of equation (3.14). In figure 2, we depict the resulting marginal stability curves varying the parameters *α*_{μ} and *α*_{H}.

We denote by *γ*_{cr} the critical value of *γ*, i.e. the minimum value of the marginal stability curve obtained by fixing *α*_{H} and *α*_{μ}. We denote by the critical wavenumber, namely the value of for which the marginal stability curve has a minimum. All the critical values of the marginal stability curves have been found by using Newton’s method with the software Mathematica 11.0 (Wolfram Research, Champaign, IL, USA).

In figure 3, we plot the critical values *γ*_{cr} and when varying *α*_{H} and *α*_{μ}. We find that *γ*_{cr} strongly depends on *α*_{μ} and *α*_{H}. In figure 3*a*, we find that if we increase the parameter *α*_{μ} the critical value *γ*_{cr} also increases, so that high values of *α*_{μ} have a stabilizing effect. On the contrary, in figure 3*c* we find that if we increase *α*_{H} the critical value *γ*_{cr} decreases. We highlight that, if *α*_{H} tends to zero, we obtain that *γ*_{cr}≃6.22, which is the single-layer limit discussed before. The same limit is found for *α*_{μ} tending to zero, as it represents the case in which the bottom layer in figure 1*a* becomes infinitely soft. The critical wavelength is always of the same order as the body thickness, and is more influenced by the parameter *α*_{H} if *α*_{μ}>1, as we can note from figure 3*b*,*d*.

The case in which *α*_{ρ}=1 is of particular interest in the applications because it is reproducible in experiments using hydrogels. In fact, these soft materials are mainly composed of water, thus they have a density which is of the order of 10^{3} kg m^{−3}. Nonetheless, by small variations in the cross-link concentration, it is possible to obtain a shear modulus *μ* ranging from 100 Pa to 10 kPa. For example, if we consider two hydrogel layers with *H*_{a}=*H*_{b} (figure 2*b*), where the clamped one has *μ*_{a}=300 Pa and the other *μ*_{b}=600 Pa, we find that *γ*_{cr}≃4.4366. Accordingly, an instability would appear at *H*_{a}≥*μ*_{a}*γ*_{cr}/(*ρg*)≃13.57 cm.

The general case in which *α*_{ρ}≠1 is a bit more complex; in fact, equation (3.14) can be written in the following compact form:
3.15where the coefficients *c*_{1}, *c*_{2} and *c*_{3} depend on *α*_{H}, *α*_{μ}, *α*_{ρ} and , as reported in appendix B.

Even if their expressions are very cumbersome, we can still make some general observations. In fact, we observe that *c*_{1} does not depend on *α*_{ρ}, whereas, if we fix the other variables, *c*_{3} has a different sign if *α*_{ρ}>1 or if 0<*α*_{ρ}<1. Hence, one of the two real roots of equation (3.15) changes sign if we consider *α*_{ρ}>1 or 0<*α*_{ρ}<1.

Thus, we make a distinction in the following between these two cases, which physically correspond to the two configurations depicted in figure 1.

#### (i) Case (a): free surface instability (*γ*>0)

The configuration shown in figure 1*a* undergoes a morphological transition if equation (3.15) possesses at least a positive root for *γ*, because we assume *g*>0.

In figures 4 and 5, we depict the marginal stability curves when we vary the parameters *α*_{H}, *α*_{ρ} and *α*_{μ} in equation (3.15). In this case, we find that the instability is localized at the free boundary of the slab, i.e. at *Y* =*H*.

We can observe that we have the same behaviour as discussed for the case *α*_{ρ}=1: if we increase the parameter *α*_{μ} we obtain a stabilization of the system (i.e. *γ*_{cr} increases) whereas if we decrease *α*_{H} we have instability for lower values of *γ*.

#### (ii) Case (b): interfacial instability (*γ*<0)

Conversely, the configuration shown in figure 1*b* undergoes a morphological transition if equation (3.15) possesses a negative root for *γ*, as we assume *g*<0. As previously discussed, this happens only if *α*_{ρ}>1, meaning that the top layer is heavier than the bottom one. Thus, this case is the elastic analogue of the Rayleigh–Taylor instability. As found in fluids, the instability is concentrated at the interface between the layers and decays away from it.

In this configuration, we define the critical value *γ*_{cr} as the maximum of the marginal stability curve for fixed *α*_{H},*α*_{μ} and *α*_{ρ}.

In figure 6, we set *α*_{ρ}=2 and we plot the marginal stability curves for several values of the parameters *α*_{H} and *α*_{μ}. Also in this case, we highlight that increasing the parameter *α*_{μ} stabilizes the system, whereas an increase in the parameter *α*_{H} favours the onset of the interfacial instability.

In the next section, based on the results of the linear stability analysis, we build the simulation tools for studying the fully nonlinear morphological transition.

## 4. Post-buckling analysis

In this section, we numerically implement the fully nonlinear problem given by equations (2.2)–(2.4). We finally report the results of numerical simulations for the two cases under consideration, highlighting the morphological evolution of the emerging patterns in the fully nonlinear regime.

### (a) Finite-element implementation

The boundary value problem is implemented by using the open source tool for solving partial differential equations, FEniCS [21]. In order to enforce the incompressibility constraint, a mixed formulation has been chosen. If the two layers have different stiffnesses or mass densities, the pressure field may present a discontinuity at the interface between the two layers, according to the basic solution given by equation (2.5). Accordingly, we used the element *P*_{2}–*P*_{0} [22] in numerical simulations.

This element discretizes the displacement with piecewise quadratic functions and the pressure field with piecewise constant functions, so that we can correctly account for a discontinuous pressure field. It is also numerically stable in linear elasticity [22] and it has been successfully used in several nonlinear applications [14,15].

We use a rectangular mesh whose height is *H*=1 and whose length is the critical wavelength , where is the critical value arising from the previous linear stability analysis and depending on *α*_{ρ}, *α*_{H} and *α*_{μ}. We set ** u**=

**0**at

*Y*=0 and we impose periodic boundary conditions at

*X*=0 and

*X*=

*λ*. The number of elements used depends on the length of the mesh; the maximum number of elements used is 30 000.

In order to investigate the post-buckling regime, we impose a sinusoidal imperfection at the top boundary of the mesh with a wavenumber *k*_{cr} and an amplitude *h*=10^{−4}*H*, as done in [23,24].

The solution is found by using an incremental iterative Newton–Raphson method increasing (or decreasing in the fluid analogue case) the control parameter *γ*. In each iteration, the calculation is performed by using the linear algebra back-end PETSc (Portable, Extensible Toolkit for Scientific Computation) and the linear system is solved through a LU (lower–upper) decomposition. The code automatically adjusts the increment of the control parameter if *γ* is near the critical value *γ*_{cr} or when the Newton–Raphson method fails to converge.

Since secondary bifurcations may appear in such a dispersive problem, due to sub-harmonic resonance phenomena in the fully nonlinear regime, we performed further simulations using the approach proposed in [25]. Accordingly, we looked for period-doubling and period-tripling secondary bifurcations by using as the computational domain the sets [0,2 *mπ*/*k*_{cr}]×[0,1] with *m*=2 and *m*=3, respectively. However, we did not find any further bifurcation in the parameter range considered in this article, in agreement with the experimental observations performed in the single-layer case [12].

### (b) Numerical results

In the following, we report the results of the numerical simulation for the two cases under consideration.

#### (i) Case (a): free surface instability (*γ*>0)

We first implement the case described in figure 1*a*, setting *α*_{ρ}=1 in order to mimic the behaviour of a slab made of two hydrogel layers. In figure 7, we depict the results of the numerical simulations for two different values of *α*_{μ}. In particular, we highlight that the deformation is localized at the free boundary of the body, and it evolves towards the formation of stable hanging digitations. Let *Δh* be the maximum vertical distance of the points on the free surface and *Δl* be the horizontal distance between the points which have initial coordinates (*λ*/4,*H*) and (3/4 *λ*,*H*), so that *Δl*/*λ*=0.5 if *γ*<*γ*_{cr}. Thus, we employ *Δh* and *Δl* to study the nonlinear evolution of the finger morphology.

As shown in figure 8*a*, we find that the fingering height *Δh* continuously increases as the control parameter *γ* goes beyond its critical value. When performing a cyclic variation of the order parameter, where we first incremented *γ* until a value *γ*_{max}>*γ*_{cr} and later decreased it back to the initial value, we found that both *Δh* and *Δl* did not encounter any discontinuity, always following the same curve in both directions. Moreover, in the weakly nonlinear regime *Δh* increases as the square root of the distance to threshold of the order parameter, thus highlighting the presence of a supercritical pitchfork bifurcation.

The shape and the thickness of these fingers strongly depend on the stiffness and the thickness of the two layers. As shown in figure 8*b*, the fingers become thicker as we increase *α*_{μ}.

We remark that the maximum diameter *h* of the mesh elements is chosen as the maximum value such that the resulting *Δh* and *Δl*/*λ* differ by less than 10^{−3} from the corresponding values obtained using a refined mesh with *h*/2.

#### (ii) Case (b): interfacial instability (*γ*<0)

We now focus on the elastic analogue of the Rayleigh–Taylor instability, which occurs in the configuration depicted in figure 1*b*. Here we set *α*_{ρ}=2, so that the top layer is heavier than the bottom one.

The simulation results are depicted in figure 9 for two different values of *α*_{H}. In particular, we find a behaviour similar to the Rayleigh–Taylor instability in fluids: the displacement is concentrated at the interface of the two layers, forming a marginally stable undulation.

Let *Δh* denote here the maximum vertical distance of the points on the interface between the two layers and let *Δl* be the horizontal distance between the points which have initial coordinates (*λ*/4,*H*_{a}) and (3/4 *λ*,*H*_{a}), so that *Δl*/*λ*=0.5 if *γ*>*γ*_{cr}. In figure 10, we show the nonlinear evolution of such morphological parameters as a function of *γ*. As in the previous case, we measured the quantity *Δh* decreasing the parameter *γ*, and found a continuous increase in the height of the undulation, as reported in figure 10*a*. We highlight that the normalized thickness *Δl*/*λ* strongly depends on the parameter *α*_{H}, as we can see from figure 10*b*.

In fact, for thin soft layers the undulation decreases its width while decreasing *γ* beyond its critical value, thus forming a digitation. Conversely, the undulation width increases for thick top layers, thus forming a stable wrinkle. In summary, *Δl*/*λ* increases if the top layer is sufficiently thin, while it decreases if the top layer is above a critical thickness. In both cases, the resulting morphology is perfectly reversible after cyclic variations of the order parameters, highlighting the presence of a supercritical pitchfork bifurcation.

## 5. Discussion and concluding remarks

In this work, we used theoretical and computational tools to investigate the stability of a soft elastic bilayer subjected only to the bulk gravity force.

Assuming that both layers are made of incompressible neo-Hookean materials, we first formulated the boundary value problem in nonlinear elasticity, considering that the slab attached at one end to a rigid substrate and was traction-free at the other end. Considering the two configurations depicted in figure 1, we have identified their basic undeformed solutions given by equation (2.5), which are characterized by a hydrostatic pressure linearly varying in the thickness direction, possibly encountering a discontinuity across the material interface.

Secondly, we have studied the linear stability by means of the method of incremental deformations superposed on the basic elastic solution. We found that both configurations can undergo a morphological transition governed by the order parameter *γ*, representing the ratio between potential and elastic energies of the top layer. In particular, its critical value depends on three dimensionless parameters: *α*_{H}, *α*_{μ} and *α*_{ρ}, representing the thickness, the shear moduli and the density ratios between the layers, respectively.

Thirdly, we have implemented a finite-element code to solve the boundary value problem in the fully nonlinear instability regime. Other than validating the predictions of the linear stability analysis, the simulations have highlighted the nonlinear evolution of the characteristic patterns.

Compared with the classic Rayleigh–Taylor hydrodynamic instability, not surprisingly we have found that elastic effects tend to stabilize the dynamics of the surface undulations forming beyond the linear stability threshold. Nonetheless, we obtained a rich morphological diagram with respect to both the geometric and elastic parameters. In the following, we briefly discuss the main results for the two cases under consideration.

If the body hangs below a rigid wall, as depicted in figure 1*a*, we find that there always exists a critical value for the order parameter, driving a morphological transition localized at the free surface. Such a shape instability is favoured if the bottom layer is softer and thicker than the top one, having a critical horizontal wavelength of the same order as the body thickness. In the nonlinear regime, this critical undulation evolves towards forming a digitation whose characteristic penetration length continuously increases beyond the linear stability threshold, highlighting the existence of a supercritical pitchfork bifurcation.

If the body is attached to a rigid substrate at the bottom surface, as depicted in figure 1*b*, a morphological transition can occur if and only if the top layer has a higher density than the bottom one. Similarly to the previous case, the onset of an elastic bifurcation is favoured by a softer and thicker bottom layer compared with the top one, with a critical wavelength of the same order as the body thickness. However, an important difference is that the shape instability is localized at the interface between the two layers, displaying two characteristic nonlinear patterns. If the top layer is thinner than the bottom one, the undulation evolves towards forming finger-like protrusions, while in the opposite geometrical limit a stable wrinkling occurs.

In summary, we have characterized the shape instabilities occurring in a soft elastic bilayer subjected only to the action of the gravity bulk force. Unlike the Rayleigh–Taylor instability in fluids, we have demonstrated that the nonlinear elastic effects saturate the dynamic instability of the bifurcated solutions, displaying a rich morphological diagram where both digitations and stable wrinkling can emerge. The results of this work provide important guidelines for the design of novel soft systems with tunable shapes. In fact, the possibility to use external stimuli to control both the geometric and the elastic properties in smart materials, such as hydrogels or dielectric elastomers [26], can be used to provoke morphological transitions on demand [27]. Morphological changes in such soft devices may be employed, for example, to selectively change the surface roughness (e.g. to perform drag reduction in fluid–structure interactions [28]) or to fabricate tailor-made patterns (e.g. to design adaptive material scaffolds [29]).

## Author' contributions

D.R. carried out the stability analysis and the numerical simulations. P.C. conceived of and designed the study. Both authors drafted, read and approved the manuscript.

## Competing interests

The authors declare that they have no competing interests

## Funding

This work has been partially supported by Progetto Giovani GNFM 2016 funded by the National Group of Mathematical Physics (GNFM - INdAM) and by MFAG AIRC grant no. 17412.

## Appendix A. Structure of the matrix *M*

We report the matrix *M* that we used in equation (3.14). We split it into 16 blocks
where **0** is the null 2×2 matrix and

## Appendix B. Expressions of the coefficients *c*_{j}

The coefficient *c*_{1} in equation (3.15) is given by
whereas *c*_{2} is
and the expression of *c*_{3} is

## Footnotes

One contribution of 13 to a theme issue ‘Patterning through instabilities in complex media: theory and applications.’

- Accepted January 4, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.