2 Problem formulation and governing equations

Figure 2. Schematic illustration showing the centreline cross-section of the examined configuration. A thin viscous fluid layer of initial thickness $\tilde{h}_{0}$ is confined between a rigid surface and a pre-stretched elastic sheet supported at its boundaries. Non-uniform and time-varying slip velocity $\tilde{\boldsymbol{u}}_{slip}(\tilde{x},{\tilde{y}},\tilde{t})$ and external pressure $\tilde{p}_{e}(\tilde{x},{\tilde{y}},\tilde{t})$ drive the viscous–elastic interaction and create the deformation field $\tilde{d}(\tilde{x},{\tilde{y}},\tilde{t})$ .

We study the viscous–elastic dynamics of a viscous fluid of density $\tilde{\unicode[STIX]{x1D70C}}$ and viscosity $\tilde{\unicode[STIX]{x1D707}}$ confined between a flat rigid surface and a pre-stretched elastic sheet of length $\tilde{l}_{m}$ , width $\tilde{w}_{m}$ , thickness $\tilde{h}_{m}$ , Young’s modulus $\tilde{E}_{Y}$ and Poisson’s ratio $\unicode[STIX]{x1D708}$ , as shown in figure 2. We denote dimensional variables by tildes, normalized variables without tildes and characteristic values by an asterisk superscript. We employ a Cartesian coordinate system $(\tilde{x},{\tilde{y}},\tilde{z})$ , and adopt the $\Vert$ and $\bot$ subscripts to denote parallel and perpendicular directions to the $\tilde{x}-{\tilde{y}}$ plane, respectively.

The fluid velocity is $\tilde{\boldsymbol{u}}=(\tilde{\boldsymbol{u}}_{\Vert },\tilde{u} _{\bot })=(\tilde{u} ,\tilde{v},\tilde{w})$ and fluid pressure is $\tilde{p}$ . The total gap between the plates is $\tilde{h}(\tilde{x},{\tilde{y}},\tilde{t})=\tilde{h}_{0}+\tilde{d}(\tilde{x},{\tilde{y}},\tilde{t})$ , where $\tilde{t}$ is time and $\tilde{h}_{0}$ is the initial gap. The deformation field $\tilde{d}(\tilde{x},{\tilde{y}},\tilde{t})$ is induced either due to an external pressure $\tilde{p}_{e}(\tilde{x},{\tilde{y}},\tilde{t})$ , acting directly on the elastic sheet, or due to an internal pressure formed by a non-uniform slip velocity $\tilde{\boldsymbol{u}}_{slip}(\tilde{x},{\tilde{y}},\tilde{t})$ on the rigid surface, which varies over a characteristic length scale $\tilde{l}^{\ast }$ in the $\tilde{x}-{\tilde{y}}$ plane. The characteristic velocities in the $\tilde{x}-{\tilde{y}}$ plane and $\hat{\boldsymbol{z}}$ direction are respectively $\tilde{u} ^{\ast }$ and $\tilde{w}^{\ast }$ , and the characteristic pressure, deformation and time are respectively denoted as $\tilde{p}^{\ast }$ , $\tilde{d}^{\ast }$ , and $\tilde{t}^{\ast }$ .

The most general description for the dynamics of a thin elastic sheet is given by the nonlinear Föppl–von Kármán equations (Timoshenko & Woinkowsky-Krieger Reference Timoshenko and Woinkowsky-Krieger1987; Howell, Kozyreff & Ockendon Reference Howell, Kozyreff and Ockendon2009), which account for bending and tension forces, external traction, as well as the solid’s inertia. In this work, we assume that the solid’s inertia is negligible and focus on the case of a strongly pre-stretched elastic sheet with isotropic tension $\tilde{T}$ , assumed to be much larger than any internal tension, $\tilde{T}_{in}$ , forming in the system during actuation of the sheet. From scaling of the Föppl–von Kármán equation that couples the spatial variations in internal tension to Gaussian curvature (Howell et al. Reference Howell, Kozyreff and Ockendon2009, p. 176), we find $\tilde{T}_{in}\sim (\tilde{d}^{\ast }/\tilde{l}^{\ast })^{2}\tilde{E}_{Y}\tilde{h}_{m}$ , and define

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\frac{\tilde{T}_{in}}{\tilde{T}}=\left(\frac{\tilde{d}^{\ast }}{\tilde{l}^{\ast }}\right)^{2}\frac{\tilde{E}_{Y}\tilde{h}_{m}}{\tilde{T}}\ll 1. & & \displaystyle\end{eqnarray}$$

Expanding the deformation field and the tension in powers of $\unicode[STIX]{x1D6FC}$ and considering the leading order, the nonlinear Föppl–von Kármán equations reduce to a single linear plate equation containing bending and pre-stretching terms (see Timoshenko & Woinkowsky-Krieger Reference Timoshenko and Woinkowsky-Krieger1987; Howell et al. Reference Howell, Kozyreff and Ockendon2009)

(2.2) $$\begin{eqnarray}\displaystyle \tilde{p}=\tilde{B}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{4}\tilde{d}-\tilde{T}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{2}\tilde{d}+\tilde{p}_{e}, & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D735}}_{\Vert }=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\tilde{x},\unicode[STIX]{x2202}/\unicode[STIX]{x2202}{\tilde{y}})$ is the two-dimensional gradient and $\tilde{B}=\tilde{E}_{Y}\tilde{h}_{m}^{3}/12(1-\unicode[STIX]{x1D708}^{2})$ is the bending stiffness.

We restrict our analysis to shallow geometries, to negligible fluidic inertia, represented by small Womersley and reduced Reynolds numbers, and to small elastic deformations,

(2.3a-d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}=\frac{\tilde{h}_{0}}{\tilde{l}^{\ast }}\ll 1,\quad Wo=\frac{\tilde{\unicode[STIX]{x1D70C}}\tilde{h}_{0}^{2}}{\tilde{\unicode[STIX]{x1D707}}\tilde{t}^{\ast }}\ll 1,\quad \unicode[STIX]{x1D716}Re=\unicode[STIX]{x1D716}\frac{\tilde{\unicode[STIX]{x1D70C}}\tilde{u} ^{\ast }\tilde{h}_{0}}{\tilde{\unicode[STIX]{x1D707}}}\ll 1,\quad \unicode[STIX]{x1D6FD}=\frac{\tilde{d}^{\ast }}{\tilde{h}_{0}}\ll 1. & & \displaystyle\end{eqnarray}$$

The first three assumptions (2.3a–c ) allow us to apply the lubrication approximation and to obtain a nonlinear Reynolds equation (Leal Reference Leal2007), and the last assumption, $\unicode[STIX]{x1D6FD}\ll 1$ , enables its linearization (see e.g. Tulchinsky & Gat Reference Tulchinsky and Gat2016; Kodio et al. Reference Kodio, Griffiths and Vella2017).

In § 6 we relax the assumptions (2.1) and (2.3d ), and explore the effects of nonlinearity of the Reynolds and Föppl–von Kármán equations on the deformation and tension fields, considering an axisymmetric geometry.

Based on (2.3), the fluid motion is governed by the lubrication equations (Leal Reference Leal2007)

(2.4a-c ) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0,\quad \tilde{\unicode[STIX]{x1D735}}_{\Vert }\tilde{p}=\tilde{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}^{2}\tilde{\boldsymbol{u}}_{\Vert }}{\unicode[STIX]{x2202}\tilde{z}^{2}},\quad \frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}\tilde{z}}=0. & & \displaystyle\end{eqnarray}$$

We assume that the fluid is subject to the slip and the no-penetration boundary conditions on the bottom surface. Horizontal motion of the elastic sheet is negligible in the small-deformation limit, implying the no-slip and the no-penetration boundary conditions along it,

(2.5a,b ) $$\begin{eqnarray}\displaystyle (\tilde{\boldsymbol{u}}_{\Vert },\tilde{u} _{\bot })|_{\tilde{z}=0}=(\tilde{\boldsymbol{u}}_{slip}(\tilde{x},{\tilde{y}},\tilde{t}),0),\quad (\tilde{\boldsymbol{u}}_{\Vert },\tilde{u} _{\bot })|_{\tilde{z}=\tilde{h}}=\left(\mathbf{0},\frac{\unicode[STIX]{x2202}\tilde{d}}{\unicode[STIX]{x2202}\tilde{t}}\right). & & \displaystyle\end{eqnarray}$$

In this work, we consider a non-uniform slip velocity and specifically focus on the electro-osmotic slip which arises over electrically charged surfaces due to interaction of an externally applied electric field with the excess of net charge in the electric double layer. In the thin-double-layer limit, such interaction results in bulk fluid motion outside the outer edge of the electric double layer according to the Helmholtz–Smoluchowski equation (Hunter Reference Hunter2000),

