3.1. Motivation and local coordinate system

In § 2, we derived the general net balance equations for a thin sheet of viscous fluid. Now, we show how to supplement these equations with effective boundary conditions that apply at a free edge of the sheet. Near such an edge, there is a boundary layer in which the in-plane and transverse dimensions of the sheet become comparable, as illustrated in figure 2(a). We note that the solution for the flow in this inner region was found numerically by Munro & Lister (Reference Munro and Lister2018), but we show that the effective boundary conditions for the bulk flow can be obtained just using asymptotic matching. In this derivation, we consider the general situation where the edge of the sheet can be arbitrarily curved, though, for simplicity, we assume that it remains approximately planar. We use intrinsic curvilinear coordinates embedded in the sheet edge; a similar derivation is presented by O’Kiely (Reference O’Kiely2017), though without the inclusion of surface tension. Since the problem is quasi-steady, we can focus on determining the instantaneous boundary conditions and, for the moment, suppress the dependence on time $t$ .

Figure 2. (a) Sketch of the inner region at the edge of a thin sheet. (b) Sketch of the curvilinear coordinate system employed at the edge of the sheet.

An edge of the sheet is identified as a curve on which $h=0$ . As illustrated in figure 2(b), we parametrise the projection of this curve onto the $\boldsymbol {\tilde {x}}=(\tilde {x},\tilde {y})$ -plane using arc-length $s$ , and denote the corresponding planar tangent vector as $\boldsymbol {\tilde {t}}(s)$ . We fix the orientation such that the planar normal pointing outwards from the sheet edge is given by $\boldsymbol {\tilde {n}}=\boldsymbol {k}\times \boldsymbol {\tilde {t}}$ , where we recall that $\boldsymbol {k}$ denotes the unit vector in the $z$ -direction. The normal and tangent vectors are related by the Serret–Frenet formulae (Kreyszig Reference Kreyszig1959)

(3.1a,b) \begin{align} \frac {\mathrm {d}\boldsymbol {\tilde {t}}}{\mathrm {d}s}&=\kappa \boldsymbol {\tilde {n}}, & \frac {\mathrm {d}\boldsymbol {\tilde {n}}}{\mathrm {d}s}&=-\kappa \boldsymbol {\tilde {t}}, \end{align}

where $\kappa (s)$ is the curvature of the edge (projected onto the $(\tilde{x},\tilde{y})$ -plane). The position of any point in the sheet can be expressed in the form

(3.2) \begin{equation} \boldsymbol {r}(s,n,z)=\boldsymbol {\tilde {x}}+z\boldsymbol {k}=\int _0^s \boldsymbol {\tilde {t}}(s^\prime )\, \mathrm {d}s^\prime + n\boldsymbol {\tilde {n}}+z\boldsymbol {k}, \end{equation}

where $n\lt 0$ and $H^-(s,n)\lt z\lt H^+(s,n)$ . The edge of the sheet is defined to be at $n=0$ , where we have $H^-(s,0)=H^+(s,0)=H(s,0)$ .

Now, our strategy is to express the integrated governing equations (2.9), (2.11) and (2.13) using the local coordinates $(s,n)$ . Then, at the edge of the sheet, since $h(s,0)=0$ , we seemingly have five boundary conditions

(3.3) \begin{equation} T_{nn}=T_{sn}={N_n}=M_{nn}=M_{sn}=0 \quad \text {at }n=0{,} \end{equation}

where subscripts denote components of the tensor or vector. However, this is one too many boundary conditions for the outer problem. This issue was first addressed in the context of thin elastic plates (see, for example, Love Reference Love1927; Timoshenko & Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1959). We resolve the difficulty by rescaling into the boundary layer at the edge and thus deriving the appropriate effective boundary conditions to apply to the outer problem.

3.2. Edge boundary layer

We examine the boundary layer by defining

(3.4a–d) \begin{align} \kern-4em n=\epsilon \hat {n}, \qquad T_{ss}=\hat {T}_{ss}, \qquad T_{sn}=\epsilon \hat {T}_{sn}, \qquad T_{nn}=\epsilon \hat {T}_{nn}, \end{align}

(3.4e–i) \begin{align} M_{ss}=\hat {M}_{ss}, \qquad M_{sn}=\hat {M}_{sn}, \qquad M_{nn}=\epsilon \hat {M}_{nn}, \qquad N_s = \frac {\hat {N}_s}{\epsilon }, \qquad N_n = \hat {N}_n, \end{align}

where we denote variables in the boundary layer by hats (not to be confused with the transverse stress components as in § 2). The different scalings of the tensions, shears and bending moments are made to obtain non-trivial balances in the dimensionless integrated Stokes equations, (2.9), (2.11) and (2.13), which become (see, for example, van de Fliert, Howell & Ockendon Reference van de Fliert, Howell and Ockendon1995)

(3.5a) \begin{align} \frac {\partial \hat {T}_{ss}}{\partial s}+\frac {\partial }{\partial \hat {n}}(\hat {\ell } \hat {T}_{sn})-\epsilon \kappa \hat {T}_{sn}&= 0, \end{align}

(3.5b) \begin{align} \epsilon \frac {\partial \hat {T}_{sn}}{\partial s}+\frac {\partial }{\partial \hat {n}}(\hat {\ell } \hat {T}_{nn})+\kappa \hat {T}_{ss}&= 0, \end{align}

(3.5c) \begin{align} {\frac {\partial \hat {N}_{s}}{\partial s}+\frac {\partial }{\partial \hat {n}}(\hat {\ell } \hat {N}_{n})} &= 0, \end{align}

(3.5d) \begin{align} \kern 20pt\epsilon \frac {\partial }{\partial s}\left (\hat {M}_{ss}+\hat {H} \hat {T}_{ss}\right )+\frac {\partial }{\partial \hat {n}}\left [\hat {\ell }\left (\hat {M}_{sn}+\epsilon \hat {H}\hat {T}_{sn}\right )\right ]-\epsilon \kappa \left (\hat {M}_{sn}+\epsilon\!\hat {H}\hat {T}_{sn}\right ) & = \hat {\ell } {\hat {N}_{s}}, \end{align}

(3.5e) \begin{align} \frac {\partial }{\partial s}\left (\hat {M}_{sn}+\epsilon \hat {H} \hat {T}_{sn}\right )+\frac {\partial }{\partial \hat {n}}\left [\hat {\ell }\left (\hat {M}_{nn}+\hat {H}\hat {T}_{nn}\right )\right ]+\kappa \left (\hat {M}_{ss}+\hat {H}\hat {T}_{ss}\right ) & = \hat {\ell } {\hat {N}_{n}}, \end{align}

where $\hat {\ell }=1-\epsilon \kappa \hat {n}$ is the metric coefficient. The boundary conditions (3.3) at the edge of the sheet are transformed to

(3.6) \begin{equation} \hat {T}_{nn}=\hat {T}_{sn}={\hat {N}_{n}}=\hat {M}_{nn}=\hat {M}_{sn}=0 \quad \text {at }\hat {n}=0. \end{equation}

We now expand our variables as asymptotic series in powers of $\epsilon$ , i.e. $\hat {T}_{ss}\sim \hat {T}_{ss0}+\epsilon \hat {T}_{ss1}+\cdots$ as $\epsilon \to 0$ . Note that the scalings (3.4) already assume the leading-order matching conditions

(3.7) \begin{align} T_{sn0},\, T_{nn0},\, M_{nn0}\to 0 \quad \textrm {as } n\to 0,&&&{\hat {N}_{s0}}\to 0 \quad \textrm {as }\hat {n}\to -\infty . \end{align}

As anticipated above and suggested by the sketch in figure 2(a), we also assume that, although the sheet thickness $h$ varies significantly in the edge layer, the centre-surface $H$ does not, so that $\hat {H}(s,n)\sim \hat {H}_0(s)+O(\epsilon )$ .

At leading order, we find from (3.5d ) that

(3.8) \begin{equation} {\hat {N}_{s0}}=\frac {\partial \hat {M}_{sn0}}{\partial \hat {n}}. \end{equation}

Substituting this result into (3.5c ) gives, at leading order,

