## Abstract

We consider an inclined rectangular duct of constant cross-section conveying viscous fluid and covered by an elastic plate. The fluid is described by the Stokes equations and the plate by the Föppl–von Kármán equations. The equations admit an equilibrium solution in which the plate is flat and fluid flows underneath due to gravity. This base flow induces a varying traction across the plate, which can lead to out-of-plane buckling due to the associated in-plane shear. Linear stability analysis demonstrates that buckling occurs for sufficiently thin plates on steep slopes and deep channels. The most unstable modes take the form of either symmetric or antisymmetric downslope-directed chevron patterns that travel downstream at a fraction of the average speed of the base flow. An analogous analysis shows that similar buckling also occurs if the elastic plate is replaced by a thin skin of very viscous fluid. Our description provides a simple model for the formation of ropy pahoehoe lava.

## 1. Introduction

The interaction of a fluid with an overlying elastic ‘skin’ (or, indeed, any rheologically distinct superficial layer) can generate a variety of buckling patterns. Common examples include the flow-induced wrinkling of the thin skins formed atop hot milk, wax or crème anglais. Likewise, it has been suggested that motion of the molten interior of a lava flow can buckle the overlying solidified crust, creating a characteristic ‘ropy’ appearance (Fink & Fletcher 1978; figure 1). In these examples, the fluid flow plays an essential role in forming the wrinkles of the skin by exerting the stresses that promote buckling. However, the role played can be more passive, such as in semiconductor manufacturing in which an elastic plate floating on a viscous layer is compressed by external forces, with the fluid merely controlling the time scale of buckling (Huang & Suo 2002).

The goal of the current article is to explore a model problem in which fluid flow directly induces buckling. More specifically, we consider the gravity-driven flow of a viscous fluid through a rectangular duct overlain by a thin elastic plate. Because the flow is driven by gravity alone, without a pressure gradient, an equilibrium state can be established in which the surface skin remains flat. This equilibrium and its linear perturbations can be described with relative mathematical simplicity, in contrast to pressure-driven flows that generally form conduits that taper in the streamwise direction. The fluid traction acting on the base of the skin varies across the duct and generates shear within it, potentially inducing out-of-plane buckling when tractions are sufficiently strong. This configuration is one of the simplest idealizations of the particular fluid–structure interaction problem in question, and bears similarities to the buckling of an isolated elastic plate under gravity or externally imposed shear (Balmforth *et al*. 2008).

To model the configuration, we treat the viscous fluid with the Stokes approximation (valid for relatively slow flows) and represent the elastic plate as a thin Hookean isotropic elastic solid satisfying the Föppl–von Kármán equations (Love 1944). The latter are the simplest plate equations that incorporate both the bending and compressional terms necessary to account for out-of-plane buckling. The choice of this model for the skin is not essential, and we supplement our study of the elastically plated duct with an exploration of a related problem in which the skin is composed of a yet-more viscous, thin immiscible fluid (see also Teichman 2002).

Despite the common occurrence of flow-induced buckling of a superficial skin, relatively few preceding studies exist on the problem. Perhaps the closest earlier work to our current study is by Luo & Pozrikidis (2006, 2007). They considered the buckling of an arbitrarily shaped elastic section in an otherwise rigid plate suspended in flowing fluid. Compression and shear are induced in the elastic section by the flow, resulting in a complicated buckling instability. However, although the elastic portion of the plate deforms out of its plane on buckling, these authors do not incorporate the resulting feedback on the fluid dynamics. This simplifies the analysis, but alters the nature of the instability since fluid pressure variations cannot provide a restoring force.

The paper is structured as follows. In §2, we formulate the problem and note the assumptions required. In §3, we present the flat base state profile about which we perturb. We analyse the linear stability numerically in §4. A brief comparison of this theory with a qualitative experiment is presented in §5. In §6, we consider the particular geological application of ropy pahoehoe lava. Finally, in §7, we summarize our most significant findings. Appendix A furnishes the derivation of the Föppl–von Kármán equations in the present context and outlines their limits of validity. Appendix B briefly discusses the supplementary problem in which the elastic skin is replaced by a very viscous fluid one.

