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.

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 [26].In that setup, 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.
In the above mentioned study, Oshri [26] analyzed 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 setup will open new avenues for designing advanced technological devices, such as micro-mechanical switches [27][28][29] and microfluidic mixing devices [30][31][32].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 required for the sheet to release its stored energy, thereby increasing, for example, the timescales 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 [5,28,33], 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 to the lateral dimension of the chamber, ∆, 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, λ.We show that for fixed dimensions of the chamber, L y and ∆, the system exhibits two asymptotic solutions as a function of λ.The sheet's inertia dominates the dynamics when λ ≫ 1, and is therefore referred to below as the "solid-dominated" region, while the dynamics is governed by the fluid's inertia when λ ≪ 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, σ, and the lowest frequency of vibration, ω, around these initial states.The two solutions exhibit similar behaviour as a function of λ, namely, they converge to a constant in the solid-dominated region, while they exhibit the scaling λ 1/2 in the fluiddominated 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., ∆ ≪ 1 [34].
At moderate times, the 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 analyze 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 σt p ≃ ln c∆ 1/2 /v du (0) , where σ 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 a kinetic energy of the sheet if λ ≫ 1, and into a hydrodynamic energy if λ ≪ 1.We show that at t = t p relatively large spike of pressure drop is applied on the sheet.
The paper is organized as follows.In § II, 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 § III, 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 § IV, 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 § V, we discuss a possible experimental realization of the system, and in § VI we draw conclusions, and propose a direction for future study.