(3.9) \begin{equation} \frac {\partial }{\partial \hat {n}}\left ({\hat {N}_{n0}}+\frac {\partial \hat {M}_{sn0}}{\partial s}\right )=0. \end{equation}

By applying the boundary conditions (3.6), we deduce that

(3.10) \begin{equation} {\hat {N}_{n0}}+\frac {\partial \hat {M}_{sn0}}{\partial s}=0, \end{equation}

and, by matching to the outer region, we deduce the leading-order effective boundary condition

(3.11) \begin{equation} {\hat {N}_{n0}}+\frac {\partial M_{sn0}}{\partial s}=0\quad \textrm {at }n=0. \end{equation}

However, by combining (3.5b ) and (3.5e ) at leading order, we obtain

(3.12) \begin{equation} {\hat {N}_{n0}}=\frac {\partial \hat {M}_{sn0}}{\partial s} +\frac {\partial \hat {M}_{nn0}}{\partial \hat {n}} +\kappa \hat {M}_{ss0}, \end{equation}

which can be used to eliminate the shear stress and express (3.11) purely in terms of the bending-moment tensor. In summary, we can express the leading-order effective boundary conditions for the outer problem as

(3.13) \begin{equation} T_{sn}=T_{nn}=M_{nn}=2\frac {\partial M_{sn}}{\partial s}+\frac {\partial M_{nn}}{\partial n}+\kappa M_{ss}=0\qquad \textrm {at }n=0. \end{equation}

A benefit of this method, when compared with similar derivations carried out by Howell, Kozyreff & Ockendon (Reference Howell, Kozyreff and Ockendon2009)and O’Kiely (Reference O’Kiely2017), for example, is that we did not need to calculate any velocity components of the fluid; instead, we worked with the tensions and bending moments. Moreover, incorporating surface tension contributions into the definitions of the net tensions and bending moments made it straightforward to generalise the boundary conditions found by O’Kiely (Reference O’Kiely2017) to include surface tension effects. We note that an alternative derivation of the effective boundary conditions based on a virtual work argument is presented by Srinivasan et al. (Reference Srinivasan, Wei and Mahadevan2017), though there appear to be some sign inconsistencies in their formulation.

Armed with the boundary conditions (3.13), we are ready to tackle the outer governing equations (2.7), (2.9), (2.11) and (2.13). As noted in § 2, we must first derive constitutive relations for the integrated tensions and bending moments by analysing the asymptotic limit as $\epsilon \rightarrow 0$ . In doing so, we choose to focus on a model geometrical set-up in which a disc of viscous fluid retracts under surface tension, and then examine the response of the system to small transverse perturbations.

4.1. Leading-order solution

Now, we invoke the dimensionless Newtonian constitutive relations, namely

(4.1a–c) \begin{align} \tilde {\boldsymbol {\sigma }}&=-p\tilde {\boldsymbol {I}}+\tilde {\boldsymbol {\nabla }}\tilde {\boldsymbol {u}}+{\tilde {\boldsymbol {\nabla }}\tilde {\boldsymbol {u}}}^t, & \epsilon ^2\sigma _{zz}&=-p-2\tilde {\boldsymbol {\nabla }}\cdot \tilde {\boldsymbol {u}}, & \epsilon ^2\hat {\boldsymbol {\sigma }}&= \frac {\partial \tilde {\boldsymbol {u}}}{\partial z}+ \epsilon ^2\tilde {\boldsymbol {\nabla }}w, \end{align}

where the pressure $p$ has been made dimensionless with $\gamma /\epsilon L$ , the same scaling as the in-plane stress. Here, we have assumed that the viscosity $\eta$ is constant; the theory developed below is generalised to include small viscosity variations in Appendix A. When we express the dependent variables as asymptotic expansions of the form $\tilde {\boldsymbol {u}}\sim \tilde {\boldsymbol {u}}_0+\epsilon ^2\tilde {\boldsymbol {u}}_1+\cdots$ , we immediately see from (4.1c ) that the flow is extensional to leading order, with the in-plane velocity $\tilde {\boldsymbol {u}}_0$ independent of $z$ , i.e.

(4.2) \begin{equation} \tilde {\boldsymbol {u}}_0 = \tilde {\boldsymbol {u}}_0\left (\tilde {\boldsymbol {x}},t\right ). \end{equation}

The net mass-conservation equation (2.7) thus reduces to

(4.3) \begin{equation} \frac {\partial h_0}{\partial t}+ \tilde {\boldsymbol {\nabla }}\cdot \left (h_0\tilde {\boldsymbol {u}}_0\right )=0. \end{equation}

Next, we use the constitutive relation (4.1a ) to evaluate the leading-order in-plane stress $\tilde {\boldsymbol {\sigma }}_0$ and thus from (2.10), the in-plane tension tensor, namely

(4.4) \begin{equation} \boldsymbol{T}_0\,{=}\,\left (2+2h_0\tilde {\boldsymbol {\nabla }}\cdot \tilde {\boldsymbol {u}}_0\right )\tilde {\boldsymbol {I}}+h_0\left (\tilde {\boldsymbol {\nabla }}\tilde {\boldsymbol {u}}_0+\tilde {\boldsymbol {\nabla }}\tilde {\boldsymbol {u}}_0^t\right ). \end{equation}

In this expression, the first factor of $2$ is the contribution due to surface tension, and the remaining terms (proportional to $h_0$ ) are the viscous contributions. Let us denote the region of the $\tilde {\boldsymbol {x}}$ -plane occupied by the sheet by $\Omega$ , with boundary $\partial \Omega$ . Then, the governing equation and boundary condition for $\boldsymbol{T}_0$ , namely

(4.5a,b) \begin{align} \tilde {\boldsymbol {\nabla }}\cdot \boldsymbol{T}_0&=\boldsymbol {0} \quad \text {in }\Omega, & \boldsymbol{T}_0\cdot \tilde {\boldsymbol {n}}&=\boldsymbol {0} \quad \text {on }\partial \Omega , \end{align}

follow from (2.9) and (3.13), respectively. In principle, given $h_0$ , the boundary-value problem (4.4) and (4.5) determines both $\boldsymbol{T}_0$ and $\tilde {\boldsymbol {u}}_0$ (up to an irrelevant rigid-body motion), and then $h_0$ can be stepped forward in time using (4.3).

In this paper, we focus on the behaviour of a sheet whose thickness is spatially uniform to leading order, i.e. for which

(4.6) \begin{equation} h_0\left (\tilde {\boldsymbol {x}},t\right )=\psi (t). \end{equation}

In this case, the problem (4.4) and (4.5) implies that

(4.7) \begin{equation} \boldsymbol{T}_0\left (\tilde {\boldsymbol {x}},t\right )=\boldsymbol {0}. \end{equation}

Although the flow is extensional at leading order, the viscous and surface tension terms in (2.10) exactly balance, so the leading-order tension in the sheet is identically zero. Up to an arbitrary rigid-body translation and rotation, the corresponding leading-order velocity is found from (4.4) to be given by

(4.8) \begin{equation} \tilde {\boldsymbol {u}}_0\left (\tilde {\boldsymbol {x}},t\right )=-\frac {\tilde {\boldsymbol {x}}}{3\psi (t)}. \end{equation}

Then, the mass-conservation equation (4.3) reduces to $\dot {\psi }-2/3=0$ (with the dot denoting differentiation) and, therefore,

(4.9) \begin{equation} h_0(\tilde {\boldsymbol {x}},t)=\psi (t)=1+\frac {2t}{3}. \end{equation}

In this leading-order solution, the initially uniform sheet thickness remains uniform and grows linearly with $t$ , as the sheet retracts under surface tension. If we define in-plane Lagrangian variables $\tilde {\boldsymbol {X}}$ by

(4.10) \begin{equation} \tilde {\boldsymbol {x}}=\frac {\tilde {\boldsymbol {X}}}{\sqrt {\psi (t)}}, \end{equation}

then, with respect to $\tilde {\boldsymbol {X}}$ , the sheet domain, which we will now denote by $\Omega _X$ , remains fixed for all time. Of course, this result is subject to the caveat that the aspect ratio of the sheet must remain small, which requires that $\psi (t)\ll \epsilon ^{-2/3}$ .