## 2. Formulation

As shown in figure 2, we consider a duct of infinite length and rectangular cross section (having width 2*y*_{0} and depth 2*αy*_{0}, where *α* is the aspect ratio) inclined at angle *θ* to the horizontal. The duct is described by a Cartesian coordinate system with directed downslope, across the slope and perpendicular to the slope; the origin is located on the centre line. The bottom, , and sides, , are rigid; the top, , is an isotropic Hookean elastic plate of thickness 2*d*≪2*y*_{0}, Young's modulus *E* and Poisson ratio *ν*, where is the out-of-plane displacement of the initially flat plate. An incompressible, inertialess Newtonian fluid of density *ρ* and viscosity *μ* fills the duct and flows downslope under gravity.

The fluid motion is described by Stokes' equations(2.1)where , is the stress tensor, *g* is the acceleration due to gravity, is the pressure, is the velocity, is the identity tensor, and † denotes the transpose. On the rigid lower and side surfaces of the duct, we impose no-slip,(2.2)

The strains in the elastic plate are assumed to be small so that the plate may be described by the Föppl–von Kármán equations in Eulerian coordinates. In appendix A, we give a discussion of the underlying assumptions that are required for these equations to apply. We also assume that the weight of the plate is negligible compared with fluid tractions acting upon it and that elastic waves are much faster than the adjustment time scale of the coupled elastic-fluid system so that the plate is in instantaneous equilibrium. Thus, the in-plane and out-of-plane force balances are given by(2.3)respectively, where the subscript *h* denotes in-plane ( and ) components, , and is the fluid traction acting on the base of the plate. Here, the in-plane stresses (integrated over the thickness) and strains are(2.4)respectively, where is the in-plane displacement and tr(.) is the trace. The terms in the second equation of (2.3) represent a resistance to bending, out-of-plane forcing by the fluid and a buckling term coupling in-plane stresses with out-of-plane displacement. The plate is clamped along its lateral edges to the sides of the duct,(2.5)

The plate displacement and fluid velocities are connected through the kinematic boundary condition,(2.6)where is time. We have neglected the and time-derivative terms in the in-plane equations because the in-plane displacements must be small relative to the out-of-plane displacements (essential scalings for deriving the Föppl–von Kármán equations, as given in (A 4) in appendix A) whereas the three velocity components are all of similar order. The traction on the plate is exerted by fluid shear stresses as follows:(2.7)where is the normal to the plate.

### (a) Non-dimensionalization

We non-dimensionalize the governing equations based on the half-width of the channel and the gravitational driving force as follows:(2.8)Substituting (2.8) into (2.1)–(2.7), we obtain(2.9)subject to(2.10)(2.11)This flow field is coupled to the plate through(2.12)with(2.13)subject to(2.14)The dimensionless parameters are the inclination of the duct *θ*, the aspect ratio *α*, the Poisson ratio of the plate *ν*,(2.15)measuring the relative importance of gravity-induced fluid forcing to bending stiffness and . is large when the plate is relatively flimsy, and it is small when the plate is relatively rigid. *G* measures how much strain is introduced in the base state by the in-plane stresses. Because linear instability is controlled by these stresses, rather than the strains, *G* will not feature in the stability analysis.

## 3. Base state

Our base state (denoted by subscript 0) is driven by gravity alone; the plate is flat and the flow is steady and unidirectional. The flow profile is (Rosenhead 1963)(3.1)with *v*_{0}=*w*_{0}=0 and .

The fluid flow imposes a traction on the plate, causing downstream displacement according to(3.2)with *η*_{0}=*ζ*_{0}=0. For convenience, we define the scaled shear stress,(3.3)For shallow ducts (small *α*), the flow becomes independent of *y* away from the edges; the corresponding shear profile is . For deep ducts (large *α*), the flow close to the plate becomes independent of the depth of the duct and the corresponding shear profile becomes independent of *α*. Sample velocity profiles and displacements for the base state are shown in figure 3.