(2.6) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{u}}_{slip}=-\frac{\tilde{\unicode[STIX]{x1D700}}\tilde{\unicode[STIX]{x1D701}}(\tilde{x},{\tilde{y}})}{\tilde{\unicode[STIX]{x1D707}}}\tilde{\boldsymbol{E}}, & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D700}}$ is the fluid permittivity, $\tilde{\unicode[STIX]{x1D701}}$ is the zeta potential on the bottom surface and $\tilde{\boldsymbol{E}}$ is the tangential imposed electric field.

2.1 Scaling analysis and non-dimensionalization

Scaling by the characteristic dimensions, we define the following normalized quantities: $(x,y,z)=(\tilde{x}/\tilde{l}^{\ast },{\tilde{y}}/\tilde{l}^{\ast },\tilde{z}/\tilde{h}_{0})$ , $(u,v,w)=(\tilde{u} /\tilde{u} ^{\ast },\tilde{v}/\tilde{u} ^{\ast },\tilde{w}/\tilde{w}^{\ast })$ , $p=\tilde{p}/\tilde{p}^{\ast }$ , $t=\tilde{t}/\tilde{t}^{\ast }$ , $d=\tilde{d}/\tilde{d}^{\ast }$ and $h=\tilde{h}/\tilde{h}_{0}$ . As noted by Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), the bending-tension length scale $\tilde{l}_{BT}=\sqrt{\tilde{B}/\tilde{T}}$ determines the relative importance of bending and tension forces in the elastic response of the sheet; for $\tilde{l}^{\ast }\ll \tilde{l}_{BT}$ bending forces dominate, whereas for $\tilde{l}^{\ast }\gg \tilde{l}_{BT}$ tension forces dominate. A convenient dimensionless number when scaling (2.2) is $\unicode[STIX]{x1D706}$ , defined as

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\equiv \frac{\text{Bending}}{\text{Stretching}}=\left(\frac{\tilde{l}_{BT}}{\tilde{l}^{\ast }}\right)^{2}=\frac{\tilde{B}}{\tilde{T}\tilde{l}^{\ast 2}}=\frac{1}{12(1-\unicode[STIX]{x1D708}^{2})}\frac{\tilde{E}_{Y}\tilde{h}_{m}}{\tilde{T}}\left(\frac{\tilde{h}_{m}}{\tilde{l}^{\ast }}\right)^{2}. & & \displaystyle\end{eqnarray}$$

In this study, our main focus is on deformations which are comparable to, or larger than the sheet thickness, and therefore the appropriate scaling for deformation is based on tension and not on bending

(2.8) $$\begin{eqnarray}\displaystyle \tilde{d}^{\ast }=\frac{\tilde{p}^{\ast }\tilde{l}^{\ast 2}}{\tilde{T}}, & & \displaystyle\end{eqnarray}$$

with which the elastic equation (2.2) takes the form

(2.9) $$\begin{eqnarray}\displaystyle p=\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FB}_{\Vert }^{4}d-\unicode[STIX]{x1D6FB}_{\Vert }^{2}d+p_{e}. & & \displaystyle\end{eqnarray}$$

From order-of-magnitude analysis of the continuity and in-plane momentum equations (2.4a,b ), it follows that the perpendicular velocity and the pressure scale as $\tilde{w}^{\ast }\sim \unicode[STIX]{x1D716}\tilde{u} ^{\ast }$ and $\tilde{p}^{\ast }\sim 12\tilde{\unicode[STIX]{x1D707}}\tilde{u} ^{\ast }/\unicode[STIX]{x1D716}^{2}\tilde{l}^{\ast }$ , respectively. Using the kinematic boundary condition (2.5b ) we obtain the scaling for the viscous–elastic time scale $\tilde{t}^{\ast }\sim \tilde{d}^{\ast }/\tilde{w}^{\ast }\sim \tilde{d}^{\ast }/\unicode[STIX]{x1D716}\tilde{u} ^{\ast }$ . Owing to a linear relation between $\tilde{d}^{\ast }$ and $\tilde{p}^{\ast }$ as well as between $\tilde{p}^{\ast }$ and $\tilde{u} ^{\ast }$ , the ratio $\tilde{d}^{\ast }/\tilde{u} ^{\ast }$ is not dependent on actuation force, nor on $\tilde{\boldsymbol{u}}_{slip}$ or $\tilde{p}_{e}$ . Thus, the viscous–elastic time scale solely depends on the properties of the fluid and the elastic medium, specifically the ratio $\tilde{\unicode[STIX]{x1D707}}/\tilde{T}$ and the geometry:

(2.10) $$\begin{eqnarray}\displaystyle \tilde{t}^{\ast }=\frac{12\tilde{\unicode[STIX]{x1D707}}\tilde{l}^{\ast 4}}{\tilde{T}\tilde{h}_{0}^{3}}=\frac{12\tilde{\unicode[STIX]{x1D707}}\tilde{l}^{\ast }}{\unicode[STIX]{x1D716}^{3}\tilde{T}}. & & \displaystyle\end{eqnarray}$$

Based on this scaling, and following standard lubrication theory, from (2.4), (2.5) and (2.10) we obtain the normalized Reynolds equation (Leal Reference Leal2007, p. 313)

(2.11) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }[h^{3}\unicode[STIX]{x1D735}_{\Vert }p]=-\frac{1}{2}\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }[h\boldsymbol{u}_{slip}], & & \displaystyle\end{eqnarray}$$

where using the definition (2.3d ), the fluid thickness $h$ can be expressed as $h=1+\unicode[STIX]{x1D6FD}d$ .

2.2 Viscous–elastic governing equations for a pre-stretched elastic sheet

Substituting (2.9) into (2.11), yields the nonlinear governing equation

(2.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }[h^{3}\unicode[STIX]{x1D735}_{\Vert }(\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FB}_{\Vert }^{4}d-\unicode[STIX]{x1D6FB}_{\Vert }^{2}d)]=-\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }[h\boldsymbol{f}_{F}]+\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }[h^{3}\unicode[STIX]{x1D735}_{\Vert }f_{E}], & & \displaystyle\end{eqnarray}$$

where we define $\boldsymbol{f}_{F}=f_{Fx}\hat{\boldsymbol{x}}+f_{Fy}\hat{\boldsymbol{y}}=\boldsymbol{u}_{slip}/2$ and $f_{E}=p_{e}$ as the actuation mechanisms. The subscript $F$ refers to driving force applied to the sheet due to the non-uniform slip velocity in the fluid and the subscript $E$ refers to driving force due to the pressure applied directly to the elastic sheet.

Assuming small elastic deformations, $\tilde{d}^{\ast }\ll \tilde{h}_{0}$ , yields the linearized viscous–elastic governing equation in terms of deformation accounting for both bending and pre-stretching (Kodio et al. Reference Kodio, Griffiths and Vella2017), which contains a source term that depends on the external forces

(2.13) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FB}_{\Vert }^{6}d+\unicode[STIX]{x1D6FB}_{\Vert }^{4}d=-\unicode[STIX]{x1D735}_{\Vert }\boldsymbol{\cdot }\boldsymbol{f}_{F}+\unicode[STIX]{x1D6FB}_{\Vert }^{2}f_{E}. & & \displaystyle\end{eqnarray}$$

For completeness, the corresponding dimensional equation is given by

(2.14) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\tilde{d}}{\unicode[STIX]{x2202}\tilde{t}}-\frac{\tilde{h}_{0}^{3}}{12\tilde{\unicode[STIX]{x1D707}}}\tilde{\unicode[STIX]{x1D735}}_{\Vert }\boldsymbol{\cdot }[\tilde{\unicode[STIX]{x1D735}}_{\Vert }(\tilde{B}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{4}\tilde{d}-\tilde{T}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{2}\tilde{d})]=-\tilde{h}_{0}\tilde{\unicode[STIX]{x1D735}}_{\Vert }\boldsymbol{\cdot }\tilde{\boldsymbol{f}}_{F}+\frac{\tilde{h}_{0}^{3}}{12\tilde{\unicode[STIX]{x1D707}}}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{2}\tilde{f}_{E}. & & \displaystyle\end{eqnarray}$$

The governing equations (2.12)–(2.14) are supplemented by the boundary conditions