II. Formulation of the problem
We consider an inextensible thin sheet of total length L, bending modulus B, thickness h, and density ρ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 Lx , Ly , and W , respectively.A Cartesian coordinate system is located on the left edge of the sheet.A cross-section of the chamber on the xỹ plane is placed at 0 ≤ x ≤ Lx and − Ly /2 ≤ ỹ ≤ Ly /2.When t < 0, the valve connecting the two parts of the chamber is closed, and the volumes above and below the sheet, ṽu ( t) and ṽd ( t), are filled with an incompressible, inviscid fluid of density ρℓ .Hereafter, we denote quantities related to the upper and lower parts of the chamber by the subscripts 'u' and 'd', respectively.At t ≥ 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, L, and we normalize time by the inertial time-scale of the sheet t⋆ = L2 (ρ sh h/ B) 1/2 , i.e., 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 are driven by the sheet's motion, we chose the sheet's inertial timescale 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 tilde over the symbol, and the corresponding nondimensional 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 to the total volume of the chamber, i.e., h L/( Lx Ly ) ≪ 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 = ±L y /2.Fourthly, we assume that the vertical dimension of the chamber, L y , is larger than the typical length scale, ℓ, 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 1. Schematic overview of the system.A thin sheet of total length L, bending modulus B, density ρsh , and thickness h divides a closed rectangular chamber of dimensions Lx × Ly into two parts.The excess length of the sheet compared to the lateral dimension of the chamber is given by ∆ = L − Lx (not shown in the figure).The volumes of the chamber above and below the sheet, ṽi(t) (i =u,d), are filled with an inviscid and irrotational fluid of density ρℓ .At t ≥ 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 setup that corresponds to the above model involves a valve-controlled channel that connects the two compartments of the chamber.
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 ϕ i (x, y, t), where i=u,d, from which we can determine the velocity profile of the fluid as v i = ∇ϕ i , where ∇ 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 ϕ i = φi (ρ sh h/ B) 1/2 , and the pressures, by p i = pi L3 / B. The evolution of these hydrodynamic fields, in space and over time, is determined by the continuity equation and Bernoulli's equation: 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 [35].In addition, in Eqs.(2b) and ( 5), we define the dimensionless parameters: The structure-to-fluid mass ratio, λ, 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 [36][37][38][39].The parameter ∆ 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 ∆ = L − Lx .For a given system, the parameters λ and ∆ remain constant throughout the dynamic evolution.
To solve the continuity equation, Eq. (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 = ±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: , t , and In addition to the periodic boundary conditions at y = ±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, Eq. (2a), and the periodic boundary conditions, it becomes apparent that this requirement is satisfied when c d (t) = c u (t) ≡ 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 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 ∈ [0, 1].With this change of reference frame, we define the position vector to a point on the sheet as 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 θ(s, t); see figure 1.These three elastic fields, i.e., x sh (s, t), y sh (s, t) and θ(s, t), are not independent, since they are related by the geometric constraints: By using these definitions, the kinematic boundary conditions on the sheet-fluid interfaces are given by: where D/Dt = ∂/∂t + v i • ∇ is the two-dimensional convective derivative.The balance of moments and forces on the sheet is given by: where 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 F = F L2 / B. In addition, nd = (− sin θ, cos θ) is a unit normal vector on the sheet that points outwards from the lower part of the chamber, and the hydrodynamic pressures in Eq. (8b) are calculated on their respective sides of the sheet-fluid interfaces.Note that in Eq. (8a) we neglect the rotational inertia term.This is justified in the limit of a thin and inextensible sheet, as assumed in this analysis [40][41][42].Equations (8) are supplemented by the following boundary conditions on the sheet's edges: where we assume hinged boundary conditions in Eq. (9c).This completes the formulation of the problem.In summary, given the excess length ∆, the vertical dimension of the chamber L y , the parameter λ, 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)- (9).In the analysis that follows, we will always assume that the sheet and the fluid are initially at rest, i.e., ∂x sh ∂t (s, 0) = 0 and ϕ 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., ∆ ≪ 1.For this reason, in the next section, we reduce our model to this so-called small-amplitude approximation [34] 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 p sh (t) and E k sh (t), respectively, and E f (t) accounts for the kinetic energy of the fluid.Therefore, the total energy is given by: where | .| corresponds to the norm of the enclosed vector, and our normalization implies that E = Ẽ L/ B. Since our system starts from rest, the total energy of the system equals the initial potential energy of the sheet, E = E p sh (0), and this energy is conserved throughout the system's evolution.

A. The small-amplitude approximation
The assumption that the amplitude of the sheet remains small, or equivalently that ∆ ≪ 1, implies that the geometric relations, Eq. ( 6), reduce to ∂y sh /∂s ≃ θ and ∂x sh /∂s ≃ 1 − 1 2 (∂y sh /∂s) 2 .The non-linearity 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, Eq. (9a).Indeed, in the small-amplitude approximation, this constraint is given by: Here, we replace the arclength coordinate of the sheet with the Eulerian coordinate of the fluid, s ≃ x, according to our level of approximation.Correspondingly, the balance of forces and moments on the sheet, Eq. ( 8), reduces to: 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 Eq. (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) ∝ ∆, and that the total energy of the system is conserved, the energy of the fluid is, at most, proportional to E f ∼ ∆.Therefore, if we approximate the energy of the fluid by E f ∼ |v| 2 ℓ, where |v| is the typical velocity in the chamber and ℓ is the decay length of the disturbances in the flow, we obtain |v| ∼ ∆/ℓ.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, Eq. (2b), and the kinematic boundary conditions by: These approximations will further be verified a posteriori in § IV, where we analyze the nonlinear dynamics of the system.In particular, we will compare the results of this approximation with the numerical solution of the nonlinear model, Eqs. ( 2)- (9).Note that since the continuity equations, Eq. (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, Eqs.(2a) and ( 11)-( 13), supplemented by the linearized form of the boundary conditions, Eqs. ( 5), ( 4), (9b) and (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 S = T 0 Ldt, where with respect to the elastic fields y sh (x, t) and F x (t) and the hydrodynamic fields ϕ 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.

Modal expansion
The continuity equations, Eq. (2a), and their corresponding boundary conditions on the fluid-chamber interfaces, Eq. ( 4), are satisfied when the potential functions are given by: where a m (t) (m = 0, 1, 2, ...) and c m (t) (m = 1, 2, 3, ...) are unknown time-dependent coefficients, and the ± 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: where the functions sin(πnx) automatically satisfy the boundary conditions on the sheet edges, Eqs.(9b) and (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, Eq. ( 12), Bernoulli's equation, Eq. (13a), the kinematic boundary conditions, Eq. (13b), and the geometric constraint, Eq. (11).Equation ( 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, Eq. ( 14).To this end, we follow the analysis in Appendix C and substitute the potential functions, Eq. ( 15), and the height function, Eq. ( 16), into the Lagrangian, Eq. ( 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: where Einstein's summation rule is implied for repeated indices, and we define the following symmetric matrices: where δ nk is the Kronecker delta, and for n ̸ = 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/λ.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 [43][44][45].
The second comment is related to the potential functions ϕ i (x, y, t) that result from the minimization.Using Eq. ( C2b), we substitute c m (t) = −a m (t) in the potential functions, Eq. ( 15), to obtain This solution implies that at y = ±L y /2 the velocity of the fluid is oriented only in the y-direction, i.e., ∂ϕ i /∂x(x, ±L y /2, t) = 0.It also implies that ϕ i (x, ±L y /2, t) = 0, which, given Eqs.( 5) and (13a), yields a constant zero pressure along the inlet and the outlet walls of the chamber, p i (x, ±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 ≫ ℓ, where ℓ ≃ 1/(πm), for the smallest nonzero 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., (dA n /dt)(0) = 0, we can determine the dynamic evolution of the system from the minimization of Eq. ( 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: Once A n (t) are determined from the solution of Eq. ( 20), the position of the sheet in time and in space is given by Eq. ( 16), and the hydrodynamic potentials and the pressure fields are determined from Eqs. ( 15), (C2) and (13a).In the next section, we utilize this formulation to investigate the early time evolution.Then, in § IV we use it to analyze the dynamics at later times.

III. 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.

A. Recap of the quasi-static solution
Following the analysis in the study of Oshri [26], 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 ≤ v du (0) ≤ v cr du , where ∆ 1/2 , the system is governed by the asymmetric branch.In this branch, the lateral compression is constant, F x (0) = 4π 2 , and the height functions are given by: 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: Note, that when v du (0) → 0 we know from Eq. (21b) that the pressure difference vanishes, p ud (x, y, 0) → 0, and the elastic configuration converges to the second, asymmetric, mode of buckling, y sh (x, 0) → ∆/π 2 sin(2πx).The total energy of the system in this configuration is given by E as = 4π 2 ∆.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, Eq. ( 21), around the horizontal axis.
On the other hand, when the volume difference is set as v cr du ≤ v du (0) < 2∆/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 2∆/3.The height functions in this branch are given by the parametric solution: where u = F x (0)/2 is a function of the lateral compression.Given the initial volume difference, v du (0), and the excess length, ∆, we can determine the lateral compression, F x (0), from Eq. (23c), and then substitute this solution into Eqs.(23a) and (23b) to obtain the height profile.When p ud (x, y, 0) → 0, the height function converges to the first, symmetric, mode of buckling, which is given by y sh (x, 0) → 4∆/π 2 sin(πx) and v du (0) = 8∆ 1/2 /π 2 .The total elastic energy of this shape is given by E s = π 2 ∆.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 ≤ v du (0) < 2∆/3, is plotted in figure 2. In the following sections, we will use these height functions, Eqs.(21a) and (23a), as the base solutions for our perturbative time-dependent expansion.
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 Eqs.(21b) and (23b), 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 Eqs. (21a) and (23a) in the limit p ud (x, y, 0) → 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.

B. 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 x sh (x,0) Evolution of the static solution in the small-amplitude approximation, where ∆ = 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 Eq. (21b), while in the symmetric branch it is given by Eqs. the base solutions, i.e., A n (t) = A n (0) + ϵ Ān e σt and F x (t) = F x (0) + ϵ Fx e σt , where Ān and Fx are unknown constants, σ is the growth rate, and ϵ ≪ 1 is an arbitrarily small parameter.Then, we substitute these perturbed functions into the equations of motion, Eq. ( 20), and expand them up to linear order in ϵ.The leading order of this expansion, order ϵ 0 , is given by: and the subleading order, order ϵ, is given by The above equations in the subleading order always have the trivial solution Ān = 0 and Fx = 0, unless their corresponding determinant vanishes.This condition gives the growth rate, σ.Once σ is determined, its corresponding eigenfunction is obtained from the solution of Eq. ( 25).The hydrodynamic fields related to this eigenfunction are determined from Eqs. ( 15) and (C2).

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 Eq. (21a) in the limit p ud (x, y, 0) → 0. This solution reads y sh (x, 0) = ∆/π 2 sin(2πx) and F x (0) = 4π 2 .A projection of this configuration on the normal mode expansion, Eq. ( 16), gives A n (0) = ∆/π 2 δ 2n .As expected, this initial state exactly satisfies the equilibrium equations at order ϵ 0 , Eq. (24).At the next order, i.e., order ϵ, we find that Fx = 0 and Ān = 0 for all the even perturbations, i.e., n = 2, 4, 6... .Consequently, Eq. ( 25) yields linear and homogeneous equations that involve only the odd perturbations.These equations always has the trivial solution Ān = 0, except when its 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: In figure 3, we plot this analytical approximation for the growth rate as a function of λ, and compare it with the numerical solution of Eq. ( 25) for the case where N = 8.In addition, we compare this analytical solution with the growth rate obtained from the linearization of Eqs. ( 2)-( 9), i.e., where ∆ is assumed to be finite; see Appendix D for the details of this solution.2)-( 9), and solid and dashed lines correspond to the growth rates obtained from the small-amplitude approximation.When λ ≫ 1, the growth rate converges to the constant σ ≃ √ 3π 2 , whereas when λ ≪ 1 the growth rate is given by σ ≃ 3π 6 /8(λ/Ly) 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 ∆.
Equation ( 26) 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 Eq. ( 26) depends on the parameter λ/L y = ρsh h/(ρ ℓ Ly ), this result is not general but depends on the order of approximation.Had we solved Eq. ( 25) with N ≥ 3, the two parameters λ 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 Eq. ( 26) well describes the growth rate in the limit ∆ ≪ 1.
Thirdly, while the solution of the small-amplitude approximation is independent of the excess length, ∆, the more general solution for finite values of ∆ does depend on this parameter; see figure 3 for a comparison.In particular, for a fixed value of λ, the growth rate increases with an increase in ∆.
Fourthly, the analytical solution of the growth rate, Eq. ( 26), exhibits two different regions as a function of λ.When λ/L y ≫ 1, the growth rate is a constant, σ ≃ √ 3π 2 , that coincides with the growth rate of a sheet that is uncoupled from an external fluid.However, when λ/L y ≪ 1, the growth rate exhibits the scaling σ ≃ 3π 6 /8(λ/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 σ ∝ (λ/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 /(π 2 λ), in the denominator of the growth rate.When λ is large, the added mass approaches zero and the sheet's inertia is not affected by the fluid's motion.In contrast, when λ 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 Eq. ( 21).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 Eqs. ( 25) and ( 18) become diagonal in the solid-dominated region.Therefore, up to small corrections of the order 1/λ, Ā1 alone is excited at the instability.Indeed, in figure 4(a), in which we plot the eigenfunction for the case where λ = 100 for both N = 2 and N = 8, we find that the two solutions are almost identical.In contrast, in the fluid-dominated region, λ ≪ 1, and the matrices in Eq. ( 25) have nonzero off-diagonal terms.Therefore, all the odd modes become coupled and are excited at the instability.Nonetheless, our numerical investigation indicates that Ān / Ā1 ≪ 1 for all n ≥ 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, Eq. (21a), when we expand Eq. (21a) 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  Sixthly, in Fig. 5(a), we plot the flow field obtained from the linear stability analysis at finite ∆.Note that the maximum velocity of the fluid is obtained at the sheet's center, 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 ≥ 3, are included.This is because each coefficient a i m (t) in the fluid's potential functions, Eq. ( 15), depends on all the excited modes; see Eq. (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 fluid-dominated region.
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 β ≡ dv du /dp ud , where pud (t) is the average pressure drop on the sheet.In the small-amplitude approximation the average pressure drop on the sheet is given by pud (t) = 1 0 [p u (x, 0, t) − p d (x, 0, t)] dx and the volume difference in the chamber is given by v du (t) = 2 1 0 y sh (x, t)dx [26].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)/ Ā1 = (4/π)e σt .In addition, the average pressure difference at this order is given by pud (t)/ Ā1 = 2L y σ 2 e σt /(πλ).Consequently, at the onset of the instability, the compressibility is a constant, which is given by: In figure 6, we plot this result as a function of λ and compare the lowest order approximation with the solution at finite values of ∆.On the one hand, since the growth rate is constant in the solid-dominated region, we find that the compressibility scales linearly with λ, i.e., β → 2λ/(3π 4 L y ).On the other hand, since in the fluid-dominated region, the growth rate scales as (λ/L y ) 1/2 , the compressibility converges to the constant β → 16/(3π 6 ) ≃ 5.54 × 10 −3 .As expected, this result for the fluid-dominated region is very close to the compressibility obtained in the quasi-static solution, Eq. (21b), which gives  6. Log-log plot of the compressibility as a function of λ close to the onset of the instability.Symbols correspond to the compressibility at finite values of ∆, and the solid black line corresponds to the analytical solution obtained from the small-slope approximation.While in the solid-dominated region β ∝ λ, in the fluid-dominated region, β converges to a constant.

Linear stability around the first mode of buckling
This section analyzes 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 σ ≡ iω, where i = √ −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 Eq. ( 23) that the base solution is given by y(x, 0) = ( 2√ ∆/π) sin(πx) and F x (0) = π 2 .A projection of this configuration on the normal mode expansion, Eq. ( 16), gives A n (0) = ( 2√ ∆/π)δ 1n .This initial state, as expected, satisfies the leading order of our perturbative expansion, Eq. (24).Substituting this leading order in Eq. ( 25) and solving for the unknown constants, we find that Fx = 0 and Ān = 0 for all the odd modes (n = 1, 3, 5...).Consequently, Eq. ( 25) 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 [46].This gives, In figure 7(a), we plot this solution and compare it with the solution of Eq. ( 25) for the case of N = 8.Note that our solution, Eq. ( 28) with N = 2, approximates well the small-amplitude limit, ∆ ≪ 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 Eq. ( 28) with the eigenvalues obtained from the linearization of Eqs. ( 2)-( 9), where the small-amplitude approximation is relaxed.We observe that as ∆ increases, the oscillation frequency decreases with increasing excess length.However, the overall trend of the dependence of ω on λ remains consistent with the small-amplitude solution.
We also note that the relatively large decrease in ω when λ ≫ 1 is a result of our chosen hinged boundary conditions.In contrast, systems with clamped boundary conditions display a much milder dependence on ∆, as reported in previous studies [42,47].x sh (x,t) y sh (x,t) FIG. 7. The lowest oscillation frequency and the sheet's eigenfunction obtained from the linear stability analysis around the first buckling mode.In both panels Ly = 2. (a) Log-log plot of the oscillation frequency as a function of the parameter λ.Solid and dashed lines correspond to the solution of Eq. ( 28) when N = 2 and N = 8, respectively.All symbols correspond to the linear stability analysis of Eqs. ( 2)-( 9) at finite ∆.(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 ŷsh (1/4) = 1; see Appendix D).Although only one mode is excited when λ ≫ 1 and infinite modes are excited when λ ≪ 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 ω exemplifies the two limiting scenarios that we encountered for the growth rate in Eq. ( 26).On the one hand, in the solid-dominated region, where λ ≫ 1, the frequency converges to the constant, ω ≃ 2 √ 3π 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 λ ≪ 1, we find the scaling, ω ≃ 27π 7 /32 [λ/ tanh(L y π/2)] 1/2 , i.e., ω ∝ λ 1/2 .
Alternatively, since the added mass is given by 128 9π 3 tanh(πLy/2) λ , when λ ≫ 1 the sheet's inertia is almost unaffected by the fluid motion, but when λ ≪ 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, Eq. ( 25) becomes diagonal, up to corrections of the order of 1/λ, and essentially only the second mode, i.e., Ā2 , is excited at the instability, while in the fluid-dominated region Eq. ( 25) has nonzero off-diagonal terms, and all the even modes are excited.Nonetheless, our numerical investigation of the solution of Eq. ( 25) in the fluiddominated region indicates that the ratio Ān / Ā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 ω, Eq. ( 7), 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.
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 7(b)).Since only even modes are excited at the instability, the average pressure difference on the sheet, pud (t), vanishes.Therefore, in this case, there is no analogue to the compressibility calculated in § III B 1.

IV. The evolution at moderate times
In this section, we relax the assumption that t ≪ 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) ≪ v cr du , where ∆ 1/2 (see § III A), 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 pud (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 § II A. This formulation yields the simplified set of nonlinear equations, Eq. ( 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 dA k dt dAn dt + V nk A k A n , that corresponds to the total energy of the system, Eq. ( 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: We keep in mind that A 1 (t) and A 2 (t) are related through the constraint of the excess length, Eq. (20b).
To obtain some insight regarding the validity of this two-mode approximation, we solve Eq. ( 20) numerically with N = 2 and compare the results with the solution of the nonlinear model, i.e., the numerical solution of Eqs. ( 2)-( 9).In our investigation, we set ∆ = 0.01 and L y = 2, and consider two different values for the parameter λ.The initial configuration is given by Eq. ( 21), where v du (0) = 0.01v cr du [48].The results of these numerical solutions are presented in figures 9(a) and 9(b), where we follow the time-dependent behaviour of the mid-point 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 (λ = 100), the two-mode approximation holds over one period of motion, while in the fluid-dominated region (λ = 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 § III B 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 Eq.( 20) with N = 3 instead of N = 2, we observe a convergence towards the numerical data for longer times; see the dashed-gray line in figure 9(b).
Nonetheless, in both cases, i.e., the solid-and fluid-dominated regions, the agreement between the numerical solution of Eqs. ( 2)-( 9) and the analytical approximation, Eq. ( 20) with N = 2, holds up to the first maximum.Similar results are also obtained when we perturb the system's parameters, ∆ and L y , and the initial configuration.Therefore, in the following analysis, we will 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.x sh (x,t) 9.The sheet's mid-point as a function of time in (a) the solid-dominated (λ = 100) and (b) the fluid-dominated (λ = 0.1) regions.In both panels, ∆ = 0.01, Ly = 2, and the growth rate, σ, is approximated by Eq. ( 26).The solid black line denotes the two-mode approximation, i.e., Eq. ( 20) with N = 2, and the open blue circles denote the solid black line denotes the solution of the nonlinear model, Eqs. ( 2)-( 9).In panel (b), the dashed gray line denotes the solution of Eq. ( 20) with N = 3.The initial configuration of the sheet is given by Eq. ( 21) with v du (0) = 0.01v cr du .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 panel (b); see the corresponding markers in panel (b).Between 1 ○- 3  ○ the sheet releases potential energy as it transforms from the second mode of buckling to the first mode of buckling.After 3 ○, the height of the sheet's mid-point decreases and the sheet gains back potential energy, as seen in 4 ○.

A. The elasto-hydrodynamic 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 Eq. ( 22).In the limit v du (0) ≪ v cr du , this energy is approximated as E as ≃ 4π 2 ∆, 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 = π 2 ∆.Therefore, at most, our system can convert δE = E as − E s ≃ 3π 2 ∆ 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 k sh (t) and E p sh (t), and the energy of the fluid, E f (t), are given by the first, second, and third terms, respectively, in the right-hand side of Eq. ( 10).In the small-amplitude approximation, these energies reduce to 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 Eq. ( 20) with N = 2 in the solid-and the fluid-dominated regions, λ = 10 and λ = 0.1 respectively, where the initial configuration is given by Eq. ( 21) with v du (0) = 0.01v cr du .In both cases, we find that, after some initial delay, the potential energy of the sheet drops from E ≃ 4π 2 ∆, i.e., the energy of the second mode of buckling, Eq. ( 22), to E ≃ π 2 ∆, i.e., the energy of the first mode of buckling.The energy released in this process, δE ≃ 3π 2 ∆, 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 10(a)) the kinetic energy of the sheet becomes much larger than the energy of the fluid, while as λ decreases, the opposite picture emerges, i.e., the energy of the fluid becomes much larger than the kinetic energy of the sheet (figure 10(b)).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).
Our two-mode approximation allows us to quantify these findings and to further estimate the maximum values of E k sh (t p ) and E f (t p ) as a function of λ, where t p denotes the time at which the potential energy of the sheet reaches its minimum value.Indeed, at E p sh (t p ) = π 2 ∆, the sheet is close in shape to the first mode of buckling, i.e., A 1 (t p ) ≃ 2 √ ∆/π.The constraint on the excess length, Eq. (20b), then implies that A 2 (t p ) ≃ 0 and that (dA 1 /dt)(t p ) ≃ 0. In addition, the derivative of the second mode at that moment, (dA 2 /dt)(t p ), is obtained from Eq. ( 29) when we substitute E = 4π 2 ∆ for the total energy.Taken together, these approximations yield the following maximum energies at t = t p : In figure 10(c), we compare these maximum energies with the numerical solution of Eqs. ( 2)-( 9).Overall, we find a good fit between the analytical approximation and the numerical solution over the entire range of λ.Using Eq. ( 30), we find that in the solid-dominated region, λ ≫ 1, most of the initial energy is converted into the kinetic energy of the sheet, i.e., E k sh (t p ) ≃ 3π 2 ∆ − δ 1 and E f (t p ) ≃ δ 1 , where δ 1 = 128∆ tanh(πL y /2)/(3πλ), while in the fluid-dominated region, λ ≪ 1, most of the energy is converted into the energy of the fluid, i.e., E k sh (t p ) ≃ δ 2 and E f (t p ) ≃ 3π 2 ∆ − δ 2 , where δ 2 = 27π 5 ∆λ/ [128 tanh(πL y /2)].Furthermore, if we estimate the average instantaneous velocity of the fluid, v(t p ), by using E f (t p ) ∝ v(t p ) 2 ℓ/λ, we find that in the solid-dominated region the average velocity is proportional to v(t p ) ∝ (∆/ℓ) 1/2 , while in the fluid-dominated region it is proportional to v(t p ) ∝ (λ∆/ℓ) 1/2 ; namely, while E f (t p , λ ≪ 1) ≫ E f (t p , λ ≫ 1), their corresponding velocities present the opposite relationship, i.e., v(t p , λ ≪ 1) ≪ v(t p , λ ≫ 1).This is because the momentum of the fluid, p f ∝ v(t p )/λ -but not its velocity -increases in the fluid-dominated region.λ Eq. (4.2b) Eq. (4.2a) FIG. 10.The energetic interplay between the three components of the total energy in (a) the solid-dominated (λ = 10) and (b) the fluid-dominated (λ = 0.1) regions of the system.In these two panels, we solve Eq. ( 20) for the case where N = 2, ∆ = 0.01, Ly = 2, and σ is given by Eq. ( 26).The initial configuration of the sheet is given by Eq. ( 21), where v du (0) = 0.01v cr du .After some initial delay, the potential energy of the sheet drops from E p sh (t ≪ 1) ≃ 4π 2 ∆ to E p sh (tp) ≃ π 2 ∆.Open blue circles represent the potential energy obtained from the solution of the nonlinear equations, Eqs. ( 2)-( 9).The energy released from the sheet is divided between the kinetic energy of the sheet, E k sh (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 k sh (tp) ≫ E f (tp), while in the fluid-dominated region we find the opposite relation, The kinetic energy of the sheet and the fluid at t = tp as a function of λ, i.e., Ei(tp) = E k sh (tp), E f (tp).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 color symbols correspond to the numerical solution of Eqs. ( 2)-( 9), where the initial conditions are similar to those used in panels (a) and (b).

B. 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., Eq. ( 21) with v du (0) ≪ v cr du , we aim to find the time t = t p at which the potential energy of the sheet first drops to the minimum value, E p sh (t p ) ≃ π 2 ∆.To do so, we first use Eqs.(20b) and ( 26) to eliminate A 2 (t) and λ in favor of A 1 (t) and σ, respectively, in Eq. (29).Then, we substitute the energy of the initial configuration, Eq. ( 22), into Eq.( 29) and integrate it between t ∈ [0, t p ].This gives: where A 1 (t p ) ≃ 2∆ 1/2 /π is approximately the amplitude of the sheet at time t p , and A 1 (0) = 2 1 0 y sh (x, 0) sin(πx)dx = 32v du (0)/(3π + π 3 ) is the projection of the initial configuration, Eq. ( 21), on the first mode of the sheet.
In figure 11(a), we fix the excess length and the vertical dimension of the chamber at ∆ = 0.01 and L y = 2, respectively, and plot σt p for λ ∈ [10 −3 , ).We find that at a given initial volume difference, v du (0), the peak time, σt p , changes by less than five percent over more than six orders of magnitude in λ.While σt p is almost independent of λ, it does depend strongly on the initial configuration of the sheet, v du (0), and the excess length, ∆.An analytical approximation of this dependence can be extracted from Eq. ( 31), if we assume that the system is in the solid-dominated region, where σ ≃ √ 3π 2 .Under this assumption, we can integrate the right-hand side of this equation and take the limit v du (0) ≪ 1 of the resulting expression [49].This gives: where c ≃ 2.9.While the scaling in Eq. ( 32) 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 ≃ 2.0; see figure 11(b).We note that the independence of L y in the right-hand side of Eq. ( 32) is a result of our assumption that λ ≫ 1. Had we derived this scaling using the fluid-dominated region, we would have found that the integral in Eq. ( 31) does depend on the vertical dimension of the chamber.However, this dependence diminishes to zero when L y > ∼ 1, as is required by our fourth assumption in § II.
We also note that the logarithmic divergence of σt p in the limit v du (0) → 0 is similar to the divergence of the period of a pendulum near the separatrix [50].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.31), and symbols with corresponding colors represent the numerical solution of Eqs. ( 2)- (9).For a given v du (0), the time σtp changes by less than 5% over six orders of magnitude of the parameter λ.(b) Comparison between the analytical scaling, Eq. (32), 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 ≃ 2.9 slightly overestimates the numerical prediction.

C. The pud -v du relation
In this section, we investigate the behaviour of the pud -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, Eq. (13a), and the normal mode expansion, Eqs. ( 15) and (C2), we have that the average pressure difference on the sheet is given by pud (t) = 2Ly πλ d 2 A1 dt 2 .In addition, using Eq. ( 16), we find the volume difference as a function of time, v du (t) = 2 1 0 y sh (x, t)dx = 4A 1 (t)/π.Thereafter, we solve Eq. ( 20) numerically with N = 2 and plot the parametric solution (p ud (t),v du (t)) in the range t ∈ [0, t p ].
The results of these solutions are plotted in figures 12(a) and 12(b), for the solid-(λ = 100) and fluid-(λ = 0.01) dominated regions, respectively.Qualitatively, the two regions exhibit similar behavior.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 pud (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 12(b)), than in the solid-dominated region (figure 12(a)).
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, pud = 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, pud (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 ) ≃ 2∆ 1/2 /π and A 2 (t p ) ≃ 0, and that E ≃ 4π 2 ∆, we find from Eqs. (20b) and ( 29) an expression for d 2 A1 dt 2 (t p ) as a function of the system's parameters.This gives the maximum backward pressure, This analytical approximation compares well with the numerical solution of Eqs. ( 2)-( 9); see figure 12(c).Therefore, in the solid-dominated region, the maximum backward pressure decays to zero as pud (t p , λ ≫ 1) ≃ −48π 2 L y ∆ 1/2 /λ, whereas in the fluid-dominated region it is independent of λ, i.e., pud (t p , λ ≪ 1) ≃ −(27π 5 /8)L y ∆ 1/2 / tanh (πL y /2).Note also that the pud -v du relation approximately follows the static solution, Eqs.(21b) and (23b), in the fluiddominated region; see the dashed lines in figure 12(b).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.In both panels, the solid black line corresponds to the solution of Eq. ( 20) with N = 2, ∆ = 0.01, Ly = 2, and v du (0) = 0.01v cr du .The blue points correspond to times t = 0 and t = tp.The dashed gray lines correspond to the static solution, Eqs.(21b) and (23b).The insets show pud (t) as a function of time, where σ is given by Eq. ( 26).In the fluid-dominated region, the maximum backward pressure (p ud (tp) ≃ −200) is much larger than that in the solid-dominated region (p ud (tp) ≃ −1).In addition, the evolution of the pud -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.

V. Discussion on experimental consequences
To fully define the system, eight physical parameters are needed: Four parameters specify the properties of the sheet { L, ρsh , B, h}, three parameters define the dimensions of the chamber { Lx , Ly , W }, and an additional parameter characterizes the fluid ρℓ .The control parameter is the initial volume difference in the chamber ṽdu (0) = ṽd (0)− ṽ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 σ, which is defined in Eq.( 26) and characterizes the early stage of the evolution, and the time t p , which is given in Eq.( 32) and characterizes the behaviour at moderate times.In dimensional form, these quantities are given by: ∆/ L ≪ 1 and ṽdu (0)/ṽ cr du ≪ 1 : where ṽcr ∆1/2 L3/2 W is the volume difference at the asymmetric-to-symmetric transition in the quasistatic solution ( § III A), and we keep in mind that the normalized excess length is assumed small in our analysis, i.e., ∆/ L = ( L − Lx )/ L ≪ 1.
The second prediction of our theory corresponds to the frequency of oscillations ω around the first buckling mode, Eq. ( 28).In dimensional form the frequency is given by: ω where in the small-amplitude approximation the first buckling mode is obtained when ṽdu (0) = 8 π 2 ∆1/2 L3/2 W .It is important to note that this analytical prediction is valid only for very small excess lengths ∆/ L < ∼ 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 Lx = Ly = W = 5 mm, that is filled with water (ρ ℓ ≃ 10 3 kg/m 3 ).In addition, let us assume that the sheet is made of polyethylene terephthalate with Young's modulus Ẽ ≃ 1 GPa, Poisson's ratio ν ≃ 0.3, and thickness h ≃ 0.1 mm, such that the bending modulus is B = Ẽh3 /[12(1 − ν 2 )] ≃ 9 × 10 −5 J [10].The density of the sheet is approximately ρsh ≃ 1500 kg/m 3 , and its total length is L = 5.5 mm ( ∆/ L = 0.09).Under these conditions, the inertial timescale of the sheet is t⋆ = (ρ sh h L4 / B) 1/2 ≃ 10 −3 s, and the structure-to-fluid mass ratio is given by λ ≃ 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 ∼ ρℓ ṽ( tp) L μ , where μ ≃ 10 −3 Pa • s is the dynamic viscosity, and ṽ( tp ) ∼ (λ ∆/ Ly ) 1/2 ( L/ t⋆ ) is our scaling for the fluid's velocity at time tp in the fluid-dominated region (see § IV A) we find that Re ∼ 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 Eqs. (34a) and ( 35) that the growth rate and the frequency of the oscillations are σ ≃ 3.2/ t⋆ and ω ≃ 8.5/ t⋆ , respectively.In addition, given an initial volume difference, the peak time is obtained from Eq. (34b).This gives tp ≃ 3.1 t⋆ if ṽdu (0)/ṽ cr du = 10 −4 , and tp ≃ 1.7 t⋆ if ṽdu (0)/ṽ cr du = 10 −2 .At that moment, the average backward pressure difference on the sheet is given by pud (t p ) ≃ −160 KPa; see Eq. (33).
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 = ±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 pud (t p ), it will have a negligible effect on the dynamics.

VI. 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 analyze the dynamics at t ≪ 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, ∆ ≪ 1, we obtained analytical solutions for the highest growth rate and the smallest oscillation frequency, Eqs. ( 26) and ( 28) respectively, which agree well with the numerical solutions.Yet, experimental data is needed to validate these central analytical predictions.To facilitate comparisons with experiments, we repeated these results in dimensional form in § V.
Given the chamber's dimensions, we showed that both σ and ω converge to constants in the solid-dominated region and exhibit λ 1/2 scaling in the fluid-dominated region.This scaling highlights the effect of the fluid on the sheet's motion.When λ ≫ 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 λ ≪ 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/λ, 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 ≫ 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 analyze the energetic interplay between the sheet and the fluid; see Eq. ( 30) 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 the order λ −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, Eq. ( 32), constitutes another central result of our theory, which can be verified experimentally.In particular, we showed that σt p is almost independent of the parameter λ, and in the limit v du (0) ≪ 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, pud , 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 elasto-hydrodynamic 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.
A. Derivation of Eq. ( 10) Equations ( 2)-( 9) have a conserved first integral that corresponds to the total energy in the system.To derive this conserved quantity, we multiply Eq. (8a) by ∂θ/∂t, and Eq.(8b) by ∂x sh /∂t, and subtract the second equation from the first.Then, we integrate the resulting equation between s ∈ [0, 1], and use integration by parts and the geometric constraints, Eq. ( 6), to simplify the result.This gives, 2 ds are readily identified as the kinetic and potential energies of the sheet, respectively.Given the boundary conditions on the sheet's edges, Eq. ( 9), the first term in the right hand-side of Eq. (A1) vanishes.Therefore, to complete the derivation, it remains to show that the second term in the right-hand side of Eq. (A1) equals −dE f (t)/dt.Following Ref. [35], the kinetic energy of an incompressible fluid is given by, where δv i (t) are the perimeters of the upper or the lower parts of the chamber, ds is an infinitesimal element on δv i (t) (on the sheet ds = ds), and ni (s, t) are the corresponding local unit normal vectors on δv i (t).Since (∇ϕ i • ni ) x=0,1 = 0 on the sidewalls of the chamber, in accordance with Eq. (4b), and since we have periodic boundary conditions on the upper and lower walls, the right hand side of Eq. (A2) reduces to an integral over the configuration of the sheet.When summed over the two parts of the chamber this gives Substituting Eq. (A4) into Eq.(A1) and integrating once with respect to time completes the derivation.

B. Minimization of the action
In this Appendix, we show that the minimization of the action, S = T 0 Ldt where L is given by Eq. ( 14), yields the complete set of equilibrium equations in the small-amplitude approximation, Eqs. ( 11)- (13).To do so, we minimize the action with respect to the elastic fields, y sh (x, t) and F x (t), and the hydrodynamic fields, ϕ u (x, y, t) and ϕ d (x, y, t), in the standard way.We consider a small perturbation in each of these variables, for example, y sh → y sh + δy sh , and then expand the action to linear order in the perturbation, δy sh (x, t).This procedure gives, after integration by parts, the variation, where v i are the volumes of the chamber above and below the sheet in the small-amplitude approximation, δv i are the perimeters of the upper and lower volumes, ni are the unit normal vectors on δv i , and ds is an infinitesimal line element on δv i (on the sheet ds = ds).
we obtain a closed set of equations for the linearization of the system in this, more general, case and explain the direction we take to obtain the numerical solution.
To linearize Eqs. ( 2)-( 9), we first expand the elastic and the hydrodynamic fields around their base solutions; for example, y sh (s, t) = y sh (s, 0) + ϵe σt ŷsh (s), where ŷsh (s) is a yet-to-be-determined eigenfunction, and ϵ is an arbitrary small parameter.Similarly, we define the eigenfunctions {x sh (s), θ(s), Fx (s), Fy (s)} for the elastic sheet, and { φi (x, y), pi (x, y)} for the fluid.We keep in mind that the fluid starts from rest, and therefore the base solutions for the hydrodynamic fields are equal to zero.Thereafter, we substitute these expansions in the continuity and Bernoulli's equations, Eq. ( 2), and expand them to a linear order in ϵ.This expansion reads, where in the last equation we determine the constant c i (t) such that Eq. ( 5) is satisfied.Similarly, an expansion of the geometric constrains, Eq. ( 6), and the force balance equations on the sheet, Eq. ( 8), gives, To solve this set of equations for given L y , ∆, and λ, we first obtain numerically the base solution for the position of the sheet, i.e., x sh (s, 0) and θ(s, 0).Then, we substitute this solution into the linearized equations and discretize them.The discrete equations are solved using a finite-difference scheme for the elastic sheet and a finite-element scheme for the solution of Eq. (D1) in the bulk of the fluid.
FIG.2.Evolution of the static solution in the small-amplitude approximation, where ∆ = 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 Eq. (21b), while in the symmetric branch it is given by Eqs.(23b) and (23c).The volume difference at the asymmetricto-symmetric transition, v du (0) = v cr du , is labeled by 3○.The pressure difference in the chamber vanishes when the sheet accommodates either the second or the first mode of buckling, labels 1○ and 4 ○.As the volume difference approaches its limiting value v du (0) → (2∆/3) 1/2 , the pressure difference diverges.(b) Evolution of the sheet's profile as the volume difference increases; see the corresponding labeled numbers in panel (a).Note, that despite the relatively large change in the pressure difference in the symmetric branch, the elastic configurations remain almost unchanged.
FIG. 8. (a) The flow fields and (b) the hydrodynamic pressure fields obtained from the linear stability analysis of Eqs.(2)-(9) when ∆ = 0.01, λ = 0.1, and Ly = 2.The normalization of the eigenfunctions is as in figure 7(b).This normalization implies [High, low] = [5.3,26.5] in the flow fields, and [High, low] = [−560, 560] in the pressure fields.In panel (a), arrows represent the streamlines, and colors represent the relative magnitudes of the velocity.

FIG. 11 .
FIG.11.The peak time, σtp, as a function of the system's parameters.(a) The peak time as a function of λ, where ∆ = 0.01 and Ly = 2, for two different values of the initial volume difference.Dashed lines correspond to the numerical integration of Eq. (31), and symbols with corresponding colors represent the numerical solution of Eqs.(2)-(9).For a given v du (0), the time σtp changes by less than 5% over six orders of magnitude of the parameter λ.(b) Comparison between the analytical scaling, Eq. (32), 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 ≃ 2.9 slightly overestimates the numerical prediction.

FIG. 12 .
FIG.12.The pud -v du relation and the maximum backward pressure.In all three panels, open green triangles correspond to the numerical solution of Eqs.(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 (λ = 100), and (b) the fluid-dominated (λ = 0.01) regions.In both panels, the solid black line corresponds to the solution of Eq. (20) with N = 2, ∆ = 0.01, Ly = 2, and v du (0) = 0.01v cr du .The blue points correspond to times t = 0 and t = tp.The dashed gray lines correspond to the static solution, Eqs.(21b) and (23b).The insets show pud (t) as a function of time, where σ is given by Eq. (26).In the fluid-dominated region, the maximum backward pressure (p ud (tp) ≃ −200) is much larger than that in the solid-dominated region (p ud (tp) ≃ −1).In addition, the evolution of the pud -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.

(
d dt E k sh (t) + E p sh (t) = − p u − p d ) ∂x sh ∂t • nd ds,

p 1 0(
i ∇ϕ i • ni ds = − f (t) = i=u,d 1 2λ vi(t) |∇ϕ i |2 dxdy is the energy of the fluid, and in the second equality we used the kinematic boundary condition, Eq. (7), to replace the normal velocity of the fluid with the normal velocity of the sheet.Finally, note that on the sheet the normal vectors are related by, nu (s, t) = −n d (s, t).Using this relation and Eq.(A3), we obtain that,dE f (t) dt = p u − p d ) ∂x sh ∂t • nd ds.(A4)
Log-log plot of the highest growth rate as a function of λ where Ly = 2. Symbols correspond to the linear stability analysis of Eqs. (