2.1. Problem set-up

Figure 1. $(a)$ An elastic sheet (radius $\tilde {R}$ , density $\tilde {\rho }_s$ , bending rigidity $\tilde {B}$ , Poisson’s ratio $\nu$ , thickness $\tilde {e}$ ) immersed in a fluid (ambient density $\tilde {\rho }_a$ , ambient pressure $\tilde {p}_a$ , dynamic viscosity $\tilde {\mu }$ ) and forced periodically at its centre (force $\tilde {F}_a$ , angular frequency $\tilde {\omega }$ , radius $\tilde {\ell }$ ) is placed below a solid substrate with gravity pointing downward. $\boldsymbol{x}_\perp =(x,y)$ represent the horizontal coordinates. The dynamics is characterised by dimensionless numbers defined in table 1: Reynolds number $ \textit{Re}_{{bv}}$ , squeeze number ${Sq}_{{bv}}$ , solid inertia $\mathcal{I}_{{bv}}$ , weight $\mathcal{G}$ and forcing strength $\alpha$ . $(b)$ When $\mathcal{G}$ is small enough, the sheet hovers around an equilibrium position $h_{{eq}}$ (illustrated here for $\alpha =5$ , $\mathcal{G}=2.25$ and a purely viscous dynamics, $\mathcal{I}_{{bv}}={Re}_{{bv}}={Sq}_{{bv}}=0$ ). Above a critical weight, the sheet cannot adhere to the substrate (shown for $\mathcal{G}=2.425$ ). Thin lines represent the gap thickness at the centre of the sheet $h(x=0,t)$ and thick lines the time-averaged $\langle h \rangle (x=0,t)$ .

We consider the system shown in figure 1: an elastic sheet of radius $\tilde {R}$ , thickness $\tilde {e}$ , density $\tilde {\rho }_s$ , Poisson ratio $\nu$ , Young’s modulus $\tilde {E}$ and bending modulus $\tilde {B}= \tilde {E} \tilde {e}^3/12(1-\nu ^2)$ placed near a solid substrate. The surrounding fluid is Newtonian with ambient density $\tilde {\rho }_a$ , ambient pressure $\tilde {p}_a$ and a constant viscosity $\tilde {\mu }$ . The sheet is forced at its centre by a harmonic active load with angular frequency $\tilde {\omega }$ , radius $\tilde {\ell }$ and force magnitude $\tilde {F}_a$ . We aim to establish equilibrium conditions under which the elastic sheet hovers at a finite, stable time-averaged distance from the wall despite the gravitational pull.

An important characteristic scale of the system is the elastohydrodynamic height (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025)

(2.1a) \begin{align} \tilde {H}_{{bv}}={\tilde {R}^2}{\left (\frac {\tilde {\mu } \tilde {\omega }}{\tilde {B}}\right )^{1/3}} \!, \end{align}

for which viscous stresses scaling as $\tilde {\mu } \tilde {\omega } \tilde {R}^2 / \tilde {H}_{{bv}}^2$ and bending stresses scaling as $\tilde {B} \tilde {H}_{{bv}} / \tilde {R}^4$ balance. Associated with $\tilde {H}_{{bv}}$ is the aspect ratio $\varepsilon _{{bv}}=\tilde {H}_{{bv}}/\tilde {R}$ and an elastohydrodynamic force scale

(2.1b) \begin{align} \tilde {F}_{{bv}} = (\tilde {\mu } \tilde {\omega } \tilde {B}^2 )^{1/3}. \end{align}

The scale for elastohydrodynamic stresses is then $\tilde {F}_{{bv}}/\tilde {R}^2$ . In the remainder of the article, we mostly work with dimensionless quantities (written throughout without a tilde, in contrast to dimensional quantities):

(2.2) \begin{align} \begin{split} t = \tilde {t}\tilde {\omega }, \quad \boldsymbol{x}_\perp =\frac {\tilde {\boldsymbol{x}}_\perp }{\tilde {R}}, \quad h(\boldsymbol{x}_\perp ,t)=\frac {\tilde {h}(\tilde {\boldsymbol{x}}_\perp ,\tilde {t})}{\tilde {H}_{{bv}}}, \\ p(\boldsymbol{x}_\perp ,t) = \frac {\tilde {p}(\tilde {\boldsymbol{x}}_\perp ,\tilde {t})}{\tilde {F}_{{bv}}/\tilde {R}^2},\quad \rho (\boldsymbol{x}_\perp ,t)=\frac {\tilde {\rho }(\tilde {\boldsymbol{x}}_\perp ,\tilde {t})}{\tilde {\rho }_a}, \end{split} \end{align}

with $\boldsymbol{x}_\perp =(x,y)$ the horizontal coordinates, $p(\boldsymbol{x}_\perp ,t)$ the fluid pressure relative to the ambient pressure and $\rho (\boldsymbol{x}_\perp ,t)$ the fluid density.

We consider small aspect ratios, $\varepsilon _{{bv}} \ll 1$ , for which the viscous flow in the thin layer separating the sheet and the wall dominates the dynamics. Two other dimensionless numbers describe the flow. Inertial effects are characterised by the film Reynolds number $ \textit{Re}_{{bv}}=\tilde {\rho } \tilde {\omega } \tilde {H}_{{bv}}^2/\tilde {\mu }$ that compares the inertial pressure $\tilde {\rho } \tilde {R}^2 \tilde {\omega }^2$ to the viscous stress $\tilde {F}_{{bv}}/\tilde {R}^2$ . The Reynolds number can also be written as ${\textit{Re}}_{{bv}}=(\tilde {H}_{{bv}}/\tilde {\delta })^2$ , with $\tilde {\delta }=(\tilde {\mu }/\tilde {\rho } \tilde {\omega })^{1/2}$ the viscous penetration length, the length scale for diffusion of vorticity. It may also be interpreted as a Womersley number, and we refer to Appendix A.1 for a detailed discussion of the scaling of inertial effects. We only emphasise here that ${\textit{Re}}_{{bv}}$ is constructed based on the vertical length and velocity scales, $\tilde {H}_{{bv}}$ and $\tilde {\omega } \tilde {H}_{{bv}}$ , as appropriate in the lubrication limit $\varepsilon _{{bv}}\ll 1$ (Batchelor Reference Batchelor1967). Compressible effects are characterised by the squeeze number ${Sq}_{{bv}} = \tilde {F}_{{bv}}/\tilde {p}_a \tilde {R}^2$ that compares the viscous stress to the ambient pressure $\tilde {p}_a$ (Taylor & Saffman Reference Taylor and Saffman1957). With (2.1), the inertial and compressible effects in the fluid are quantified respectively by

(2.3) \begin{align} \textit{Re}_{{bv}}=\frac {\tilde {\rho }_a \tilde {\omega }^2 \tilde {R}^4 }{\big (\tilde {\mu } \tilde {\omega } \tilde {B}^2\big )^{1/3}}, \qquad {Sq}_{{bv}}=\frac {\big (\tilde {\mu } \tilde {\omega } \tilde {B}^2\big )^{1/3}}{\tilde {p}_a \tilde {R}^2}. \end{align}

2.2. Inertial lubrication

Let us find a depth-integrated description of the flow in the thin gap between the sheet and the wall ( $\varepsilon _{{bv}}\ll 1$ ) which includes the first-order effects of inertia at $\mathcal{O}({{Re}_{{bv}}})$ and compressibility at $\mathcal{O}({{Sq}_{{bv}}})$ . Mass conservation yields

(2.4a) \begin{align} {\frac {\partial \left (\rho h\right )}{\partial t}} + \boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\left (\rho \boldsymbol{q}\right ) &= 0, \\[-12pt] \nonumber \end{align}

(2.4b) \begin{align} \rho &= 1 + {Sq}_{{bv}} p, \end{align}

with $\boldsymbol{q}=\int _0^h \boldsymbol{v}_\perp\, {\mathrm{d}}z=\tilde {\boldsymbol{q}} / \tilde {\omega } \tilde {R}$ the horizontal volumetric fluid flux, $\boldsymbol{v}_\perp$ the horizontal fluid velocity profile and $\boldsymbol{\nabla }_\perp =(\partial /\partial x, \partial /\partial y)$ the horizontal gradient operator. Equation (2.4b ) is the dimensionless ideal gas law assuming isothermal conditions, $\tilde {\rho }/\tilde {\rho }_a=1+\tilde {p}/\tilde {p}_a$ , an assumption we discuss later in § 5. The volumetric flux $\boldsymbol{q}$ appearing in (2.4) is found from the Navier–Stokes equations. Rojas et al. (Reference Rojas, Argentina, Cerda and Tirapegui2010) describe a procedure to consider the first-order inertial corrections to lubrication theory for free surface flows. We adapt their derivation to a fluid layer bounded by two solid walls (Appendix A), which yields the depth-integrated horizontal momentum balance, to $\mathcal{O}(\varepsilon _{{bv}}^2,{Re}_{{bv}},\varepsilon _{{bv}}^2{Re}_{{bv}},{Re}_{{bv}}{Sq}_{{bv}},{Sq}_{{bv}})$ :