4.2. Small thickness perturbations

We have seen that the leading-order tension in the sheet is identically zero when the sheet thickness is spatially uniform. We now introduce small thickness perturbations of order $\epsilon ^2$ which, as we will demonstrate, are sufficient to induce regions of compression and thus the possibility of buckling. To simplify the analysis, we make the change of variables from $ (\tilde {\boldsymbol {x}},z,t )$ to $ (\tilde {\boldsymbol {X}},z,t )$ , where $\tilde {\boldsymbol {X}}$ are the Lagrangian in-plane variables introduced in (4.10). We then perturb about the above leading-order solution as follows:

(4.11a) \begin{align} h\left (\tilde {\boldsymbol {X}},t\right )&\sim \psi (t) +\epsilon ^2h_1\left (\tilde {\boldsymbol {X}},t\right ) +O\bigl (\epsilon ^4\bigr ), \end{align}

(4.11b) \begin{align} \bar {\boldsymbol {u}}\left (\tilde {\boldsymbol {X}},t\right )&\sim -\frac {\tilde {\boldsymbol {X}}}{3\psi (t)^{3/2}} +\epsilon ^2\bar {\boldsymbol {u}}_1\left (\tilde {\boldsymbol {X}},t\right ) +O\bigl (\epsilon ^4\bigr ), \end{align}

where the initial thickness perturbation $h_1 (\tilde {\boldsymbol {X}},0 )$ is assumed to be specified. We impose the constraint

(4.12) \begin{equation} \iint _{{\Omega _X}} h_1\left (\tilde {\boldsymbol {X}},0\right )\,\mathrm {d}\tilde {\boldsymbol {X}}=0, \end{equation}

so that the mass of the sheet is accounted for entirely by the leading-order solution.

We also make small perturbations to the centre-surface $H$ , so that

(4.13) \begin{equation} H\left (\tilde {\boldsymbol {X}},t\right )\sim \delta H_1\left (\tilde {\boldsymbol {X}},t\right ), \end{equation}

where $0\lt \delta \ll 1$ . The initial centre-surface displacement $\delta H_1 (\tilde {\boldsymbol {X}},0 )$ is again assumed to be specified and small. The restriction to small centre-surface perturbations allows us to linearise about the base state $H=0$ , and the size of $\delta$ in relation to $\epsilon$ is irrelevant. The resulting theory models the onset of buckling, should it occur, and remains valid so long as $H_1$ remains smaller than $O (\delta ^{-1} )$ .

We recall that the in-plane tension tensor $\boldsymbol{T}$ is zero at leading order, and its asymptotic expansion thus takes the form

(4.14) \begin{equation} \boldsymbol{T}\left (\tilde {\boldsymbol {X}},t\right )\sim \epsilon ^2\boldsymbol{T}_1\left (\tilde {\boldsymbol {X}},t\right )+O\left (\epsilon ^4,\epsilon ^2\delta ^2\right ). \end{equation}

The first-order in-plane stress $\tilde {\boldsymbol {\sigma }}_1$ is found by substituting the expansions (4.11)–(4.13) into the governing equations (2.3)–(2.5) and constitutive relations (4.1). The first non-zero contribution $\boldsymbol{T}_1$ to the tension is then found from (2.10), which produces

(4.15) \begin{equation} \boldsymbol{T}_1=2\left (\psi ^{3/2}\tilde {\boldsymbol {\nabla }}\cdot \bar {\boldsymbol {u}}_1-\frac {h_1}{\psi }\right )\tilde {\boldsymbol {I}}+ \psi ^{3/2}\left (\tilde {\boldsymbol {\nabla }}\bar {\boldsymbol {u}}_1+\tilde {\boldsymbol {\nabla }}\bar {\boldsymbol {u}}_1^t\right ), \end{equation}

where, now, the gradient operator $\tilde {\boldsymbol {\nabla }}$ is performed with respect to the new in-plane variables $\tilde {\boldsymbol {X}}$ .

The first-order tension satisfies a boundary-value problem analogous to (4.5), that is,

(4.16a,b) \begin{align} \tilde {\boldsymbol {\nabla }}\cdot \boldsymbol{T}_1&=\boldsymbol {0} \quad \text {in }{\Omega _X}, & \boldsymbol{T}_1\cdot \tilde {\boldsymbol {n}}&=\boldsymbol {0} \quad \text {on }\partial {\Omega _X} \end{align}

(with no contributions due to perturbations in $\partial {\Omega _X}$ because $\boldsymbol{T}_0$ is identically zero). As in § 4.1, if $h_1$ is known, then the problem (4.16) and constitutive relation (4.15) in principle determine both $\boldsymbol{T}_1$ and $\bar {\boldsymbol {u}}_1$ , up to an arbitrary rigid-body motion. The evolution of $h_1$ is then determined from the first-order mass conservation equation (2.7), namely

(4.17) \begin{equation} \frac {\partial h_1}{\partial t}-\frac {2h_1}{3\psi }+\psi ^{3/2}\tilde {\boldsymbol {\nabla }}\cdot \bar {\boldsymbol {u}}_1=0. \end{equation}

We can simplify the problem (4.15)–(4.17) by introducing a scaled Airy stress function $\mathcal {A} (\tilde {\boldsymbol {X}},t )$ such that

(4.18) \begin{equation} \boldsymbol{T}_1=\psi ^{-3/4} \mathfrak {H}^{\text {c}}[\mathcal {A}] =\psi ^{-3/4}{\begin{pmatrix} \frac {\partial ^2\mathcal {A}}{\partial Y^2} & -\frac {\partial ^2\mathcal {A}}{\partial X\partial Y} \\[4pt] -\frac {\partial ^2\mathcal {A}}{\partial X\partial Y} & \frac {\partial ^2\mathcal {A}}{\partial X^2} \end{pmatrix}}, \end{equation}

which satisfies (4.16a ) identically. Here, we have introduced the notation $\mathfrak {H}[\cdot ]$ for the two-dimensional Hessian matrix and $\mathfrak {H}^{\text {c}}$ for the corresponding cofactor matrix. By eliminating $\bar {\boldsymbol {u}}_1$ from (4.15), we find that $\mathcal {A}$ satisfies the forced biharmonic equation

(4.19) \begin{equation} \tilde {\nabla }^4\mathcal {A}+\psi ^{-1/4}\tilde {\nabla }^2h_1=0, \end{equation}

and the mass-conservation equation (4.17) can be expressed as

(4.20) \begin{equation} 6\frac {\partial h_1}{\partial t}+\psi ^{-3/4}\tilde {\nabla }^2\mathcal {A}=0. \end{equation}

By eliminating $\mathcal {A}$ from (4.19) and (4.20), we find that $h_1$ satisfies

(4.21) \begin{equation} \frac {\partial \tilde {\nabla }^2h_1}{\partial t} =\frac {\dot {\psi }}{4\psi }\,\tilde {\nabla }^2h_1, \end{equation}

and hence

(4.22) \begin{equation} \tilde {\nabla }^2h_1\left (\tilde {\boldsymbol {X}},t\right ) =\psi (t)^{1/4}\tilde {\nabla }^2h_1\left (\tilde {\boldsymbol {X}},0\right ). \end{equation}

The first-order tension in the sheet is thus given by (4.18), where $\mathcal {A}$ satisfies

(4.23a) \begin{equation} \tilde {\nabla }^4\mathcal {A}+\tilde {\nabla }^2h_1\left (\tilde {\boldsymbol {X}},0\right )=0 \quad \text {in }{\Omega _X} \end{equation}

and (from (4.16b ))

(4.23b) \begin{equation} \mathcal {A}=\frac {\partial \mathcal {A}}{\partial n}=0 \quad \text {on }\partial {\Omega _X}. \end{equation}

Since $\Omega _X$ is fixed with respect to the Lagrangian variables $\tilde {\boldsymbol {X}}$ , the scaled stress function $\mathcal {A}$ is independent of $t$ , and determined once and for all by the boundary-value problem (4.23). Thus, the spatial form of the stress field (4.18) is likewise fixed, and it simply scales with $\psi (t)^{-3/4}$ as time increases. The evolution of the thickness perturbations is then given by