(2.15a-c ) $$\begin{eqnarray}\displaystyle d=0,\quad \frac{\unicode[STIX]{x2202}^{2}d}{\unicode[STIX]{x2202}x^{2}}=0,\quad \frac{\unicode[STIX]{x2202}^{4}d}{\unicode[STIX]{x2202}x^{4}}=0\quad \text{at }x=0,l_{m}\text{ and }0\leqslant y\leqslant w_{m}, & & \displaystyle\end{eqnarray}$$

(2.16a-c ) $$\begin{eqnarray}\displaystyle d=0,\quad \frac{\unicode[STIX]{x2202}^{2}d}{\unicode[STIX]{x2202}y^{2}}=0,\quad \frac{\unicode[STIX]{x2202}^{4}d}{\unicode[STIX]{x2202}y^{4}}=0\quad \text{at }y=0,w_{m}\text{ and }0\leqslant x\leqslant l_{m}, & & \displaystyle\end{eqnarray}$$

and the initial condition $d(x,y,t=0)=0$ . The first two boundary conditions (2.15a,b ) and (2.16a,b ) correspond to no deflection and no moment at the boundaries, whereas the last conditions (2.15c ) and (2.16c ) are obtained from (2.2) by further assuming that the fluidic and external pressures are both zero at the boundaries, $p=p_{e}=0$ (see, e.g. Kodio et al. Reference Kodio, Griffiths and Vella2017). We note that the condition $p=0$ is motivated by the experimental set-up shown in figure 1, where the fabricated chamber has open boundaries, which are well represented by zero gauge pressure. In addition, for the case of a tension-dominant regime, which is the main focus of this work, setting $\unicode[STIX]{x1D706}=0$ in (2.12)–(2.14) yields fourth-order governing equations. The boundary conditions (2.15a,b ) and (2.16a,b ) hold for this case, with the second derivative requirement arising directly from the zero pressure condition (rather than from the no-moment requirement in the sixth-order case).

The governing equation (2.13) can be solved by several methods. For the boundary conditions under consideration (2.15)–(2.16), eigenfunction expansions or a Green’s functions approach are particularly suitable. We refer the reader to Prosperetti (Reference Prosperetti2011) for both methods of solution, and note that in practice these calculations may be cumbersome. The general solution based on a Green’s function approach for the case presented here is provided in appendix B of the supplementary material. However, we stress that both eigenfunctions and Green’s functions approaches would be far more complex to implement if boundary conditions other than those presented here are adopted. Furthermore, analytical treatment would be greatly complicated, if at all possible, in non-rectangular and non-circular domains of interest.

3 Deformations due to a square-shaped actuation region

We now consider the non-uniform slip velocity (2.6) as a driving force and specifically focus on a square-shaped actuation of the form

(3.1) $$\begin{eqnarray}\displaystyle f_{F_{x}}(x,y,t)=\left\{\begin{array}{@{}ll@{}}\displaystyle E(t) & |x-c_{x}|\leqslant L\text{ and }|y-c_{y}|\leqslant L\\ \displaystyle 0 & \text{otherwise},\end{array}\right. & & \displaystyle\end{eqnarray}$$

corresponding to a square-shaped zeta-potential distribution, where $E(t)$ is a spatially homogeneous and time-dependent electric field along the $\hat{\boldsymbol{x}}$ axis. Here $2L$ is the side length of the actuation square, whereas $c_{x}$ and $c_{y}$ indicate the $x$ - and $y$ -coordinates of its centre. For convenience, hereafter we set $l_{m}=w_{m}=\unicode[STIX]{x03C0}$ .

For a suddenly applied actuation, $E(t)=H(t)$ , where $H$ is the Heaviside step function, using closed-form solutions derived in appendix B of the supplementary material, we obtain the deformation field resulting from (3.1)

(3.2a ) $$\begin{eqnarray}\displaystyle d(x,y,t)=\frac{16}{\unicode[STIX]{x03C0}^{2}}\mathop{\sum }_{m,n=1}^{\infty }{\mathcal{A}}(m,n;L,c_{x},c_{y})\frac{\sin (mx)\sin (ny)}{nF(m,n;\unicode[STIX]{x1D706})}(1-\text{e}^{-F(m,n;\unicode[STIX]{x1D706})t}), & & \displaystyle\end{eqnarray}$$

(3.2b,c ) $$\begin{eqnarray}\displaystyle {\mathcal{A}}=\cos (mc_{x})\sin (nc_{y})\sin (mL)\sin (nL),\quad F=(m^{2}+n^{2})^{2}[\unicode[STIX]{x1D706}(m^{2}+n^{2})+1].\qquad & & \displaystyle\end{eqnarray}$$

In appendix C of the supplementary material, we use numerical computations to verify the analytical solution (3.2) for the case of $\unicode[STIX]{x1D706}=0$ , showing excellent agreement.

Figure 3. Investigation of the deformation field resulting from two square-shaped regions with opposite signs of zeta potential subjected to a constant electric field suddenly applied at $t=0$ . (a) The deformation field (colour map) at $t=0.1$ , superposed on a schematic illustration of the two oppositely directed electro-osmotic actuation regions. (b) The scaled steady-state deformation $d/d_{max}$ along the $\hat{\boldsymbol{x}}$ axis, showing the evolution of the deformation for $\unicode[STIX]{x1D706}=10^{-3}$ , and a comparison between the steady-state deformation at $\unicode[STIX]{x1D706}=10^{-3}$ (tension dominant, with $d_{max}=3.7\times 10^{-2}$ ) and $\unicode[STIX]{x1D706}=10$ (bending dominant, with $d_{max}=1.2\times 10^{-3}$ ). (c,d) The evolution of the maximum deformation and pressure as a function of time. (e) The maximum deformation $\tilde{d}_{max}$ at steady state, scaled by $\tilde{h}_{0}$ , as a function of $\unicode[STIX]{x1D706}$ , for $\tilde{E}_{Y}=1~\text{MPa}$ and $\tilde{E}_{Y}=1~\text{GPa}$ . The dashed and solid grey lines indicate the validity boundaries $\tilde{d}_{max}/\tilde{h}_{0}=0.1$ and $\unicode[STIX]{x1D6FC}=0.1$ , respectively, where the greyed-out region is beyond the validity of the model. All calculations were performed using $\tilde{l}^{\ast }=5~\text{mm}$ , $\tilde{h}_{0}=100~\unicode[STIX]{x03BC}\text{m}$ , $\tilde{p}^{\ast }\approx 1~\text{Pa}$ , $\tilde{h}_{m}=10~\unicode[STIX]{x03BC}\text{m}$ , $\unicode[STIX]{x1D708}=0.5$ , $c_{x,1}=2\unicode[STIX]{x03C0}/5$ , $c_{x,2}=3\unicode[STIX]{x03C0}/5$ , $c_{y,1}=c_{y,2}=\unicode[STIX]{x03C0}/2$ and $L=\unicode[STIX]{x03C0}/10$ .

Figure 3(a) presents the deformation field (3.2) at $t=0.1$ , resulting from two square-shaped regions with opposite signs of zeta potential, subjected to an electric field suddenly applied at $t=0$ (see also supplementary movie 1). The opposing flows result in a positive internal pressure at the interface between the two regions, and negative pressure at the far edges of the squares. Figure 3(b) presents the deformation field (along the $\hat{\boldsymbol{x}}$ axis) as a function of time, for a pre-stretching-dominant regime $\unicode[STIX]{x1D706}=10^{-3}$ . The steady-state deformation ( $t=1$ ) with $d_{max}=3.7\times 10^{-2}$ is compared with that obtained for a bending-dominant regime ( $\unicode[STIX]{x1D706}=10$ ) having $d_{max}=1.2\times 10^{-3}$ , showing that the higher order of the dominant term in the bending regime effectively results in additional averaging of the fluidic pressure and thus reduces the spatial resolution and the magnitude of the deformation field. Figure 3(c,d) presents the development of the maximum pressure and deformation as a function of time towards steady state, where $p_{max}=3.4\times 10^{-1}$ and $d_{max}=3.7\times 10^{-2}$ . Importantly, while for a rigid configuration actuation of the electric field would result in an instantaneous jump in pressure (on an acoustic time scale), in the elastic case the evolution of the pressure is coupled to that of the deformation, as given by (2.11). Nevertheless, for the case presented, the rise time to 50 % of the final pressure is $t=10^{-3}$ , nearly two orders of magnitude shorter than the viscous–elastic time scale $t=0.7\times 10^{-1}$ in which the deformation reaches 50 % of its steady-state value. This time delay between pressure and deformation is also evident in periodic actuations (see supplementary movie 2) where the deformation phase lags behind that of the pressure. We note that at extremely short time scales, inertial effects must be included to properly describe the deformation dynamics. Figure 3(e) presents the maximum deformation (scaled by $\tilde{h}_{0}$ ) at steady state as a function of $\unicode[STIX]{x1D706}$ , for a fixed sheet thickness and two different elastic moduli. Consistent with the scaling (2.8) and non-dimensional solution (3.2), as $\unicode[STIX]{x1D706}$ increases ( $\tilde{T}$ decreases), the magnitude of the deformation increases until it reaches a constant value when the problem is dominated entirely by bending. The dashed and solid grey lines respectively indicate the range of validity of the model, where the deformation is no longer small, $\tilde{d}_{max}/\tilde{h}_{0}=0.1$ , and where internal stretching is non-negligible, $\unicode[STIX]{x1D6FC}=0.1$ . For clarity, the regions of the graph which are beyond the validity of the model are greyed out.

4 Effect of pre-stretching on the resolution, magnitude and time scale of the deformation field

As we are aiming to utilize fluidic actuation as a mechanism to create desired deformation in the elastic sheet, it is of interest to examine the trade-off between the magnitude and time scale of deformation’s and the attainable resolution, for a given amplitude of the actuation force. Any desired deformation can be written as a Fourier series on a finite domain, where the resulting spectrum of frequencies would depend on the desired shape and on the deformation’s position relative to the boundaries. Consider a steady-state deformation of the form

(4.1) $$\begin{eqnarray}\displaystyle \tilde{d}(\tilde{x},{\tilde{y}})=\tilde{d}_{0}\sin \left(k_{x}\unicode[STIX]{x03C0}\frac{\tilde{x}}{\tilde{l}_{m}}\right)\sin \left(k_{y}\unicode[STIX]{x03C0}\frac{{\tilde{y}}}{\tilde{w}_{m}}\right), & & \displaystyle\end{eqnarray}$$

created due to a non-uniform Helmholtz–Smoluchowski slip velocity, (2.6), acting in the $\hat{\boldsymbol{x}}$ direction. Here $k_{x}$ and $k_{y}$ are wavenumbers in the $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{y}}$ directions, respectively, and $\tilde{d}_{0}$ is the amplitude yet to be determined. Using (2.6) and (2.14), the zeta potential required for generating the deformation (4.1) is given by