(2.5) \begin{align} 12 {\boldsymbol{q}}+ { h^3}{\boldsymbol{\nabla }}_\perp p + {Re}_{{bv}} h^3 \left [ \frac 65 {\frac {\partial }{\partial t}}\left (\frac {{\boldsymbol{q}}}{ h}\right ) +\frac {54}{35} \frac {{\boldsymbol{q}} }{ h} \boldsymbol{\cdot }\boldsymbol{\nabla }_\perp \left (\frac {{\boldsymbol{q}}}{ h}\right ) - \frac 6{35} \frac {{\boldsymbol{q}}}{ h^2}{\frac {\partial h}{\partial t}} \right ]&=0. \end{align}

The first two terms in (2.5) reduce to the Reynolds equation from inertialess lubrication theory, valid for vanishing Reynolds number (Batchelor Reference Batchelor1967). The term proportional to $ \textit{Re}_{{bv}}$ corresponds to the first-order correction due to both the unsteady and convective inertia of the fluid. It corrects the parabolic Poiseuille velocity profile of lubrication theory and predicts a sextic profile, discussed in § 4. Rojas et al. (Reference Rojas, Argentina, Cerda and Tirapegui2010) found good agreement between experiments and this extended lubrication theory for Reynolds numbers of order one. Previous theoretical works from Ishizawa (Reference Ishizawa1966) and Jones & Wilson (Reference Jones and Wilson1975), albeit limited to non-deformable boundaries, suggest that this correction may even be valid as long as the Reynolds number is less than $100$ . Equation (2.5) recovers the inertial corrections derived in these studies and generalises them to an arbitrary height $h(\boldsymbol{x}_\perp , t)$ , suggesting that (2.5) may be valid for finite, relatively large values of $ \textit{Re}_{{bv}}$ , say $\mathcal{O}({10})$ .

2.3. Elastic deformations

To close the system formed by (2.4) and (2.5), the fluid pressure $p$ must be linked to elastic stresses by considering a vertical momentum balance of the elastic sheet. We adopt the Kirchhoff–Love model (Timoshenko & Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1959; Landau & Lifshitz Reference Landau and Lifshitz1986) to describe the deformation of the sheet. We therefore neglect any in-plane stretching, which is justified for a sheet undergoing cylindrical bending, i.e. deforming in a single direction and preserving a zero Gaussian curvature (in which case, the model reduces to the Euler–Bernoulli beam) and for two-dimensional deformations whose amplitude remains small compared with the thickness of the sheet. We also neglect the tension induced by the shear stress from the flow (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025). The normal force per unit area from the aerodynamic interactions in the thin gap is $p + \mathcal{O}({\varepsilon _{{bv}}^2})$ , so that the force balance in the thin-film limit reads

(2.6) \begin{align} \begin{split} \mathcal{I}_{{bv}} {\frac {\partial ^2h}{\partial t^2}} &= p + \boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\left ( \boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\boldsymbol{M}\right ) + f_a(\boldsymbol{x}_\perp ,t)+\mathcal{G} + f_{w}(h), \\ \boldsymbol{M} &=-\left [\left (1-\nu \right ) \boldsymbol \kappa + \nu \operatorname {tr}(\boldsymbol \kappa ) \mathbf I \right ]. \end{split} \end{align}

The dimensionless number $\mathcal{I}_{{bv}}=\tilde {\rho }_s \tilde {e} {\tilde {H}}_{{bv}} \tilde {\omega }^2 \tilde {R}^2/ \tilde {F}_{{bv}}$ compares solid inertia to the elastohydrodynamic stress scale. The right-hand side of (2.6) corresponds respectively to the stress from the fluid, the bending stress (with $\boldsymbol \kappa$ the Hessian of $h$ – the curvatures of the sheet), the periodic active stress $f_a$ , the sheet’s areal weight with $\mathcal{G}=\tilde {\rho }_s \tilde {e} \tilde {R}^2 \tilde {g}/\tilde {F}_{{bv}}\gt 0$ and a term preventing collision with the wall. We also note that in the case of a constant bending rigidity considered here, the bending stresses simplify to $\boldsymbol{\nabla }_\perp \boldsymbol{\cdot }( \boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\boldsymbol{M} )=-\boldsymbol{\nabla} ^4_\perp h$ . We emphasise that the sheet is assumed to have a uniform weight: $\mathcal{G}$ is constant. This differs from the experiments of Colasante (Reference Colasante2016) and Weston-Dawkes et al. (Reference Weston-Dawkes, Adibnazari, Hu, Everman, Gravish and Tolley2021), where it is heavier at its centre and may be submitted to an additional pulling force. We have made this choice to isolate the role of elastohydrodynamics without introducing additional complications. Nonetheless, this model is versatile, and non-uniform loading could be studied by prescribing a spatially varying mass distribution $\mathcal{G}(\boldsymbol{x}_\perp )$ , leading also to a spatially varying inertia $\mathcal{I}_{{bv}}(\boldsymbol{x}_\perp )$ , and/or by adding a localised force to (2.6).

The active forcing behind $f_a$ is harmonic in time and is distributed at the centre of the sheet with a dimensionless radius $\ell =\tilde {l}/\tilde {R}$ (figure 1). Letting $\mathbb H$ the Heaviside function ( $\mathbb H(x)=1$ if $x\geqslant 0$ , $=0$ otherwise) and $\alpha =\tilde {F}_a/\tilde {F}_{{bv}}$ the dimensionless strength of the forcing, we have

(2.7) \begin{align} f_a(\boldsymbol{x}_\perp ,t)=\alpha \cos (t) \frac {1-\mathbb H(\lvert \boldsymbol{x}_\perp \rvert -\ell )}{\ell }. \end{align}

With (2.1), the dimensionless numbers in (2.6) read

(2.8) \begin{align} \begin{split} \mathcal{I}_{{bv}}&=\frac {\tilde {\rho }_s \tilde {e} \tilde {\omega }^2 \tilde {R}^4}{\tilde {B}} , \quad \mathcal{G} = \frac {\tilde {F}_g}{\tilde {F}_{{bv}}} = \frac {\tilde {\rho }_s \tilde {e} \tilde {g} \tilde {R}^2}{\big (\tilde {\mu } \tilde {\omega } \tilde {B}^{2}\big )^{1/3}}, \quad \alpha =\frac {\tilde {F}_a}{\tilde {F}_{{bv}}}=\frac {\tilde {F}_a}{\big (\tilde {\mu } \tilde {\omega } \tilde {B}^{2}\big )^{1/3}}. \end{split} \end{align}

Finally, as discussed later in § 3.2, we observe that when the forcing and weight both exceed a critical value, the edges of the sheet may contact the wall. To handle this numerically, we introduce a local repulsive force $f_w(h)= (A/h )^n$ , $n \gt 1$ , modelling an elastic collision with a contact that effectively occurs at the dimensionless height $A\ll 1$ . Indeed, $f_w$ is conservative, since it derives from the potential $A^n h^{1-n}/(n-1)$ , and is only significant for $h\lesssim A$ . The introduction of a short-ranged repulsion is similar to the penalty method commonly employed in contact mechanics (Wriggers Reference Wriggers2006) and in fluid–structure interaction problems (e.g. Glowinski et al. Reference Glowinski, Pan, Hesla, Joseph and Periaux2001).

2.4. Time scales

To discuss the different physical effects that contribute to the lubrication dynamics, it will prove useful to define characteristic time scales associated with each mechanism. From the viscous lubrication equations (2.4) and (2.5) with ${\textit{Re}}_{{bv}}={{{Sq}}}_{{bv}}=0$ , $12\tilde {\mu } \partial \tilde {h}/\partial \tilde {t} = \tilde {\boldsymbol{\nabla }}_\perp \boldsymbol{\cdot }(\tilde {h}^3 \tilde {\boldsymbol{\nabla }}_\perp \tilde {p})$ , a time $\tilde {T} = \tilde {\mu } \tilde {R}^2/ (\tilde {H}^2 \tilde {P})$ can be defined with $\tilde {H}$ a height and $\tilde {P}$ the characteristic pressure associated with the mechanism of interest. Considering each effect separately, this defines the following time scales:

(2.9a) \begin{align} \mathrm{gravitational}:\,\tilde {T}_g(\tilde {h})=\frac {\tilde {\mu } \tilde {R}^2}{\tilde {\rho }_s \tilde {g} \tilde {e} \tilde {h}^2}; \end{align}

active, $\tilde {T}_a(\tilde {h}) = {\tilde {\mu } \tilde {R}^4}/{(\tilde {F}_a \tilde {h}^2)}$ ; elastohydrodynamic, $\tilde {T}_{{bv}}(\tilde {h}) = {\tilde {\mu } \tilde {R}^6}/{(\tilde {B} \tilde {h}^3)}$ ; compressible, $\tilde {T}_c(\tilde {h})={\tilde {\mu } \tilde {R}^2}/({\tilde {p}_a \tilde {h}^2})$ ; and inertial, $\tilde {T}_i(\tilde {h})={\tilde {\rho }_a \tilde {h}^2}/{\tilde {\mu }}$ . The latter is characteristic of momentum diffusion and is not defined by the same lubrication scaling. While gravity is always present, the effects of elasticity, compressibility and inertia are only relevant when the sheet is forced into motion by the active force. In these cases, we expect the time-averaged and long-term dynamic of the sheet to only depend on even powers of the active force $\tilde {F}_a$ , since a change of sign of the force magnitude $\tilde {F}_a$ is equivalent to a phase shift. We therefore define visco-active time scales as the square of the active forcing time scale $\tilde {T}_a$ divided by the response time scale of the mechanism of interest:

(2.9b) \begin{align} \begin{split} \text{elastohydrodynamic,}\;\tilde {T}_{{a},\textit{bv}}(\tilde {h}) &= \frac {\tilde {T}_a^2(\tilde {h})}{\tilde {T}_{{bv}}(\tilde {h})}=\frac {\tilde {\mu } \tilde {R}^2 \tilde {B}}{\tilde {F}_a^2 \tilde {h}}; \\[-3pt] \text{inertial,}\; \tilde {T}_{{a,i}}(\tilde {h})=\frac {\tilde {T}_a^2(\tilde {h})}{\tilde {T}_{i}(\tilde {h})}=\frac {\tilde {\mu }^3 \tilde {R}^8}{\tilde {F}_a^2 \tilde {\rho }_a \tilde {h}^6}; \quad & \text{compressible,}\;\tilde {T}_{{a,c}}(\tilde {h}) = \frac {\tilde {T}_a^2(\tilde {h})}{\tilde {T}_{c}(\tilde {h})}=\frac {\tilde {\mu } \tilde {R}^6 \tilde {p}_a }{\tilde {F}_a^2 \tilde {h}^2}. \end{split} \end{align}

2.5. Boundary conditions

The edges of the elastic sheet are free of bending moment, twisting moment and shear force. The appropriate boundary conditions are (Naghdi Reference Naghdi1973, p. 586)

(2.10) \begin{align} \boldsymbol{M} \boldsymbol{e_r} \boldsymbol{\cdot }\boldsymbol{e_r}= 0, \quad \left (\boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\boldsymbol{M} \right )\boldsymbol{\cdot }\boldsymbol{e_r} + \boldsymbol{\nabla }_\perp \left ( \boldsymbol{M} \boldsymbol{e_r} \boldsymbol{\cdot }\boldsymbol{e_\theta }\right )\boldsymbol{\cdot }\boldsymbol{e_\theta }=0, \end{align}

where $\boldsymbol{e_r}$ is the outward unit normal at the edge of the sheet and $\boldsymbol{e_\theta }$ is the unit tangent. As the third and final boundary condition, we set the pressure at the edges of the sheet as

(2.11) \begin{align} p= \begin{cases} 0\,&\text{if} \quad \boldsymbol{q} \boldsymbol{\cdot }\boldsymbol{e}_r \gt 0\ \text{(outflow)}, \\ - \dfrac {k}{2} {Re}_{{bv}}\left (\dfrac {\boldsymbol{q} \boldsymbol{\cdot }\boldsymbol{e}_r }{h}\right )^2\,&\text{if} \quad \boldsymbol{q} \boldsymbol{\cdot }\boldsymbol{e}_r \lt 0 \ \text{(inflow)} ,\end{cases} \end{align}

with $k=1/2$ , a loss coefficient. This boundary condition models the pressure loss when fluid enters the thin gap; it is justified and discussed in more detail in Appendix B. In the inertialess case ${\textit{Re}}_{{bv}}=0$ , (2.11) simplifies to the classical condition $p=0$ , which imposes the ambient pressure at the edges.

2.6. Numerical model

We have formulated the governing equations (2.4a ), (2.5), (2.6), and boundary conditions (2.10) and (2.11) for a two-dimensional (2-D) sheet for completeness. For purely viscous flows, we have shown (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025) that there are no qualitative differences between the one-dimensional (1-D) and 2-D axisymmetric situations. Accordingly, we restrict the present study to 1-D for simplicity: $\boldsymbol{x}_\perp \rightarrow x$ and $\partial /\partial y=0$ . We solve the governing equations in conservative form

(2.12a) \begin{align} {\frac {\partial \boldsymbol{S}}{\partial t}} + {\frac {\partial \boldsymbol{F}}{\partial x}} = \boldsymbol{Q}, \end{align}

with

(2.12b) \begin{align} \begin{array}{lllccccccr} \boldsymbol{S}(\boldsymbol{U}) = & \Big[\rho h, & \mathcal{I}_{{bv}} v, & 0, & \dfrac {6{Re}_{{bv}}}{5}u, & h \Big], \\[9pt] \boldsymbol{F}(\boldsymbol{U}) = & \Big[ 0, & -\dfrac {\partial m}{\partial x}, & -\dfrac {\partial h}{\partial x}, & \dfrac {27}{35}{Re}_{{bv}} u^2 + p, & hu \Big], \\[9pt] \boldsymbol{Q}(\boldsymbol{U}) = & \Big[ \rho v, & p + f_a(x,t)+\mathcal{G} + f_w(h), & m, & -\dfrac {12u}{h^2} + \dfrac {6}{35}{Re}_{{bv}}\dfrac {uv}{h}, & 0 \Big], \end{array} \end{align}

and the five primary unknowns $\boldsymbol{U}= [v, m, h, u, p]$ , representing the vertical sheet velocity $v={\partial h}/{\partial t}$ , bending moment $m=- \partial ^2 h/\partial x^2$ , height $h$ , average horizontal fluid velocity $u=q/h$ and fluid pressure $p$ . We note that $\rho =1+{Sq}_{{bv}} p$ .

We consider a domain $x \in [0,1]$ with symmetric boundary conditions at $x=0$ and the appropriate boundary conditions at $x=1$ :

(2.12c) \begin{align} m = 0, \quad {\frac {\partial m}{\partial x}} = 0, \quad p&= \begin{cases} 0\,&\text{if}\,u \geqslant 0, \\ - k{Re}_{{bv}} u^2/2\,&\text{if}\,u \lt 0. \end{cases} \end{align}

We solve (2.12) using the DuMu $^x$ library (Koch, Weishaupt & Gläser Reference Koch, Weishaupt and Gläser2021). We discretise space using a staggered finite volume scheme, where pressure unknowns are located at cell centres and other unknowns are located at vertices. This set-up avoids checkerboard oscillations in the fluid pressure. The flux term in $u^2$ is treated with a first-order upwind scheme. The equations are advanced in time using a diagonally implicit third-order Runge–Kutta scheme (Alexander Reference Alexander1977, Thm. 5) and the nonlinear system at each Runge–Kutta stage is solved with Newton’s method. Lower-order methods either yielded unsatisfactory accuracy or required excessively small time step sizes. We used time step sizes $10^{-4} \leqslant \Delta t \leqslant 0.2$ and spatial step sizes $0.002\leqslant \Delta x \leqslant 0.02$ . The numerical simulations have been run to a time-averaged steady state – which could take from $t=\mathcal{O}({10})$ up to $t=\mathcal{O}({10^6})$ depending on the parameters – or until the height diverged. We systematically ensured that any divergence of the numerical solution was independent of the numerical parameters and therefore corresponded to adhesion failure. For initial conditions, we considered a flat sheet, $\boldsymbol{U}(x,t=0)=[0, 0, h(x,t=0), 0, 0]$ , with $h(x,t=0)$ a constant. For large values of $\alpha$ and $\mathcal{G}$ , this initialisation sometimes leads to divergence, even though a time-averaged steady state exists. In these cases, we initialised the simulation using the steady-state solution from a run with the same $\alpha$ but smaller $\mathcal{G}$ (numerical continuation). The repulsion force was either turned off ( $A=0$ ) or, when needed, chosen as $f_w(h)= (A/h )^n$ with $A=10^{-5}$ and $n=5$ . We have verified that this choice does not significantly affect the results as long as $A$ is small and $n$ is large.

2.7. Choice of dimensionless parameters

Equations (2.4b ), (2.5) and (2.6) depend on five dimensionless numbers defined in (2.1) and (2.8), with $\varepsilon _{{bv}}$ not appearing in the governing equations and $\ell$ defined in (2.7). They are summarised in table 1. Since dimensional quantities such as the sheet’s bending rigidity $\tilde {B}$ or the excitation frequency $\tilde {\omega }$ enter multiple dimensionless groups, independently varying them in experiments is not feasible. Nevertheless, numerical simulations enable us to disentangle the respective roles of the dimensionless numbers in the dynamics and to clarify the underlying physical mechanisms they influence. We consider a uniform sheet with a uniformly distributed weight: $\mathcal{I}_{{bv}}$ and $\mathcal{G}$ are constants. We also consider the limit where the forcing is localised at a single point, $\ell \rightarrow 0$ ; in practice, we set $\ell =0.05$ in our numerical simulations. The effect of a finite $\ell$ and a sheet locally rigid at its centre has already been discussed in prior work (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025).

Table 1. Characteristic scales and dimensionless parameters.

To guide our study, we consider the experiments of Weston-Dawkes et al. (Reference Weston-Dawkes, Adibnazari, Hu, Everman, Gravish and Tolley2021), who report a time-averaged equilibrium height (figure 1 $b$ ) $\tilde {h}_{{eq}} \approx 600\,\unicode{x03BC}\mathrm{m}$ for a sheet of thickness $\tilde {e} \approx 300\,\unicode{x03BC}\mathrm{m}$ , Young’s modulus $\tilde {E} \approx 3\,\mathrm{GPa}$ , radius $\tilde {R} \simeq 10\,\mathrm{cm}$ , vibrating at frequency $\tilde {\omega } = 2\pi \times 200\,\mathrm{Hz}$ and supporting a weight $\tilde {W} \approx 5\,\mathrm{N}$ . The density ratio between the sheet and air is $\tilde {\rho }_s/\tilde {\rho }_a\approx 10^3$ . The vibrations are generated by an eccentric rotating mass motor, with an estimated mass $\tilde {m} \approx 0.6\,\mathrm{g}$ and gyration radius $\tilde {r} \approx 1\,\mathrm{mm}$ , yielding a driving force $\tilde {F}_a=\tilde {m} \tilde {r} \tilde {\omega }^2$ . This gives a dimensionless forcing strength $\alpha =\tilde {F}_a/\tilde {F}_{{bv}}\approx 90$ . While this is overestimated due to the rigid central support (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025), it nevertheless suggests that the regime $\alpha =\mathcal{O}({10})$ is relevant experimentally.