## 4. Linear stability

To examine the stability of the base state, we perturb each of the variables about that equilibrium, *f*=*f*_{0}+*f*_{1}, where *f* denotes any of *u*, *v*, *w*, *p* and *ζ*, and the subscript 1 refers to the perturbation. We then linearize in the perturbation amplitudes *f*_{1} and calculate the growth rate, , of infinitesimal normal mode perturbations with downstream wavenumber *k*: . This furnishes the eigenvalue problem(4.1)(4.2)(4.3)(4.4)with boundary conditions(4.5)(4.6)where(4.7)and(4.8)Note that the in-plane displacement equations are no longer needed, and that we have used the in-plane force balance (the first equation of (2.12)) to replace with (1/2)∂/∂*y* in the kinematic condition (the first equation of (4.6)).

The out-of-plane plate equation (4.7) regulates instability: the first three terms describe the stabilizing effects of bending stiffness, fluid pressure and hydrostatic restoring force, respectively, whereas the final term describes the destabilizing effect of differential fluid traction on the base of the plate.

We solved the systems (4.1)–(4.8) numerically using a Chebyshev collocation scheme (Trefethen 2000) with the boundary conditions imposed explicitly (Weideman & Reddy 2000) and *u*_{1} and *p*_{1} eliminated. The resulting generalized eigenvalue problem was solved using the QZ algorithm. We found that employing an equal number of collocation points (between 24 and 32) in each direction yielded reasonable convergence rates for most values of *α*.

The dispersion relations for the four least-stable modes for a square duct on a 45° slope (*α*=cot *θ*=1) and are plotted in figure 4; associated contour plots of out-of-plane displacements for a selection of wavenumbers are provided in figure 5. For *k*→0, the dominant mode *E*_{2} corresponds to over- or under-filling the channel uniformly in *x* and travels faster than the base state flow. At intermediate and large wavenumbers, an even *E*_{1} and odd *O*_{1} pair of modes dominate, and for a range of intermediate wavenumbers both are unstable. These modes are destabilized by the traction exerted by the fluid on the base of the plate. They are stabilized at large wavenumber by bending stiffness and at small wavenumber by a combination of cross-stream bending stiffness and the restoring pressure force of the fluid. The out-of-plane displacements for these modes take the form of patterns of downslope-directed chevrons with crests having an appreciable angle to the axis of the channel, despite the base displacement being small (cf. the shear-induced buckles of an isolated elastic plate; Mansfield 1964; Wong & Pellegrino 2006; Balmforth *et al*. 2008). At smaller *k*, the odd mode (figure 5*b*) with fewer cross-slope oscillations is favoured, while at larger *k* the effect of shear is increased and the even mode (figure 5*a*), with more oscillations in *y* (augmenting the destabilizing shear), is favoured. At large wavenumber, the eigenfunctions become concentrated near the sidewalls where there is maximum shear. Little deflection occurs near the centre-line, with the result that the flows to either side of the channel do not interact and the growth rates of the dominant odd and even pair become indistinguishable. These unstable waves propagate downslope at a fraction of the average base state velocity. A multitude of other modes exist with decreasing growth rates, as illustrated by the additional odd mode *O*_{2} in figures 4 and 5*d*. These modes are similar to the dominant modes, but with an increasing number of oscillations in *y*.