(4.24) \begin{equation} h_1\left (\tilde {\boldsymbol {X}},t\right ) =\left (1-\psi (t)^{1/4}\right )\tilde {\nabla }^2\mathcal {A}\left (\tilde {\boldsymbol {X}}\right )+h_1\left (\tilde {\boldsymbol {X}},0\right ). \end{equation}

Note that the mass constraint (4.12) on the initial thickness perturbation holds for all time, i.e.

(4.25) \begin{equation} \iint _{{\Omega _X}} h_1\left (\tilde {\boldsymbol {X}},t\right )\,\mathrm {d}\tilde {\boldsymbol {X}}=0 \end{equation}

for all $t$ .

4.3. Evolution of the centre-surface

For consistency with (4.13), we find that the bending moment tensor scales with

(4.26) \begin{equation} \boldsymbol{M}\left (\tilde {\boldsymbol {X}},t\right )\sim \epsilon ^2\delta \boldsymbol{M}_1\left (\tilde {\boldsymbol {X}},t\right ), \end{equation}

where

(4.27) \begin{equation} \boldsymbol{M}_1=-\frac {\psi ^4}{6}\frac {\partial }{\partial t}\left (\mathfrak {H}[H_1]+(\tilde {\nabla }^2H_1)\tilde {\boldsymbol {I}}\right ) -\frac {\psi ^3}{18}\left (4\mathfrak {H}[H_1]+(\tilde {\nabla }^2H_1)\tilde {\boldsymbol {I}}\right ). \end{equation}

By using (2.11) to eliminate $\boldsymbol{N}$ from (2.13), we thus obtain the moment balance equation in the form

(4.28) \begin{equation} \frac {\psi ^{15/4}}{3}\left (\psi \frac {\partial \tilde {\nabla }^4H_1}{\partial t}+\frac {5}{6}\tilde {\nabla }^4H_1\right )=\mathfrak {H}^{\text {c}}[\mathcal {A}]:\mathfrak {H}[H_1]. \end{equation}

We can slightly simplify this equation by defining the function

(4.29) \begin{equation} J\left (\tilde {\boldsymbol {X}},t\right )=\psi (t)^{5/4}H_1\left (\tilde {\boldsymbol {X}},t\right ), \end{equation}

which satisfies

(4.30) \begin{equation} \frac {\partial \tilde {\nabla }^4J}{\partial t}=3\psi (t)^{-19/4}\mathfrak {H}^{\text {c}}[\mathcal {A}]:\mathfrak {H}[J]. \end{equation}

The effective boundary conditions (3.13) may also be expressed in terms of $J$ in the forms

(4.31a) \begin{align} \frac {\partial }{\partial t}\left ( \frac {\partial ^2J}{\partial n^2}+\tilde {\nabla }^2J\right ) +\frac {1}{2\psi (t)}\left (\frac {\partial ^2J}{\partial n^2}-\tilde {\nabla }^2J\right )&=0 &\text {on }&\partial {\Omega _X}, \end{align}

(4.31b) \begin{align} \frac {\partial }{\partial t}\left ( \frac {\partial ^3J}{\partial n^3}-3\frac {\partial \tilde {\nabla }^2J}{\partial n}+3\kappa _0\tilde {\nabla }^2J\right ) +\frac {1}{2\psi (t)}\left (\frac {\partial ^3J}{\partial n^3}-\frac {\partial \tilde {\nabla }^2J}{\partial n}-\kappa _0\tilde {\nabla }^2J\right )&=0 &\text {on }&\partial {\Omega _X}. \end{align}

We emphasise that these boundary conditions are again expressed in the Lagrangian frame, in which $\Omega _X$ is a fixed domain, with known boundary $\partial\Omega_{X}$ , whose curvature $\kappa _0(\tilde {\boldsymbol {X}})$ is thus independent of time. The curvature $\kappa$ in the Eulerian domain can be recovered using $\kappa (\tilde {\boldsymbol {x}},t)=\sqrt {\psi (t)}\kappa _0(\tilde {\boldsymbol {x}}\sqrt {\psi (t)})$ .

To summarise, given the initial thickness perturbation $h_1(\tilde {\boldsymbol {X}},0)$ , the scaled Airy stress function $\mathcal {A}(\tilde {\boldsymbol {X}})$ is fully determined by the boundary-value problem (4.23). The evolution of the sheet centre-surface is then governed by the partial differential equation (4.30), subject to the boundary conditions (4.31) and the initial condition

(4.32) \begin{equation} J(\tilde {\boldsymbol {X}},0)=H_1(\tilde {\boldsymbol {X}},0). \end{equation}

Of particular interest is whether certain choices of initial data $h_1(\tilde {\boldsymbol {X}},0)$ and $H_1(\tilde {\boldsymbol {X}},0)$ can give rise to temporal growth in the centre-surface displacement $H_1(\tilde {\boldsymbol {X}},t)$ .

From (4.18), we see that the sum of the principal stresses is given by $\operatorname {Tr}(\boldsymbol{T}_1)=\psi ^{-3/4}\tilde {\nabla }^2\mathcal {A}$ , and the boundary conditions (4.23b ) thus imply that

(4.33) \begin{equation} \iint _{{\Omega _X}} \operatorname {Tr}(\boldsymbol{T}_1)\,\mathrm {d}\tilde {\boldsymbol {X}}=0. \end{equation}

It follows that, except for the trivial case where $\tilde {\nabla }^2h_1(\tilde {\boldsymbol {X}},0)=0$ and so $\boldsymbol{T}_1$ is identically zero, there must be a subset of $\Omega _X$ in which $\operatorname {Tr}(\boldsymbol{T}_1)\lt 0$ , i.e. where at least one of the principal stresses is negative and the sheet is thus locally under compression. In the next section, we will show that these compressive zones can indeed give rise to transient growth in the sheet centre-surface by focusing on the relatively simple special case where $\Omega _X$ is a disc.

4.4. Model for a retracting viscous disc

Now let us apply the general theory developed thus far to the particular case where $\Omega _X$ is a disc subject to axisymmetric thickness perturbations. The disc is defined by $0\leq \zeta \lt 1$ , where $\zeta$ is the radial Lagrangian variable, related to the usual plane polar variable $r$ by $\zeta =r\sqrt {\psi (t)}$ . The sheet thickness perturbations are given by $h_1(\zeta ,t)$ , for which the net mass conservation condition (4.25) reduces to

(4.34) \begin{equation} \int _0^1\zeta h_1(\zeta ,t)\,\mathrm {d}\zeta =0. \end{equation}

Given this constraint, we measure the size of the thickness perturbations using a scalar amplitude $A$ , defined by

(4.35) \begin{equation} A=\left [\int _0^1\zeta h_1(\zeta ,{0})^2\,\mathrm {d}\zeta \right ]^{1/2}. \end{equation}

From (4.22) with the assumption of axisymmetry, we have

(4.36) \begin{equation} \frac {1}{\zeta }\,\frac {\mathrm {d}}{\mathrm {d}\zeta } \left (\zeta \frac {\mathrm {d}}{\mathrm {d}\zeta }\right ) \bigl [h_1(\zeta ,t)-\psi (t)^{1/4}h_1(\zeta ,0)\bigr ]=0. \end{equation}

Imposing boundedness at the origin and the mass constraint (4.34), we deduce that

(4.37) \begin{equation} h_1(\zeta ,t)=\psi (t)^{1/4}h_1(\zeta ,0). \end{equation}

Similarly, (4.23a ) can be integrated directly in this case to give

(4.38) \begin{equation} \tilde {\nabla }^2\mathcal {A}= \frac {1}{\zeta }\,\frac {\mathrm {d}}{\mathrm {d}\zeta } \left (\zeta \frac {\mathrm {d}\mathcal {A}}{\mathrm {d}\zeta }\right ) =-h_1\left (\zeta ,0\right ). \end{equation}

The in-plane tension is given by