(4.2) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D701}}(\tilde{x},{\tilde{y}})=\frac{\tilde{h}_{0}^{2}}{6\tilde{\unicode[STIX]{x1D700}}\tilde{E}}\int (-\tilde{B}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{6}\tilde{d}+\tilde{T}\tilde{\unicode[STIX]{x1D6FB}}_{\Vert }^{4}\tilde{d})\,\text{d}\tilde{x}=\tilde{\unicode[STIX]{x1D701}}_{0}\cos \left(k_{x}\unicode[STIX]{x03C0}\frac{\tilde{x}}{\tilde{l}_{m}}\right)\sin \left(k_{y}\unicode[STIX]{x03C0}\frac{{\tilde{y}}}{\tilde{w}_{m}}\right), & & \displaystyle\end{eqnarray}$$

enabling us to explicitly express the amplitude $\tilde{d}_{0}$ in terms of relevant physical quantities,

(4.3) $$\begin{eqnarray}\displaystyle \tilde{d}_{0}=-\frac{\tilde{6\unicode[STIX]{x1D700}}\tilde{\unicode[STIX]{x1D701}}_{0}\tilde{E}}{\unicode[STIX]{x03C0}^{3}\tilde{h}_{0}^{2}}\frac{k_{x}/\tilde{l}_{m}}{((k_{x}/\tilde{l}_{m})^{2}+(k_{y}/\tilde{w}_{m})^{2})^{2}(\unicode[STIX]{x03C0}^{2}\tilde{B}[(k_{x}/\tilde{l}_{m})^{2}+(k_{y}/\tilde{w}_{m})^{2}]+\tilde{T})}. & & \displaystyle\end{eqnarray}$$

Substituting (4.2) into (2.14), we obtain the corresponding time-dependent solution

(4.4) $$\begin{eqnarray}\displaystyle \tilde{d}(\tilde{x},{\tilde{y}},\tilde{t})=\tilde{d}_{0}\sin \left(k_{x}\unicode[STIX]{x03C0}\frac{\tilde{x}}{\tilde{l}_{m}}\right)\sin \left(k_{y}\unicode[STIX]{x03C0}\frac{{\tilde{y}}}{\tilde{w}_{m}}\right)(1-\text{e}^{-\tilde{t}/\tilde{\unicode[STIX]{x1D70F}}_{0}}), & & \displaystyle\end{eqnarray}$$

that evolves towards steady-state deformation $\tilde{d}(\tilde{x},{\tilde{y}})$ , (4.1), and provides the viscous–elastic time scale $\tilde{\unicode[STIX]{x1D70F}}_{0}$ required to achieve this steady state,

(4.5) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D70F}}_{0}=\frac{12\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x03C0}^{4}\tilde{h}_{0}^{3}}\frac{1}{((k_{x}/\tilde{l}_{m})^{2}+(k_{y}/\tilde{w}_{m})^{2})^{2}(\unicode[STIX]{x03C0}^{2}\tilde{B}[(k_{x}/\tilde{l}_{m})^{2}+(k_{y}/\tilde{w}_{m})^{2}]+\tilde{T})}. & & \displaystyle\end{eqnarray}$$

Figure 4. The effect of wavenumber on the magnitude of the deformation and the time scale of viscous–elastic interaction. (a) Deformation in dimensional form along the $\hat{\boldsymbol{x}}$ axis for various values of $k_{x}$ , with $k_{y}=1$ and $\tilde{h}_{m}=10~\unicode[STIX]{x03BC}\text{m}$ . (b,c) The maximum (dimensional) deformation $\tilde{d}_{0}$ and time scale $\tilde{\unicode[STIX]{x1D70F}}_{0}$ , respectively, as a function of wavenumber $k_{x}$ , for $k_{y}=1$ and different values of the elastic sheet thickness, $\tilde{h}_{m}$ . All calculations were performed using non-uniform electro-osmotic actuation as the driving mechanism, with $\tilde{h}_{0}=100~\unicode[STIX]{x03BC}\text{m}$ , $\tilde{l}_{m}=\tilde{w}_{m}=4~\text{cm}$ , $\tilde{E}_{Y}=0.3~\text{GPa}$ , $\unicode[STIX]{x1D708}=0.5$ , $\tilde{T}=15~\text{Pa}~\text{m}$ , $\tilde{\unicode[STIX]{x1D701}}_{0}=-70~\text{mV}$ , $\tilde{E}=100~\text{V}~\text{cm}^{-1}$ and $\tilde{\unicode[STIX]{x1D700}}=7.08\times 10^{-10}~\text{F}~\text{m}^{-1}$ .

Using the physical values noted in its caption, figure 4(a) presents the resulting steady-state deformation (4.1) along the $\hat{\boldsymbol{x}}$ axis for various values of $k_{x}$ , with $k_{y}=1$ and $\tilde{h}_{m}=10~\unicode[STIX]{x03BC}\text{m}$ , and clearly shows the reduction in deformation magnitude as the wavenumber increases. Furthermore, (4.3) and (4.5) indicate that the scaling of the amplitude $\tilde{d}_{0}$ and time scale $\tilde{\unicode[STIX]{x1D70F}}_{0}$ with the wavenumber depends on the relative contribution of the bending and tension terms in the denominator. When $\unicode[STIX]{x03C0}^{2}\tilde{B}k_{x}^{2}\ll \tilde{T}\tilde{l}_{m}^{2}$ , pre-stretching is dominant over bending forces and the deformation $\tilde{d}_{0}$ and time scale $\tilde{\unicode[STIX]{x1D70F}}_{0}$ scale as $k_{x}^{-3}$ and $k_{x}^{-4}$ , respectively. However, when $\unicode[STIX]{x03C0}^{2}\tilde{B}k_{x}^{2}\gg \tilde{T}\tilde{l}_{m}^{2}$ , bending is dominant, and the deformation and time scale decrease as $k_{x}^{-5}$ and $k_{x}^{-6}$ , respectively. From the definition of $\tilde{B}$ , the condition for pre-stretching dominance can be expressed as $\unicode[STIX]{x03C0}^{2}\tilde{E}_{Y}\tilde{h}_{m}^{3}k_{x}^{2}/12(1-\unicode[STIX]{x1D708}^{2})\tilde{T}\tilde{l}_{m}^{2}\ll 1$ . Figure 4(b,c) presents the maximum deformation and the time scale as a function of $k_{x}$ , respectively, for three values of the membrane thickness, $\tilde{h}_{m}$ , showing that the deformation of a $10~\unicode[STIX]{x03BC}\text{m}$ thick sheet scales as $k_{x}^{-3}$ throughout the investigated range, whereas the time scale scales as $k_{x}^{-4}$ . As the membrane thickness increases, the bending effect becomes apparent and for $\tilde{h}_{m}=100~\unicode[STIX]{x1D707}\text{m}$ the amplitude and the time scale decrease as $k_{x}^{-3}$ and $k_{x}^{-4}$ for sufficiently low $k_{x}$ values, but settle on a $k_{x}^{-5}$ and $k_{x}^{-6}$ dependence for high wavenumbers. For $\tilde{h}_{m}=1~\text{mm}$ , the bending effect is dominant and thus the amplitude and the time scale decrease as $k_{x}^{-5}$ and $k_{x}^{-6}$ even for low wavenumbers. We note that for the set of parameters chosen here, even the first modes corresponding to the largest deformation satisfy the small-deformation requirement of the model.