The growth rates and phase speeds of the most unstable mode are shown on the (*k*, )-plane for the same square duct (with *α*=cot *θ*=1) in figure 6. As increases, the resistance of the plate to bending decreases and the flat base state becomes unstable above a critical value over an expanding window of wavenumbers. The fourth-order form of the bending stiffness ensures that this term is large even for moderate wavenumbers and cross-stream scales of variation, so the plate must be very flimsy ( must be very large) to ensure instability. Beyond , the most unstable modes are found at increasingly shorter wavelengths and have slowly decreasing wavespeeds. A selection of neutral stability curves for different slope angles and aspect ratios is shown in figure 7. For steeper slopes, onset occurs for smaller driving forces and longer waves, because the component of gravity driving the underlying shear flow is increased and, to a lesser extent, because the stabilizing hydrostatic pressure is reduced. As the channel becomes deeper, larger velocities can be generated with associated larger shears on the base of the plate. Onset thus occurs at smaller driving force and for longer waves. However, the behaviour becomes independent of the aspect ratio once *α* exceeds two or so.

The behaviour of the critical value as a function of slope angle and aspect ratio is shown in figure 8. Figure 8*a* also displays contours of growth rate and wavenumber for the most unstable mode above onset with *α*=1.

Finally, figure 9 presents mode profiles for the most unstable mode (here *E*_{1}) and the most unstable wavenumbers for the -values marked in figure 6. The characteristics of the out-of-plane displacement are identical to those described in figure 7. The figure also highlights the increasing concentration of the eigenfunctions to the regions near the sidewalls with increasing . At onset, the perturbation to the velocity field is appreciable throughout the depth of the channel, but when the growth rates are larger only the uppermost fluid layer responds to the plate motion.

### (a) Local analysis

A crude short-wavelength-style analysis can be used to complement the numerical stability theory. We assume that the perturbations vary more rapidly in *y* than the background state and decompose the perturbations as . We then integrate the fluid equations through the depth of the flow and impose the boundary conditions at *z*=±*α* to obtain a ‘dispersion relation’,(4.9)where and . Because is a function of *y*, this relation can only be viewed as a short-wavelength approximation, suitable when the length scale *m*^{−1} is much shorter than the channel width.

To estimate the neutral stability curves, we look for the smallest value of that gives *ω*_{i}=0 over all possible *m* and *y*. This minimization requires the value of , which is plotted in figure 10. Taking the value at the sidewalls (although the short-wavelength form does not satisfy the boundary conditions there) provides a crude lower bound on the neutral stability curves: solutions of this algebraic problem are included in figures 6 and 7 for comparison with the numerical data. It performs best at the lowest and the highest wavenumbers where the eigenfunction possesses short-wavelength structure in *y* (as demanded by the local theory). For small *k* and large , *m*∼*k*^{−1} and the approximate long-wave cut-off isreflecting an out-of-plane balance of forces between cross-stream bending stiffness, pressure and shear. At large *k*, *m*∼*k* and the short-wave cut-off is given byvia a balance between bending stiffness and shear. The critical value of is poorly predicted by the local analysis because *m* is not large there.

For *k*≫1, the local analysis also predicts that the maximum growth rate occurs for *m*=−*k* sign (i.e. wavy perturbations with crests aligned at 45° and orientated to take advantage of the destabilizing effect of the local shear), implying a most unstable wavenumber,which is also included in figure 6. Overall, the local analysis offers a useful guide to the stability characteristics, but is only quantitatively accurate at large *k*.

## 5. Qualitative experiments

To provide qualitative verification of the theoretical predictions, we performed a suite of experiments. A 7.6 cm wide, 2.6 cm deep and 1.2 m long channel was mounted on an inclinable table (the achievable fluid flux effectively limited us to an operating range of (0°, 35°)). ‘Thera-Band’ latex exercise bands were placed on top of the channel to form an elastic skin. Four bands of different thickness were used and were clamped into place using slats. An elongation test indicated that Young's modulus for all bands was in the range 2±0.4 GPa. The sheets used were colour coded: tan (2*d*=0.12±0.02 mm), yellow (2*d*=0.15±0.02 mm), red (2*d*=0.20±0.02 mm) and green (2*d*=0.25±0.02 mm). At the upper end of the channel, a reservoir and a solid top plate were inserted (the latter to help generate the desired base state). Golden syrup (*ρ*=1.4 g cm^{−3}; *μ*=2.6±0.2 or 3.4±0.2 Pa s at 20°C) was supplied to the reservoir, and the flux through the duct was controlled manually and made as constant as possible using a valve. For *θ*≲10°, experiments could be maintained for several minutes with the help of a pump; for *θ*>10°, experiments could be maintained for at least 30 s. Still images of a deflected laser beam at 1/3 s intervals permitted measurements of the wavelength and wavespeed. The experimental set-up is shown in figure 11*a*.