(4.39) \begin{align} \boldsymbol{T}_1(\zeta ,t)=\operatorname {diag}\left [T_{1rr},T_{1\theta \theta }\right ] =\psi (t)^{-3/4}\operatorname {diag}\left [\frac {1}{\zeta }\,\frac {\mathrm {d}\mathcal {A}}{\mathrm {d}\zeta },\frac {\mathrm {d}^2\mathcal {A}}{\mathrm {d}\zeta ^2}\right ]. \end{align}

By integrating (4.38), we thus obtain

(4.40a) \begin{align} T_{1rr} &= -A\psi (t)^{-3/4}\,\frac {F(\zeta )}{\zeta ^2}, \end{align}

(4.40b) \begin{align} T_{1\theta \theta } &= A\psi (t)^{-3/4} \,\frac {F(\zeta )-\zeta F'(\zeta )}{\zeta ^2}, \end{align}

where we have defined the function $F$ such that

(4.41) \begin{equation} AF(\zeta )=\int _0^\zeta sh_1(s,0)\,\mathrm {d}s. \end{equation}

By including the factor $A$ in (4.41), we ensure that $F$ satisfies the normalisation condition

(4.42) \begin{equation} \int _0^1\frac {F'(\zeta )^2}{\zeta }\,\mathrm {d}\zeta =1, \end{equation}

along with the boundary conditions

(4.43) \begin{equation} F(0)=F'(0)=F(1)=0. \end{equation}

Otherwise, $F$ may be chosen freely by varying the initial thickness perturbation $h_1(\zeta ,0)$ .

It follows from (4.40) that

(4.44) \begin{equation} T_{1rr}+T_{1\theta \theta }=-A\psi (t)^{-3/4}\,\frac {F'(\zeta )}{\zeta } =-\psi (t)^{-3/4}h_1(\zeta ,0) \end{equation}

and hence, as pointed out in § 4.3, for any non-trivial initial centre-surface perturbation, there must always be regions of the disc where $T_{1rr}+T_{1\theta \theta }\lt 0$ so the sheet is locally under compression.

Although we have restricted to axisymmetric thickness perturbations, it is possible for the azimuthal tension $T_{1\theta \theta }$ to be negative. We therefore make no such restriction to the sheet centre-surface displacement, which may well be unstable to non-axisymmetric perturbations. As the problem for $H_1$ is linear, we can write the solution as a sum over azimuthal modes, that is,

(4.45) \begin{equation} H_1(\zeta ,\theta ,t)=\psi (t)^{-5/4}J(\zeta ,\theta ,t)=b(t) + c(t) \zeta \mathrm {e}^{\mathrm {i}\theta }+\psi (t)^{-5/4}\sum _{m=0}^\infty J^{(m)}(\zeta ,t)\mathrm {e}^{\mathrm {i}m\theta } \end{equation}

(real part assumed). The two scalars $b$ and $c$ are included to account for arbitrary rigid-body motions. They are chosen such that

(4.46a) \begin{align} \int _0^{2\pi }\int _0^1H_1(\zeta ,\theta ,t)\zeta \,\mathrm {d}\zeta \mathrm {d}\theta &=0, \end{align}

(4.46b) \begin{align} \int _0^{2\pi }\int _0^1H_1(\zeta ,\theta ,t)\mathrm {e}^{-\mathrm {i}\theta }\zeta ^2\,\mathrm {d}\zeta \mathrm {d}\theta &=0, \end{align}

which eliminate the net transverse displacement and rotation of the sheet, respectively. We assume that the coordinates are oriented such that the constraints (4.46) are satisfied at $t=0$ .

The centre-surface equation (4.30) becomes

(4.47) \begin{equation} \frac {\partial \triangle _m^2J^{(m)}}{\partial t} +3A\psi (t)^{-19/4}\left \{ \frac {1}{\zeta }\,\frac {\partial }{\partial \zeta }\left (\frac {F(\zeta )}{\zeta }\,\frac {\partial J^{(m)}}{\partial \zeta }\right ) -\frac {m^2}{\zeta ^2}\,\frac {\mathrm {d}}{\mathrm {d}\zeta }\left (\frac {F(\zeta )}{\zeta }\right )J^{(m)} \right \}=0, \end{equation}

where

(4.48) \begin{equation} \triangle _m:=\frac {\partial ^2}{\partial \zeta ^2}+\frac {1}{\zeta }\,\frac {\partial }{\partial \zeta }-\frac {m^2}{\zeta ^2} \end{equation}

is the Laplace operator for mode $m$ . The operator $\triangle _m$ is of Cauchy–Euler form and singular at $\zeta =0$ , and the appropriate conditions to impose on $J^{(m)}(\zeta ,t)$ as $\zeta \rightarrow 0$ depend somewhat on the value of $m$ . For $m\gt 2$ , bounded solutions for $J^{(m)}(\zeta ,t)$ are proportional to $\zeta ^m$ or $\zeta ^{m+2}$ as $\zeta \rightarrow 0$ . For $m=2$ , the value of $J^{(2)}(\zeta ,0)$ must be set to zero to ensure that $\triangle _2J^{(2)}$ is bounded. For $m=1$ , the value of $\partial J^{(1)}/\partial \zeta (\zeta ,0)$ is indeterminate and, without loss of generality, may be set to zero by choosing $c(t)$ appropriately in (4.45). Similarly, no generality is lost by setting $J^{(0)}(0,t)$ to zero, by adjusting the function $b(t)$ . Thus, for all mode numbers $m$ , we can select a unique solution for $J{(m)}$ by imposing the boundary conditions

(4.49) \begin{align} J^{(m)}(0,t)=\frac {\partial J^{(m)}}{\partial \zeta }(0,t)=0. \end{align}

The parameters $b$ and $c$ are then given by

(4.50a) \begin{align} H_1(0,\theta ,t)=b(t)&=-2\psi (t)^{-5/4}\int _0^1J^{(0)}(\zeta ,t)\zeta \,\mathrm {d}\zeta , \end{align}

(4.50b) \begin{align} \mathrm {e}^{-\mathrm {i}\theta }\frac {\partial H_1}{\partial \zeta } (0,\theta ,t)=c(t)&=-3\psi (t)^{-5/4}\int _0^1J^{(1)}(\zeta ,t)\zeta ^2\,\mathrm {d}\zeta . \end{align}

The boundary conditions (4.31) at the disc edge are transformed to

(4.51a) \begin{align} \frac {\partial }{\partial t}\left [2\frac {\partial ^2 J^{(m)}}{\partial \zeta ^2}+\frac {\partial J^{(m)}}{\partial \zeta }-m^2J^{(m)}\right ]+\frac {1}{2\psi (t)}\left (m^2J^{(m)}-\frac {\partial J^{(m)}}{\partial \zeta }\right )&=0, \end{align}

(4.51b) \begin{align} \frac {\partial }{\partial t} \left [ 2\frac {\partial ^3 J^{(m)}}{\partial \zeta ^3} - 3(m^2 + 1)\frac {\partial J^{(m)}}{\partial \zeta } + 6m^2J^{(m)} \right ] + \frac {1}{2\psi (t)}\left (1 - m^2\right )\frac {\partial J^{(m)}}{\partial \zeta }&=0\\ & \text {at }\zeta =1, \nonumber \end{align}

and the initial condition for $J^{(m)}$ is given by

(4.52) \begin{equation} J^{(m)}(\zeta ,0)=\frac {1}{2\pi }\int _0^{2\pi }H_1(\zeta ,\theta ,0)\text {e}^{-\mathrm {i} m \theta }\,\textrm {d} \theta . \end{equation}

6.1. Axisymmetric eigenvalue problem

We have seen in § 5 that it is typical for the centre-surface to grow transiently, then decay for large time. We also see that certain modes can be selected, with either $m=0$ or $m=2$ appearing to be dominant for most parameter values. We now show that this behaviour can be quantified by making some approximations to the boundary conditions (4.51) at the edge of the disc. For simplicity, we begin by considering an axisymmetric centre-surface, before generalising to a non-axisymmetric centre-surface to understand the mode selection.

Seeking a separable solution to the axisymmetric centre-surface equation (5.1), we make the ansatz

(6.1) \begin{equation} H_1(\zeta ,t)=\psi (t)^{-5/4} J^{(0)}(\zeta ,t) =\psi (t)^{-5/4}\exp {\left [\frac {6A}{5\lambda }\left (\psi (t)^{-15/4}-1\right )\right ]}g(\zeta ), \end{equation}

where $\lambda$ is an eigenvalue. Then the axisymmetric centre-surface equation and boundary conditions (5.1) and (5.2) become

(6.2) \begin{equation} \zeta g^{\prime \prime \prime }(\zeta )+g^{\prime \prime }(\zeta )-\frac {1}{\zeta }\,g^\prime (\zeta )=\lambda \,\frac {F(\zeta )}{\zeta }\,g^\prime (\zeta ), \end{equation}

and

(6.3) \begin{align} g(0)=g^\prime (0)&=0, \end{align}

(6.4) \begin{align} 2g^{\prime \prime }(1)+g^\prime (1)+\frac {\lambda }{6A}\psi (t)^{15/4}g^\prime (1)&=0. \end{align}

We see that, due to the final term in (6.4), the problem does not accept a fully separable solution. However, in the limit of large thickness perturbations where $A\gg 1$ , the boundary condition (6.4) may be approximated by

(6.5) \begin{equation} 2 g^{\prime \prime }(1)+g^\prime (1)=0. \end{equation}

This approximation breaks down for large times where $t=\textit {O} (A^{4/15} )$ , but allows us to capture the early dynamics where buckling may occur, even if it is transient. For given $F(\zeta )$ , the eigenvalue problem (6.2) with boundary conditions (6.3) and (6.5) may be solved numerically by shooting, with asymptotic behaviour $g(\zeta )\sim \zeta ^2$ as $\zeta \to 0$ and $\lambda$ determined as a shooting parameter by imposing the boundary condition (6.5).

Given $F(\zeta )$ (satisfying the conditions (4.42) and (4.43)), (6.2), along with the boundary conditions (6.3) and (6.5), constitutes an eigenvalue problem for $g$ and $\lambda$ . The eigenfunctions, $g_k$ , satisfy an orthogonality condition, given by

(6.6) \begin{equation} \langle g_j,g_k \rangle =\int _0^1\frac {F(\zeta )}{\zeta }\,g_j^\prime (\zeta ) g_k^\prime (\zeta )\,\textrm {d} \zeta =0\quad \text {for }j\neq k. \end{equation}

(We note that $F$ need not be positive on $(0,1)$ , in which case, $\langle \cdot ,\cdot \rangle$ does not formally define an inner product.) Having computed all the eigenvalues $\lambda _k$ and eigenfunctions $g_k$ , we can reconstruct the solution for the centre-surface as an eigenfunction expansion, namely

(6.7) \begin{equation} H_1(\zeta ,t)=\psi (t)^{-5/4}\sum _{k} \frac {\langle H_1(\zeta ,0),g_k\rangle }{\langle g_k,g_k \rangle }\exp \left [\frac {6A}{5\lambda _k}\left (\psi (t)^{-15/4}-1\right )\right ]g_k(\zeta ). \end{equation}

We check the validity of using the approximate boundary condition (6.5) instead of (6.4) with a thickness perturbation $h_1(\zeta ,0)=10\sin (2 \pi \zeta )/\zeta$ , corresponding to $A\approx 12.48$ , and the initial centre-surface given by $H_1(\zeta ,0)=\zeta ^2(15-6\zeta )/9$ . In figure 11, we show the evolution of the centre-surface displacement $H_1(1,t)$ at the edge of the disc, predicted by the full numerical solution described in § 5 and by the approximate solution (6.7). We see that there is very good agreement in the early-time behaviour and good qualitative agreement between the two solutions for all times, with the eigenfunction expansion (6.7) capturing well the growth and decay of the full solution. Furthermore, we will demonstrate below that the approximate solution (6.7) provides good estimates of both the duration and the amplitude of the transient growth.

Figure 11. The evolution of the centre-surface displacement at the edge of the disc, $H_1(1,t)$ , calculated from the full centre-surface boundary-value problem (5.1) and (5.2) in red, and via the eigenvalue approximation (6.7) in black. The initial thickness and centre-surface perturbations are given by $h_1(\zeta ,0)=10\sin (2 \pi \zeta )/\zeta$ and $H_1(\zeta ,0)=\zeta ^2(15-6\zeta )/9$ .

6.2. Quantifying the centre-surface deviation

It is difficult to make much analytical progress with the full expansion (6.7), so let us consider for now the case where the centre-surface perturbation is exactly an eigenfunction. Then, we only have a contribution from one term in the series, say

(6.8) \begin{equation} H_1(\zeta ,t)=c\psi (t)^{-5/4} \exp {\left [\frac {6A}{5 \lambda }\left (\psi (t)^{-15/4}-1\right )\right ]}g(\zeta ). \end{equation}

In this instance, assuming we have again normalised $H_1$ such that $d(0)=1$ , we explicitly calculate the maximum difference (5.9) to be given by

(6.9) \begin{equation} \textrm{d}(t)=\psi (t)^{-5/4}\exp \left [\frac {6A}{5\lambda }\left (\psi (t)^{-15/4}-1\right )\right ]. \end{equation}

It is thus possible for the solution to grow only if $\lambda$ is negative. However, by taking the inner product of the eigenvalue equation (6.2) with $g'$ , we find that the eigenvalues are given by

(6.10) \begin{equation} \frac {\lambda }{A}=\frac {g^\prime (1)^2/2+\int _0^1\left [\zeta {g^{\prime \prime }}(\zeta )^2+{g^\prime }(\zeta )^2/\zeta \right ]\,\textrm {d} \zeta }{\int _0^1\zeta T_{1rr}(\zeta ,0){g^\prime }(\zeta )^2\,\textrm {d} \zeta }. \end{equation}

Here, the numerator is non-negative and the denominator depends on the initial radial tension. As we would expect (e.g. Filippov & Zheng Reference Filippov and Zheng2010), if the radial tension is positive everywhere, then all of the eigenvalues $\lambda$ are positive and transient growth is impossible. However, if the radial tension is negative everywhere, then the eigenvalues are negative and transient buckling is possible; if $T_{1rr}$ changes sign, then we can have both positive and negative eigenvalues.

Assuming that $\lambda$ is negative, we find that the stationary point of $d^\prime (t)=0$ occurs at $t=t_*$ , where

(6.11) \begin{equation} \psi (t_*)=\frac {2t_*}{3}+1=\left (\frac {-18A}{5 \lambda }\right )^{4/15}. \end{equation}

To have $d(t)$ initially increasing, we need

(6.12) \begin{equation} A\gt -\frac {5\lambda }{18}\gt 0, \end{equation}

i.e. we need both for the problem (6.2)–(6.5) to admit a negative eigenvalue $\lambda$ and for $A$ to be sufficiently large. As seen by Ryan et al. (Reference Ryan, Breward, Griffiths and Howell2024), there is a threshold for the amplitude of the thickness perturbation, above which there is transient buckling and below which there is not. At the stationary point $t=t_*$ , we calculate the maximum centre-surface deformation

(6.13) \begin{equation} d_*=d(t_*)=\left (\frac {-5\lambda }{18A}\right )^{1/3}\exp \left (-\frac {1}{3}-\frac {6A}{5\lambda }\right ). \end{equation}

We infer that thickness perturbations of amplitude $\epsilon ^2A$ where $A=O\bigl (1/\log {(1/\delta )\bigr )}$ can cause $H_1$ to grow by an order of magnitude in $\delta$ , thus invalidating the neglect of nonlinear terms in § 4.2. We note also that, in the full eigenfunction expansion (6.7), the term corresponding to the largest negative eigenvalue $\lambda _*$ (i.e. the negative eigenvalue of smallest amplitude) will dominate the solution when $t\sim t_*$ , so we can continue to use the approximations (6.11) and (6.13) for general centre-surface profiles comprising a mix of eigenfunctions. We thus predict that the maximal time and centre-surface deformation amplitude should satisfy $\psi (t_*)=O (A^{4/15} )$ and $\textrm{d}(t_*)=O (A^{-1/3}\mathrm {e}^{kA} )$ as $A\rightarrow \infty$ , for some constant $k$ .

We now compare these predicted relationships to numerical results calculated using the full centre-surface equation (5.1) and boundary conditions (5.2). We take the initial centre-surface profile $H_1(\zeta ,0)=\zeta ^2$ and Gaussian thickness perturbation (5.7). In figure 12(a), we observe, as expected, a threshold value of $A$ for transient growth, above which there is a clear $4/15$ power law. We see that there is excellent agreement between the predicted relationships, (6.11) and (6.13), and the numerical solution (figure 12) once the threshold for transient growth has been reached. Moreover, the asymptotically straight curves seen using log-linear axes in figure 12(b) are consistent with the predicted exponential dependence of $\textrm{d}_*$ on $A$ .

We note that the maximal time $t_*\sim A^{4/15}$ occurs precisely when the approximation (6.5) breaks down. Nevertheless, we conclude from the excellent agreement observed in figure 12 that the asymptotic predictions (6.11) and (6.13) correctly capture the power-law behaviour for large $A$ , though not necessarily the prefactors.

Figure 12. Plot of (a) $\psi (t_*)$ and (b) $\textrm{d}(t_*)A^{1/3}$ versus the thickness perturbation amplitude $A$ , as calculated from the full boundary-value problem (5.1) and (5.2). We use a thickness perturbation given by (5.7), with $\mu _h=0.3,0.4,0.5$ , and initial centre-surface displacement $H_1(\zeta ,0)=\zeta ^2$ .

6.3. Maximising axisymmetric buckling

We recall from figure 8 that the magnitude of the transient growth strongly depends on both the initial centre-surface and the thickness profiles. Now we pose the question of which combination of thickness and centre-surface perturbations gives rise to the largest transient growth. The above analysis suggests the following related problem: which function $F(\zeta )$ , satisfying the normalisation condition (4.42) and boundary conditions (4.43), gives rise to the smallest possible (in magnitude) negative eigenvalue $\lambda _*$ of the problem (6.2)–(6.5)? We then maximise over centre-surface perturbations by choosing $H_1(\zeta ,0)$ to be proportional to the eigenfunction $g_*(\zeta )$ corresponding to the extremal eigenvalue $\lambda _*$ .

Mathematically, our problem is then to find

(6.14) \begin{equation} \lambda _*=\min _{F(\zeta )}\left \{|\lambda |:\lambda \lt 0\right \}, \end{equation}

subject to $F$ satisfying the the constraint (4.42) and boundary conditions (4.43), and $\{g,\lambda \}$ solving the eigenvalue problem (6.2), (6.3) and (6.5). We perturb around the extremal solutions by setting $g^\prime \mapsto g_*^\prime +\chi$ , $F\mapsto F_*+\phi$ , while $\lambda =\lambda _*$ remains stationary. Then, substituting into (6.2), (6.3) and (6.5), we get

(6.15a) \begin{gather} \left (\zeta \chi (\zeta )^\prime \right )^\prime -\frac {1+\lambda _* F_*(\zeta )}{\zeta }\,\chi (\zeta )=\lambda _*\,\frac {g_*^\prime (\zeta ) \phi (\zeta )}{\zeta }, \end{gather}

(6.15b) \begin{gather} \chi (0)=2\chi ^\prime (1)+\chi (1)=0. \end{gather}

This problem for $\chi$ is self-adjoint, with the homogeneous problem satisfied by $g_*^\prime (\zeta )$ . By the Fredholm Alternative Theorem, we obtain the solvability condition

(6.16) \begin{equation} \int _0^1\frac {g_*^\prime (\zeta )^2}{\zeta }\,\phi (\zeta )\,\textrm {d} \zeta =0. \end{equation}

Meanwhile, by perturbing the conditions (4.42) and (4.43) on $F$ , we find

(6.17a) \begin{gather} \int _0^1\frac {\mathrm {d}}{\mathrm {d}\zeta }\left (\frac {F_*'(\zeta )}{\zeta }\right )\phi (\zeta )\,\textrm {d} \zeta =0, \end{gather}

(6.17b) \begin{gather} \phi (0)=\phi '(0)=\phi (1)=0. \end{gather}

From (6.15) and (6.17), we deduce that the extremal functions $g_*$ and $F_*$ satisfy the boundary-value problems

(6.18a) \begin{align} g_*^{\prime \prime \prime }(\zeta )+\frac {g_*^{\prime \prime }(\zeta )}{\zeta }-\frac {1+\lambda _* F_*(\zeta )}{\zeta ^2}\,g_*^\prime (\zeta ) & = 0, \end{align}

(6.18b) \begin{align} g_*(0)=g_*^\prime (0)=2g_*^{\prime \prime }(1)+g_*^\prime (1)&=0, \end{align}

(6.18c) \begin{align} F_*^{\prime \prime }(\zeta )-\frac {F_*^\prime (\zeta )}{\zeta }-\mu g_*^\prime (\zeta )^2 & = 0, \end{align}

(6.18d) \begin{align} F_*(0)=F_*'(0)=F_*(1)&=0. \end{align}

The extremal eigenvalue $\lambda _*$ is determined as part of the solution, while the additional eigenvalue $\mu$ is associated with the constraint (4.42) and may be set to $\pm 1$ by scaling $g_*$ appropriately. We solve the problem (6.18) by shooting from $\zeta =0$ , with the asymptotic behaviour $g_*(\zeta )\sim \zeta ^2$ and $F_*(\zeta )\sim c\zeta ^2$ as $\zeta \rightarrow 0$ , where $c$ and $\lambda$ are determined as shooting parameters by imposing the boundary conditions at $\zeta =1$ .

To validate the results of the above approach, we also calculate the extremal kernel function $F_*$ and the corresponding extremal eigenvalue $\lambda _*$ and eigenfunction $g_*$ numerically using the Rayleigh–Ritz method (see, for example, Collins Reference Collins2006). We write (6.10) in the form $\lambda =I[g]/K[g]$ , where

(6.19a,b) \begin{align} I[g]&=\frac {g^{\prime \prime }(1)^2}{2}+\int _0^1\left (\zeta g^{\prime \prime }(\zeta )^2+\frac {g^\prime (\zeta )^2}{\zeta }\right )\,\textrm {d} \zeta , & K[g]&=\int _0^1\frac {F(\zeta )g^\prime (\zeta )^2}{\zeta }\,\textrm {d} \zeta . \end{align}

We approximate $g(\zeta )$ and $F(\zeta )$ by truncated power series in $\zeta$ , with the coefficients chosen to satisfy the boundary conditions (6.18b ) and (6.18d ), as well as the normalisation conditions (4.42) and $K[g]=1$ . The remaining coefficients are then varied to minimise $I[g]$ .

For this exercise, we fix three degrees of freedom (DoFs) in $g$ (which is therefore approximated by a polynomial of degree 6) while taking one, two or three DoFs in $F$ (which is approximated by a polynomial of degree 4, 5 or 6). The approximate values thus obtained for the smallest negative eigenvalue are given in table 1. We see that this sequence of eigenvalues approaches a limit as the number of DoFs is increased, and that the limiting value agrees with the value of $\lambda _*$ computed from the “optimal” boundary-value problem (6.18). This extremal value of $\lambda$ tells us about the absolute maximum axisymmetric transient growth that can be observed for a given (large) perturbation amplitude $A$ .

Table 1. Value of the smallest negative eigenvalue $\lambda _*$ , computed using the Rayleigh–Ritz method with three degrees of freedom (DoFs) in $g$ and varying DoFs in $F$ . The ‘optimal’ value is obtained by solving the boundary-value problem (6.18).

We plot the calculated thickness perturbation profiles in figure 13(a) and indeed see that three DoFs in both $g$ and $F$ are sufficient to give an excellent polynomial approximation to the thickness perturbation that maximises axisymmetric transient growth. This extremal perturbation corresponds to the sheet being slightly thicker at the centre and thinner towards the edge, and indeed these kinds of perturbations were also found to promote axisymmetric buckling in the numerical experiments performed in § 5. The corresponding optimal initial centre-surface displacement is shown in figure 13(b). This characteristic bowl-like shape is very similar to the maximal axisymmetric displacement shown in figure 9(c), illustrating again how the results of this section can help us to understand what kinds of centre-surface profiles are likely to be selected by the dynamics.

Figure 13. (a) Plot of the extremal thickness perturbation, $h_1(\zeta ,0)=F_*'(\zeta )/\zeta$ , versus $\zeta$ . The solid curves are obtained using the Rayleigh–Ritz approximation with three degrees of freedom (DoFs) in $g$ and varying DoFs in $F$ . The dashed curve is the ‘optimal’ perturbation, given by the solution of (6.18). (b) Plot of the optimal initial centre-surface profile, $H_1(\zeta ,0)=g_*(\zeta )$ versus $\zeta$ .

6.4. Non-axisymmetric eigenvalue problem

In the limit of large $A$ , the dynamics can be approximately described by an eigenvalue problem also in the non-axisymmetric case. Now, when we make the ansatz

(6.20) \begin{equation} H_1(\zeta ,\theta ,t)=\psi (t)^{-5/4}J^{(m)}(\zeta ,t)\text {e}^{\mathrm {i}m\theta }=\psi (t)^{-5/4}\exp {\left [\frac {6A}{5\lambda ^{(m)}}\left (\psi (t)^{-15/4}-1\right )\right ]}g(\zeta )\text {e}^{\mathrm {i}m\theta }, \end{equation}

the centre-surface equation (4.47) is transformed to

(6.21) \begin{equation} \triangle _m^2\,g(\zeta ) =\lambda ^{(m)}\left [ \frac {1}{\zeta }\,\frac {\partial }{\partial \zeta }\left (\frac {F(\zeta )}{\zeta }\,\frac {\partial J^{(m)}}{\partial \zeta }\right ) -\frac {m^2}{\zeta ^2}\,\frac {\mathrm {d}}{\mathrm {d}\zeta }\left (\frac {F(\zeta )}{\zeta }\right )J^{(m)} \right ], \end{equation}

with the boundary conditions

(6.22a) \begin{align} g(0)&=0, \end{align}

(6.22b) \begin{align} g^\prime (0)&=0, \end{align}

(6.22c) \begin{align} 2g^{\prime \prime }(1)+g^\prime (1)-m^2g(1)&=0, \end{align}

(6.22d) \begin{align} 2g^{\prime \prime \prime }(1)-3(m^2+1)g^\prime (1)+6m^2g(1)&=0. \end{align}

As in § 6.1, terms of order $\psi (t)^{15/4}/A$ have been neglected in the boundary conditions (6.22c ) and (6.22d ), so this approximation breaks down for sufficiently large $t$ .

By taking the inner product of (6.21) with $g$ , we find that the eigenvalue $\lambda ^{(m)}$ can be expressed as

(6.23) \begin{equation} \frac {\lambda ^{(m)}}{A} = \frac {\begin{split}\displaystyle 2\int _0^1\left [\zeta g^{\prime \prime }(\zeta )^2 + \left (1+2m^2\right )\frac {g^\prime (\zeta )^2}{\zeta }+ m^2\left (m^2- 4\right )\frac {g(\zeta )^2}{\zeta ^3}\right ] \,\textrm {d} \zeta\\ +g^\prime (1)^2 - 2m^2g(1)g^\prime (1)-5m^2g(1)^2\end{split}}{\displaystyle 2\int _0^1\left [\zeta T_{1rr}(\zeta ,0)g^\prime (\zeta )^2+m^2\frac {T_{1\theta \theta }(\zeta ,0)g(\zeta )^2}{\zeta }\right ]\,\textrm {d} \zeta }. \end{equation}

In the limit as $m\rightarrow \infty$ , (6.23) becomes

(6.24) \begin{equation} \frac {\lambda ^{(m)}}{A} \sim \frac {\displaystyle m^2\int _0^1g(\zeta )^2/\zeta ^3\,\mathrm {d}\zeta }{\displaystyle \int _0^1 T_{1\theta \theta }(\zeta ,0)g(\zeta )^2/\zeta \,\textrm {d} \zeta }. \end{equation}

Therefore, negative eigenvalues can exist, implying that non-axisymmetric buckling is possible, whenever the hoop tension $T_{1\theta \theta }$ is negative. However, we note that the eigenvalues grow like $m^2$ for large $m$ , so that the magnitude of any transient growth will decrease exponentially for larger mode numbers.

We now use the eigenvalue approximation to explain the results concerning mode selection found in figure 8. Assuming that the behaviour of the centre-surface is dominated by the smallest (in magnitude) negative eigenvalue, it follows that the mode with the largest deformation amplitude, $d_*$ , will be that with the smallest negative eigenvalue. Figure 8 suggests that only modes $m=0,1,2$ can be dominant. Motivated by this observation, we calculate the smallest negative eigenvalue for modes $m=0,1,2$ by solving (6.21)–(6.22) numerically, for the Gaussian thickness perturbation given by (5.7) with varying $\mu _h$ . The results are shown in figure 14, where we see that the axisymmetric mode dominates (i.e. $\lambda ^{(0)}$ is closest to zero) for $0\leq \mu _h \lesssim 0.6$ , while the $m=2$ mode dominates for $\mu _h\gtrsim 0.6$ . The point of intersection at $\mu _h\approx 0.6$ corresponds to the region in the contour plot in figure 8 where the dominant mode switches between $m=0$ and $m=2$ , as $\mu _h$ varies. The locations of the maxima in $\lambda ^{(0)}$ and $\lambda ^{(2)}$ (indicated by dashed lines) are also encouragingly consistent with the values of $\mu _h$ that locally maximise $d_*$ in figure 8. The maximum value of $\lambda ^{(2)}$ is closer to zero than the maximum in $\lambda ^{(0)}$ , which explains why larger values of $d_*$ are attained with $m=2$ than with $m=0$ .

We recall that figure 8 shows small regions of parameter values where the $m=1$ mode dominates, which appears to contradict figure 14. In these regions, the initial centre-surface displacement is approximately orthogonal to the dominant eigenfunction, allowing other subdominant modes to play a role in the dynamics.

Figure 14. A plot of the smallest (in magnitude) negative eigenvalue, $\lambda ^{(m)}$ , satisfying the eigenvalue problem (6.21)–(6.22), where the thickness perturbation is given by (5.7) with $m=0$ (blue), $m=1$ (red) and $m=2$ (black). The local maxima are indicated by dashed lines for $m=0,2$ .

6.5. Maximising non-axisymmetric buckling

We now ask what thickness perturbation leads to the smallest negative eigenvalue in (6.21) for a non-axisymmetric centre-surface. We use a Rayleigh–Ritz approximation, as in § 6.3, to calculate the permissible functions $F$ and $g$ that give the smallest (in magnitude) eigenvalue $\lambda _*^{(m)}$ for each mode number $m$ , using (6.23). The results are presented in figure 15, in which the square root modulus of each extremal eigenvalue is plotted versus $m$ , clearly showing that the eigenvalues grow with $m^2$ for large $m$ , in agreement with (6.24). We see that the closest eigenvalue to zero is $\lambda _*^{(2)}\approx -4.38$ , with $|\lambda _*^{(0)}|$ and $|\lambda _*^{(3)}|$ being the next smallest. There is also the special case $m=1$ , where the minimum eigenvalue is approximately the same as for $m=6$ . We conclude that $m=2$ is the easiest mode to excite, in that it can undergo transient growth at smaller values of the amplitude $A$ than any other mode. The corresponding extremal exigenvalue $\lambda _*^{(2)}\approx -4.38$ gives a bound on the transient growth that can be observed for any initial thickness and centre-surface perturbations. For $m\geq 3$ , we calculate that $0\gt \lambda ^{(2)}\gt \lambda _*^{(m)}$ , meaning that, even for the thickness perturbation that is optimal for a given $m\geq 3$ , the mode $m=2$ will be more dominant. This result explains why $m=0$ and $m=2$ were shown to be dominant in § 5.

Figure 15. The square root modulus of the extremal eigenvalues of (6.21) and (6.22) versus mode number $m$ .