5 Effect of discretized actuation on the deformation field

Equation (2.13) may be solved to obtain the actuation field required to achieve pre-defined deformation patterns of elastic plates, which may be desired in engineering applications. However, implementation of the resulting continuous actuation distribution may be challenging in practice. We here consider a discrete actuation profile, consisting of $D\times D$ individual squares, and seek to minimize the error between the resulting and desired deformation, for a given level of discretization. Specifically, actuating square $i$ with amplitude $a_{i}$ yields a deflection $a_{i}d_{i}(x,y)$ given by (3.2a ), such that the resulting deformation is $\sum _{i=1}^{D^{2}}a_{i}d_{i}(x,y)$ . For a given desired deformation $d(x,y)$ , we follow a least-squares method (Bjorck Reference Bjorck1996) and seek to minimize the error

(5.1) $$\begin{eqnarray}\displaystyle \int _{0}^{l_{m}}\int _{0}^{w_{m}}\left(\mathop{\sum }_{i,j=1}^{D^{2}}a_{i}d_{i}-d\right)^{2}\text{d}x\,\text{d}y=\mathop{\sum }_{i,j=1}^{D^{2}}a_{i}a_{j}\unicode[STIX]{x1D63C}_{i,j}-2\mathop{\sum }_{i=1}^{D^{2}}a_{i}b_{i}+c, & & \displaystyle\end{eqnarray}$$

where

(5.2a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}_{i,j}=\int _{0}^{l_{m}}\int _{0}^{w_{m}}d_{i}d_{j}\,\text{d}x\,\text{d}y,\quad b_{i}=\int _{0}^{l_{m}}\int _{0}^{w_{m}}d_{i}d\,\text{d}x\,\text{d}y,\quad c=\int _{0}^{l_{m}}\int _{0}^{w_{m}}d^{2}\,\text{d}x\,\text{d}y.\quad \qquad & & \displaystyle\end{eqnarray}$$

Differentiating (5.1) with respect to each $a_{i}$ and equating to zero yields the optimality condition

(5.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}\boldsymbol{a}=\boldsymbol{b}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{a}=[a_{1},a_{2},\ldots ,a_{D^{2}}]^{\text{T}}$ , $\boldsymbol{b}=[b_{1},b_{2},\ldots ,b_{D^{2}}]^{\text{T}}$ and $\unicode[STIX]{x1D63C}$ is a symmetric $D^{2}\times D^{2}$ matrix with rank $D^{2}-D$ and thus is not invertible (singular). The singularity of $\unicode[STIX]{x1D63C}$ stems from the fact that the source term in (2.13) depends on gradients of the actuation field, thus allowing an associated gauge freedom in the choice of actuation (coefficients $a_{i}$ ) without modifying the resulting deformation. Specifically, for the case of a driving force $f_{F_{x}}$ acting in $\hat{\boldsymbol{x}}$ the direction, it follows that adding an arbitrary function $f_{0}(y)$ to the driving force, which has $D$ values in discrete form, will not modify the resulting deformation.

Figure 5. The effect of actuation discretization on the shape and the magnitude of the deformation. (a) Deformation formed by the continuous driving force $f_{F_{x}}$ shown in (c), acting in the $\hat{\boldsymbol{x}}$ direction. (b) Deformation obtained from $4\times 4$ square-shaped actuation regions, each having a uniform value which is found by solving the least-squares problem (5.3) and shown in (d). (e,f) Deformation along the $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{y}}$ axes, respectively, for different levels of discretization.

Figure 6. Demonstration of the use of a least-squares discretization method for creation of a complex deformation. (a) Top view of a microfluidic configuration, in which black lines represent the confining walls, characterized by step-like behaviour. (b) The deformation field (colour map) obtained from $28\times 28$ square-shaped actuation regions, each having a uniform value found by solving the least-squares problem (5.3). (c) Deformation along the two cross-sections, depicted in (a) and (b), showing that discretization results in undesired oscillations of the deformation field. Dashed and solid lines correspond to the chosen and obtained deformation, respectively.

To determine the coefficients $a_{i}$ , we first reduced the matrix $\unicode[STIX]{x1D63C}$ and vector $\boldsymbol{b}$ only to rows with corresponding non-zero eigenvalues of $\unicode[STIX]{x1D63C}$ , obtaining a $(D^{2}-D)\times (D^{2})$ matrix and a $(D^{2}-D)\times (1)$ vector, respectively. We then solved the reduced system (5.3) and found the vector $\boldsymbol{a}$ using MATLAB’s routine lsqminnorm (release R2017b, Mathworks, USA), which computes the minimum least-squares solution of the system. As an illustrative example, we consider a localized deformation,

(5.4) $$\begin{eqnarray}\displaystyle d(x,y)=(2/\unicode[STIX]{x03C0})^{12}x^{3}(\unicode[STIX]{x03C0}-x)^{3}y^{3}(\unicode[STIX]{x03C0}-y)^{3}, & & \displaystyle\end{eqnarray}$$

shown in figure 5(a), and for simplicity assume a membrane regime (dominant pre-stretching, $\unicode[STIX]{x1D706}=0$ ). Using (2.13), the continuous driving force $f_{F_{x}}$ required to create this deformation (at steady state) is given by

(5.5) $$\begin{eqnarray}\displaystyle f_{F_{x}}(x,y)=-\int _{\unicode[STIX]{x03C0}/2}^{x}\unicode[STIX]{x1D6FB}_{\Vert }^{4}d(x^{\prime },y)\,\text{d}x^{\prime }, & & \displaystyle\end{eqnarray}$$

and shown in figure 5(c). Figure 5(d) shows the discretization of this field into $4\times 4$ square-shaped regions, each assigned with a uniform value obtained from the least-squares method solution of (5.3). Figure 5(b) presents the resulting deformation obtained by superposition of solutions for individual squares (3.2), using the values for forcing shown in figure 5(d). Figure 5(e,f) presents the deformation along the $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{y}}$ axes, respectively, for different levels of discretization. While, as expected, discretization results in undesired oscillations of the deformation field due to high wavenumbers which are unbalanced, even a relatively coarse discretization with $D=4$ allows maintaining of a fairly localized deformation, with amplitude only larger by $4\,\%$ than the desired one. For higher levels of discretization ( $D=6$ and $D=8$ ), the difference between the resulting and desired deformations is almost indistinguishable.

To challenge our method of solution, we consider the microfluidic configuration shown in figure 6(a), in which the deformation magnitude is dictated strictly by 0 and 1, represented in the figure by white and black colours, respectively. Clearly, the exact description of this deformation field requires an infinite number of wavenumbers which cannot be represented by a finite grid discretization. Nevertheless, the deformation obtained using $28\times 28$ square-shaped actuation regions is able to reproduce the main features of the geometry as shown in figure 6(b), although with undesired oscillations. For further clarification, figure 6(c) presents the deformation along the two cross-sections, illustrated in figure 6(a,b), showing that the maximal and minimal deformations are approximately by $11\,\%$ and $18\,\%$ larger than the desired one, respectively.

6 Investigation of nonlinear effects

In the previous sections we restricted our analysis to a strongly pre-stretched elastic sheet, $\unicode[STIX]{x1D6FC}\ll 1$ , and small elastic deformations, $\unicode[STIX]{x1D6FD}\ll 1$ , and solved the linear governing equation (2.13) for the deformation field. In this section, we relax these restrictions and explore the nonlinear effects arising from the hydrodynamic nonlinearity in the Reynolds equation, $h^{3}$ , and the nonlinear coupling between the tension and deformation in the Föppl–von Kármán equations. However, we limit our nonlinear analysis to positive deformations, $d>0$ . We first study theoretically the weakly nonlinear effect of internal tension on the deformation (§ 6.2) and then investigate numerically the nonlinear effects considering finite $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ (§ 6.3). To this end, we consider an axisymmetric geometry which allows reduction to a single spatial variable and thus greatly simplifies theoretical and numerical investigation. This simple case is useful in providing physical insight on the nonlinear effects of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ on the deformation and tension field. We also note the analysis of the axisymmetric configurations is relevant to adaptive optics applications, in which the liquid within a circular microchamber, covered with a thin elastic membrane, is pressurized to deform the membrane forming a plano-convex microlens (Chronis et al. Reference Chronis, Liu, Jeong and Lee2003).