With this arrangement, we successfully observed wrinkling of the latex sheet due to the underlying fluid flow: a sample pattern is shown in figure 11*b*. Moreover, the wrinkles propagated downstream at speeds that were a small fraction of the average fluid speed. Overall, the observed properties of the wrinkling patterns and the trends on varying the inclination and the thickness of the latex sheet were consistent with the theoretical predictions: for experiments performed at angles below a critical value, the sheet appeared flat; at higher inclination, buckles were observed. As expected from theory, the critical value was of the order of tens of degrees and increased with increasing thickness of the elastic sheet. Similarly, the observed wavelengths were of the order of centimetres and decreased with increasing slope and with decreasing thickness of the elastic sheet, whereas the wavespeed increased with slope. Distinct chevron patterns much like the linear eigenfunctions were observed, although the maxima were generally located near the edges and rarely along the centre line. Frequently, the maxima on either side were not in phase suggesting a competition between even and odd modes of instability, exactly as in the theory.

Despite this qualitative success, we were unable to provide a quantitative validation of the theory because we found significant problems regarding reproducibility: experiments with apparently identical set-up had variations as much as 10° in onset angle. We believe that the problem originates from some initial compression or tension of the sheet, which is almost impossible to avoid and is unquantifiable. A similar, although less serious, issue arose in our previous experiments with an isolated sheared elastic sheet (Balmforth *et al*. 2008). There, theory predicted that displacements of the plate in *y* by fractions of a millimetre could appreciably shift onset. We therefore anticipate that the presence of tiny amounts of compression or tension can significantly alter the shear instability, with compression accentuating it and tension eliminating it.

## 6. Geological application

The surface of ‘ropy pahoehoe’ lava flows has a characteristic corrugated appearance consisting of folds a few centimetres in amplitude and a few centimetres to tens of centimetres in wavelength (Fink & Fletcher 1978), resulting from an imposition of stresses on the cooled, rheologically distinct upper crust. Based on field studies of solidified flows and movies of active ones, Fink & Fletcher (1978) suggest that wrinkles are formed where a flow encounters a constriction or a sudden change in slope, and the faster moving lava advects upstream crust into crust near the obstacle to compress and buckle it. Beyond the generation point, the flow rotates the wrinkles into their characteristic parabolic shape.

A possible alternative mechanism is suggested by the current study: the corrugations are a manifestation of shear-induced wrinkling. Although our model is a crude approximation of a real lava flow, it does capture the most fundamental features. Regular wrinkles are particularly formed on channelized lava flows bordered by ‘levees’ (Fink & Fletcher 1978; Garry *et al*. 2006), providing some justification for our geometry. From a rheological perspective, although lava is a heterogeneous fluid with significant non-Newtonian, temperature-dependent rheology (Griffiths 2000), a number of observations indicate that our approximation is reasonable at leading order. First, the formation of the upper boundary layer (solid or otherwise) insulates the underlying flow and the temperature in the interior is essentially constant and uniform (Hon *et al*. 1994). Second, ropes form preferentially relatively close to the vent where the bubble and crystal contents are comparatively low, and so the interior fluid may be approximately Newtonian. The overlying crust is observed to be a visco-elastic shell, with a brittle upper casing and a high-viscosity lower cushion (Hon *et al*. 1994). Our elastic plate provides one possible idealization of this crust (cf. Iverson 1990); we also provide a corresponding analysis for a very viscous fluid plate in appendix B (cf. Biot 1961; Fink & Fletcher 1978). The two models share many common features and we concentrate on the former here. We note that the most immediate shortcoming of the model is that it enforces conservation of crust material, whereas stretching and fracture in conjunction with solidification of freshly exposed lava can generate new crust in response to shear. Fracture at the plate edges also allows the underlying flow to advect the crust, modifying the basic flow and the lateral boundary conditions on the plate. This advection modifies the details of the base flow and the lateral boundary conditions on the plate; our model is not directly applicable to this case.

