Fluttering-induced flow in a closed chamber

Abstract

We study the emergence of fluid flow in a closed chamber that is driven by dynamical deformations of an elastic sheet. The sheet is compressed between the sidewalls of the chamber and partitions it into two separate parts, each of which is initially filled with an inviscid fluid. When fluid exchange is allowed between the two compartments of the chamber, the sheet becomes unstable, and its motion displaces the fluid from rest. We derive an analytical model that accounts for the coupled, two-way, fluid–sheet interaction. We show that the system depends on four dimensionless parameters: the normalized excess length of the sheet compared with the lateral dimension of the chamber, $\varDelta$; the normalized vertical dimension of the chamber; the normalized initial volume difference between the two parts of the chamber, $v_{du}(0)$; and the structure-to-fluid mass ratio, $\lambda$. We investigate the dynamics at the early times of the system's evolution and then at moderate times. We obtain the growth rates and the frequency of vibrations around the second and the first buckling modes, respectively. Analytical solutions are derived for these linear stability characteristics within the limit of the small-amplitude approximation. At moderate times, we investigate how the sheet escapes from the second mode. Given the chamber's dimensions, we show that the initial energy of the sheet is mostly converted into hydrodynamic energy of the fluid if $\lambda \ll 1$ and into kinetic energy of the sheet if $\lambda \gg 1$. In both cases, most of the initial potential energy is released at time $t_{p}\simeq \ln [c \varDelta ^{1/2}/v_{du}(0)]/\sigma$, where $\sigma$ is the growth rate and $c$ is a constant.

1. Introduction