6.1 Coupled Reynolds and nonlinear Föppl–von Kärmán equations

We consider a circular elastic sheet of radius $\tilde{R}_{m}$ subjected to axisymmetric driving forces acting either through the pressure applied directly to the elastic sheet or through the non-uniform slip velocity in the fluid. For simplicity, we neglect the effect of bending and consider a membrane (tension-dominant) regime, $\unicode[STIX]{x1D706}\ll 1$ . Based on these assumptions, the deformations of the elastic membrane are described by the axisymmetric Föppl–von Kármán equations for the membrane (see Landau & Lifshitz Reference Landau and Lifshitz1959; Lister et al. Reference Lister, Peng and Neufeld2013; Zheng, Griffiths & Stone Reference Zheng, Griffiths and Stone2015)

(6.1) $$\begin{eqnarray}\displaystyle \tilde{p}=-\tilde{\unicode[STIX]{x1D735}}_{r}\boldsymbol{\cdot }(\tilde{T}_{r}\tilde{\unicode[STIX]{x1D735}}_{r}\tilde{d})+\tilde{p}_{e}, & & \displaystyle\end{eqnarray}$$

and

(6.2) $$\begin{eqnarray}\displaystyle \frac{1}{\tilde{r}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\tilde{r}}\left(\tilde{r}^{3}\frac{\unicode[STIX]{x2202}\tilde{T}_{r}}{\unicode[STIX]{x2202}\tilde{r}}\right)=-\frac{\tilde{E}_{Y}\tilde{h}_{m}}{2}\left(\frac{\unicode[STIX]{x2202}\tilde{d}}{\unicode[STIX]{x2202}\tilde{r}}\right)^{2}, & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D735}}_{r}=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\tilde{r})\hat{\boldsymbol{r}}$ is the gradient in polar coordinates, $\tilde{r}=\sqrt{\tilde{x}^{2}+{\tilde{y}}^{2}}$ is the radial coordinate and $\tilde{T}_{r}(\tilde{r},\tilde{t})$ is the radial tension resulting from both external and internal tension formed in the membrane.

We define the normalized radial coordinate $r=\tilde{r}/\tilde{l}$ , radial tension $T_{r}=\tilde{T}_{r}/\tilde{T}$ and dimensionless size of the membrane $R_{m}=\tilde{R}_{m}/\tilde{l}$ . Substituting the normalized variables into (6.1)–(6.2) and using the results of scaling analysis in § 2.1 we obtain the non-dimensional Föppl–von Kármán equations

(6.3) $$\begin{eqnarray}\displaystyle & \displaystyle p=-\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }(T_{r}\unicode[STIX]{x1D735}_{r}d)+p_{e}, & \displaystyle\end{eqnarray}$$

(6.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r^{3}\frac{\unicode[STIX]{x2202}T_{r}}{\unicode[STIX]{x2202}r}\right)=-\frac{1}{2}\unicode[STIX]{x1D6FC}\left(\frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}r}\right)^{2}. & \displaystyle\end{eqnarray}$$

Substituting (6.3) into the nonlinear Reynolds equation (2.11), we obtain

(6.5) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }[h^{3}\unicode[STIX]{x1D735}_{r}[\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }(T_{r}\unicode[STIX]{x1D735}_{r}d)]]=-\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }[h\boldsymbol{f}_{F}]+\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }[h^{3}\unicode[STIX]{x1D735}_{r}f_{E}]. & & \displaystyle\end{eqnarray}$$

The evolution equation (6.5) and the second Föppl–von Kármán equation (6.4) are two-way coupled nonlinear equations, governing the fluid–structure interaction, that should be solved at once to obtain both deformation and tension fields.

The governing equations (6.4) and (6.5) are supplemented by six boundary conditions. At the centre of the membrane, we require regularity of $d(r,t)$ and $T_{r}(r,t)$

(6.6a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}r}=0\quad \text{and}\quad \frac{\unicode[STIX]{x2202}T_{r}}{\unicode[STIX]{x2202}r}=0\quad \text{at }r=0. & & \displaystyle\end{eqnarray}$$

We also assume a zero flux $q$ at the origin, given by $q=r(-h^{3}\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}r+(1/2)hu_{slip})$ , that under assumptions of finite $u_{slip}$ and $\unicode[STIX]{x2202}p_{e}/\unicode[STIX]{x2202}r$ at $r=0$ from (6.1) implies

(6.7) $$\begin{eqnarray}\displaystyle \lim _{r\rightarrow 0}rh^{3}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left[\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(rT_{r}\frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}r}\right)\right]=0. & & \displaystyle\end{eqnarray}$$

At the edge of the membrane, we consider the following boundary conditions:

(6.8a-c ) $$\begin{eqnarray}\displaystyle d=r\frac{\unicode[STIX]{x2202}T_{r}}{\unicode[STIX]{x2202}r}+(1-\unicode[STIX]{x1D708})(T_{r}-1)=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(rT_{r}\frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}r}\right)=0\quad \text{at }r=R_{m}. & & \displaystyle\end{eqnarray}$$

The first two conditions (6.8a,b ) correspond to no transverse displacement, and fixed horizontal displacement at an outer circular frame holding the membrane, respectively. The last condition (6.8c ) is obtained, as previously, from the elastic balance (6.3) by assuming the fluidic and external pressures are both zero at the boundaries, $p=p_{e}=0$ .

6.2 Asymptotic analysis for weakly nonlinear effects due to induced tension

In this section, we use asymptotic analysis to study the weakly nonlinear effects arising from internal tension formed in the elastic sheet during the deflection. We apply asymptotic expansions to decouple the two-way coupled nonlinear equations (6.5) and (6.4), assuming small deformation and strong pre-stretching of the elastic membrane, and obtain a correction for the effect of induced internal tension.

Assuming small elastic deformations, $\unicode[STIX]{x1D6FD}=\tilde{d}^{\ast }/\tilde{h}_{0}\ll 1$ , we eliminate the hydrodynamic nonlinearity in the Reynolds equation and obtain

(6.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D6FB}_{r}^{2}[\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }(T_{r}\unicode[STIX]{x1D735}_{r}d)]=-\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }\boldsymbol{f}_{F}+\unicode[STIX]{x1D6FB}_{r}^{2}f_{E}. & & \displaystyle\end{eqnarray}$$

For the case of strongly pre-stretched elastic membrane, $\unicode[STIX]{x1D6FC}=\tilde{T}_{in}/\tilde{T}\ll 1$ , we expand the deformation and the tension in powers of $\unicode[STIX]{x1D6FC}$

(6.10a,b ) $$\begin{eqnarray}\displaystyle d=d^{(0)}+\unicode[STIX]{x1D6FC}d^{(1)}+O(\unicode[STIX]{x1D6FC}^{2}),\quad T_{r}=1+\unicode[STIX]{x1D6FC}T_{r}^{(1)}+O(\unicode[STIX]{x1D6FC}^{2}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is a small parameter satisfying $\text{max}(\unicode[STIX]{x1D716}Re,Wo,\unicode[STIX]{x1D716}^{2},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D706})\ll \unicode[STIX]{x1D6FC}\ll 1$ . Substituting (6.10) into the coupled governing equations (6.4) and (6.9), results in three one-way coupled linear equations for the leading-order deformation and the first-order correction for the deformation and tension fields. In appendix D of the supplementary material, we present a detailed derivation of the governing equations and the appropriate boundary conditions at each order, and provide closed-form solutions for the leading-order deformation and the first-order correction for the tension.

Similarly to § 3, as an illustrative example we consider the non-uniform electro-osmotic slip velocity as a driving force and specifically focus on a spatially non-uniform actuation of the form

(6.11) $$\begin{eqnarray}\displaystyle \boldsymbol{f}_{F}(r,t)=V_{r}\text{J}_{1}\left(\frac{\unicode[STIX]{x1D712}_{n}r}{R_{m}}\right)H(t)\hat{\boldsymbol{r}}, & & \displaystyle\end{eqnarray}$$