For the lava flows observed by Fink & Fletcher (1978), *θ*=5°, 2*y*_{0}≳2.5 m, 2*d*≲5 cm and the wavelength is of the order of 10 cm. Taking representative values *α*=0.5 (Calvari *et al*. 1994), *E*=95 GPa, *ν*=0.27 (Schilling *et al*. 2003) and *ρ*=3000 kg m^{−3}, we predict that a least-stable wavelength of 10 cm requires a skin thickness of order 0.2 mm. This prediction is appreciably thinner than estimated for the actual flow, but nevertheless is geologically reasonable (Hon *et al*. 1994). Thus, buckling by shear alone is a possible explanation for ropy pahoehoe formation, although the theoretically predicted requirements are at the limits of what is observed. The mode shapes are also in qualitative agreement with the field observations (compare figures 1 and 9).

It is also worth noting that there is a considerable experimental literature on laboratory analogue flows using cooling PEG wax (Griffiths *et al*. (2003) and Garry *et al*. (2006) and references therein). Qualitatively, the features produced by these extruded, cooling flows of wax mimic those of real lava flows: they form levees, become channelized, and once a solidified skin forms may develop surface texture and regular folding. A second laboratory analogue flow is the intrusion of a salty gravity current into a less dense surfactant solution. At the interface between the two fluids, a micellar gel is formed that is contorted by the dynamics of the underlying flowing fluid in a manner resembling an active lava flow (Clayton 2004). Our predicted mode shapes are in qualitative agreement with both types of laboratory flows.

## 7. Conclusions

In this article, we have demonstrated that an elastic plate can be buckled through traction exerted by an underlying fluid shear flow. More precisely, we have explored linear buckling instabilities in an elastic plate clamped over an inclined duct filled with fluid flowing under gravity alone. The varying fluid traction across the top of the duct leads to a destabilizing in-plane shear in the plate, and we have mapped out the conditions under which this effect drives buckling and classified its character. Instability is most readily observed for thin sheets on steep slopes and deep channels, and the most unstable modes take the form of downslope-directed chevron patterns. These patterns resemble the structures seen on the crusts of pahoehoe lava flows, leading us to explore whether crustal shear due to underlying lava flow could be responsible for their formation.

We also successfully conducted a suite of simple experiments to confirm the theory qualitatively. Unfortunately, the experiments were plagued by technical complications. In particular, an uncontrollable amount of slack or tension was probably introduced in the elastic plate when it was initially clamped into place. The resulting lateral compression swamps the shear-induced instability and prevents any quantitative comparison with theory. To surmount this issue, a more sophisticated experiment is needed, which we leave for future work.

Although our main focus has been the buckling of an elastic plate, we have also given a brief discussion of the corresponding problem when the fluid is coated by a thin film of much more viscous fluid. Such a skin can also be buckled by an underlying fluid shear flow, giving another example of G. I. Taylor's analogy between an elastic plate and a viscous sheet.

## Acknowledgments

We thank Mike Hoog for generously providing a sample of Thera-Band, Mark Martinez for generous use of his laboratory facilities and Ron Schott for kind permission to use the photograph in figure 1. R.V.C. thanks the Departments of Mathematics at the Universities of Alberta and British Columbia for their hospitality while this work was undertaken. A.C.S. thanks the Killam Foundation for support. This work was begun at the 2004 Geophysical Fluid Dynamics Summer Program (Woods Hole Oceanographic Institution), which is supported by the National Science Foundation and the Office of Naval Research; we thank the participants for useful discussions.

## Footnotes

- Received April 5, 2008.
- Accepted August 15, 2008.

- © 2008 The Royal Society