From the above-mentioned parameters, the height scale is $\tilde {H}_{{bv}} \approx 14\,\mathrm{mm}$ , while the dimensionless equilibrium height is only a small fraction of this value: $h_{{eq}}=\tilde {h}_{{eq}}/\tilde {H}_{{bv}} \approx 0.04$ . Our theoretical and numerical analysis will recover this observation. However, this shows that the dimensionless numbers based on $\tilde {H}_{{bv}}$ , such as $ \textit{Re}_{{bv}}$ , ${Sq}_{{bv}}$ and $\mathcal{I}_{{bv}}$ , which allow for a compact theoretical description of the system dynamics, are not accurate indicators of the relative effect of fluid inertia, fluid compressibility or solid inertia when the system is at equilibrium height. Hence, we additionally define dimensionless numbers using $\tilde {h}_{{eq}}$ as the vertical scale:

(2.13) \begin{align} {Re}_{{eq}} = \frac {\tilde {\rho }_a \tilde {\omega } \tilde {h}_{{eq}}^2}{\tilde {\mu }}=h_{{eq}}^2 {Re}_{{bv}}, \quad {{{Sq}}}_{{eq}}=\frac {\tilde {\mu } \tilde {\omega } \tilde {R}^2}{\tilde {h}_{{eq}}^2 \tilde {p}_a} = h_{{eq}}^{-2} {Sq}_{{bv}}, \quad \mathcal{I}_{{eq}}=\frac {\tilde {\rho }_s\tilde {e} \tilde {h}_{{eq}}^3 \tilde {\omega }}{\tilde {\mu } \tilde {R}^2}=h_{{eq}}^{3}\mathcal{I}_{{bv}}. \end{align}

Unlike the original dimensionless groups, these quantities cannot be computed a priori since $\tilde {h}_{{eq}}$ is selected by the system and is initially unknown. The experimental values yield $\alpha =\mathcal{O}({10})$ and we systematically vary $\alpha$ in § 3. The Reynolds number is $ \textit{Re}_{{eq}} \approx 30$ , indicating that fluid inertia may play a significant role, studied in § 4. We also note that in some of the experiments of Weston-Dawkes et al. (Reference Weston-Dawkes, Adibnazari, Hu, Everman, Gravish and Tolley2021), when the sheet is steadily pulled instead of equilibrating under a constant load, the central height can reach up to 3 mm, corresponding to Reynolds numbers of several hundred. This further highlights the importance of understanding inertial effects from the surrounding fluid. The parameter controlling solid inertia, however, is smaller, $\mathcal{I}_{{eq}} \approx 0.6$ , suggesting a weaker yet potentially non-negligible effect of solid inertia. In this study, we focus on the influence of fluid effects and leave a systematic investigation of solid inertia for future work. A related investigation into the coupling between solid inertia and elastohydrodynamics can be found in Ramanarayanan & Sánchez (Reference Ramanarayanan and Sánchez2024). Finally, although ${{{Sq}}}_{{eq}}\approx 0.006$ , we show in § 5 that compressibility could still noticeably affect the dynamics even for small squeeze numbers.

3.1. Weak active forcing ( $\alpha \lesssim 1$ )

Figure 2. Asymptotic results for $\alpha \lesssim 1$ , $\mathcal{I}_{{bv}}={Re}_{{bv}}={Sq}_{{bv}}=0$ , adapted from Poulain et al. (Reference Poulain, Koch, Mahadevan and Carlson2025). $(a)$ Schematic illustration of the link between the active force direction and the sheet’s convexity. $(b)$ Equilibrium height $h_{{eq}}$ as a function of the the rescaled dimensionless weight $\mathcal{G}/\alpha ^2$ . Symbols are results from numerical simulations, the lines are the prediction of (3.2) obtained by numerical continuation (with a cutoff $N=5$ ). For $\mathcal{G}/\alpha ^2\gt \mathcal{G}_{\textit{max}}/\alpha ^2\simeq 0.137$ , no equilibrium is possible and the sheet always detaches from the substrate (greyed area).

We first consider the regime for which inertia and compressibility are negligible, $\mathcal{I}_{{bv}}={Re}_{{bv}}={Sq}_{{bv}}=0$ , so that the dynamics is solely governed by viscous elastohydrodynamic interactions. Then, the governing equations (2.4) and (2.5) simplify to the Reynolds equation

(3.1) \begin{align} 12{\frac {\partial h}{\partial t}}-\boldsymbol{\nabla }_\perp \boldsymbol{\cdot }\big (h^3 \boldsymbol{\nabla }_\perp p\big )=0. \end{align}