corresponding to a zeta-potential distribution of $\unicode[STIX]{x1D701}(r)=r\text{J}_{1}(\unicode[STIX]{x1D712}_{1}r/R_{m})$ , subjected to a suddenly applied electric field $V_{r}/r$ in the $\hat{\boldsymbol{r}}$ direction, where $\text{J}_{1}(\unicode[STIX]{x1D712}_{n}r/R_{m})$ is the Bessel function of the first kind and of the first order, and $\unicode[STIX]{x1D712}_{n}$ is the $n$ th root of $\text{J}_{0}(\unicode[STIX]{x1D712}_{n})=0$ . In our analysis, we hereafter focus on the first root, $\unicode[STIX]{x1D712}_{1}=2.4048$ .

Using the governing equations and closed-form solutions derived in appendix D of the supplementary material, we can obtain analytical solutions for $d^{(0)}(r,t)$ and $T_{r}^{(1)}(r,t)$ resulting from (6.11). The leading-order deformation field is

(6.12) $$\begin{eqnarray}\displaystyle d^{(0)}(r,t)=-\frac{V_{r}R_{m}^{3}}{\unicode[STIX]{x1D712}_{n}^{3}}\text{J}_{0}\left(\frac{\unicode[STIX]{x1D712}_{n}r}{R_{m}}\right)\left[1-\exp \left(-\frac{\unicode[STIX]{x1D712}_{n}^{4}}{R_{m}^{4}}t\right)\right], & & \displaystyle\end{eqnarray}$$

while the first-order correction for tension distribution reads

(6.13) $$\begin{eqnarray}\displaystyle T_{r}^{(1)}(r,t)=\frac{V_{r}^{2}R_{m}^{4}}{4\unicode[STIX]{x1D712}_{n}^{4}}\left[\frac{\unicode[STIX]{x1D708}\text{J}_{1}(\unicode[STIX]{x1D712}_{n})^{2}}{1-\unicode[STIX]{x1D708}}+\text{J}_{0}(\bar{r})^{2}+\text{J}_{1}(\bar{r})^{2}-\frac{\text{J}_{0}(\bar{r})\text{J}_{1}(\bar{r})}{\bar{r}}\right][1-\text{e}^{-(\unicode[STIX]{x1D712}_{n}^{4}/R_{m}^{4})t}]^{2}, & & \displaystyle\end{eqnarray}$$

where $\bar{r}=\unicode[STIX]{x1D712}_{n}r/R_{m}$ . The maximum value $T_{r,max}^{(1)}$ is obtained at the centre of the membrane

(6.14) $$\begin{eqnarray}\displaystyle T_{r,max}^{(1)}=\frac{V_{r}^{2}R_{m}^{4}}{4\unicode[STIX]{x1D712}_{n}^{4}}\left[\frac{\unicode[STIX]{x1D708}\text{J}_{1}(\unicode[STIX]{x1D712}_{n})^{2}}{1-\unicode[STIX]{x1D708}}+\frac{1}{2}\right]. & & \displaystyle\end{eqnarray}$$

While it is difficult to obtain a closed-form solution for transient first-order deformation $d^{(1)}(r,t)$ , the corresponding steady-state deformation depends solely on the spatial coordinate and admits a closed-form solution

(6.15) $$\begin{eqnarray}\displaystyle d_{s}^{(1)}(r)=\int _{r}^{R_{m}}T_{r,s}^{(1)}\frac{\text{d}d_{s}^{(0)}}{\text{d}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}, & & \displaystyle\end{eqnarray}$$

where the subscript $s$ denotes the steady state.

Figure 7. The effect of weak nonlinearity due to induced tension. (a,b) Comparison between asymptotic and numerical results for the maximum deformation and tension at steady state as a function of $\unicode[STIX]{x1D6FC}$ . Dashed black and grey lines correspond to the leading- and first-order asymptotic solutions, respectively, while black dots correspond to the numerical solution. (c,d) Comparison between asymptotic (dashed lines) and numerical (solid lines) solutions for the steady-state deformation and tension, respectively, for $\unicode[STIX]{x1D6FC}=0,0.5$ and 1. All calculations were performed using $\unicode[STIX]{x1D6FD}=0$ , $V_{r}=-2$ , $\unicode[STIX]{x1D708}=0.5$ and $R_{m}=2$ .

Figure 7 summarizes the effect of weak nonlinearity due to induced tension on the deformation and tension resulting from the forcing (6.11). Figure 7(a,b) presents the maximum deformation and tension at steady state as a function of $\unicode[STIX]{x1D6FC}$ . Dashed black and grey lines correspond to the leading- and first-order asymptotic solutions obtained from (6.12), (6.14) and (6.15), whereas black dots correspond to the numerical solution. While the leading-order deformation and tension are independent of $\unicode[STIX]{x1D6FC}$ , the first-order correction and numerical solution clearly show the reduction in deformation magnitude and the increase in the resulting tension, respectively, as the parameter $\unicode[STIX]{x1D6FC}$ increases. In figure 7(c,d) we compare the asymptotic (dashed lines) and numerical (solid lines) solutions for the steady-state deformation and tension fields in the case of $\unicode[STIX]{x1D6FC}=0,0.5$ and 1, showing good agreement. As can be inferred from the results of figure 7, for $\unicode[STIX]{x1D6FC}\ll 1$ the first-order asymptotic solutions (6.13) and (6.15) accurately describe the behaviour of the deformation and tension fields, but as $\unicode[STIX]{x1D6FC}$ approaches and passes $O(1)$ the asymptotic solutions overpredict them. Nevertheless, $d_{max}$ and $T_{r,max}$ continue to decrease/increase monotonically with $\unicode[STIX]{x1D6FC}$ and for $\unicode[STIX]{x1D6FC}\gg 1$ scale like $\unicode[STIX]{x1D6FC}^{-1/3}$ and $\unicode[STIX]{x1D6FC}^{1/3}$ , respectively, as shown by the nonlinear investigation in figure 9.

Figure 8. Investigation of nonlinear effects on the transient behaviour and the magnitude of deformation. (a,b) The time evolution of the deformation field for $\unicode[STIX]{x1D6FD}=0$ (a) and $\unicode[STIX]{x1D6FD}=0.1$ (b), in the case of a strongly pre-stretched elastic membrane, $\unicode[STIX]{x1D6FC}=0$ . Grey solid lines represent the numerical results and black dashed lines represent the leading-order asymptotic solution (6.12). (c) The maximum deformation, obtained at the centre of the membrane, as a function of time for various values of $\unicode[STIX]{x1D6FD}$ , with $\unicode[STIX]{x1D6FC}=0$ . (d) The maximum deformation, in dimensional form, scaled by the initial fluid thickness, $\tilde{d}_{max}/\tilde{h}_{0}=\unicode[STIX]{x1D6FD}d_{max}$ as a function of the parameter $\unicode[STIX]{x1D6FD}$ , for various values of $\unicode[STIX]{x1D6FC}$ . In (c,d) solid lines and dots represent the numerical results, respectively, and the dashed black curve represents the asymptotic solution (6.12). All calculations were performed using $V_{r}=-2$ , $\unicode[STIX]{x1D708}=0.5$ and $R_{m}=2$ .

6.3 Numerical investigation of nonlinear effects

To investigate the nonlinear effects of finite $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ on the deformation and tension field, we proceed with a numerical analysis of the viscous–elastic governing equations (6.4)–(6.5). To study the effect of hydrodynamic nonlinearity ( $h^{3}$ ) on the transient behaviour and on the magnitude of deformation, we consider for simplicity the case of a strongly pre-stretched elastic membrane with $T_{r}=1$ ( $\unicode[STIX]{x1D6FC}=0$ ) and solve numerically the nonlinear evolution equation (6.5) for the deformation using finite differences. In addition, to explore the combined effect of internal tension and hydrodynamic nonlinearity on the steady-state behaviour, we solve numerically the corresponding steady-state boundary value problem (6.4)–(6.5) subject to the six boundary conditions (6.6)–(6.8a-c ) using MATLAB’s routine bvp4c. Further details of the numerical procedures and their cross-validation are discussed in appendix E of the supplementary material.