Many natural processes and technological applications rest on fluid–structure interactions to maintain their regular functionality. Of particular interest are the mutual interactions between slender elastic objects and a fluid medium that trigger elastohydrodynamic instabilities. Such instabilities are vital for the control, for example, of the passage of air through the lungs (Ishizaka & Flanagan 1972; Grotberg & Jensen 2004), the directionality of blood flow (Pedley, Brook & Seymour 1996) and the blood pressure of tall animals (Pedley et al. 1996). Moreover, bending deformations of slender objects in viscous or inertial fluids have been manipulated for applications in soft robotics (Kim, Laschi & Trimmer 2013; Matia & Gat 2015; Rothemund et al. 2018), the fabrication of microfluidic soft actuators (Thorsen, Maerkl & Quake 2002; Hosoi & Mahadevan 2004; Holmes et al. 2013; Fargette, Neukirch & Antkowiak 2014; Gomez, Moulton & Vella 2017; Christov et al. 2018; Boyko et al. 2019; Jiao & Liu 2021), the manufacture of semiconductors (King 1989) and the design of soft and active matter through catalytic reactions (Laskar et al. 2022; Manna et al. 2022) and dynamical wrinkles (Chopin, Dasgupta & Kudrolli 2017; Kodio, Griffiths & Vella 2017; Box et al. 2019; Pocivavsek et al. 2019; O'Kiely et al. 2020; Diamant 2021; Guan et al. 2022, 2023).

Despite recent achievements, novel designs of small-scale devices still call for a deeper understanding of elastohydrodynamic couplings. One such design was recently introduced by Oshri (2021). In that set-up, a thin sheet is compressed between the two sides of a closed chamber and divides it into two separate parts that are connected by a valve (figure 1). At time $t< 0$, the valve is closed, and each part of the chamber is filled with an incompressible fluid. In the absence of fluids, the sheet would have accommodated its minimum energetic state, i.e. the lowest mode of buckling, but in the presence of fluids, the sheet is forced to accommodate a higher energetic state. The additional energy can be exploited to displace the fluid from rest, if, for example, the valve is opened to allow the transfer of fluids between the two compartments of the chamber.

Figure 1. Schematic overview of the system. A thin sheet of total length $\tilde {L}$, bending modulus $\tilde {B}$, density $\tilde {\rho }_{sh}$ and thickness $\tilde {h}$ divides a closed rectangular chamber of dimensions $\tilde {L}_x\times \tilde {L}_y$ into two parts. The excess length of the sheet compared with the lateral dimension of the chamber is given by $\tilde {\varDelta }=\tilde {L}-\tilde {L}_x$ (not shown in the figure). The volumes of the chamber above and below the sheet, $\tilde {v}_{i}(t)$ ($i=u,d$), are filled with an inviscid and irrotational fluid of density $\tilde {\rho }_{\ell }$. At $\tilde {t}\geq 0$, fluid is allowed to exchange freely between the two compartments of the chamber. In our formulation, the fluid exchange occurs through the upper and lower walls of the chamber (represented by dashed-dotted blue lines). To model this exchange, we apply periodic boundary conditions along these walls. One possible experimental set-up that corresponds to the above model involves a valve-controlled channel that connects the two compartments of the chamber.

In the above mentioned study, Oshri (2021) analysed the quasi-static evolution of the system, wherein the volume of fluid exchanged between the two parts of the chamber is the control parameter. In contrast, the present work focuses on the dynamical evolution of the system, wherein the fluid is driven by the spontaneous relaxation of the sheet from higher to lower energetic states. We believe that the dynamic analysis of this set-up will open new avenues for designing advanced technological devices, such as micromechanical switches (Krylov et al. 2008; Zhang et al. 2014; Preston et al. 2019) and microfluidic mixing devices (Stroock et al. 2002; Liu, Kim & Sung 2004; Lee et al. 2011). Indeed, the additional coupling between the sheet and the surrounding fluid confers increased flexibility in the design of such switches. Different fluids with different viscosities can be used to manipulate the time that is required for the sheet to release its stored energy, thereby increasing, for example, the time scales over which the switches operate. In addition, when the two parts of the chamber are filled with different fluids, the elastic energy released from the sheet can be exploited for mixing: the pressure field induced in the chamber can be utilized to inject the fluid from one side of the chamber into the fluid on the other side, thereby inducing mixing of the two fluids. Typically, such devices function in conditions of low Reynolds numbers, where the effects of viscosity are significant. However, our system can also be applied in the design of pneumatic time-delay switches and soft pneumatic actuators (Rothemund et al. 2018; Preston et al. 2019; Drotman et al. 2021), which typically operate in the opposing limit of high Reynolds numbers.

While successful implementation of these applications is in itself a challenging task (which we plan to pursue in future research), in this work, we aim to answer more fundamental questions related to the underlying physical behaviour of the system. For example, how much of the initial elastic energy is subsequently transferred from the sheet to the fluid? How is the velocity of the fluid that is induced in the chamber related to the elastic properties of the sheet? What is the maximum pressure difference that the sheet induces in the chamber?

As a first step to answering these questions, we derive an analytical model that encompasses the elasticity of thin sheets and the hydrodynamics of inviscid fluids. Our model reveals that the system depends on four dimensionless parameters: the normalized excess length of the sheet compared with the lateral dimension of the chamber, $\varDelta$, where the total length of the sheet is used to normalize all lengths; the normalized vertical dimension of the chamber, $L_y$; the normalized initial volume difference in the chamber, $v_{du}(0)$; and the structure-to-fluid mass ratio, $\lambda$. We show that for fixed dimensions of the chamber, $L_y$ and $\varDelta$, the system exhibits two asymptotic solutions as a function of $\lambda$. The sheet's inertia dominates the dynamics when $\lambda \gg 1$, and is therefore referred to below as the ‘solid-dominated’ region, while the dynamics is governed by the fluid's inertia when $\lambda \ll 1$, and is therefore referred to as the ‘fluid-dominated’ region.

We investigate the system's behaviour both in the early stages of its evolution and at moderate times during which nonlinear effects control the dynamics. For the early stages, we employ linear stability analysis around the (unstable) second buckling mode and the (stable) first buckling mode. We obtain the highest growth rate, $\sigma$, and the lowest frequency of vibration, $\omega$, around these initial states. The two solutions exhibit similar behaviour as a function of $\lambda$, namely they converge to a constant in the solid-dominated region, while they exhibit the scaling $\lambda ^{1/2}$ in the fluid-dominated region. Furthermore, we show that in the solid-dominated region only one mode of the sheet is essentially excited at the instability, while an infinite number of modes are excited in the fluid-dominated region. Analytical approximations are derived for each of these cases under the assumption that the amplitude of the sheet remains small, i.e. $\varDelta \ll 1$ (Landau & Lifshitz 1986).

At moderate times, a weakly nonlinear analysis is performed around the second buckling mode. Given a small initial volume difference between the upper and lower parts of the chamber, we analyse the dynamic evolution of the system up to the peak time $t_{p}$, at which the sheet releases most of its initial potential energy. We show that, after some initial delay, the sheet rapidly escapes from the unstable state. We derive the approximation $\sigma t_{p}\simeq \ln [c\varDelta ^{1/2}/v_{du}(0)]$, where $\sigma$ is the growth rate of the linear instability and $c$ is a constant, and show that it agrees well with the numerical results. At $t_{p}$, most of the initial potential energy is converted into kinetic energy of the sheet if $\lambda \gg 1$ and into hydrodynamic energy if $\lambda \ll 1$. We show that at $t=t_{p}$ a relatively large spike of pressure drop is applied on the sheet.

The paper is organized as follows. In § 2, we first formulate the problem for finite excess lengths. Then, we reduce this formulation to the small-amplitude approximation and introduce the modal expansion of the solution. In § 3, we investigate the early stages of the evolution. After recalling the static solution, we employ a linear stability analysis around the second and the first modes of buckling. In § 4, we investigate the system's evolution at moderate times. In particular, we examine the energetic interplay between the sheet and the fluid, derive the scaling for the peak time, $t_{p}$, and explore the relation between the volume difference and the pressure drop on the sheet. Finally, in § 5, we discuss a possible experimental realization of the system and, in § 6, we draw conclusions, and propose a direction for future study.

2. Formulation of the problem

We consider an inextensible thin sheet of total length $\tilde {L}$, bending modulus $\tilde {B}$, thickness $\tilde {h}$ and density $\tilde {\rho }_{sh}$. The sheet divides a rectangular closed chamber into two parts, which are connected by a valve (figure 1). The lateral, the vertical and the width dimensions of the chamber are denoted by $\tilde {L}_x$, $\tilde {L}_y$ and $\tilde {W}$, respectively. A Cartesian coordinate system is located on the left edge of the sheet. A cross-section of the chamber on the $\tilde {x}\tilde {y}$ plane is placed at $0\leq \tilde {x}\leq \tilde {L}_x$ and $-\tilde {L}_y/2 \leq \tilde {y}\leq \tilde {L}_y/2$. When $\tilde {t}<0$, the valve connecting the two parts of the chamber is closed, and the volumes above and below the sheet, $\tilde {v}_{u}(\tilde {t})$ and $\tilde {v}_{d}(\tilde {t})$, are filled with an incompressible, inviscid fluid of density $\tilde {\rho }_{\ell }$. Hereafter, we denote quantities related to the upper and lower parts of the chamber by the subscripts ‘$u$’ and ‘$d$’, respectively. At $\tilde {t}\geq 0$, the valve is opened, and free exchange of fluid is allowed in the chamber.

In the analysis that follows, we normalize all lengths by the total length of the sheet, $\tilde {L}$, and we normalize time by the inertial time scale of the sheet $\tilde {t}_{\star }=\tilde {L}^2(\tilde {\rho }_{sh} \tilde {h}/\tilde {B})^{1/2}$, i.e.

(2.1a–d)\begin{equation} t=\tilde{t}/\tilde{t}_{{\star}},\quad x=\tilde{x}/\tilde{L},\quad L_x=\tilde{L}_x/\tilde{L},\quad v_{d}(t)=\tilde{v}_{d}(\tilde{t})/\tilde{L}^3,\quad \text{etc.} \end{equation}

We choose this normalization because we anticipate that the wavelengths on the sheet will scale with the sheet's total length. In addition, since the dynamics in the system is driven by the sheet's motion, we chose the sheet's inertial time scale for the normalization. Note that our normalization with respect to lengths and time implies the normalization of the hydrodynamic fields and of the elastic fields, as will be emphasized further during the formulation. Hereafter, we denote all dimensional quantities with tildes over the symbols, and the corresponding non-dimensional quantities without a tilde.

Our model is based on the following six assumptions. Firstly, we assume that the system remains uniform along the width dimension of the chamber. Therefore, we set $W=1$ and consider a two-dimensional system. Secondly, we assume that the volume occupied by the elastic sheet is negligible compared with the total volume of the chamber, i.e. $\tilde {h}\tilde {L}/(\tilde {L}_x \tilde {L}_y) \ll 1$, and as a result $v_{u}(t)+v_{d}(t)=L_x L_y$. Thirdly, we assume that the fluid exchange between the two parts of the chamber occurs through the upper and lower walls, i.e. the walls located at $y=\pm L_y/2$. Fourthly, we assume that the vertical dimension of the chamber, $L_y$, is larger than the typical length scale, $\ell$, over which the disturbances in the flow caused by the sheet's motion decay to zero. In addition, we assume that there is no contact between the sheet and the sidewalls of the chamber, or of the sheet with itself, at any time during the system's evolution. Lastly, we assume that at $t=0$ the system is at rest and that the sheet accommodates a configuration that is dictated by the volume difference $v_{du}(0)=v_{d}(0)-v_{u}(0)$.

For an inviscid and irrotational fluid, the state of the flow is determined by four fields. Two of these are the fluid's potential functions $\phi _i(x,y,t)$, where $i=u,d$, from which we can determine the velocity profile of the fluid as ${\boldsymbol {v}}_i=\boldsymbol {\nabla } \phi _i$, where $\boldsymbol {\nabla }$ is the two-dimensional gradient operator. The other two fields that characterize the flow are the pressures $p_i(x,y,t)$ in each side of the chamber. Using our normalization convention, we find that the potential functions are normalized by $\phi _i=\tilde {\phi }_i(\tilde {\rho }_{sh} \tilde {h}/\tilde {B})^{1/2}$ and the pressures by $p_i=\tilde {p}_i\tilde {L}^3/\tilde {B}$. The evolution of these hydrodynamic fields, in space and over time, is determined by the continuity equation and Bernoulli's equation:

(2.2a)\begin{gather} \nabla^2 \phi_i=0, \end{gather}

(2.2b)\begin{gather}\lambda p_{i}+\frac{\partial \phi_i}{\partial t}+\frac{1}{2}|\boldsymbol{\nabla} \phi_i|^2 =c_i(t), \end{gather}

where $c_i(t)$ are arbitrary functions that depend on time. Throughout the system's development, these functions are employed to maintain a constant pressure at a point within each part of the chamber (Lamb 1945). In addition, in (2.2b) and (2.5), we define the dimensionless parameters:

(2.3a,b)\begin{equation} \lambda=\frac{\tilde{\rho}_{sh}\tilde{h}}{\tilde{\rho}_{\ell}\tilde{L}},\quad \varDelta=1-L_x. \end{equation}

The structure-to-fluid mass ratio, $\lambda$, accounts for the ratio between the densities of the sheet and the fluid and the slenderness of the sheet. This dimensionless parameter plays a role, for example, in the problem of a flag flapping under a uniform axial flow (Argentina & Mahadevan 2005; Connel & Yue 2007; Alben 2008; Alben & Shelley 2008). The parameter $\varDelta$ accounts for the difference between the total length of the sheet and the lateral dimension of the chamber. In dimensional form, it may be expressed as $\tilde {\varDelta }=\tilde {L}-\tilde {L}_x$. For a given system, the parameters $\lambda$ and $\varDelta$ remain constant throughout the dynamic evolution.

To solve the continuity equation (2.2a), we must first specify the boundary conditions on the chamber's walls and the fluid–sheet interfaces. Since the fluid that exits the upper wall of the chamber enters through the lower wall, we set periodic boundary conditions through $y=\pm L_y/2$. Thereafter, we ensure that there is no penetration of fluid through the sidewalls of the chamber. These restrictions give the boundary conditions:

(2.4a)\begin{gather} \phi_{u}\left(x,\frac{L_y}{2},t\right)=\phi_{d}\left(x,-\frac{L_y}{2},t\right) \quad \text{and}\quad \frac{\partial\phi_{u}}{\partial y}\left(x,\frac{L_y}{2},t\right)= \frac{\partial\phi_{d}}{\partial y}\left(x,-\frac{L_y}{2},t\right), \end{gather}

(2.4b)\begin{gather}\frac{\partial \phi_i}{\partial x}(0,y,t)=\frac{\partial \phi_i}{\partial x}(1-\varDelta,y,t)=0. \end{gather}

In addition to the periodic boundary conditions at $y=\pm L_y/2$, it is necessary to ensure that $p_{u}(x,L_y/2,t)=p_{d}(x,-L_y/2,t)$ along these walls. By utilizing Bernoulli's equation (2.2a), and the periodic boundary conditions, it becomes apparent that this requirement is satisfied when $c_{d}(t)=c_{u}(t)\equiv c(t)$. Consequently, we can determine the function $c(t)$ by fixing the pressure at a specific point in the lower part of the chamber. In the following analysis we choose

(2.5)\begin{equation} p_{d}\left(\frac{1-\varDelta}{2},-\frac{L_y}{2},t\right)=0. \end{equation}

Two sets of equations model the contact between the sheet and the fluid. The first set corresponds to the kinematic boundary conditions that ensure continuous contact between the sheet and the fluid. The second set corresponds to the force balance equations on the sheet that ensure proper transfer of the momentum between the solid and the fluid. To obtain these two sets of equations, we first define the elastic fields that describe the position of the sheet on the $xy$ plane. It is important to note that, as we assumed the sheet to be inextensible, the elastic model accounts only for bending deformations and does not include stretching deformations. In contrast to the Eulerian description of the fluid, it is convenient to adopt a Lagrangian description for the sheet and to define the elastic fields as functions of the normalized arclength parameter on the sheet, $s\in [0,1]$. With this change of reference frame, we define the position vector to a point on the sheet as ${\boldsymbol {x}}_{sh}(s,t)=(x_{sh}(s,t),y_{sh}(s,t))$ and the angle between the tangent to the sheet and the $x$ axis as $\theta (s,t)$ (see figure 1). These three elastic fields, i.e. $x_{sh}(s,t)$, $y_{sh}(s,t)$ and $\theta (s,t)$, are not independent, since they are related by the geometric constraints

(2.6a)\begin{gather} \frac{\partial x_{sh}}{\partial s}=\cos\theta, \end{gather}

(2.6b)\begin{gather}\frac{\partial y_{sh}}{\partial s}=\sin\theta. \end{gather}

By using these definitions, the kinematic boundary conditions on the sheet–fluid interfaces are given by

(2.7a,b)\begin{equation} y=y_{sh}(x(s),t),\quad \frac{{\rm D} y_{sh}}{{\rm D} t}= \frac{\partial \phi_i}{\partial y}, \end{equation}

where ${\rm D}/{\rm D}t=\partial /\partial t +{\boldsymbol {v}}_i \boldsymbol {\cdot } \boldsymbol {\nabla }$ is the two-dimensional convective derivative. The balance of moments and forces on the sheet is given by

(2.8a)\begin{gather} \frac{\partial^2\theta}{\partial s^2}={-}F_x \sin\theta+F_y \cos\theta, \end{gather}

(2.8b)\begin{gather}\frac{\partial^2 {\boldsymbol{x}}_{sh}}{\partial t^2}={-}\frac{\partial {\boldsymbol{F}}}{\partial s} +\left[p_{d}(x_{sh},y_{sh},t)-p_{u}(x_{sh},y_{sh},t)\right]{\hat{\boldsymbol{n}}}_{d}, \end{gather}

where ${\boldsymbol {F}}=(F_x(s,t),F_y(s,t))$ is the vector of reaction forces per unit length at a cross-section of the sheet and our normalization implies that ${\boldsymbol {F}}=\tilde {{\boldsymbol {F}}}\tilde {L}^2/\tilde {B}$. In addition, ${\hat {\boldsymbol {n}}}_{d}=(-\sin \theta,\cos \theta )$ is a unit normal vector on the sheet that points outwards from the lower part of the chamber, and the hydrodynamic pressures in (2.8b) are calculated on their respective sides of the sheet–fluid interfaces. Note that in (2.8a) we neglect the rotational inertia term. This is justified in the limit of a thin and inextensible sheet, as assumed in this analysis (Neukirch et al. 2012; Goriely 2017; Kodio, Goriely & Vella 2020). Equations (2.8) are supplemented by the following boundary conditions on the sheet's edges:

(2.9a)\begin{gather} x_{sh}(0,t)=0,\quad x_{sh}(1,t)= 1-\varDelta, \end{gather}

(2.9b)\begin{gather}y_{sh}(0,t)=0,\quad y_{sh}(1,t)=0, \end{gather}

(2.9c)\begin{gather}\frac{\partial \theta}{\partial s}(0,t)=0,\quad \frac{\partial \theta}{\partial s}(1,t)=0, \end{gather}

where we assume hinged boundary conditions in (2.9c).

This completes the formulation of the problem. In summary, given the excess length $\varDelta$, the vertical dimension of the chamber $L_y$, the parameter $\lambda$ and the initial volume difference in the chamber $v_{du}(0)$, the dynamic evolution of the system is determined from the solution of the coupled equations (2.2)–(2.9). In the analysis that follows, we always assume that the sheet and the fluid are initially at rest, i.e. $({\partial {\boldsymbol {x}}_{sh}}/{\partial t}) (s,0)=0$ and $\phi _i(x,y,0)=0$.

While solutions to our set of nonlinear equations can, in practice, be sought only numerically, some analytical progress that sheds light on the underlying physics of the system can be achieved under the assumption that the excess length remains small, i.e. $\varDelta \ll 1$. For this reason, in the next section, we reduce our model to this so-called small-amplitude approximation (Landau & Lifshitz 1986) and exploit this formulation to study the time-dependent behaviour of the system.

However, before we proceed to the next section, we should add a comment regarding the system's energy. Since we assumed an ideal fluid, i.e. one without viscous dissipation, and since we consider an elastic model, our equations have a conserved first integral that corresponds to the system's total energy. In accordance with Appendix A, it can be shown that the total energy of the system is given by the sum of the energies of the sheet and the fluid, $E=E_{sh}(t)+E_{f}(t)$, where $E_{sh}(t)$ accounts for the sum of the kinetic and the potential energies of the sheet, which are designated $E_{sh}^{p}(t)$ and $E_{sh}^{k}(t)$, respectively, and $E_{f}(t)$ accounts for the kinetic energy of the fluid. Therefore, the total energy is given by

(2.10)\begin{equation} E=\frac{1}{2}\int_0^1 \left[\left|\frac{\partial{\boldsymbol{x}}_{sh}}{\partial t}\right|^2+ \left(\frac{\partial \theta}{\partial s}\right)^2\right]{\rm d}s+ \sum_{i={u,d}}\frac{1}{2\lambda}\iint_{v_i(t)}|\boldsymbol{\nabla}\phi_i|^2\,{\rm d}\kern0.7pt x\,{{\rm d}y}, \end{equation}

where $|\cdot |$ corresponds to the norm of the enclosed vector and our normalization implies that $E=\tilde {E}\tilde {L}/\tilde {B}$. Since our system starts from rest, the total energy of the system equals the initial potential energy of the sheet, $E=E_{sh}^{p}(0)$, and this energy is conserved throughout the system's evolution.

2.1. The small-amplitude approximation

The assumption that the amplitude of the sheet remains small, or equivalently that $\varDelta \ll 1$, implies that the geometric relations (2.6) reduce to $\partial y_{sh}/\partial s\simeq \theta$ and $\partial x_{sh}/\partial s\simeq 1-\frac {1}{2}(\partial y_{sh}/\partial s)^2$. The nonlinearity in the derivative of $x_{sh}(s,t)$ is retained in the leading order of the theory so as to satisfy the constraint on the excess length (2.9a). Indeed, in the small-amplitude approximation, this constraint is given by

(2.11)\begin{equation} \varDelta=\frac{1}{2}\int_0^1\left(\frac{\partial y_{sh}}{\partial x}\right)^2{\rm d}\kern0.7pt x. \end{equation}

Here, we replace the arclength coordinate of the sheet with the Eulerian coordinate of the fluid, $s\simeq x$, according to our level of approximation. Correspondingly, the balance of forces and moments on the sheet (2.8) reduces to

(2.12)\begin{equation} \frac{\partial^2 y_{sh}}{\partial t^2}+\frac{\partial^4 y_{sh}}{\partial x^4}+F_x(t) \frac{\partial^2 y_{sh}}{\partial x^2}+\left[p_{u}(x,0,t)-p_{d}(x,0,t)\right]=0, \end{equation}

where the lateral compression, $F_x(t)$, is a function that depends solely of time. In addition, the pressure difference that the fluid exerts on the sheet, i.e. the last term in (2.12), is calculated at the sheet–fluid interface, $y=0$.

Thus far, we have approximated only the elastic part of the model. To further simplify the hydrodynamic part, we need to estimate the order of its corresponding fields. Given that the initial energy of the sheet scales linearly with the excess length, $E_{sh}(0)\propto \varDelta$, and that the total energy of the system is conserved, the energy of the fluid is, at most, proportional to $E_{f}\sim \varDelta$. Therefore, if we approximate the energy of the fluid by $E_{f}\sim |{\boldsymbol {v}}|^2 \ell$, where $|{\boldsymbol {v}}|$ is the typical velocity in the chamber and $\ell$ is the decay length of the disturbances in the flow, we obtain $|{\boldsymbol {v}}|\sim \sqrt {\varDelta /\ell }$. Furthermore, if we assume that the order of approximation of a derivative over the potential function, with respect to either a spatial dimension or time, does not change, then we can approximate Bernoulli's equation (2.2b), and the kinematic boundary conditions by

(2.13a)\begin{gather} \lambda p_i+\frac{\partial\phi_i}{\partial t}=c(t), \end{gather}

(2.13b)\begin{gather}\frac{\partial y_{sh}}{\partial t}=\left(\frac{\partial \phi_i}{\partial y}\right)_{y=0}. \end{gather}

These approximations will further be verified a posteriori in § 4, where we analyse the nonlinear dynamics of the system. In particular, we compare the results of this approximation with the numerical solution of the nonlinear model (2.2)–(2.9). Note that since the continuity equations (2.2a) are already linear in the potential functions, they remain unchanged in our approximated model.

This completes the reduction of our model to the small-amplitude limit. In summary, (2.2a) and (2.11)–(2.13), supplemented by the linearized form of the boundary conditions (2.5), (2.4), (2.9b) and (2.9c), form a closure and describe the coupled dynamics of the sheet and the fluid in the small-amplitude approximation. A comment is necessary regarding this simplified formulation. In accordance with the derivation in Appendix B, it can be shown that the reduced model emanates from the minimization of the action $\mathcal {S}=\int _0^{\rm T}\mathcal {L}\,{\rm d}t$, where

(2.14)\begin{align} \mathcal{L} &= \int_0^1\left[\frac{1}{2}\left(\frac{\partial y_{sh}}{\partial t}\right)^2-\frac{1}{2} \left(\frac{\partial^2y_{sh}}{\partial x^2}\right)^2+F_x(t) \left(\frac{1}{2}\left(\frac{\partial y_{sh}}{\partial x}\right)^2-\varDelta\right)\right. \nonumber\\ &\quad \left.\vphantom{\left(\frac{\partial y_{sh}}{\partial t}\right)^2} +\frac{1}{\lambda}\left[\phi_{d}(x,0,t)-\phi_{u}(x,0,t)\right] \frac{\partial y_{sh}}{\partial t}\right]{{\rm d}\kern0.7pt x} \nonumber\\ &\quad -\frac{1}{2\lambda}\int_{0}^{{L_y}/{2}} \int_0^1|\boldsymbol{\nabla}\phi_{u}|^2\,{\rm d}\kern0.7pt x\,{\rm d}y- \frac{1}{2\lambda}\int_{-{L_y}/{2}}^0 \int_0^1|\boldsymbol{\nabla}\phi_{d}|^2\,{\rm d}\kern0.7pt x\,{\rm d}y, \end{align}

with respect to the elastic fields $y_{sh}(x,t)$ and $F_x(t)$ and the hydrodynamic fields $\phi _{i}(x,y,t)$. In the next section, we employ a modal expansion of these fields and combine it with the Lagrangian formulation to derive a simplified set of equations that are dependent solely on time.

2.1.1. Modal expansion

The continuity equations (2.2a), and their corresponding boundary conditions on the fluid–chamber interfaces (2.4), are satisfied when the potential functions are given by

(2.15)\begin{align} \phi_i(x,y,t) &= a_0(t)\left(y\pm \frac{L_y}{2}\right)+\sum_{m=1}^{\infty}\cos ({\rm \pi} m x)\left[a_m(t)\exp\left({{\rm \pi} m\left(y\pm \frac{L_y}{2}\right)}\right) \right. \nonumber\\ &\quad \left.+\,c_m(t)\exp\left({-{\rm \pi} m \left(y\pm\frac{L_y}{2}\right)}\right)\right], \end{align}

where $a_m(t)$ ($m=0,1,2,\ldots$) and $c_m(t)$ ($m=1,2,3,\ldots$) are unknown time-dependent coefficients and the $\pm$ signs correspond to the solutions of the potential functions in the lower and the upper parts of the chamber, respectively. Similarly, we expand the solution of the sheet's height function in the following normal modes:

(2.16)\begin{equation} y_{sh}(x,t)=\sum_{n=1}^{\infty}A_n(t)\sin({\rm \pi} n x), \end{equation}

where the functions $\sin ({\rm \pi} n x)$ automatically satisfy the boundary conditions on the sheet edges (2.9b) and (2.9c), and $A_n(t)$ are as-yet unknown coefficients.

With these expansions, the solution to our problem reduces to finding the unknown coefficients, $a_m(t)$, $c_m(t)$ and $A_n(t)$, and the compression force $F_x(t)$, from the solution of the force balance equation (2.12), Bernoulli's equation (2.13a), the kinematic boundary conditions (2.13b) and the geometric constraint (2.11). Equation (2.16) involves infinite summation over the modes of the height function, but, in practice, we will truncate this series at $n=N$. A closed system of equations is then obtained when the coefficients of $a_m(t)$ and $c_m(t)$ are truncated at $N-1$.

However, instead of directly using these equations, we take a different – yet equivalent – approach, by utilizing the Lagrangian formulation (2.14). To this end, we follow the analysis in Appendix C and substitute the potential functions (2.15), and the height function (2.16), into the Lagrangian (2.14). We then integrate over the spatial coordinates. Thereafter, we minimize the Lagrangian with respect to $a_m(t)$ and $c_m(t)$ and express these coefficients in terms of $A_n(t)$. Substituting $a_m(t)$ and $c_m(t)$ back into the Lagrangian gives

(2.17)\begin{equation} \mathcal{L}\left[A_1,\ldots,A_N,F_x\right]= T_{nk}\frac{{\rm d}A_k}{{\rm d}t}\frac{{\rm d}A_n}{{\rm d}t}-V_{nk}A_k A_n+F_x(t)C_{nk}A_k A_n-\varDelta F_x(t), \end{equation}

where Einstein's summation rule is implied for repeated indices, and we define the following symmetric matrices:

(2.18a)\begin{gather} T_{nk}=\frac{1}{4}\delta_{nk}+\frac{L_y}{2\lambda}W(n,0)W(k,0)+\sum_{m=1}^{N-1} \frac{2}{{\rm \pi} m \lambda}\tanh\left(\frac{{\rm \pi} m L_y}{2}\right)W(k,m)W(n,m), \end{gather}

(2.18b)\begin{gather}V_{nk}=\frac{{\rm \pi}^4}{4}n^2k^2\delta_{nk},\quad C_{kn}=\frac{{\rm \pi}^2}{4}nk\delta_{nk}, \end{gather}

where $\delta _{nk}$ is the Kronecker delta, and $W(n,m)= ({n}/{{\rm \pi} })({(1-(-1)^{n+m})}/{(n^2-m^2)})$ for $n\neq m$ and zero otherwise.

Two comments are in order regarding this Lagrangian. First, since the matrix $T_{nk}$ is coupled to the kinetic terms in the Lagrangian, it takes on the role of a mass matrix in this formulation. This mass matrix has contributions from both the inertia of the sheet, i.e. the first term in $T_{nk}$, and the hydrodynamics of the fluid, i.e. the terms proportional to $1/\lambda$. The latter hydrodynamic terms are frequently referred to as added mass or virtual mass, since they describe an additional mass that the sheet appears to acquire when it accelerates in the fluid (Munk 1924; Lighthill 1960; Coene 1992).

The second comment is related to the potential functions $\phi _i(x,y,t)$ that result from the minimization. Using (C2b), we substitute $c_m(t)=-a_m(t)$ in the potential functions (2.15) to obtain

(2.19)\begin{equation} \phi_i(x,y,t)=a_0(t)\left(y\pm \frac{L_y}{2}\right)+2\sum_{m=1}^{N-1}a_m(t)\cos({\rm \pi} m x)\sinh\left[{\rm \pi} m\left(y\pm \frac{L_y}{2}\right)\right]. \end{equation}

This solution implies that at $y=\pm L_y/2$ the velocity of the fluid is oriented only in the $y$ direction, i.e. $\partial \phi _i/\partial x(x,\pm L_y/2,t)=0$. It also implies that $\phi _i(x,\pm L_y/2,t)=0$, which, given (2.5) and (2.13a), yields a constant zero pressure along the inlet and the outlet walls of the chamber, $p_i(x,\pm L_y/2,t)=0$. We anticipate that these conditions will occur only when the disturbances that the sheet induces in the flow decay to zero. Therefore, the small-amplitude model holds strictly when $L_y\gg \ell$, where $\ell \simeq 1/({\rm \pi} m)$, for the smallest non-zero mode, is now explicitly identified as the decay length of the hydrodynamic disturbances.

Keeping in mind these limitations of the small-amplitude model, we go back to derive the equations for the coefficients $A_n(t)$. Given an initial volume difference, which, in turn, corresponds to an initial configuration of the sheet, i.e. a set of initial conditions for the coefficients $A_n(0)$, and keeping in mind that the system starts from rest, i.e. $({\rm d}A_n/{\rm d}t)(0)=0$, we can determine the dynamic evolution of the system from the minimization of (2.17) with respect to $A_n(t)$ and $F_x(t)$. This minimization yields $N+1$ algebraic differential equations that, in our matrix notation, read

(2.20a)\begin{gather} T_{nk}\frac{{\rm d}^2A_k}{{\rm d} t^2}+(V_{nk}-F_x(t)C_{nk}) A_k=0, \end{gather}

(2.20b)\begin{gather}C_{nk}A_kA_n=\varDelta. \end{gather}

Once $A_n(t)$ are determined from the solution of (2.20), the position of the sheet in time and in space is given by (2.16), and the hydrodynamic potentials and the pressure fields are determined from (2.15), (C2) and (2.13a). In the next section, we utilize this formulation to investigate the early time evolution. Then, in § 4 we use it to analyse the dynamics at later times.

3. The early time evolution

In this section, we investigate the system's stability close to an initial equilibrium state. The section is divided into two parts. In the first, we recall the static solutions of the system in the small-amplitude approximation. In the second, we employ a linear stability analysis around the first two buckling modes to extract the growth rates and the flow fields of the perturbation around these modes.

3.1. Recap of the quasi-static solution

Following the analysis in the study of Oshri (2021), the quasi-static evolution of the system is governed by two different branches of solutions, which we call ‘asymmetric’ and ‘symmetric’. Here, we recall the height functions in these branches.

On the one hand, when the initial volume difference is set as $0\leq v_{du}(0)\leq v_{du}^{cr}$, where $v_{du}^{cr}= ({2(3+{\rm \pi} ^2)}/{{\rm \pi} \sqrt {3(15+2{\rm \pi} ^2)}})\varDelta ^{1/2}$, the system is governed by the asymmetric branch. In this branch, the lateral compression is constant, $F_x(0)=4{\rm \pi} ^2$, and the height functions are given by

(3.1a)\begin{align} y_{sh}(x,0) &= \frac{p_{ud}(x,y,0)}{16{\rm \pi}^4}\left[2{\rm \pi}^2(1-x)x+1-\cos\left(2{\rm \pi} x\right)\right] \nonumber\\ &\quad +\frac{1}{\rm \pi}\sqrt{\varDelta-\frac{15+2{\rm \pi}^2}{768{\rm \pi}^6}p_{ud}(x,y,0)^2}\sin\left(2{\rm \pi} x\right), \end{align}

(3.1b)\begin{gather} p_{du}(x,y,0)=\frac{24{\rm \pi}^4}{3+{\rm \pi}^2}v_{du}(0), \end{gather}

where $p_{ud}(x,y,0)=p_{u}(x,y,0)-p_{d}(x,y,0)$ is the pressure difference between the upper and lower parts of the chamber. The potential energy of the sheet in this branch is given by

(3.2)\begin{equation} E_{as}=4{\rm \pi}^2\varDelta-\frac{6{\rm \pi}^4}{3+{\rm \pi}^2}v_{du}(0)^2. \end{equation}

Note that when $v_{du}(0)\rightarrow 0$ we know from (3.1b) that the pressure difference vanishes, $p_{ud}(x,y,0)\rightarrow 0$, and the elastic configuration converges to the second, asymmetric, mode of buckling, $y_{sh}(x,0)\rightarrow \sqrt {\varDelta /{\rm \pi} ^2}\sin (2{\rm \pi} x)$. The total energy of the system in this configuration is given by $E_{as}=4{\rm \pi} ^2\varDelta$. Note also that we considered solutions with an initial volume difference that is greater than zero. This is because the static solution has mirror symmetry around the $x$ axis. Solutions with $v_{du}(0)<0$ (and $p_{ud}(x,y,t)<0$) are obtained by a reflection of the height functions (3.1) around the horizontal axis.

On the other hand, when the volume difference is set as $v_{du}^{cr}\leq v_{du}(0)< \sqrt {2\varDelta /3}$, the system is governed by the symmetric branch. In this case, the inextensibility of the sheet implies an upper limit on the volume difference. In the case of a hinged sheet, this limit is given by $\sqrt {2\varDelta /3}$. The height functions in this branch are given by the parametric solution:

(3.3a)\begin{gather} y_{sh}(x,0)=\frac{p_{ud}(x,y,0)}{8u^2}(1-x)x+ \frac{p_{ud}(x,y,0)}{16u^4}\left[1-\frac{\cos\left[2u(x-1/2)\right]}{\cos u}\right], \end{gather}

(3.3b)\begin{gather}p_{ud}(x,y,0)={-}\frac{16\sqrt{6}u^{7/2}\cos u}{\sqrt{6u+4u(6+u^2)\cos^2 u-15\sin (2u)}}\varDelta^{1/2}, \end{gather}

(3.3c)\begin{gather}v_{du}(0)={-}\frac{2\sqrt{2}\left[u(3+u^2)-3\tan u\right] \cos u}{\sqrt{3}u^{3/2}\sqrt{6u+4u(6+u^2)\cos^2 u-15\sin(2u)}}\varDelta^{1/2}, \end{gather}

where $u=\sqrt {F_x(0)}/2$ is a function of the lateral compression. Given the initial volume difference, $v_{du}(0)$, and the excess length, $\varDelta$, we can determine the lateral compression, $F_x(0)$, from (3.3c), and then substitute this solution into (3.3a) and (3.3b) to obtain the height profile. When $p_{ud}(x,y,0)\rightarrow 0$, the height function converges to the first, symmetric, mode of buckling, which is given by $y_{sh}(x,0)\rightarrow \sqrt {4\varDelta /{\rm \pi} ^2}\sin ({\rm \pi} x)$ and $v_{du}(0)=8\varDelta ^{1/2}/{\rm \pi} ^2$. The total elastic energy of this shape is given by $E_{s}={\rm \pi} ^2\varDelta$. An example of the evolution of the sheet and the $p_{ud}(x,y,0)$–$v_{du}(0)$ relation in this static solution, where $0\leq v_{du}(0)< \sqrt {2\varDelta /3}$, is plotted in figure 2. In the following sections, we use these height functions (3.1a) and (3.3a) as the base solutions for our perturbative time-dependent expansion.

Figure 2. Evolution of the static solution in the small-amplitude approximation, where $\varDelta =0.01$. (a) The volume difference, $v_{du}(0)$, as a function of the pressure difference, $p_{ud}(x,y,0)$. In the asymmetric branch the $p_{ud}(x,y,0)$–$v_{du}(0)$ relation is given by (3.1b), while in the symmetric branch it is given by (3.3b) and (3.3c). The volume difference at the asymmetric-to-symmetric transition, $v_{du}(0)=v_{du}^{cr}$, is labelled by ③. The pressure difference in the chamber vanishes when the sheet accommodates either the second or the first mode of buckling, labels ① and ④. As the volume difference approaches its limiting value $v_{du}(0)\rightarrow (2\varDelta /3)^{1/2}$, the pressure difference diverges. (b) Evolution of the sheet's profile as the volume difference increases; see the corresponding labelled numbers in (a). Note, that despite the relatively large change in the pressure difference in the symmetric branch, the elastic configurations remain almost unchanged.

Before we proceed, we emphasize that while the sheet's configuration evolves continuously from the moment that we open the valve, the static pressure difference, given in (3.1b) and (3.3b), changes instantaneously at $t=0$. This is because we assumed an incompressible fluid, in which the speed of sound is infinite. Nonetheless, the asymmetric and the symmetric modes of buckling, obtained respectively from (3.1a) and (3.3a) in the limit $p_{ud}(x,y,0)\rightarrow 0$, are exceptions. These configurations remain in static equilibrium, which can nevertheless be unstable, when the valve is opened. For this reason, in the next section, we investigate the linear stability of the system around these two limiting initial states.

3.2. Linear stability

To derive the linear stability around the two limiting scenarios, i.e. the second and first buckling modes, we assume that the sheet's height function is given by the static solution, up to a small perturbation that grows exponentially with time. Correspondingly, we first perturb the amplitudes of the normal modes and the lateral compression around the base solutions, i.e. $A_n(t)=A_n(0)+\epsilon \bar {A}_n\,{\rm e}^{\sigma t}$ and $F_x(t)=F_x(0)+\epsilon \bar {F}_x\,{\rm e}^{\sigma t}$, where $\bar {A}_n$ and $\bar {F}_x$ are unknown constants, $\sigma$ is the growth rate and $\epsilon \ll 1$ is an arbitrarily small parameter. Then, we substitute these perturbed functions into the equations of motion (2.20), and expand them up to linear order in $\epsilon$. The leading order of this expansion, order $\epsilon ^0$, is given by

(3.4a)\begin{gather} V_{nk}A_k(0)-F_x(0)C_{nk}A_k(0)=0, \end{gather}

(3.4b)\begin{gather}C_{nk}A_k(0)A_n(0)=\varDelta, \end{gather}

and the subleading order, order $\epsilon$, is given by

(3.5a)\begin{gather} \left[\sigma^2 T_{nk}+V_{nk}-F_x(0)C_{nk}\right]\bar{A}_k-C_{nk}A_k(0)\bar{F}_x=0, \end{gather}

(3.5b)\begin{gather}C_{nk}A_k(0)\bar{A}_n=0. \end{gather}

The above equations in the subleading order always have the trivial solution $\bar {A}_n=0$ and $\bar {F}_x=0$, unless their corresponding determinant vanishes. This condition gives the growth rate, $\sigma$. Once $\sigma$ is determined, its corresponding eigenfunction is obtained from the solution of (3.5). The hydrodynamic fields related to this eigenfunction are determined from (2.15) and (C2).

3.2.1. Linear stability around the second mode of buckling

When the initial configuration of the sheet is given by the second mode of buckling, the base solution is derived from (3.1a) in the limit $p_{ud}(x,y,0)\rightarrow 0$. This solution reads $y_{sh}(x,0)=\sqrt {\varDelta /{\rm \pi} ^2}\sin (2{\rm \pi} x)$ and $F_x(0)=4{\rm \pi} ^2$. A projection of this configuration on the normal mode expansion (2.16) gives $A_n(0)=\sqrt {\varDelta /{\rm \pi} ^2}\delta _{2n}$. As expected, this initial state exactly satisfies the equilibrium equations at order $\epsilon ^0$ (3.4). At the next order, i.e. order $\epsilon$, we find that $\bar {F}_x=0$ and $\bar {A}_n=0$ for all the even perturbations, i.e. $n=2,4,6,\dots$. Consequently, (3.5) yields linear and homogeneous equations that involve only the odd perturbations. These equations always have the trivial solution $\bar {A}_n=0$, except when the corresponding determinant vanishes. A tractable solution to this condition, which also gives a good approximation to the highest growth rate, is obtained at the lowest order when $N=2$. This solution reads

(3.6)\begin{equation} \sigma=\frac{\sqrt{3}{\rm \pi}^2}{\sqrt{1+\dfrac{8 L_y}{{\rm \pi}^2\lambda}}}. \end{equation}

In figure 3, we plot this analytical approximation for the growth rate as a function of $\lambda$, and compare it with the numerical solution of (3.5) for the case where $N=8$. In addition, we compare this analytical solution with the growth rate obtained from the linearization of (2.2)–(2.9), i.e. where $\varDelta$ is assumed to be finite. See Appendix D for the details of this solution.

Figure 3. Log–log plot of the highest growth rate as a function of $\lambda$ where $L_y=2$. Symbols correspond to the linear stability analysis of (2.2)–(2.9), and solid and dashed lines correspond to the growth rates obtained from the small-amplitude approximation. When $\lambda \gg 1$, the growth rate converges to the constant $\sigma \simeq \sqrt {3}{\rm \pi} ^2$, whereas when $\lambda \ll 1$ the growth rate is given by $\sigma \simeq \sqrt {3{\rm \pi} ^6/8}(\lambda /L_y)^{1/2}$. Note that the differences between the solid ($N=2$) and dashed ($N=8$) lines are almost not visible in the figure. While the growth rate is independent of the excess length in the small-amplitude approximation, the more general solution (symbols) shows that the growth rate increases with $\varDelta$.

Equation (3.6) is one of the central results in this paper. Several comments are in order regarding this solution. Firstly, note that the highest growth rate is always real and positive, i.e. the second mode of buckling is always an unstable state of the system.

Secondly, while (3.6) depends on the parameter $\lambda /L_y=\tilde {\rho }_{sh} \tilde {h}/(\tilde {\rho }_{\ell }\tilde {L}_y)$, this result is not general but depends on the order of approximation. Had we solved (3.5) with $N\geq 3$, the two parameters $\lambda$ and $L_y$ would have appeared independently in the solution. Nonetheless, comparing the lowest order solution, i.e. the solution with $N=2$, with that of, say, $N=8$, we find that the growth rate remains almost unchanged (compare the solid and dashed lines in figure 3). For this reason, we conclude that (3.6) well describes the growth rate in the limit $\varDelta \ll 1$.

Thirdly, while the solution of the small-amplitude approximation is independent of the excess length, $\varDelta$, the more general solution for finite values of $\varDelta$ does depend on this parameter. See figure 3 for a comparison. In particular, for a fixed value of $\lambda$, the growth rate increases with an increase in $\varDelta$.

Fourthly, the analytical solution of the growth rate (3.6), exhibits two different regions as a function of $\lambda$. When $\lambda /L_y\gg 1$, the growth rate is a constant, $\sigma \simeq \sqrt {3}{\rm \pi} ^2$, that coincides with the growth rate of a sheet that is uncoupled from an external fluid. However, when $\lambda /L_y\ll 1$, the growth rate exhibits the scaling $\sigma \simeq \sqrt {3{\rm \pi} ^6/8}(\lambda /L_y)^{1/2}$. These asymptotic solutions define two limiting behaviours of the system. The former scenario represents a ‘solid-dominated’ region, in which the pressure difference exerted by the fluid on the sheet is negligible compared with the inertia of the sheet. The latter scenario, where $\sigma \propto (\lambda /L_y)^{1/2}$, represents the opposite limit of a ‘fluid-dominated’ region, in which the inertia of the sheet is negligible compared with the pressure difference exerted by the fluid on the sheet. Similarly, these two regions are also reflected in the added mass term, $8L_y/({\rm \pi} ^2\lambda )$, in the denominator of the growth rate. When $\lambda$ is large, the added mass approaches zero and the sheet's inertia is not affected by the fluid's motion. In contrast, when $\lambda$ is small, the effective mass of the sheet increases, resulting in slower dynamics. It is worth noting that since the dynamics of the system becomes very slow in the fluid-dominated region, we would expect some aspects of this solution to align with our earlier, quasi-static solution in the asymmetric branch, as shown in (3.1). This is because the quasi-static solution describes a slow spontaneous relaxation of the system, where the inertia of the sheet is negligible. This convergence to the quasi-static solution is further demonstrated in the following analysis.

Fifthly, note that if we fix $L_y$, the matrices in (3.5) and (2.18) become diagonal in the solid-dominated region. Therefore, up to small corrections of the order $1/\lambda$, $\bar {A}_1$ alone is excited at the instability. Indeed, in figure 4(a), in which we plot the eigenfunction for the case where $\lambda =100$ for both $N=2$ and $N=8$, we find that the two solutions are almost identical. In contrast, in the fluid-dominated region, $\lambda \ll 1$ and the matrices in (3.5) have non-zero off-diagonal terms. Therefore, all the odd modes become coupled and are excited at the instability. Nonetheless, our numerical investigation indicates that $\bar {A}_n/\bar {A}_1\ll 1$ for all $n\geq 3$. Therefore, while we would expect the leading-order solution, i.e. $N=2$, to approximate the eigenfunction well, it will not coincide with the higher-order solution. Indeed, in figure 4(b), we compare the eigenfunctions obtained from the lowest order and the high orders of approximations and find finite differences between them. We note that the eigenfunction in the fluid-dominated region emerges from the quasi-static configuration (3.1a), when we expand (3.1a) in powers of $p_{ud}(x,y,0)$, extract the linear order of this expansion and normalize it in accordance with our convention. The agreement between this quasi-static profile and the eigenfunction obtained from the linear stability analysis is shown in figure 4(b) (open squares).

Figure 4. Eigenfunctions of the sheet's height function in both the solid- and fluid-dominated regions. In both panels, $L_y=2$ and open blue circles correspond to the linear stability analysis at finite $\varDelta$, i.e. derived from (2.2)–(2.9). The eigenfunctions are normalized such that $[y_{sh}(1/2,t)-y_{sh}(1/2,0)]/{\rm e}^{\sigma t}=1$ (numerically we choose $\hat {y}_{sh}(1/2)=1$; see Appendix D). (a) In the solid-dominated region $(\lambda =100)$, only one mode is excited, i.e. the low ($N=2$) and the high ($N=8$) mode approximations coincide. (b) In the fluid-dominated region ($\lambda =0.1$), all odd modes are excited. Therefore, the lowest approximation, $N=2$, does not coincide with that obtained with higher modes, $N=8$. Nonetheless, since $\bar {A}_n/\bar {A}_1\ll 1$, the differences between the low and the high orders of the approximations are still small. Open squares represent the quasi-static approximation obtained from (3.1).

Sixthly, in figure 5(a), we plot the flow field obtained from the linear stability analysis at finite $\varDelta$. Note that the maximum velocity of the fluid is obtained at the sheet's centre, where the eigenfunction of the sheet's height is maximized (see figure 4). Unlike the growth rate, for which the small-amplitude approximation provided us with a good estimation already at $N=2$, the flow field converges at a slower pace. Convergence to the spatially dependent solution, presented in figure 5(b), is obtained only when higher modes, say $N\geq 3$, are included. This is because each coefficient $a_m^i(t)$ in the fluid's potential functions (2.15) depends on all the excited modes; see (C2). In addition, in figure 5(b), we plot the eigenfunctions of the pressure fields. Note that the sheet moves upwards towards the higher-pressure field. This is because the sheet's motion drives the flow, and the pressure drop in the fluid acts to slow down the onset of the elastic instability. We note that there are no qualitative differences in the flow fields and the pressure distributions between the solid- and fluid-dominated regions. For this reason, the plots in figure 5 refer only to the fluid-dominated region.

Figure 5. (a) The flow field and (b) the hydrodynamic pressures obtained from the linear stability analysis of (2.2)–(2.9) at the highest growth rate. In both panels, $\varDelta =0.01$, $\lambda =0.1$ and $L_y=2$. The eigenfunctions are normalized as indicated in figure 4. This gives $[\text {Low}, \text {High}]=[0.98,2.94]$ in the colour bar of the flow field and $[\text {Low}, \text {High}]=[-78,78]$ in the colour bar of the hydrodynamic pressures. The solid black line corresponds to the initial configuration of the sheet, i.e. the asymmetric second mode of buckling. In (a), arrows represent the streamlines and colours represent the relative magnitudes of the velocity.

Lastly, in addition to the growth rate, the flow field and the hydrodynamic pressures, another experimentally measurable quantity is the system's ‘compressibility’, i.e. the change in the volume difference relative to the change in the pressure difference. Since the pressure difference varies in space, we define the compressibility as $\beta \equiv {\rm d}v_{du}/{\rm d}\bar {p}_{ud}$, where $\bar {p}_{ud}(t)$ is the average pressure drop on the sheet. In the small-amplitude approximation the average pressure drop on the sheet is given by $\bar {p}_{ud}(t)=\int _0^1 [p_{u}(x,0,t)-p_{d}(x,0,t)]\,{{\rm d}\kern0.7pt x}$ and the volume difference in the chamber is given by $v_{du}(t)=2\int _0^1 y_{sh}(x,t)\,{{\rm d}\kern0.7pt x}$ (Oshri 2021). Keeping in mind that the base solution is asymmetric, and therefore does not contribute to the above integral, we find that the leading order ($N=2$) is $v_{du}(t)/\bar {A}_1=(4/{\rm \pi} )\,{\rm e}^{\sigma t}$. In addition, the average pressure difference at this order is given by $\bar {p}_{ud}(t)/\bar {A}_1=2L_y\sigma ^2 \,{\rm e}^{\sigma t}/({\rm \pi} \lambda )$. Consequently, at the onset of the instability, the compressibility is a constant, which is given by

(3.7)\begin{equation} \beta=\frac{2\lambda}{\sigma^2 L_y}. \end{equation}

In figure 6, we plot this result as a function of $\lambda$ and compare the lowest-order approximation with the solution at finite values of $\varDelta$. On the one hand, since the growth rate is constant in the solid-dominated region, we find that the compressibility scales linearly with $\lambda$, i.e. $\beta \rightarrow 2\lambda /(3{\rm \pi} ^4 L_y)$. On the other hand, since in the fluid-dominated region the growth rate scales as $(\lambda /L_y)^{1/2}$, the compressibility converges to the constant $\beta \rightarrow 16/(3{\rm \pi} ^6)\simeq 5.54\times 10^{-3}$. As expected, this result for the fluid-dominated region is very close to the compressibility obtained in the quasi-static solution (3.1b), which gives $\beta =(3+{\rm \pi} ^2)/(24{\rm \pi} ^4)\simeq 5.50\times 10^{-3}$.

Figure 6. Log–log plot of the compressibility as a function of $\lambda$ close to the onset of the instability. Symbols correspond to the compressibility at finite values of $\varDelta$ and the solid black line corresponds to the analytical solution obtained from the small-slope approximation. While in the solid-dominated region $\beta \propto \lambda$, in the fluid-dominated region $\beta$ converges to a constant.

3.2.2. Linear stability around the first mode of buckling

This section analyses the system's stability around the first mode of buckling. In contrast to the second mode, which is always unstable, the first mode represents the minimum of the elastic potential energy. Therefore, the first mode is expected to remain stable and to yield oscillatory motion under a small dynamical perturbation. Since a purely imaginary growth rate represents this oscillatory motion, we set $\sigma \equiv i\omega$, where $i=\sqrt {-1}$, and look for the lowest frequency of oscillation at the onset of the instability.

When the initial state of the sheet is given by the first mode of buckling, we have from (3.3) that the base solution is given by $y(x,0)=(2\sqrt {\varDelta }/{\rm \pi} )\sin ({\rm \pi} x)$ and $F_x(0)={\rm \pi} ^2$. A projection of this configuration on the normal mode expansion (2.16) gives $A_n(0)=(2\sqrt {\varDelta }/{\rm \pi} )\delta _{1n}$. This initial state, as expected, satisfies the leading order of our perturbative expansion (3.4). Substituting this leading order in (3.5) and solving for the unknown constants, we find that $\bar {F}_x=0$ and $\bar {A}_n=0$ for all the odd modes ($n=1,3,5,\dots$). Consequently, (3.5) yields linear and homogeneous equations that involve only the even modes of the height function. As in the previous case that we considered, a tractable solution to the vanishing determinant condition is obtained when we cut the normal mode expansion at the smallest value, $N=2$. (We note that when $L_y\gg 1$, the solution obtained from the two-mode approximation ($N=2$) is pre-empted by a different branch of solutions. The new branch emanates from a higher-order correction in the modal expansion, and its details are beyond the scope of the present study.) This gives

(3.8)\begin{equation} \omega= \frac{2\sqrt{3}{\rm \pi}^2}{\sqrt{1+\dfrac{128}{9{\rm \pi}^3}\dfrac{\tanh({\rm \pi} L_y/2)}{\lambda}}}. \end{equation}

In figure 7(a), we plot this solution and compare it with the solution of (3.5) for the case of $N=8$. Note that our solution, (3.8) with $N=2$, approximates well the small-amplitude limit, $\varDelta \ll 1$. Higher-order corrections, e.g. with $N=8$, almost do not alter this solution; compare the solid and dashed lines in figure 7(a). In addition, we compare (3.8) with the eigenvalues obtained from the linearization of (2.2)–(2.9), where the small-amplitude approximation is relaxed. We observe that as $\varDelta$ increases, the oscillation frequency decreases with increasing excess length. However, the overall trend of the dependence of $\omega$ on $\lambda$ remains consistent with the small-amplitude solution. We also note that the relatively large decrease in $\omega$ when $\lambda \gg 1$ is a result of our chosen hinged boundary conditions. In contrast, systems with clamped boundary conditions display a much milder dependence on $\varDelta$, as reported in previous studies (Neukirch et al. 2012; Pandey et al. 2014).

Figure 7. The lowest oscillation frequency and the sheet's eigenfunction obtained from the linear stability analysis around the first buckling mode. In both panels $L_y=2$. (a) Log–log plot of the oscillation frequency as a function of the parameter $\lambda$. Solid and dashed lines correspond to the solution of (3.8) when $N=2$ and $N=8$, respectively. All symbols correspond to the linear stability analysis of (2.2)–(2.9) at finite $\varDelta$. (b) The sheet's eigenfunctions in the solid- and fluid-dominated regions. All eigenfunctions are normalized such that their height at $x=1/4$ is equal to one (numerically we choose $\hat {y}_{sh}(1/4)=1$; see Appendix D). Although only one mode is excited when $\lambda \gg 1$ and infinite modes are excited when $\lambda \ll 1$, the eigenfunctions of the two regions are almost identical. Inset: an example of the sheet's oscillations around the base solution. Dashed lines correspond to an illustration of the dynamic oscillations, and the solid line to the base solution.

The frequency $\omega$ exemplifies the two limiting scenarios that we encountered for the growth rate in (3.6). On the one hand, in the solid-dominated region, where $\lambda \gg 1$, the frequency converges to the constant $\omega \simeq 2\sqrt {3}{\rm \pi} ^2$ that coincides with the frequency of oscillation of a sheet that is uncoupled from a fluid flow. On the other hand, in the fluid-dominated region, where $\lambda \ll 1$, we find the scaling $\omega \simeq \sqrt {27{\rm \pi} ^7/32}[\lambda /\tanh (L_y{\rm \pi} /2)]^{1/2}$, i.e. $\omega \propto \lambda ^{1/2}$. Alternatively, since the added mass is given by $({128}/{9{\rm \pi} ^3})({\tanh ({\rm \pi} L_y/2)}/{\lambda })$, when $\lambda \gg 1$ the sheet's inertia is almost unaffected by the fluid motion, but when $\lambda \ll 1$, the added mass increases and slows down the dynamics.

These two scenarios are also manifested in the eigenfunction of the sheet's height function. In the solid-dominated region, (3.5) becomes diagonal, up to corrections of the order of $1/\lambda$, and essentially only the second mode, i.e. $\bar {A}_2$, is excited at the instability, while in the fluid-dominated region (3.5) has non-zero off-diagonal terms, and all the even modes are excited. Nonetheless, our numerical investigation of the solution of (3.5) in the fluid-dominated region indicates that the ratio $\bar {A}_n/\bar {A}_2$ remains small for all $n$. Therefore, in both regions, deviations of the eigenfunction from the second mode are almost not visible in figure 7(b).

The oscillations of the sheet induce rotational flow in the chamber, whose magnitude decreases monotonically as $y$ approaches the far distant walls; see figure 8(a). Indeed, the solution for $\omega$ (3.8) becomes independent of $L_y$ as the vertical dimension of the chamber increases. Furthermore, the fluid's maximum velocity is obtained close to the centre point, $x=1/2$, where the sheet's velocity equals zero. This is because the sheet moves up and down around this centre line and drives a net flux across it. The velocity of this flux is maximized at the centre of the sheet. It is important to note that there is a slight asymmetry in the flow patterns observed in the upper and lower regions of the chamber, which may be attributed to the initial non-zero volume difference in the base solution.

Figure 8. (a) The flow fields and (b) the hydrodynamic pressure fields obtained from the linear stability analysis of (2.2)–(2.9) when $\varDelta =0.01$, $\lambda =0.1$ and $L_y=2$. The normalization of the eigenfunctions is as in figure 7(b). This normalization implies $[\text {High}, \text {Low}]=[5.3,26.5]$ in the flow fields and $[\text {High}, \text {Low}]=[-560,560]$ in the pressure fields. In (a), arrows represent the streamlines and colours represent the relative magnitudes of the velocity.

The pressure fields induced by the elastic oscillations are plotted in figure 8(b) and show an asymmetric profile in correlation with the eigenfunction of the sheet (figure 7b). Since only even modes are excited at the instability, the average pressure difference on the sheet, $\bar {p}_{ud}(t)$, vanishes. Therefore, in this case, there is no analogue to the compressibility calculated in § 3.2.1.

4. The evolution at moderate times

In this section, we relax the assumption that $t\ll 1$ and extend the analysis up to moderate times. The limits of this analysis are discussed at the end of this section. In particular, we require that the initial configuration of the sheet be close in shape to the second mode of buckling, i.e. $v_{du}(0)\ll v_{du}^{cr}$, where $v_{du}^{cr}=({2(3+{\rm \pi} ^2)}/{{\rm \pi} \sqrt {3(15+2{\rm \pi} ^2)}})\varDelta ^{1/2}$ (see § 3.1), and we investigate the following questions: (i) What is the maximum amount of energy that is transferred from the sheet to the fluid? (ii) How long does it take the system to convert this maximum elastic energy into a fluid flow? (iii) What is the time-dependent behaviour of the $\bar {p}_{ud}(t)$–$v_{du}(t)$ relation. To address these questions, we assume that the amplitude of the sheet remains small during the dynamic evolution of the system and utilize the approximated formulation derived in § 2.1. This formulation yields the simplified set of nonlinear equations (2.20) that describe the coupling between the elastic and the hydrodynamic equations. As may be seen, this set of equations has a conserved first integral, $E=T_{nk}({{\rm d}A_k}/{{\rm d}t})({{\rm d}A_n}/{{\rm d}t})+V_{nk}A_k A_n$, that corresponds to the total energy of the system (2.10). This conserved energy constitutes the starting point for the discussion that follows.

When the initial configuration of the sheet is close in shape to the second mode of buckling, we expect the system's dynamics at moderate times to depend strongly on the first two modes, $A_1(t)$ and $A_2(t)$. This is because the initial shape of the sheet and its corresponding eigenfunction at the highest growth rate are described approximately by these two modes (see figure 4). Therefore, we reduce the expression for the total energy to the case in which $N=2$ and obtain the following equation:

(4.1)\begin{equation} E=\frac{1}{4}\left(1+\frac{8L_y}{{\rm \pi}^2\lambda}\right) \left(\frac{{\rm d} A_1}{{\rm d} t}\right)^2+\frac{1}{4} \left(1+\frac{128\tanh({\rm \pi} L_y/2)}{9{\rm \pi}^3\lambda}\right) \left(\frac{{\rm d} A_2}{{\rm d} t}\right)^2+\frac{{\rm \pi}^4}{4}A_1^2+4{\rm \pi}^4A_2^2. \end{equation}

We keep in mind that $A_1(t)$ and $A_2(t)$ are related through the constraint of the excess length (2.20b).

To obtain some insight regarding the validity of this two-mode approximation, we solve (2.20) numerically with $N=2$ and compare the results with the solution of the nonlinear model, i.e. the numerical solution of (2.2)–(2.9). In our investigation, we set $\varDelta =0.01$ and $L_y=2$, and consider two different values for the parameter $\lambda$. The initial configuration is given by (3.1), where $v_{du}(0)=0.01v_{du}^{cr}$. (The initial conditions are given by $A_1(0)=2\int _0^1 y_{sh}(x,0)\sin ({\rm \pi} x)\,{{\rm d}\kern0.7pt x}={32v_{du}(0)}/{(3{\rm \pi} +{\rm \pi} ^3)}$ and $A_2(0)=\sqrt {{\varDelta }/{{\rm \pi} ^2}-{256v_{du}(0)^2}/{(3{\rm \pi} +{\rm \pi} ^3)^2}}$, such that (2.20b) is satisfied. In addition, we keep in mind that the system starts from rest, $({{\rm d}A_1}/{{\rm d}t})(0)=({{\rm d}A_2}/{{\rm d}t})(0)=0$.) The results of these numerical solutions are presented in figures 9(a) and 9(b), where we follow the time-dependent behaviour of the midpoint on the sheet, $y_{sh}(1/2,t)$. The configurations of the sheet along the trajectory depicted in figure 9(b) are presented in figure 9(c). In both cases, we find that the approximated solution breaks down slightly after $y_{sh}(1/2,t)$ reaches its first maximum. In the solid-dominated region ($\lambda =100$), the two-mode approximation holds over one period of motion, while in the fluid-dominated region ($\lambda =0.1$), the approximation breaks down a little earlier. This difference is probably due to the different number of excited modes at the instability in each region of the system (see discussion in § 3.2.1). We conjecture that higher modes become active beyond moderate times in the fluid-dominated region and perturb the system's trajectory from the two-mode approximation. Indeed, by solving (2.20) with $N=3$ instead of $N=2$, we observe a convergence towards the numerical data for longer times (see the dashed grey line in figure 9b).

Figure 9. The sheet's midpoint as a function of time in (a) the solid-dominated ($\lambda =100$) and (b) the fluid-dominated ($\lambda =0.1$) regions. In both panels, $\varDelta =0.01$, $L_y=2$ and the growth rate, $\sigma$, is approximated by (3.6). The solid black line denotes the two-mode approximation, i.e. (2.20) with $N=2$, and the open blue circles denote the solution of the nonlinear model (2.2)–(2.9). In (b), the dashed grey line denotes the solution of (2.20) with $N=3$. The initial configuration of the sheet is given by (3.1) with $v_{du}(0)=0.01v_{du}^{cr}$. In both cases, the approximated solution breaks down after $y_{sh}(1/2,t)$ reaches its first maximum. In the solid-dominated region, the two-mode approximation holds for longer times compared with the fluid-dominated region. (c) The configurations of the sheet along the trajectory depicted in (b); see the corresponding markers in (b). Between ① and ③ the sheet releases potential energy as it transforms from the second mode of buckling to the first mode of buckling. After ③, the height of the sheet's midpoint decreases and the sheet gains back potential energy, as seen in ④.

Nonetheless, in both cases, i.e. the solid- and fluid-dominated regions, the agreement between the numerical solution of (2.2)–(2.9) and the analytical approximation, (2.20) with $N=2$, holds up to the first maximum. Similar results are also obtained when we perturb the system's parameters, $\varDelta$ and $L_y$, and the initial configuration. Therefore, in the following analysis, we utilize this approximation to examine the system's behaviour up to the point where the midpoint of the sheet reaches its first maximum. We refer to this stage of the system as moderate times.

4.1. The elastohydrodynamic energetic interplay

Since the initial configuration of the sheet is close in shape to the second mode of buckling and the total energy of the system is conserved, the total energy at any $t>0$ is given by (3.2). In the limit $v_{du}(0)\ll v_{du}^{cr}$, this energy is approximated as $E_{as}\simeq 4{\rm \pi} ^2 \varDelta$, i.e. the energy of the second mode of buckling. In addition, since the first mode of buckling is the global minimizer of the elastic sheet's potential energy, the potential energy of the sheet cannot fall below $E_{s}={\rm \pi} ^2\varDelta$. Therefore, at most, our system can convert $\delta E=E_{as}-E_{s}\simeq 3{\rm \pi} ^2 \varDelta$ of the initial potential energy either into kinetic energy of the sheet or into hydrodynamic energy of the fluid. We remind the reader that the kinetic and potential energies of the sheet, $E_{sh}^{k}(t)$ and $E_{sh}^{p}(t)$, and the energy of the fluid, $E_{f}(t)$, are given by the first, second and third terms, respectively, on the right-hand side of (2.10). In the small-amplitude approximation, these energies reduce to $E_{sh}^{k}(t)=(\lim _{\lambda \rightarrow \infty } T_{kn})({{\rm d}A_k}/{{\rm d}t})({{\rm d}A_n}/{{\rm d}t})$, $E_{sh}^{p}(t)=V_{kn}A_k A_n$ and $E_{f}(t)=T_{kn}({{\rm d}A_k}/{{\rm d}t})({{\rm d}A_n}/{{\rm d}t})-E_{sh}^{k}$.

The typical evolution of the three components of the energy, i.e. the kinetic energy of the sheet, the potential energy of the sheet and the energy of the fluid, is plotted in figures 10(a) and 10(b). These plots are obtained from the solution of (2.20) with $N=2$ in the solid- and the fluid-dominated regions, $\lambda =10$ and $\lambda =0.1$ respectively, where the initial configuration is given by (3.1) with $v_{du}(0)=0.01v_{du}^{cr}$. In both cases, we find that, after some initial delay, the potential energy of the sheet drops from $E\simeq 4{\rm \pi} ^2 \varDelta$, i.e. the energy of the second mode of buckling (3.2), to $E\simeq {\rm \pi}^2\varDelta$, i.e. the energy of the first mode of buckling. The energy released in this process, $\delta E\simeq 3{\rm \pi} ^2\varDelta$, is converted to the kinetic energy of the sheet and the energy of the fluid, while the total energy remains fixed. In the solid-dominated region (figure 10a) the kinetic energy of the sheet becomes much larger than the energy of the fluid, while as $\lambda$ decreases, the opposite picture emerges, i.e. the energy of the fluid becomes much larger than the kinetic energy of the sheet (figure 10b). In both scenarios, shortly after the initial peak, a portion of the kinetic energy is converted back into potential energy of the sheet. As a result, the sheet tends to return to its configuration of the second buckling mode, thereby decreasing the height of the sheet's midpoint, as shown in figure 9(c).

Figure 10. The energetic interplay between the three components of the total energy in (a) the solid-dominated ($\lambda =10$) and (b) the fluid-dominated ($\lambda =0.1$) regions of the system. In these two panels, we solve (2.20) for the case where $N=2$, $\varDelta =0.01$, $L_y=2$, and $\sigma$ is given by (3.6). The initial configuration of the sheet is given by (3.1), where $v_{du}(0)=0.01v_{du}^{cr}$. After some initial delay, the potential energy of the sheet drops from $E_{sh}^{p}(t\ll 1)\simeq 4{\rm \pi} ^2 \varDelta$ to $E_{sh}^{p}({t}_p)\simeq {\rm \pi}^2\varDelta$. Open blue circles represent the potential energy obtained from the solution of the nonlinear equations (2.2)–(2.9). The energy released from the sheet is divided between the kinetic energy of the sheet, $E_{sh}^{k}(t)$, and the hydrodynamic energy of the fluid, $E_{f}(t)$, such that the total energy, $E$, remains constant. In the solid-dominated region, $E_{sh}^{k}(t_{p})\gg E_{f}(t_{p})$, while in the fluid-dominated region we find the opposite relation, $E_{sh}^{k}(t_{p})\ll E_{f}(t_{p})$. (c) The kinetic energy of the sheet and the fluid at $t=t_{p}$ as a function of $\lambda$, i.e. $E_i(t_{p})=E_{sh}^{k}(t_{p}), E_{f}(t_{p})$. A logarithmic scale is used on the $x$ axis. The dotted and the dashed-dotted lines correspond to our analytical solution from the two-mode approximation, and the colour symbols correspond to the numerical solution of (2.2)–(2.9), where the initial conditions are similar to those used in (a,b).

Our two-mode approximation allows us to quantify these findings and to further estimate the maximum values of $E_{sh}^{k}( t_{p})$ and $E_{f}(t_{p})$ as a function of $\lambda$, where $t_{p}$ denotes the time at which the potential energy of the sheet reaches its minimum value. Indeed, at $E_{sh}^{p}(t_{p})={\rm \pi} ^2\varDelta$, the sheet is close in shape to the first mode of buckling, i.e. $A_1(t_{p})\simeq 2\sqrt {\varDelta }/{\rm \pi}$. The constraint on the excess length (2.20b) then implies that $A_2(t_{p})\simeq 0$ and that $({\rm d}A_1/{\rm d}t)(t_{p})\simeq 0$. In addition, the derivative of the second mode at that moment, $({\rm d}A_2/{\rm d}t)(t_{p})$, is obtained from (4.1) when we substitute $E=4{\rm \pi} ^2\varDelta$ for the total energy. Taken together, these approximations yield the following maximum energies at $t=t_{p}$:

(4.2a)\begin{gather} \frac{E_{sh}^{k}(t_{p})}{3{\rm \pi}^2\varDelta}= \frac{1}{1+\dfrac{128 \tanh({\rm \pi} L_y/2)}{9{\rm \pi}^3\lambda}}, \end{gather}

(4.2b)\begin{gather}\frac{E_{f}(t_{p})}{3{\rm \pi}^2\varDelta}=\frac{1}{1+\dfrac{9{\rm \pi}^3\lambda}{128 \tanh({\rm \pi} L_y/2)}}. \end{gather}

In figure 10(c), we compare these maximum energies with the numerical solution of (2.2)–(2.9). Overall, we find a good fit between the analytical approximation and the numerical solution over the entire range of $\lambda$.

Using (4.2), we find that in the solid-dominated region, $\lambda \gg 1$, most of the initial energy is converted into the kinetic energy of the sheet, i.e. $E_{sh}^{k}(t_{p})\simeq 3{\rm \pi} ^2 \varDelta -\delta _1$ and $E_{f}(t_{p})\simeq \delta _1$, where $\delta _1=128\varDelta \tanh ({\rm \pi} L_y/2)/(3{\rm \pi} \lambda )$, while in the fluid-dominated region, $\lambda \ll 1$, most of the energy is converted into the energy of the fluid, i.e. $E_{sh}^{k}(t_{p})\simeq \delta _2$ and $E_{f}(t_{p})\simeq 3{\rm \pi} ^2\varDelta -\delta _2$, where $\delta _2=27{\rm \pi} ^5\varDelta \lambda /[128\tanh ({\rm \pi} L_y/2)]$. Furthermore, if we estimate the average instantaneous velocity of the fluid, $\bar {v}(t_{p})$, by using $E_{f}(t_{p})\propto \bar {v}(t_{p})^2 \ell /\lambda$, we find that in the solid-dominated region the average velocity is proportional to $\bar {v}(t_{p})\propto (\varDelta /\ell )^{1/2}$, while in the fluid-dominated region it is proportional to $\bar {v}(t_{p})\propto (\lambda \varDelta /\ell )^{1/2}$; namely, while $E_{f}(t_{p},\lambda \ll 1)\gg E_{f}(t_{p},\lambda \gg 1)$, their corresponding velocities present the opposite relationship, i.e. $\bar {v}(t_{p},\lambda \ll 1)\ll \bar {v}(t_{p},\lambda \gg 1)$. This is because the momentum of the fluid $p_{f}\propto \bar {v}(t_{p})/\lambda$ – but not its velocity – increases in the fluid-dominated region.

4.2. The peak time

In the previous section, we demonstrated through energetic considerations that the sheet tends to release most of its stored potential energy. In this section, we investigate the time it takes for the system to release this energy. Given an initial configuration of the sheet that is close in shape to the second mode of buckling, i.e. (3.1) with $v_{du}(0)\ll v_{du}^{cr}$, we aim to find the time $t=t_{p}$ at which the potential energy of the sheet first drops to the minimum value, $E_{sh}^{p}(t_{p})\simeq {\rm \pi}^2 \varDelta$. To do so, we first use (2.20b) and (3.6) to eliminate $A_2(t)$ and $\lambda$ in favour of $A_1(t)$ and $\sigma$, respectively, in (4.1). Then, we substitute the energy of the initial configuration (3.2) into (4.1) and integrate it between $t\in [0,t_{p}]$. This gives

(4.3)\begin{equation} \sigma t_{p}=\int_{A_1(0)}^{A_1(t_{p})} \frac{\sqrt{432{\rm \pi}^3\varDelta+\left[9{\rm \pi}\left({-}12{\rm \pi}^4+\sigma^2\right)+ \dfrac{16\left(3{\rm \pi}^4-\sigma^2\right)}{L_y\coth({\rm \pi} L_y/2)}\right]A_1^2}}{6\sqrt{3}{\rm \pi}^{3/2}\sqrt{\left(A_1^2-\dfrac{8v_{du}(0)^2}{3+{\rm \pi}^2}\right) \left(4\varDelta-{\rm \pi}^2A_1^2\right)}}{\rm d} A_1, \end{equation}

where $A_1(t_{p})\simeq 2\varDelta ^{1/2}/{\rm \pi}$ is approximately the amplitude of the sheet at time $t_{p}$ and $A_1(0)=2\int _0^1 y_{sh}(x,0)\sin ({\rm \pi} x)\,{{\rm d}\kern0.7pt x}=32v_{du}(0)/(3{\rm \pi} +{\rm \pi} ^3)$ is the projection of the initial configuration (3.1) on the first mode of the sheet.

In figure 11(a), we fix the excess length and the vertical dimension of the chamber at $\varDelta =0.01$ and $L_y=2$, respectively, and plot $\sigma t_{p}$ for $\lambda \in [10^{-3},10^3]$ ($\sigma \in [0.4,\sqrt {3}{\rm \pi} ^2]$). We find that at a given initial volume difference, $v_{du}(0)$, the peak time, $\sigma t_{p}$, changes by less than 5 % over more than six orders of magnitude in $\lambda$. While $\sigma t_{p}$ is almost independent of $\lambda$, it does depend strongly on the initial configuration of the sheet, $v_{du}(0)$, and the excess length, $\varDelta$. An analytical approximation of this dependence can be extracted from (4.3), if we assume that the system is in the solid-dominated region, where $\sigma \simeq \sqrt {3}{\rm \pi} ^2$. Under this assumption, we can integrate the right-hand side of this equation and take the limit $v_{du}(0)\ll 1$ of the resulting expression. (We use Mathematica (Wolfram Research 2018) for the symbolic integration.) This gives

(4.4)\begin{equation} \sigma t_{p}\simeq \ln \left(c\frac{\varDelta^{1/2}}{v_{du}(0)}\right), \end{equation}

where $c\simeq 2.9$. While the scaling in (4.4) agrees well with our numerical solution of the nonlinear model, there is a small deviation in the numerical prefactor. The best fit to the nonlinear model gives $c\simeq 2.0$ (see figure 11b). We note that the independence of $L_y$ on the right-hand side of (4.4) is a result of our assumption that $\lambda \gg 1$. Had we derived this scaling using the fluid-dominated region, we would have found that the integral in (4.3) does depend on the vertical dimension of the chamber. However, this dependence diminishes to zero when $L_y\gtrsim 1$, as is required by our fourth assumption in § 2.

Figure 11. The peak time, $\sigma t_{p}$, as a function of the system's parameters. (a) The peak time as a function of $\lambda$, where $\varDelta =0.01$ and $L_y=2$, for two different values of the initial volume difference. Dashed lines correspond to the numerical integration of (4.3) and symbols with corresponding colours represent the numerical solution of (2.2)–(2.9). For a given $v_{du}(0)$, the time $\sigma t_{p}$ changes by less than $5\,\%$ over six orders of magnitude of the parameter $\lambda$. (b) Comparison between the analytical scaling (4.4) and the solution obtained from the nonlinear model. A logarithmic scale is used on the $x$ axis. Symbols correspond to the numerical solution of the nonlinear model. While the scaling of the analytical solution agrees well with the numerical data, the prefactor $c\simeq 2.9$ slightly overestimates the numerical prediction.

We also note that the logarithmic divergence of $\sigma t_{p}$ in the limit $v_{du}(0)\rightarrow 0$ is similar to the divergence of the period of a pendulum near the separatrix (Butikov 1999). Within this analogy between the two problems, the small initial deviation of the pendulum from the unstable vertical position is analogous to the small initial volume difference $v_{du}(0)$. The separatrix of the pendulum is analogous to the trajectory $v_{du}(0)=0$ in the phase space of our system.

4.3. The $\bar {p}_{ud}$–$v_{du}$ relation

In this section, we investigate the behaviour of the $\bar {p}_{ud}$–$v_{du}$ relation at moderate times. To this end, we first use the two-mode approximation to calculate the pressure difference in the chamber. From Bernoulli's equation (2.13a), and the normal mode expansion (2.15) and (C2), we have that the average pressure difference on the sheet is given by $\bar {p}_{ud}(t)=({2L_y}/{{\rm \pi} \lambda })({{\rm d}^2A_1}/{{\rm d}t^2})$. In addition, using (2.16), we find the volume difference as a function of time, $v_{du}(t)=2\int _0^1 y_{sh}(x,t)\,{{\rm d}\kern0.7pt x}=4A_1(t)/{\rm \pi}$. Thereafter, we solve (2.20) numerically with $N=2$ and plot the parametric solution ($\bar {p}_{ud}(t)$, $v_{du}(t)$) in the range $t\in [0,t_{p}]$.

The results of these solutions are plotted in figures 12(a) and 12(b) for the solid-dominated ($\lambda =100$) and fluid-dominated ($\lambda =0.01$) regions, respectively. Qualitatively, the two regions exhibit similar behaviour. The pressure difference in the chamber increases from almost zero up to a maximum positive value, from which it rapidly decreases and becomes negative. The backward pressure is maximized at the peak time $t_{p}$; see the time-dependent behaviour of $\bar {p}_{ud}(t)$ in the insets of these figures. Quantitatively, however, the two profiles are considerably different, because the maximum backward pressure is much larger in the fluid-dominated region (figure 12b) than in the solid-dominated region (figure 12a).

Figure 12. The $\bar {p}_{ud}$–$v_{du}$ relation and the maximum backward pressure. In all three panels, open green triangles correspond to the numerical solution of (2.2)–(2.9) and solid lines correspond to the solution of the two-mode approximation. The volume difference as a function of the average pressure difference on the sheet is plotted in (a) the solid-dominated ($\lambda =100$) and (b) the fluid-dominated ($\lambda =0.01$) regions. In both panels, the solid black line corresponds to the solution of (2.20) with $N=2$, $\varDelta =0.01$, $L_y=2$ and $v_{du}(0)=0.01v_{du}^{cr}$. The blue points correspond to times $t=0$ and $t=t_{p}$. The dashed grey lines correspond to the static solution (3.1b) and (3.3b). The insets show $\bar {p}_{ud}(t)$ as a function of time, where $\sigma$ is given by (3.6). In the fluid-dominated region, the maximum backward pressure ($\bar {p}_{ud}(t_{p})\simeq -200$) is much larger than that in the solid-dominated region ($\bar {p}_{ud}(t_{p})\simeq -1$). In addition, the evolution of the $\bar {p}_{ud}$–$v_{du}$ relation in the fluid-dominated region follows the static solution, except for some deviations close to the asymmetric-to-symmetric transition. (c) The absolute value of the maximum average pressure difference on the sheet.

The transition from positive to negative pressure differences can be explained as follows. Initially, the system exhibits a ‘negative feedback’ between the sheet and the fluid, meaning that the sheet's motion drives the fluid's dynamics, which, in turn applies a positive pressure difference and resists the sheet's motion. As the system evolves, the fluid continuously gains kinetic energy and reduces its resistance to the sheet's motion. Then, at some instant, the pressure difference vanishes, $\bar {p}_{ud}=0$, and the resistance of the fluid is almost eliminated. Beyond this moment, the pressure difference becomes negative, and the system exhibits a ‘positive feedback’, meaning that the sheet transfers energy to the fluid, which, in turn, enhances the sheet's motion. This positive feedback accelerates the system's dynamics and creates a spike of pressure drop in the chamber. The process terminates when the system meets the constraint on the maximum volume difference.

To estimate the magnitude of the pressure spike, $\bar {p}_{ud}(t_{p})$, we use the two-mode approximation. Recalling that near the peak time when the sheet is close in shape to the first mode of buckling, i.e. $A_1(t_{p})\simeq 2\varDelta ^{1/2}/{\rm \pi}$ and $A_2(t_{p})\simeq 0$, and that $E\simeq 4{\rm \pi} ^2\varDelta$, we find from (2.20b) and (4.1) an expression for $({{\rm d}^2A_1}/{{\rm d}t^2})(t_{p})$ as a function of the system's parameters. This gives the maximum backward pressure:

(4.5)\begin{equation} \bar{p}_{ud}(t_{p})\simeq{-}\frac{432{\rm \pi}^5L_y\varDelta^{1/2}}{9{\rm \pi}^3\lambda+128\tanh({\rm \pi} L_y/2)}. \end{equation}

This analytical approximation compares well with the numerical solution of (2.2)–(2.9) (see figure 12c). Therefore, in the solid-dominated region, the maximum backward pressure decays to zero as $\bar {p}_{ud}(t_{p},\lambda \gg 1)\simeq -48{\rm \pi} ^2 L_y\varDelta ^{1/2}/\lambda$, whereas in the fluid-dominated region it is independent of $\lambda$, i.e. $\bar {p}_{ud}(t_{p},\lambda \ll 1)\simeq -(27{\rm \pi} ^5/8)L_y\varDelta ^{1/2}/\tanh ({\rm \pi} L_y/2)$.

Note also that the $\bar {p}_{ud}$–$v_{du}$ relation approximately follows the static solution (3.1b) and (3.3b) in the fluid-dominated region (see the dashed lines in figure 12b). This is because the dynamics of the fluid in this region is much slower than that in the solid-dominated region. Nonetheless, deviations between the two solutions, static and dynamic, are observed close to the asymmetric-to-symmetric transition. These deviations may be attributed to inertial effects in the dynamics solution.

5. Discussion on experimental consequences

To fully define the system, eight physical parameters are needed: four parameters specify the properties of the sheet $\{ \tilde {L},\tilde {\rho }_{sh},\tilde {B},\tilde {h} \}$, three parameters define the dimensions of the chamber $\{ \tilde {L}_x,\tilde {L}_y,\tilde {W} \}$ and an additional parameter characterizes the fluid $\tilde {\rho }_{\ell }$. The control parameter is the initial volume difference in the chamber $\tilde {v}_{du}(0)=\tilde {v}_{d}(0)-\tilde {v}_{u}(0)$.

To facilitate comparisons between our theory and experimental observations, we present two of our central predictions in dimensional form. The first prediction relates to the scenario where the initial configuration of the sheet is close in shape to the second buckling mode. In this case, the experimentally measurable quantities are the growth rate $\sigma$, which is defined in (3.6) and characterizes the early stage of the evolution, and the time $t_{p}$, which is given in (4.4) and characterizes the behaviour at moderate times. In dimensional form, these quantities are given by

(5.1a)\begin{gather} \tilde{\varDelta}/\tilde{L}\ll 1 \quad \text{and} \quad \tilde{v}_{du}(0)/ \tilde{v}_{du}^{cr}\ll 1: \tilde{\sigma}=\frac{\sqrt{3}{\rm \pi}^2}{\sqrt{1+\dfrac{8}{{\rm \pi}^2} \dfrac{\tilde{\rho}_{\ell}\tilde{L}_y}{\tilde{\rho}_{sh} \tilde{h}}}} \left(\frac{\tilde{B}}{\tilde{\rho}_{sh} \tilde{h} \tilde{L}^4}\right)^{1/2}, \end{gather}

(5.1b)\begin{gather}\tilde{\sigma}\tilde{t}_{p}\simeq \ln\left(\frac{{\rm \pi}\sqrt{3(15+2{\rm \pi}^2)}}{(3+{\rm \pi}^2)}\frac{\tilde{v}_{du}^{cr}}{\tilde{v}_{du}(0)}\right), \end{gather}

where $\tilde {v}_{du}^{cr}=({2(3+{\rm \pi} ^2)}/{{\rm \pi} \sqrt {3(15+2{\rm \pi} ^2)}})\tilde {\varDelta }^{1/2}\tilde {L}^{3/2}\tilde {W}$ is the volume difference at the asymmetric-to-symmetric transition in the quasi-static solution (§ 3.1), and we keep in mind that the normalized excess length is assumed small in our analysis, i.e. $\tilde {\varDelta }/\tilde {L}=(\tilde {L}-\tilde {L}_x)/\tilde {L}\ll 1$.

The second prediction of our theory corresponds to the frequency of oscillations $\omega$ around the first buckling mode (3.8). In dimensional form the frequency is given by

(5.2)\begin{equation} \tilde{\omega}=\frac{2\sqrt{3}{\rm \pi}^2}{\sqrt{1+\dfrac{128}{9{\rm \pi}^3}\tanh\left(\dfrac{{\rm \pi} \tilde{L}_y}{2\tilde{L}}\right)\dfrac{\tilde{\rho}_{\ell}\tilde{L}}{\tilde{\rho}_{sh} \tilde{h}}}}\left(\frac{\tilde{B}}{\tilde{\rho}_{sh} \tilde{h} \tilde{L}^4}\right)^{1/2}, \end{equation}

where in the small-amplitude approximation the first buckling mode is obtained when $\tilde{v}_{du}(0)=({8}/{{\rm \pi} ^2})\tilde {\varDelta }^{1/2}\tilde {L}^{3/2}\tilde{W}$. It is important to note that this analytical prediction is valid only for very small excess lengths $\tilde {\varDelta }/\tilde {L}\lesssim 0.01$. Larger excess lengths will probably result in significant quantitative deviations from this solution.

To obtain a sense of the physical time and pressure scales that can potentially be induced in the system, let us consider a chamber, with dimensions $\tilde {L}_x=\tilde {L}_y=\tilde {W}=5$ mm, that is filled with water ($\tilde {\rho }_{\ell }\simeq 10^3\,\text {kg}\,{\rm m}^{-3}$). In addition, let us assume that the sheet is made of polyethylene terephthalate with Young's modulus $\tilde {E}\simeq 1$ GPa, Poisson's ratio $\nu \simeq 0.3$ and thickness $\tilde {h}\simeq 0.1$ mm, such that the bending modulus is $\tilde {B}=\tilde {E}\tilde {h}^3/[12(1-\nu ^2)]\simeq 9\times 10^{-5}\,$J (Gomez et al. 2017). The density of the sheet is approximately $\tilde {\rho }_{sh}\simeq 1500\,\text {kg}\,{\rm m}^{-3}$, and its total length is $\tilde {L}=5.5$ mm ($\tilde {\varDelta }/\tilde {L}=0.09$). Under these conditions, the inertial time scale of the sheet is $\tilde {t}_{\star }=(\tilde {\rho }_{sh} \tilde {h} \tilde {L}^4/\tilde {B})^{1/2}\simeq 10^{-3}$ s, and the structure-to-fluid mass ratio is given by $\lambda \simeq 0.03$, i.e. the system is in the fluid-dominated region. Since our theory is limited to inviscid fluids, it is reasonable to assume that the theory should agree with the solution of the more general, viscous equations, in the limit of high Reynolds numbers, where energy dissipation is rather small. Estimating the Reynolds number as ${Re}\sim {\tilde {\rho }_{\ell }\tilde {\bar {v}}(\tilde {t}_{p})\tilde {L}}/{\tilde {\mu }}$, where $\tilde {\mu }\simeq 10^{-3}\,\text {Pa}\,\text {s}$ is the dynamic viscosity and $\tilde {\bar {v}}(\tilde {t}_{p})\sim (\lambda \tilde {\varDelta }/\tilde {L}_y)^{1/2}(\tilde {L}/\tilde {t}_{\star })$ is our scaling for the fluid's velocity at time $\tilde {t}_{p}$ in the fluid-dominated region (see § 4.1), we find that ${Re}\sim 10^3$, i.e. it is relatively high. Yet, we are aware that higher Reynolds numbers should possibly be considered to reveal a convergence to the inviscid limit. Using these parameters, we have from (5.1a) and (5.2) that the growth rate and the frequency of the oscillations are $\tilde {\sigma } \simeq 3.2 /\tilde {t}_{\star }$ and $\tilde {\omega } \simeq 8.5 /\tilde {t}_{\star }$, respectively. In addition, given an initial volume difference, the peak time is obtained from (5.1b). This gives $\tilde {t}_{p}\simeq 3.1 \tilde {t}_{\star }$ if $\tilde {v}_{du}(0)/\tilde {v}_{du}^{cr}=10^{-4}$, and $\tilde {t}_{p}\simeq 1.7 \tilde {t}_{\star }$ if $\tilde {v}_{du}(0)/\tilde {v}_{du}^{cr}=10^{-2}$. At that moment, the average backward pressure difference on the sheet is given by $\tilde {\bar {p}}_{ud}(t_{p})\simeq -160$ kPa; see (4.5).

However, when attempting to predict the behaviour of an experimental system using our solution, it is important to keep in mind the assumptions made in the formulation. For example, we assumed that the fluid exchange between the two parts of the chamber occurs through the upper and lower walls, $y=\pm L_y/2$, but, in practice, this fluid exchange is likely to occur through a connecting channel, as shown in figure 1. In this case, the analytical analysis must take into account the geometry of the transition region and its associated pressure drop. For example, we anticipate that a narrow transition channel will slow down the dynamics and decrease the growth rate of the instability. Additionally, we expect that if the pressure drop in the channel is much smaller than $\bar {p}_{ud}(t_{p})$, it will have a negligible effect on the dynamics.

6. Concluding remarks

We investigated the dynamic interaction between a thin sheet and an inviscid fluid that are confined in a closed rectangular chamber. Our investigation focused on two different regions of the system: the early time evolution, where nonlinear effects are negligible, and the evolution at moderate times, where nonlinearity plays a crucial role in the solution. To analyse the dynamics at $t\ll 1$, we employed a linear stability analysis around the second and first buckling modes. While the second mode is always an unstable state whose highest growth rate is a positive number, the first mode is always stable and yields an imaginary growth rate, i.e. periodic oscillations. In the small-amplitude approximation, $\varDelta \ll 1$, we obtained analytical solutions for the highest growth rate and the smallest oscillation frequency (3.6) and (3.8) respectively, which agree well with the numerical solutions. Yet, experimental data are needed to validate these central analytical predictions. To facilitate comparisons with experiments, we repeated these results in dimensional form in § 5.

Given the chamber's dimensions, we showed that both $\sigma$ and $\omega$ converge to constants in the solid-dominated region and exhibit $\lambda ^{1/2}$ scaling in the fluid-dominated region. This scaling highlights the effect of the fluid on the sheet's motion. When $\lambda \gg 1$, the elastic forces are primarily balanced by the inertia of the sheet, and the fluid has minimal effect on the dynamics. On the other hand, when $\lambda \ll 1$, the elastic forces are mainly balanced by the hydrodynamic pressure difference, and the sheet's inertia has a limited effect on the dynamics. The differences between these two regions are further manifested in the eigenfunctions of the linear stability solution. In the solid-dominated region only one mode of the sheet is excited, i.e. the other modes fall to zero as $1/\lambda$, while an infinite number of modes are excited in the fluid-dominated region. While this difference did not affect the system's behaviour at moderate times, since in the leading order the dynamics is governed by the first two modes, we conjecture that it can influence the system's behaviour at $t\gg 1$; namely, beyond moderate times when higher modes have an increasing effect on the dynamics, we would expect the fluid-dominated solution to be less ordered than the solution in the solid-dominated region.

In addition, we focused on how the sheet escapes from the second buckling mode at moderate times. Key to this analysis is the two-mode approximation that allowed us to analytically analyse the energetic interplay between the sheet and the fluid; see (4.2) and figure 10. This energetic interplay is based on the fact that our model incorporates an elastic sheet and an inviscid fluid, which ensures that the system's total energy is always conserved. At each moment in time, the total energy is distributed between the kinetic energy of the sheet, the potential energy of the sheet and the energy of the fluid, in different proportions. In the solid-dominated region, the sheet's initial potential energy is converted to the sheet's kinetic energy, and only a small fraction, of order $\lambda ^{-1}$, is converted to the energy of the fluid. However, the picture is reversed in the fluid-dominated region, where most energy is used to displace the fluid.

The time at which the potential energy of the sheet reaches a minimum (4.4) constitutes another central result of our theory, which can be verified experimentally. In particular, we showed that $\sigma t_{p}$ is almost independent of the parameter $\lambda$, and in the limit $v_{du}(0)\ll 1$ it diverges logarithmically.

During the dynamic evolution at moderate times, we observed that the sheet–fluid interplay experiences a transition from negative to positive feedback. Initially, the fluid resists the sheet's motion, but at later times, the hydrodynamic pressure difference acts in the direction of the sheet's motion and promotes the sheet's dynamics. The positive feedback ends at $t=t_{p}$, when the volume difference reaches its maximum value, dictated by the inextensibility of the sheet. At that moment, the pressure difference on the sheet, $\bar {p}_{ud}$, reaches its peak value.

An important extension of the present theory is the inclusion of viscosity in the mathematical formulation of the fluid. This extension will allow us to investigate the behaviour of the sheet at both high and low Reynolds numbers and thus look for the elastohydrodynamic instabilities caused by viscous effects. It will also allow us to investigate the formation of boundary layers and examine their impact on the dynamics of the system. We will pursue this extension in a future study.