We have already studied this regime (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025) in the limit $\alpha \lesssim 1$ . More precisely, we used an asymptotic expansion to $\mathcal{O}({\alpha ^2})$ and found the results valid up to $\alpha \simeq 1$ . In short, when periodic vibrations drive an elastic sheet at its centre, it tends to adhere to a nearby surface due to the coupling between its elastic deformation and the lubrication flow in the intervening gap. As the sheet is pushed towards the surface, it adopts a convex shape that favours fluid outflow; while when it is pulled away from the surface, it adopts a concave shape that resists inflow (figure 2 $a$ ). This asymmetry in the flow response over a period of oscillation results in a net outflow, leading to a time-averaged attraction towards the surface. This symmetry breaking can be traced back to the non-time-reversible pressure distribution in a rigid squeeze film (as shown later in figure 7 $b$ ). When the soft sheet deforms under this pressure, its kinematics inherit the irreversibility, so that the dynamics is no longer constrained by the scallop theorem (Purcell Reference Purcell1977; Lauga Reference Lauga2011). The resulting rectified flow can then counteract the sheet’s weight and give rise to an equilibrium hovering height. We have used these insights to study the system analytically and present the conclusions of our analysis in the following. We refer the reader to Appendix C for further details on the underlying assumptions and to Poulain et al. (Reference Poulain, Koch, Mahadevan and Carlson2025) for the complete derivation. In short, the height $h(x,t)$ is decomposed into spatial eigenmodes (we note that similar modal expansion have recently been used in related problems (Papanicolaou & Christov Reference Papanicolaou and Christov2025)), and an asymptotic analysis allows to find an evolution equation for time-averaged height $\langle h_0 \rangle (t) = \int _t^{t+2\pi }h(x=0,t')\,{\mathrm{d}}t'$ at $\mathcal{O}({\alpha ^2})$ :

(3.2a) \begin{align} \begin{split} \frac {1}{\alpha ^2 \langle h_0\rangle ^2} \frac {{\mathrm{d}}\langle h_0\rangle }{{\mathrm{d}}t}&= \frac {1}{4}\frac {\mathcal{G}}{\alpha ^2} \langle h_0\rangle - d_0 + \sum _{i,j=1}^{\infty } d_{\textit{ij}} g_{\textit{ij}}\left (\langle h_0 \rangle \right ),\\\quad g_{\textit{ij}}(h)&= \frac {1+\left (\dfrac {h}{\sqrt {e_ie_{\kern-1pt j}}}\right )^6}{\left (1+\left (\dfrac {h}{e_i}\right )^6\right )\left (1+\left (\dfrac {h}{e_{\kern-1.5pt j}}\right )^6\right )}, \quad e_i = \frac {0.242}{i^2}. \end{split} \end{align}

In dimensional units, this can be expressed compactly using the gravitational time scale $\tilde {T}_g$ and the active elastohydrodynamic time scale $\tilde {T}_{{a,bv}}$ introduced in (2.9):

(3.2b) \begin{align} \frac {{\mathrm{d}}\langle \tilde {h}_0\rangle }{{\mathrm{d}}\tilde {t}}=\frac {\langle \tilde {h}_0 \rangle }{4\tilde {T}_g\left (\langle \tilde {h}_0 \rangle \right )} + \frac {\langle \tilde {h}_0 \rangle }{\tilde {T}_{{a,bv}}\left (\langle \tilde {h}_0 \rangle \right )}\left [ -d_0 + \sum _{i,j=1}^{\infty } d_{\textit{ij}} g_{\textit{ij}}\left (\frac {\langle \tilde {h}_0 \rangle }{\tilde {H}_{{bv}}}\right )\right ]\!. \end{align}

The first two terms on the right-hand side of (3.2) correspond respectively to the effect of gravity, with $\mathcal{G}/\alpha ^2 = \mathcal{O}({1})$ , and to the first-order effect of elastohydrodynamics, with $d_0=0.0122$ . The third term corresponds to the effects of the eigenmodes $\zeta _i$ (shown in figure 11) on the dynamics, with the $e_i=0.242/i^2$ characteristic heights below which the $i$ th spatial mode $\zeta _i$ is excited. This means that we predict higher-order excitation modes as the sheet approaches the wall. The $d_{\textit{ij}}$ are numerical coefficients of the symmetric matrix $d_{\textit{ij}}$ characterising the strength of the effect of mode $i$ ( $i=j$ ) or of the coupling between modes $i$ and $j$ on the dynamics; their values are given in Appendix C.

Figure 3. $(a)$ Illustration of the decomposition (3.3) of the sheet’s shape into a rigid-body translation $\bar h(t)-h_{{eq}}$ , a static shape $\langle h \rangle (x)$ (independent of time at the time-averaged steady state) and a time-periodic deformation $h_d(x,t)$ . $(b{,}c{,}d)$ Time-averaged shape $\langle h \rangle$ and the periodic deformation $h_d$ for $G=0.02,\,0.04,\,0.06$ , respectively, and $\alpha =1$ , $\mathcal{I}_{{bv}}={{Re}}_{{bv}}={{{Sq}}}_{{bv}}=0$ . These are obtained from numerical simulations at the time-averaged steady state. $h_d$ is shown at various times of one vibration cycle, with scale bars showing the amplitude of the deformations. As $G$ and correspondingly $h_{{eq}}$ decrease, higher-order vibration modes are excited.

The bifurcation diagram of (3.2) is shown in figure 2 $(b)$ . If the rescaled weight is too large, $\mathcal{G}\gt \mathcal{G}_{\textit{max}}=0.137 \alpha ^2$ , there is no equilibrium and the sheet fails to adhere: gravity overcomes the adhesive elastohydrodynamic effect and $\langle h_0 \rangle \rightarrow \infty$ . For $0\lt \mathcal{G}\lt \mathcal{G}_{\textit{max}}$ , the sheet finds an equilibrium near the wall, $\langle h_0 \rangle \rightarrow h_{{eq}}$ . We show in figure 2 $(b)$ that numerical simulations of the governing equations (2.6) and (3.1) agree remarkably well with the bifurcation diagram of the reduced model (3.2) for $\alpha \lesssim 1$ .

To study the sheet’s deformations numerically, we consider the time-averaged equilibrium and let

(3.3) \begin{align} h(x,t)=h_d(x,t)+\langle h\rangle (x,t) + \bar h(t)-h_{{eq}}, \end{align}

where $\langle h \rangle (x,t)=(1/2\pi )\int _t^{t+2\pi } h(x,t')\,{\mathrm{d}}t'$ is the time average (independent of time once a time-averaged steady state is reached), $\bar h (t)=(1/2) \int _{-1}^{1} h(x',t)\,{\mathrm{d}}x'$ is the spatial average, and $h_{{eq}}=\langle \bar h \rangle$ the time- and space-averaged height. The function $h_d(x,t)$ then represents the periodic deformations to the static shape $\langle h \rangle (x)$ , with $\langle h_d \rangle = \bar h_d = 0$ . The decomposition (3.3) is illustrated in figure 3, where we also show that higher-order modes of deformation are indeed excited as $\mathcal{G}/\alpha ^2$ and $h_{{eq}}$ decrease.

Finally, our analysis shows that for $\alpha \lesssim 1$ , the maximum supported weight $\mathcal{G}_{\textit{max}}$ scales as $\alpha ^2$ . Coming back to dimensional quantities gives $\tilde {W}_{\textit{max}} \sim \tilde {F}_a^2/\tilde {F}_{{bv}}$ . This predicts that for a given forcing amplitude $\tilde {F}_a$ , an increasingly soft sheet would sustain arbitrarily large weights $\tilde {W}_{\textit{max}}$ . However, as the sheet becomes softer ( $\tilde {F}_{{bv}}$ decreases), $\alpha =\tilde {F}_a/\tilde {F}_{{bv}}$ increases and cannot be considered small anymore. We study this second regime next.

3.2. Strong active forcing ( $\alpha \gg 1$ )

3.2.1. Contactless adhesion – maximum supported weight

Figure 4. Varying active forcing $\alpha$ for an inertialess and incompressible system, $\mathcal{I}_{{bv}}={{Re}}_{{bv}}={{{Sq}}}_{{bv}}=0$ . $(a{,}b)$ Equilibrium height as a function of the dimensionless weight. Open symbols represent cases where contact occurs at the edges of the sheet. $(c)$ Regime map showing the three different possibilities (adhesion with or without edge contact, and adhesion failure) as a function of the dimensionless weight and forcing. The solid line represents the prediction $\mathcal{G}_{\textit{max}}=0.137\alpha ^2$ derived for $\alpha \ll 1$ , and the dashed line is the interpolation (3.4). The inset is a zoom near the origin on a logarithmic scale. (d) Regime maps from (3.4) as $\tilde {F}_a$ varies for fixed $\tilde {F}_{{bv}}$ (left), and $\tilde {F}_{{bv}}$ varies for fixed $\tilde {F}_a$ (right).

We show in figure 4 bifurcation diagrams obtained by numerically solving the governing equations (2.6) and (3.1) for $0\lt \alpha \leq 100$ , $\mathcal{I}_{{bv}}={{Re}}_{{bv}}={{{Sq}}}_{{bv}}=0$ . The corresponding regime map of accessible dimensionless weight reveals three distinct regions: (i) a regime where adhesion is not possible and the system fails; (ii) a regime of contactless adhesion; and (iii) a regime where we predict adhesion, but with the sheet’s edges periodically touching the substrate. The latter is discussed in § 3.2.4; we first focus on contactless adhesion. For a given $\alpha$ , we define $\mathcal{G}_{\textit{max}}$ as the threshold weight above which contactless adhesion cannot take place.

We see in figure 4 $(c)$ that the relation $\mathcal{G}_{\textit{max}}=0.137 \alpha ^2$ is well verified for $\alpha \lesssim 1$ , but that $\mathcal{G}_{\textit{max}}$ reaches an approximately constant value of 7.8 for $\alpha \gtrsim 20$ . These two asymptotic behaviours are captured by the following interpolation, which also describes well the data over the whole range $0\lt \alpha \leqslant 100$ :

(3.4) \begin{align} \mathcal{G}_{\textit{max}} \simeq 0.137 \frac {\alpha ^2}{1+0.0176\alpha ^2}, \quad \tilde {W}_{\textit{max}}\simeq 0.137 \dfrac {\tilde {F}_a^2/\tilde {F}_{{bv}}}{1+0.0176 \tilde {F}_a^2/\tilde {F}_{{bv}}^2}. \end{align}

In dimensional quantities, we find $\tilde {W}_{\textit{max}}=0.137 \tilde {F}_a^2/\tilde {F}_{{bv}}$ for $\tilde {F}_a \lesssim \tilde {F}_{{bv}}$ (weak forcing), and $\tilde {W}_{\textit{max}}\simeq 7.8 \tilde {F}_{{bv}}$ for $\tilde {F}_a \gg \tilde {F}_{{bv}}$ (strong forcing). In other words, if the forcing is strong enough, the maximum supported weight becomes independent of the forcing amplitude and scales as $\tilde {F}_{{bv}}$ . Equation (3.4) is consistent with the symmetry of the system: the long-term, time-averaged behaviour can only depend on even powers of $\alpha$ or $\tilde {F}_a$ , as reversing the sign of the forcing is equivalent to a mere phase shift. It is also worth noting that the regimes of weak and strong forcing correspond to small and large deformations, respectively, of the sheet relative to $\tilde {h}_{{eq}}$ . In other soft lubrication problems, these two limits also lead to different scalings (Essink et al. Reference Essink, Pandey, Karpitschka, Venner and Snoeijer2021).

Figure 4 $(d)$ shows that (3.4) provides insights into the evolution of the maximum supported weight $\tilde {W}_{\textit{max}}$ . At fixed $\tilde {F}_{{bv}}$ (left panel), increasing the forcing strength initially enhance adhesion, but $\tilde {W}_{\textit{max}}$ saturates beyond $\tilde {F}_a/\tilde {F}_{{bv}} \approx 20$ . Conversely, at fixed $\tilde {F}_a$ (right panel), softening the sheet by decreasing $\tilde {F}_{{bv}}$ first increases the maximum weight, before weakening adhesion. This shows an optimum at $\tilde {F}_{{bv}}^\star =(\tilde {\mu } \tilde {\omega } \tilde {B}^2)^\star \approx 0.13 \tilde {F}_a$ , with an associated bending stiffness $\tilde {B}^\star \approx 0.05 \tilde {F}_a^{3/2} (\tilde {\mu } \tilde {\omega } )^{-1/2}$ . If $\tilde {F}_{{bv}}$ and $\tilde {F}_a$ can be tuned independently, the largest $\tilde {W}_{\textit{max}}$ is achieved by maximising $\tilde {F}_a$ . In practice, this means using the most powerful available actuator and adjusting the sheet’s bending rigidity accordingly. Practical design would also need to account for energy consumption and for the fact that powerful vibration motors are typically both heavy and bulky; their lateral extent also limits adhesion (Poulain et al. Reference Poulain, Koch, Mahadevan and Carlson2025). Nevertheless, our prediction resonates qualitatively with experimental observations: using relatively small motors, Weston-Dawkes et al. (Reference Weston-Dawkes, Adibnazari, Hu, Everman, Gravish and Tolley2021) demonstrated $\tilde {W}_{\textit{max}}\approx 5\,\mathrm{N}$ using thin and soft sheets, while experiments by Colasante (Reference Colasante2016) (see also experiments from the same author reported by Ramanarayanan & Sánchez (Reference Ramanarayanan and Sánchez2024)) show $\tilde {W}_{\textit{max}}=\mathcal{O}({100\mathrm{N}})$ using stronger motors attached to stiffer sheets.

3.2.2. Contactless adhesion – equilibrium height

We observe a remarkable linear relationship between $ h_{{eq}}$ and $\mathcal{G}$ for $\alpha \gg 1$ in the inset of figure 4 $(b)$ , for $\mathcal{G} \lesssim \mathcal{G}_{\textit{max}}$ . Combining this observation with the asymptotic results of the previous section, we find

(3.5a) \begin{align} h_{{eq}} \approx \begin{cases} f \left (\dfrac {\mathcal{G}}{\alpha ^2}\right ) \quad \text{if} \quad \alpha \lesssim 1 \\[20pt] 0.05 \mathcal{G} \quad \text{if} \quad \alpha \gg 1 \end{cases} \quad \text{for} \quad \mathcal{G}\lesssim \mathcal{G}_{\textit{max}}, \end{align}

or, in dimensional units,

(3.5b) \begin{align} \tilde {h}_{{eq}} \approx \begin{cases} \tilde {H}_{{bv}} \times f\! \left (\dfrac {\tilde {W} \tilde {F}_{{bv}}}{\tilde {F}_a^2}\right ) \quad \text{if} \quad \tilde {F}_a \lesssim \tilde {F}_{{bv}} \\[20pt] 0.05 \dfrac {\tilde {W} \tilde {R}^2}{\tilde {B} } \quad \text{if} \quad \tilde {F}_a \gg \tilde {F}_{{bv}} \end{cases} \quad \text{for} \quad \tilde {W}\lesssim \tilde {W}_{\textit{max}}, \end{align}

with the maximum weight given by (3.4). The function $f$ represents the stable fixed points of (3.2) and is shown in figure 2 $(a)$ ; its analytical approximations are discussed by Poulain et al. (Reference Poulain, Koch, Mahadevan and Carlson2025). Equation (3.5b ) shows that, for strong forcing $\alpha \gg 1$ , the equilibrium height results from a balance between the bending force and gravity, with viscosity setting the maximum supported weight. Remarkably, the equilibrium height at the maximum supported weight is, in both limits,

(3.6) \begin{align} \tilde {h}_{{eq}}(\tilde {W}_{\textit{max}}) \approx 0.3 \tilde {H}_{{bv}}. \end{align}

While it is not straightforward to rationalise the linear relation between $h_{{eq}}$ and $\mathcal{G}$ , we note that $\tilde {B}/\tilde {R}^2 = \tilde {F}_{{\mathrm{bv}}}/\tilde {H}_{{bv}}=\tilde {k}_{{bv}}$ , with $\tilde {k}_{{bv}}$ a natural scale for stiffness, analogous to a spring constant. Then, (3.5b ) for $\alpha \gg 1$ can be written $\tilde {h}_{{eq}} \approx 0.05 \tilde {W}/\tilde {k}_{{bv}}$ .

Figure 5. Time-averaged shape $\langle h \rangle$ and the periodic deformation $h_d$ for $\alpha =20$ , $\mathcal{I}_{{bv}}={Re}_{{bv}}={Sq}_{{bv}}=0$ , and $\mathcal{G}=1,\,2,\,4$ in panels $(a,b,c)$ , respectively.

3.2.3. Contactless adhesion – sheet’s deformations

Some of the insights gained from the weak forcing analysis presented in § 3.1 carry over to the regime $\alpha \gg 1$ . In particular, figure 5 shows that the excitation of higher-order modes of deformation as $\mathcal{G}\rightarrow 0$ , $h_{{eq}}\rightarrow 0$ remains valid. One important difference, however, is the amplitude of the sheet’s deformation. In particular, we show in figure 6 $(a{,}b)$ an illustration of the sheet’s shape for $\alpha =20$ . For $\mathcal{G}\lt \mathcal{G}_{\textit{max}}$ , panel $(a)$ , the sheet remains convex, as expected from the large gravitational pull alone. For $\mathcal{G}\gt \mathcal{G}_{\textit{max}}$ , panel $(b)$ , the sheet becomes concave when the active forcing pulls it away from the wall. This, in turn, causes the sheet’s edge to come into contact with the substrate. We confirm in figure 6 $(c)$ that turning the sheet concave during part of the vibration cycle is what leads to contact with the wall, and therefore what sets $G_{\textit{max}}$ . This change in convexity occurs when the active pull is maintained for a sufficiently long time compared with the time it takes for the sheet to undergo a significant change in shape. This process is associated with the time scale $\tilde {T}_g(\tilde {h}_{{eq}})$ defined in (2.9), with $\tilde {h}_{{eq}}\sim \tilde {W}/\tilde {k}_{bv}$ from (3.5b ). The criterion for contact is then $\tilde {\omega } \tilde {T}_g(\tilde {h}_{{eq}}) \sim (\tilde {F}_{{bv}}/\tilde {W})^3 \lesssim 1$ , and we recover the scaling (3.4), i.e. $\tilde {W}_{\textit{max}} \sim \tilde {F}_{{bv}}$ for $\alpha \gg 1$ . We note that this understanding relies heavily on the convexity of the sheet and, consequently, on the specific distribution of weight and of active force. While we have assumed a uniform weight distribution, the effects of non-uniform weight distribution remain to be explored.

Figure 6. $(a{,}b)$ Sheet’s shape and pressure field over a period of vibration for $\alpha =20$ and $(a)$ $\mathcal{G}=4$ , $(b)$ $\mathcal{G}=8$ . The arrows represent the active force periodically pushing and pulling at the centre of the sheet. In panel $(b)$ , at $t=3\pi /2$ , the edges of the sheet touch the bottom wall. $(c)$ The difference in height between the sheet’s centre $h(0,t)$ and its edge $h(1,t)$ is a measure of the sheet’s convexity: when $\mathrm{min}_t(h(x=1,t)-h(x=0,t))\gt 0$ , the sheet always remain convex as in panel $(a)$ , and the relationship $h_{{eq}}\simeq 0.05 \mathcal{G}$ between equilibrium height and weight is verified (see figure 4 $b$ ). The existence of a concave part during the vibration cycle, as shown at $t=\pi ,\,5\pi /4$ and $3\pi /2$ in panel $(b)$ , is associated with contact. Filled symbols represent the case when the sheet never touches the wall, open symbols correspond to the sheet periodically touching it.

4.1. Rigid sheets

To qualitatively understand the effect of fluid inertia, we begin with an idealised, simplified set-up: a standard squeeze film where the gap between the wall and a rigid sheet oscillates harmonically. We consider the flow in the gap and project the Navier–Stokes momentum equations along a streamline:

(4.1) \begin{align} \tilde {\rho }\left ({\frac {\partial \tilde {v}_s}{\partial \tilde {t}}}+ \tilde {v}_s{\frac {\partial \tilde {v}_s}{\partial \tilde {s}}}\right )=-{\frac {\partial \tilde {p}}{\partial \tilde {s}}} + \tilde {\mu } \tilde {{\nabla} }^2 \tilde {v}_s, \end{align}

where $\tilde {s}$ and $\tilde {v}_s$ refer to the position and velocity along the streamline, respectively. We integrate this equation along a streamline whose ends are at the top of the sheet, close to $x=0$ , and within the gap at a position $\tilde {s}$ , as shown in figure 7 $(a)$ , to find an energy balance analogous to Bernoulli’s principle:

(4.2) \begin{align} \left ( \tilde {p}\big \rvert _{\tilde {s}=0}+\frac {1}{2} \tilde {\rho } \tilde {v}_s^2\big \rvert _{\tilde {s}=0}\right ) - \left ( \tilde {p}\big \rvert _{\tilde {s}}+\frac {1}{2} \tilde {\rho } \tilde {v}_s^2\big \rvert _{\tilde {s}}\right )+\int _0^{\tilde {s}} \tilde {\rho } {\frac {\partial \tilde {v}_s}{\partial \tilde {t}}}\biggr \rvert _{\tilde {s}'}\,{\mathrm{d}} \tilde {s}'=\Delta \tilde {p}_{{visc}}, \end{align}

with $\Delta \tilde {p}_{{visc}}=\int _0^{\tilde {s}} \tilde {\mu } \tilde {\boldsymbol{\nabla} }^2 \tilde {v}_s\,{\mathrm{d}} \tilde {s}'$ the pressure loss due to viscous stresses. According to lubrication scalings, $ \tilde {v}_s\rvert _{\tilde {s}=0} = \tilde {\dot {h}} \ll \tilde {v}_s\rvert _{\tilde {s}}\simeq - \tilde {\dot {h}}\tilde {x}/\tilde {h}$ , with overdots denoting time derivatives. We use $\tilde {s} \simeq \tilde {x}$ away from $\tilde {s}=0$ since the flow in the gap is predominantly horizontal. Then, $\int _0^{\tilde {s}} (\partial { \tilde {v}_s}/\partial { \tilde {t}}){\mathrm{d}} \tilde {s}' \sim [(\tilde {\dot h})^2-\tilde {\ddot h}/\tilde {h}] \, \tilde {x}^2$ and $\Delta \tilde {p}_{{visc}} \sim -\tilde {\mu } \tilde {\dot h}\tilde {x}^2/\tilde {h}^3$ . Equation (4.2) together with these scalings gives the pressure profile $\tilde {p}(\tilde {x},\tilde {t}) \sim (\tilde {R}^2-\tilde {x}^2)[-\tilde {\mu } \tilde {\dot h}/\tilde {h}^3 - \tilde {\rho } {\tilde {\ddot h}}/{\tilde {h}}+\tilde {\rho }{\tilde {\dot h}^2}/{\tilde {h}^2} ]$ , with prefactors missing. The formal inertial lubrication analysis from (2.5) leads to the following, in dimensionless form:

(4.3) \begin{align} p(x,t)&=\big (1-x^2\big )\left (-6 \frac {\dot h}{h^3} - \frac {3{{Re}}}{5} \frac {\ddot h}{h} + \frac {51{{Re}}}{35}\frac {\dot h^2}{h^2}\right ) + p(x=1,t), \end{align}

with $p(x=1,t)$ the pressure at the edge found from the boundary condition (2.11), ${\textit{Re}} = \tilde {\rho }_a \tilde {\omega } \tilde {h}_0^2/\tilde {\mu }$ , and $\tilde {h}_0$ the characteristic gap height. The first term describes the viscous losses and is captured by standard lubrication theory. The unsteady acceleration term, proportional to $-\ddot h/h$ , is equivalent to an added mass effect as discussed later. The edge pressure tends to decrease the pressure in the gap; however, the dominant effect of inertia is to convert the kinetic energy density of the fluid drawn into (or flowing out from) the gap into pressure. This is illustrated in figure 7, where a squeeze flow is imposed by setting the height as $h(t)=1+a\cos (t)$ , $0\lt a\lt 1$ . Over one period of oscillations, the net normal force from viscous effects cancels out: $\int _{-1}^1 \langle p_{v\textit{iscous}}(x,t) \rangle \,{\mathrm{d}}x=0$ , with $p_{{ viscous}}(x,t)=-6 (1-x^2 )\dot h(t)/h^3(t)$ . However, the force from inertial effects is positive, $ \int _{-1}^1 \langle p_{\textit{inertial}}(x)\rangle \,{\mathrm{d}}x = (1-7k/24)a^2{{Re}}/1.5 + \mathcal{O}({a^4})\gt 0$ , with $p_{\textit{inertial}}=p-p_{v\textit{iscous}}$ and $k=0.5$ a loss coefficient introduced in the boundary condition (2.11), showing that inertia leads to a normal force that pushes the sheet away from the wall. This effect has long been known in the context of bearings (e.g. Kuzma Reference Kuzma1968; Tichy & Winer Reference Tichy and Winer1970; Jones & Wilson Reference Jones and Wilson1975). In addition to modifying the pressure, a non-zero Reynolds number also alters the parabolic velocity profile predicted by the purely viscous theory: the horizontal velocity profile becomes a sixth-order polynomial which allows for secondary flows (figure 7 $d$ ), a common feature of pulsatile flows (Womersley Reference Womersley1955). We finally note that the unsteady acceleration term, proportional to $\ddot {h}$ , appears as a height-dependent added mass. Indeed, integrating the pressure over the sheet’s length $2\tilde {R}$ gives rise to a force per unit length $-\tilde {m}_a \tilde {\ddot {h}}$ , with $\tilde {m}_a=4\tilde {\rho } \tilde {R}^3/5\tilde {h}$ an added mass per unit length (see also (D5)). In contrast, the added mass of a plate far from any boundary is proportional to $\tilde {\rho } \tilde {R}^2$ (Brennen Reference Brennen1982).

Figure 7. Squeeze flow of a rigid plate moving normal to a wall with a height evolving as $h(t)=1+a\cos (t)$ , shown for illustration here with $a=0.4$ . The pressure is computed from (4.3) and the fluid velocity from the calculations carried out in Appendix A. $(a)$ Streamlines of the flow associated with a rigid sheet moving towards a wall or away from it. $(b)$ Integral of pressure in space for the viscous component and inertial component. Integrated in time, the viscous component averages to zero while the inertial component gives a positive force. The inertial pressure is rescaled by the Reynolds number. $(c{,}d)$ Velocity profiles (arrows) and pressure field (colours) isolating $(c)$ the dominant viscous flow and $(d)$ the inertial corrections.

While we initially considered a sheet constrained in position for simplicity, we are more interested in the case of a weightless rigid sheet driven periodically by a dimensionless force $\alpha \cos (t)$ . In this case, we use a two-time scale asymptotic analysis to find the evolution of the height $h(t)$ in the limit $\alpha ^2{{Re}} \ll 1$ . We present the derivation in Appendix D, where we show that the sheet slowly moves away from the wall to accommodate the inertial pressure increase while maintaining a force balance. The sheet’s time-averaged position has the following evolution equation, in dimensionless and dimensional form:

(4.4) \begin{align} \begin{split} \frac {1}{\alpha ^2\langle h \rangle ^2}{\frac {{\mathrm{d}}\langle h \rangle }{\mathrm{d}t}}=c_k{{Re}} \langle h \rangle ^5, \quad {\frac {{\mathrm{d}}\langle \tilde {h} \rangle }{{\mathrm{d}}\tilde {t}}=c_k \frac {\tilde {F}_a^2 \tilde {\rho }_a \langle \tilde {h} \rangle ^7}{\tilde {\mu }^3 \tilde {R}^8}=\frac {c_k \langle \tilde {h} \rangle }{\tilde {T}_{{a,i}}\!\left (\langle \tilde {h} \rangle \right )}}, \quad c_k=\frac {1-\frac {7k}{16}}{224}. \end{split} \end{align}

The fact that ${\mathrm{d}}(\langle h \rangle )/{\mathrm{d}}t \sim \langle h\rangle ^7$ – a much stronger dependence on $\langle h\rangle$ than the effects of gravity and elastohydrodynamic, see (3.2) – follows from the expressions of the relevant time scales discussed in § 2.4. In particular, the inertial time scale shows the strongest dependence on the height, $\tilde {T}_{a,i}(\tilde {h})\sim \tilde {h}^{-6}$ . We also note that the boundary condition (2.11) appears in the prefactor $c_k$ and leads to a correction of the order of $7k/16 \simeq 22\,\%$ for $k=0.5$ as compared with the standard boundary condition $p=0$ .

4.2. Soft sheets

We expect that the destabilising influence of fluid inertia extends to soft sheets, weakening the viscous adhesion mechanism discussed in the previous section. We solve the governing equations with ${Sq}_{{bv}}=I_{{bv}}=0$ and for $\alpha =1$ and $\alpha =20$ , corresponding to the regimes of weak and strong forcing, respectively. For both values, we systematically studied the relationship between the equilibrium height $h_{\textit{eq}}$ and the weight $\mathcal{G}$ as a function of the Reynolds number $ \textit{Re}_{{bv}}$ . As discussed in § 2.7, for each set of parameters, we define the equilibrium Reynolds number ${\textit{Re}}_{{eq}}=h_{{eq}}^2 {{Re}}_{{bv}}$ that is based on the equilibrium height rather than on the height scale $\tilde {H}_{{bv}}$ and which is representative of the actual magnitude of inertial over viscous effects; hence, ${\textit{Re}}_{{eq}}$ may be small even when ${\textit{Re}}_{{bv}}$ is large.

Figure 8 summarises our results for $\alpha =1$ . There, we show examples of bifurcation diagrams and the associated regime map of adhesion. For small values of ${\textit{Re}}_{{bv}}$ , increasing the Reynolds number first marginally increases the maximum weight $\mathcal{G}_{\textit{max}}$ . Further increasing ${\textit{Re}}_{{bv}}$ then leads to a sharp decrease in $\mathcal{G}_{\textit{max}}$ , explained by the destabilising effect of fluid inertia discussed previously. Figure 8 $(d)$ shows the range of accessible equilibrium Reynolds numbers $ \textit{Re}_{{eq}}$ as a function of $ \textit{Re}_{{bv}}$ . This demonstrates that $ \textit{Re}_{{bv}}$ has little effect on the dynamics as long as $ \textit{Re}_{{eq}} \lesssim 1$ . Once $ \textit{Re}_{{eq}}=\mathcal{O}({1})$ , $\mathcal{G}_{\textit{max}}$ decreases together with $h_{{eq}}$ to ensure $ \textit{Re}_{{eq}}$ remains relatively small. In particular, we observe a sharp drop in $\mathcal{G}_{\textit{max}}$ for $ \textit{Re}_{{bv}} \gtrsim 30 \approx 1/e_1^2$ , with $e_1$ the height scale introduced in (3.2a ). As $ \textit{Re}_{{bv}}$ keeps increasing, the first branch of the bifurcation diagram (continuous black line for ${h_{{eq}}}\gt 0.1$ in figure 8 $a$ ) is eventually not accessible anymore for ${\textit{Re}}_{{bv}} \gtrsim 200 \approx 1/e_2^2$ . We observe a similar drop of $\mathcal{G}_{\textit{max}}$ for ${\textit{Re}}_{{bv}} \approx 2000 \approx 1/e_3^2$ , as the first part of the second equilibrium branch (corresponding to the excitation of the second eigenmode $\zeta _2$ for $h_{{eq}} \approx e_2$ ) is also no longer accessible.

Figure 8. $(a{,}b)$ Equilibrium height as a function of the dimensionless weight with ${Sq}_{{bv}}=\mathcal{I}_{{bv}}=0$ and $\alpha =1$ for $(a)$ $ \textit{Re}_{{bv}}\lt 200$ and $(b)$ $ \textit{Re}_{{bv}}\gt 200$ . Black lines are the stable equilibria of (3.2). $(c)$ Phase diagram showing the accessible weights as a function of $ \textit{Re}_{{bv}}$ . The first equilibrium branch corresponds to $h_{{eq}}\gt 0.1$ , the second branch to $h_{{eq}}\lt 0.1$ . $(d)$ Reynolds number based on the equilibrium height $h_{{eq}}$ , ${\textit{Re}}_{\textit{eq}}=h_{{eq}}^2 {Re}_{{bv}}$ as a function of the control parameter $ \textit{Re}_{{bv}}$ . The dashed lines represent the expected behaviour if fluid inertia did not affect the system.

Figure 9. Effect of the fluid inertia with ${{{Sq}}}_{{bv}}=\mathcal{I}_{{bv}}=0$ and ${\textit{Re}}_{{bv}}\gt 0$ for $\alpha =20$ . $(a)$ Equilibrium height as a function of the dimensionless weight for the regime of contactless adhesion. $(b)$ Regime maps and the associated $(c)$ range of Reynolds number based on equilibrium height $ \textit{Re}_{{eq}}$ . The inertial lubrication theory is not expected to be valid for ${\textit{Re}}_{\textit{eq}}\gtrsim 50.$ In panel $(b)$ , we show illustrations of $h_d$ (as defined in (3.3)) for ${\textit{Re}}_{{bv}}$ = 100, 1000 and 5000 at $t=\pi /2$ and $3\pi /2$ .

Figure 9 shows that the system’s behaviour is qualitatively similar in the regime of strong forcing ( $\alpha =20$ ): the maximum supported weight decreases as ${\textit{Re}}_{{bv}}$ increases, with sharp drops at specific Reynolds numbers, while the equilibrium height shows minor variations with ${\textit{Re}}_{{bv}}$ . Similarly to the weak forcing case, these drops are associated with higher-order modes as illustrated in panel $(b)$ . The main difference between the two regimes is that ${\textit{Re}}_{{eq}}$ can become relatively large when $\alpha \gg 1$ : we observe values up to ${\textit{Re}}_{{eq}}=\mathcal{O}({10^2})$ for $ \textit{Re}_{{bv}}=\mathcal{O}({10^4})$ . For such large values, inertial lubrication theory may no longer be valid, and higher-order corrections, or full Navier–Stokes simulations, may be necessary to accurately describe the flow and $\mathcal{G}_{\textit{max}}$ . This may be relevant to the experiments of Weston-Dawkes et al. (Reference Weston-Dawkes, Adibnazari, Hu, Everman, Gravish and Tolley2021).

5.1. Rigid sheets

To build intuition, we first consider the rigid squeeze film set-up: $\tilde {h}=\tilde {h}_0(1+a\cos (\tilde {\omega } \tilde {t}))$ . The effects of compressibility are quantified with the squeeze number ${{{Sq}}}=\tilde {\mu }\tilde {\omega }\tilde {R}^2/\tilde {h}_0^2\tilde {p}_a$ , which compares the viscous stresses with the ambient pressure $\tilde {p}_a$ . When the sheet approaches the wall, a competition arises between the outflow of fluid from the gap, resisted by viscosity, and the compression of the fluid, resisted by its pressure. For ${{{Sq}}} \ll 1$ , an incompressible description is appropriate. Conversely, for ${{{Sq}}} \gg 1$ , viscous stresses become so large that there is almost no outflow and the fluid is compressed with the sheet acting as a piston. Likewise, fluid is expanded when the surface moves away from the wall. In this limit, also examined by Salbu (Reference Salbu1964), the work per unit area to expand from $\tilde {h}_0$ to $\tilde {h}_0(1+a)$ is $\tilde {w}_{+a}=\int _{\tilde {h}_0}^{\tilde {h}_0(1+a)} \tilde {p}\,{\mathrm{d}}\tilde {h}=\tilde {p}_a \tilde {h}_0\ln (1+a)$ , and the work to compress from $\tilde {h}_0$ to $\tilde {h}_0(1-a)$ is $\tilde {w}_{-a}=\tilde {p}_a\tilde {h}_0\ln (1-a)$ . This gives $\lvert \tilde {w}_{-a}\rvert -\lvert \tilde {w}_{+a}\rvert =-\tilde {p}_a\tilde {h}_0\ln (1-a^2)=\tilde {p}_a \tilde {h}_0 \tilde {a}^2 + \mathcal{O}({a^4})\gt 0$ , i.e. it is more energetically costly to compress than to expand an ideal gas by the same volume increment. This is because the isothermal bulk modulus of the gas is its pressure, which increases during compression and decreases during expansion.

We now come back to the original setting, the displacement of the rigid sheet not constrained, but controlled by a time-periodic force $\alpha \cos (t)$ . We can anticipate that at each cycle, there will be more expansion than compression, i.e. the sheet will slowly move away from the wall. An asymptotic expansion performed in Appendix E confirms this picture and shows that the sheet’s time averaged position follows at $\mathcal{O}({\alpha ^2 {{{Sq}}}})$ , in dimensional and dimensionless form:

(5.1) \begin{align} \frac {1}{\alpha ^2\langle h \rangle ^2}{\frac {{\mathrm{d}}\langle h \rangle }{\mathrm{d}t}}=\frac {3}{8}{Sq}_{{bv}} \langle h \rangle , \quad {\frac {{\mathrm{d}}\langle \tilde {h} \rangle }{{\mathrm{d}}\tilde {t}}=\frac 38 \frac {\tilde {F}_a^2 \langle \tilde {h} \rangle ^3}{\tilde {\mu } \tilde {R}^6 \tilde {p}_a}=\frac 38 \frac {\langle \tilde {h} \rangle }{\tilde {T}_{{a,c}} (\langle \tilde {h}\rangle)}}. \end{align}

We can compare (5.1) with (3.2), where the isolated effect of gravity on a rigid sheet leads to ${\mathrm{d}}\langle h\rangle /{\mathrm{d}}t = \mathcal{G} \langle h \rangle ^3/4$ . The fact that these two different physical processes lead to the same scaling ${\mathrm{d}}\langle h\rangle /{\mathrm{d}}t \sim \langle h \rangle ^3$ originates from the expression of the compressible time scale $\tilde {T}_{{a,c}}$ and gravitational time scale $\tilde {T}_{g}$ in (2.9). It shows that the first-order effect of compressibility can be interpreted as an increase in the effective dimensionless weight $\mathcal{G}_{{eff}}$ (or dimensional weight $\tilde {W}_{{ eff}}$ ) such that

(5.2) \begin{align} \mathcal{G}_{{eff}} = \mathcal{G} + \frac 32 \alpha ^2 {Sq}_{{bv}}, \quad \tilde {W}_{{eff}} = \tilde {W} + \frac 32 \frac {\tilde {F}_a^2}{\tilde {p}_a \tilde {R}^2}. \end{align}

While we focus on small squeeze numbers, we note for completeness that the averaged normal force acting on a sheet subject to a squeeze film motion can be approximated analytically for arbitrary values of ${{{Sq}}}_{{bv}}$ (Taylor & Saffman Reference Taylor and Saffman1957; Langlois Reference Langlois1962).

5.2. Soft sheets

Figure 10. $(a{,}b)$ Effect of the fluid compressibility for $ \textit{Re}_{{bv}}=\mathcal{I}_{{bv}}=0$ and ${{{Sq}}}_{{bv}}\gt 0$ . Panels $(a{,}b)$ correspond to the weak forcing regime with $\alpha =1$ and panels $(c{,}d)$ to the strong forcing regime with $\alpha =20$ . The dashed lines in panels $(b)$ and $(d)$ are derived from (5.2) and show $\mathcal{G}_{\textit{max}}({{{Sq}}}_{{bv}})=\mathcal{G}_{\textit{max}}({{{Sq}}}_{{bv}}=0)-1.5 \alpha ^2 {{{Sq}}}_{{bv}}$ . We primarily captured the first equilibrium branch and did not systematically investigate the entire extent of the bifurcation diagrams.

We solve the governing equations (2.4)–(2.6) with boundary conditions (2.10) and (2.11) in the inertialess limit ( $ \textit{Re}_{{bv}}=\mathcal{I}_{{bv}}=0$ ); in this case, the pressure boundary condition reduces to $p=0$ . Our results for the weak forcing regime, $\alpha =1$ , are shown in figure 10 $(a{,}b)$ . They are consistent with (5.2) and the main effect of compressibility is indeed to modify the effective weight. As compressibility translates the bifurcation diagrams to lower weights, we observe the possibility of reaching an equilibrium height $h_{{eq}}\gt 0$ for $\mathcal{G}\lt 0$ , corresponding to gravity pointing towards the wall. In these cases, the repulsive effect of compressibility can balance the adhesive effects of both elastohydrodynamics and gravity. This is reminiscent of squeeze-film and near-field acoustic levitation (Shi et al. Reference Shi, Feng, Hu, Zhu and Cui2019), where objects can levitate above rapidly vibrating surfaces.

Figure 10 $(c{,}d)$ shows similar results for the strong forcing regime, $\alpha =20$ . The relation (5.2) predicts $\mathcal{G}_{\textit{max}}$ accurately up to ${{{Sq}}}_{{bv}} \approx 10^{-3}$ , but underestimates the effective weight by up to $50\,\%$ at larger squeeze numbers, likely because of the interactions between elastic deformations and compressibility are neglected by this simple scaling. The effects on the equilibrium curves are also more complex than a simple translation of the bifurcation diagram; in particular, we observe the possibility of hysteresis.