As shown in figure 8(a), we first validate our time-dependent numerical solver by comparing the numerically determined deformation field (solid lines) for $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}=0$ with the leading-order asymptotic solution (6.12) (dashed lines), showing very good agreement for all times. Figure 8(b) presents the time evolution of the deformation profile (solid lines), for $\unicode[STIX]{x1D6FD}=0.1$ and $\unicode[STIX]{x1D6FC}=0$ . Figure 8(c) shows the time evolution of the maximum deformation, obtained at the centre of the membrane, for $\unicode[STIX]{x1D6FD}=0,0.1,0.2,0.5$ and $\unicode[STIX]{x1D6FC}=0$ . Solid lines represent the numerical results, whereas dashed lines represent the leading-order asymptotic solution (6.12). It follows from figure 8(b,c) that as $\unicode[STIX]{x1D6FD}$ increases the resulting maximum deformation in dimensionless form decreases. This behaviour can be explained as follows: since the viscous resistance scales as $h^{-3}=(1+\unicode[STIX]{x1D6FD}d)^{-3}$ , the increase in $\unicode[STIX]{x1D6FD}$ leads to lower internal pressure gradient and thus lower deformation.

Figure 9. Investigation of nonlinear effect of the induced tension on the maximum displacement and tension. (a,b) Maximum displacement, $\unicode[STIX]{x1D6FD}d_{max}=h_{max}-1=\tilde{d}_{max}/\tilde{h}_{0}$ and tension, $T_{r,max}$ , as a function of $\unicode[STIX]{x1D6FC}$ , for $\unicode[STIX]{x1D6FD}=0.1,0.5$ and 1, obtained from numerical solution of (6.4)–(6.5). For $\unicode[STIX]{x1D6FC}\ll 1$ , both deformation and tension are independent of $\unicode[STIX]{x1D6FC}$ , while for $\unicode[STIX]{x1D6FC}\gg 1$ the deformation decreases as $\unicode[STIX]{x1D6FC}^{-1/3}$ and the tension increases as $\unicode[STIX]{x1D6FC}^{1/3}$ . All calculations were performed using $V_{r}=-2$ , $\unicode[STIX]{x1D708}=0.5$ and $R_{m}=2$ .

Figure 8(d) presents the scaled maximum deformation $\unicode[STIX]{x1D6FD}d_{max}$ as a function of $\unicode[STIX]{x1D6FD}$ , where the dots represent the numerical results for $\unicode[STIX]{x1D6FC}=0,2,8,$ and the dashed curve represents the asymptotic solution (6.12), corresponding to $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}=0$ . We note that when exploring the nonlinear effects on the resulting magnitude of maximum deformation it is more convenient to discuss $\unicode[STIX]{x1D6FD}d_{max}$ rather than $d_{max}$ , since the former can be expressed as $\unicode[STIX]{x1D6FD}d_{max}=h_{max}-1=\tilde{d}_{max}/\tilde{h}_{0}$ , representing the relative magnitude of $\tilde{d}$ in terms of the initial fluid thickness, $\tilde{h}_{0}$ . The numerical analysis reveals that the scaled maximum deformation $\tilde{d}_{max}/\tilde{h}_{0}$ increases nonlinearly with $\unicode[STIX]{x1D6FD}$ throughout the investigated range, showing a sub-linear behaviour, which is more pronounced as $\unicode[STIX]{x1D6FC}$ increases. Physically, this means that when $\tilde{d}$ becomes comparable to $\tilde{h}_{0}$ , the dependence of the deformation on the applied electric field is no longer linear due to internal tension and reduction in pressure, resulting in deformation that is lower than predicted by the linear response. However, in the small-deformation and strong pre-stretching limits, $\tilde{d}_{max}/\tilde{h}_{0}$ increases linearly with $\unicode[STIX]{x1D6FD}$ , consistent with the leading-order asymptotic solution (6.12) (dashed line). Furthermore, figure 8(d) shows that, as expected, (6.12) reasonably predicts the maximum deformation, yielding modest relative error ( $|(d_{max}^{num}-d_{max}^{(0)})/d_{max}^{num}|\lesssim 12\,\%$ ) for up to $\unicode[STIX]{x1D6FD}=0.1$ .

Figure 9 illustrates the nonlinear effect of induced tension on the scaled maximum deformation and tension for different values of $\unicode[STIX]{x1D6FD}$ . It is evident that at small values of $\unicode[STIX]{x1D6FC}$ , corresponding to strong pre-stretching, the resulting maximum deformation and tension are almost independent of $\unicode[STIX]{x1D6FC}$ , up to $\unicode[STIX]{x1D6FC}=O(1)$ (see also figure 7). For $\unicode[STIX]{x1D6FC}\gg 1$ and $\unicode[STIX]{x1D6FD}\lesssim 1$ , performing scaling analysis of (6.4)–(6.5) (with $r\sim 1$ ), we obtain $T_{r}\sim \unicode[STIX]{x1D6FC}d^{2}$ and $T_{r}d\sim V_{r}$ thus yielding

(6.16a,b ) $$\begin{eqnarray}\displaystyle d\sim V_{r}^{1/3}\unicode[STIX]{x1D6FC}^{-1/3}\quad \text{and}\quad T_{r}\sim V_{r}^{2/3}\unicode[STIX]{x1D6FC}^{1/3}\quad \text{for }\unicode[STIX]{x1D6FC}\gg 1,\unicode[STIX]{x1D6FD}\lesssim 1, & & \displaystyle\end{eqnarray}$$

which is consistent with the numerical results shown in figure 9. As expected, for $\unicode[STIX]{x1D6FD}\lesssim 1$ the resulting maximum tension $T_{r,max}$ indicates weak dependence on $\unicode[STIX]{x1D6FD}$ over the entire range of values of the parameter $\unicode[STIX]{x1D6FC}$ , as shown in figure 9(b).

7 Concluding remarks

In this work, we examined the effect of pre-stretching and finite boundaries on the elastohydrodynamics of an elastic sheet lying on top of a thin liquid film. Assuming strong pre-stretching and small deformations of the lubricated elastic sheet, we used the linearized Reynolds and Föppl–von Kármán equations to derive general analytical solutions describing the deformation in a finite domain due to external forces. These closed-form solutions for realistic configurations allow one to study the relationship between the magnitude, time scale and resolution of elastic deformations, as well as to examine spatially discretized actuation. The asymptotic analysis of weakly nonlinear effect due to induced tension for the case of an axisymmetric configuration showed that the first-order asymptotic solutions for the deformation and tension field accurately capture the nonlinear trend even for $\unicode[STIX]{x1D6FC}$ of $O(1)$ .

We obtained that in the small-deformation ( $\tilde{d}\ll \tilde{h}_{0}$ ) and strong pre-stretching ( $\tilde{T}_{in}\ll \tilde{T}$ ) limits, the scaling $\tilde{d}\sim \tilde{\unicode[STIX]{x1D700}}\tilde{\unicode[STIX]{x1D701}}\tilde{E}\tilde{l}^{\ast 3}/\tilde{T}\tilde{h}_{0}^{2}$ (obtained from (2.6) combined with (2.8)) is appropriate, showing a linear relation between $\tilde{d}$ and $\tilde{E}$ , and may be used to estimate the resulting deformation. However, this linear dependence of the deformation on the applied forcing breaks down as $\tilde{d}$ becomes comparable to $\tilde{h}_{0}$ , indicating a sub-linear behaviour, which is more pronounced as $\unicode[STIX]{x1D6FC}=\tilde{T}_{in}/\tilde{T}$ increases.

While our main focus was on actuations applied by the fluid (specifically by non-uniform electro-osmotic flow), the governing equation (2.12) can be readily utilized to investigate the viscous–elastic dynamics due to forces applied directly on the elastic sheet. For such actuation, as shown by the Reynolds equation (2.11), the viscous flow arises only from temporal variation of the solid deformation field. Fluid velocity and gauge pressure vanish as the deformation reaches steady state. This is in contrast to the steady state of forcing applied to the fluid, which involves both non-zero fluid velocity and gauge pressure. Furthermore, with appropriate modification of the boundary conditions (e.g. prescribing an external pressure drop), the theoretical model we presented can be readily extended to the study of fluid–structure interaction in elastic microfluidic chips, where micro-channel flow may be driven, for example, by external pressure gradients or peristaltic actuation.

Actuation at the microscale is currently implemented mostly using microelectromechanical technologies, characterized by discrete and rigid elements. The presented results lay the theoretical foundation for implementation of actuation mechanisms based on low Reynolds number fluid–structure interaction. While our analysis focused on small deformations, we believe it may be directly useful for the design of soft and continuous actuators. For example, relevant deformations in configurable optics would be of the order of a wavelength (i.e. ${<}1~\unicode[STIX]{x03BC}\text{m}$ in the visible spectrum), and thus could be well described by our model when implemented on a $10~\unicode[STIX]{x03BC}\text{m}$ liquid layer. Similarly, the field of microfluidics would highly benefit from configurable microstructures of the order of $10~\unicode[STIX]{x03BC}\text{m}$ , which may be implemented on a $100~\unicode[STIX]{x03BC}\text{m}$ thick liquid layer.