2 Problem formulation

Consider a liquid layer confined between a rigid plate and an elastic plate, in a shallow geometry limit, $h_{0}\ll L$ , where $h_{0}$ is the characteristic gap between the plates and $L$ is the typical lateral in-plane dimension. Figure 1 presents a schematic illustration of the configuration and the orientation of the Cartesian coordinate system we employ, in which the plates are parallel to the $x{-}y$ plane and $z$ is the normal coordinate to the plates. The total height of the elastic plate relative to the rigid surface is $h=h_{0}+\unicode[STIX]{x1D702}$ , where $\unicode[STIX]{x1D702}$ is the deformation of the plate.

Figure 1.

Schematic description of a Hele-Shaw configuration, consisting of a thin fluid layer confined between a rigid and an elastic plate of characteristic length $L$ and thickness $w$ , located at an initial distance of $h_{0}$ from one another. The fluid is subjected to a spatial non-uniformity in its volume and/or non-uniformity of surface properties on the rigid bottom plane, on a length scale $r_{0}$ , resulting in a small deformation $\unicode[STIX]{x1D702}\ll h_{0}$ (drawn not to scale) of the elastic plate, corresponding to the linearized lubrication model. The elastic plate spans a length scale $\ell$ . The Cartesian coordinate system is chosen such that the plates are in the $x{-}y$ plane.

We assume that the surface of the rigid plate or the volume between the plates hosts one or more regions where the fluid or wall properties differ from those in the surrounding region, and that the typical size of this region is much larger than the distance between the plates, $h_{0}$ , yet much smaller than the size of the plate,

Typical examples of surface quantities that may admit non-uniform behaviour are slip velocity over surfaces with non-uniform zeta potential or slip length over hydrophobic surfaces. Variations in bulk quantities may correspond to viscosity changes due to localized heating, or electric body forces due to application of an external electric field. In this section, we will determine the resulting flow and pressure distribution due to the presence of such non-homogeneities.

Our starting point is the Navier–Stokes equations,

accompanied by the non-compressibility condition,

where, $U_{i}$ , $F_{i}$ , $\unicode[STIX]{x1D70C}$ , $p$ and $\unicode[STIX]{x1D707}$ are the fluid velocity components, body force components, density, pressure and dynamic viscosity respectively. The indices $i$ , $j$ run over the three Cartesian coordinates $x,y,z$ , and the summation convention over repeated indices is employed. We decompose all three-dimensional vector fields such as fluid velocity, $\boldsymbol{U}$ , and body force, $\boldsymbol{F}$ , into in-plane components ( $x,y$ ) and an out-of-plane component ( $z$ ), using the convention $\boldsymbol{U}=(\boldsymbol{u},u_{z})$ and $\boldsymbol{F}=(\boldsymbol{f},f_{z})$ . The dynamics of an elastic plate can be described by the Föppl–von Karman model, which accounts for both bending and tension forces in the plate. Since the ratio of the stretching to bending free energies scales as $\unicode[STIX]{x1D702}^{2}/w^{2}$ , for small deformations ( $\unicode[STIX]{x1D702}\ll w$ ), stretching can be neglected, leading to the Kirchhoff–Love model (Landau & Lifshitz Reference Landau and Lifshitz1986). This model relates the pressure $p$ and the deformation $\unicode[STIX]{x1D702}$ according to

Here, the bending stiffness coefficient, $D$ , has been expressed in terms of the plate thickness, $w$ , Young’s modulus, $Y$ , and Poisson’s ratio, $\unicode[STIX]{x1D70E}$ , and introduces an elastic time scale $\unicode[STIX]{x1D70F}_{EL}$ according to

We focus on shallow geometries and small Womersley and reduced Reynolds numbers,

where $\mathit{Re}=\unicode[STIX]{x1D70C}u^{(0)}h_{0}/\unicode[STIX]{x1D707}_{0}$ is the Reynolds number (Panton Reference Panton2006) and $\unicode[STIX]{x1D70F}_{f}$ is the typical time scale. We assume that the size of the system, $L$ , is much larger than the typical width of the generated deformation, $\ell$ , which together with the small-deformation assumption combines to

Normalizing by typical scales,

where $u^{(0)}$ , $u_{z}$ , $\unicode[STIX]{x1D707}_{0}$ are the typical magnitude of velocity along the plane, the velocity component in the normal direction and the dynamic viscosity respectively. We can rewrite the three-dimensional Navier–Stokes equations, (2.2), in non-dimensional form as

where we have assumed that $\unicode[STIX]{x1D707}$ does not depend on the vertical coordinate $z$ . Here, the index $\unicode[STIX]{x1D6FC}$ runs over the in-plane components $x,y$ , and $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x,\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y)$ is a two-dimensional gradient. In the leading order of a small reduced Reynolds number, $\unicode[STIX]{x1D700}\mathit{Re}\ll 1$ , we can neglect the inertial terms, while in the limit of shallow geometry, $\unicode[STIX]{x1D700}\ll 1$ , which is equivalent to $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x_{\unicode[STIX]{x1D6FC}}\ll \unicode[STIX]{x2202}/\unicode[STIX]{x2202}z$ , we may drop the in-plane derivatives relative to the normal derivative. Together with the assumption that the time scale $\unicode[STIX]{x1D70F}_{f}$ is such that the Womersley number is small, we may overlook all inertial terms. Under these assumptions, together with an additional requirement for the normal body force component to vanish, $f_{z}=0$ , the three-dimensional nonlinear equations, (2.9), reduce to the following two-dimensional vector and scalar equations:

The normal component of the velocity is assumed to satisfy the no-penetration condition on the confining walls,

while the tangential components on both planes are subject to either Helmholtz–Smoluchowski slip or Navier slip boundary conditions. Helmholtz–Smoluchowski slip occurs over electrically charged surfaces under an externally applied tangential electric field. The interaction of this field with the excess of net charge in the electric double layer (EDL) results in movement of fluid outside the outer edge of the EDL according to the Helmholtz–Smoluchowski equation (Hunter Reference Hunter2001),

where $\unicode[STIX]{x1D716}$ , $\unicode[STIX]{x1D701}$ , $\unicode[STIX]{x1D707}$ and $\boldsymbol{E}$ are the dielectric constant of the liquid, the zeta potential of the confining walls, the dynamic viscosity of the liquid and the tangential imposed electric field respectively. It should be noted that throughout this work we assume a low Dukhin number, such that the ionic species concentrations in the bulk, as well as the associated electric field, are uniform outside the EDL. Otherwise, non-homogeneities of zeta potential and surface conduction introduce non-homogeneity into the resultant electric field in the bulk (Yariv Reference Yariv2004; Khair & Squires Reference Khair and Squires2008). Slip may also emerge over hydrophobic surfaces. Such slip is typically modelled by a Navier boundary condition, which states that the slip velocity near a flat surface is proportional to the local velocity gradient,

Here, the parameter $\unicode[STIX]{x1D6FD}$ is the so-called slip length, which for simple shear flow corresponds to the fictitious distance inside the solid where the fluid velocity would extrapolate to zero, and $n$ is the normal vector pointing into the fluid. Over the past few decades, experiments and molecular dynamics simulations have indeed confirmed that slip occurs in pressure-driven flows over smooth solvophobic surfaces, with slip lengths of the order of nanometres (Vinogradova Reference Vinogradova1999; Baudry et al. Reference Baudry, Charlaix, Tonck and Mazuyer2001). In this work, we investigate the effect of possible surface non-homogeneities on flow and pressure fields, and allow both types of slip velocities to vary along the surface,

Integrating (2.10a ) in the normal direction and utilizing the boundary conditions (2.14a ) and (2.14b,c ), we obtain

It is important to note that in (2.15) only one boundary condition can be used at a time, i.e. for Helmholtz–Smoluchowski slip velocity $\unicode[STIX]{x1D6FD}$ must be set to zero, while for the Navier slip condition $\boldsymbol{u}_{w}^{(0,h)}$ must be set to zero. Utilizing the kinematic condition $D(z-h(\boldsymbol{r},t))/Dt=0$ (where $D$ denotes the substantial derivative), integrating (2.3) along the normal coordinate and utilizing the no-penetration boundary condition yields

Here, the depth average is defined as $\langle \cdots \rangle =\int _{0}^{h}(\cdots )\,\text{d}z/\int _{0}^{h}\,\text{d}z$ . Similarly, integrating each of the relations (2.15) and utilizing (2.14), we obtain

where $\langle \boldsymbol{u}_{w}\rangle =(\boldsymbol{u}_{w}^{(0)}+\boldsymbol{u}_{w}^{(h)})/2$ . Substituting the mean velocity defined in (2.17) into (2.16) and subsequently substituting the pressure from (2.4) yields the nonlinear thin-film approximation,

It should be noted that in neglecting the second-order time derivatives in (2.4) we have inherently assumed that relevant processes occur over time scales much longer than $\unicode[STIX]{x1D70F}_{EL}$ . Assuming that the deformation is small, i.e. $\unicode[STIX]{x1D702}\ll h_{0}$ , yields the linearized sixth-order diffusion equation,

where $\unicode[STIX]{x1D6FC}=h_{0}^{2}/(12\unicode[STIX]{x1D707})+\unicode[STIX]{x1D6FD}h_{0}/(2\unicode[STIX]{x1D707})$ . Sixth-order diffusion equations have been employed in the study of various subjects, such as magmal intrusion under a terrestrial crust (Michaut Reference Michaut2011) and interfacial instabilities under elastic membranes (Al-Housseiny et al. Reference Al-Housseiny, Christov and Stone2013). From (2.19), the characteristic visco-elastic time scale, $\unicode[STIX]{x1D70F}_{VE}$ , is given by

We here examine configurations in which $\unicode[STIX]{x1D70F}_{EL}\ll \unicode[STIX]{x1D70F}_{VE}$ , thus focusing our analysis on visco-elastic effects rather than propagation of elastic waves. Equation (2.20) describes the time scale for the propagation of visco-elastic deformation $r_{0}$ . However, if the system is actuated on a forced time scale, $\unicode[STIX]{x1D70F}_{f}$ , its characteristic length scale of visco-elastic dynamics is determined from (2.20), $r_{0}=(\unicode[STIX]{x1D70F}_{f}\unicode[STIX]{x1D6FC}Dh_{0})^{1/6}$ . In addition, if the forced time scale is smaller than $\unicode[STIX]{x1D70F}_{VE}$ , one must ensure that the associated $\mathit{Wo}$ number remains sufficiently small.

3 Effect of local non-uniformities in a rigid Hele-Shaw cell

We seek to study the effect of a local variation in the fluid properties on the pressure distribution and flow field due to a combined effect of non-uniform viscosity in the bulk and non-uniform slip velocity on the surface, both described by the following axially symmetric and time-independent functions:

which are centred at a common point, chosen as the origin of our coordinate system. Here, $H$ stands for a Heaviside step function, $R$ is the distance from the origin at which both quantities experience a transition, and $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}^{(in)}-\unicode[STIX]{x1D707}^{(out)}$ and $\unicode[STIX]{x0394}\boldsymbol{u}^{(w)}$ stand respectively for the difference between the values of viscosity and slip velocity in the inner $(r<R)$ and outer $(r>R)$ regions. Assuming furthermore that only the inner region hosts a Helmholtz–Smoluchowski boundary condition, (2.12), the difference $\unicode[STIX]{x0394}\boldsymbol{u}^{(w)}$ between the inner and outer slip velocities is given by

While sharp transitions of surface non-homogeneities are expected to introduce deviations from the lubrication approximation, these will be confined to a narrow spatial region of order $\unicode[STIX]{x1D700}$ around the transitions (Brotherton & Davis Reference Brotherton and Davis2004). In our analysis, we overlook these and other deviations due to sharp transition of parameters in the bulk.

For the case of a static elastic plate, the continuity equation (2.16) translates into an incompressibility condition of the two-dimensional depth-averaged flow,

Taking advantage of (3.3), we introduce a stream function $\unicode[STIX]{x1D713}$ defined via $\langle \boldsymbol{u}\rangle =\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D713}\hat{z})$ and derive a higher-order differential equation for the pressure by applying the two-dimensional divergence operator to (2.17). Similarly, dividing (2.17) by $\unicode[STIX]{x1D6FC}$ and applying the normal component of the curl operator eliminates the pressure. The resultant uncoupled Poisson equations for the pressure and the stream function are given by

Table 1 in appendix A presents the emergent Poisson equation for several distinct cases.

In the particular case of constant $\unicode[STIX]{x1D6FC}$ and vanishing body force, the derived equations (3.4) reduce to two Poisson equations which govern the case of non-uniform slip velocity, derived in Boyko et al. (Reference Boyko, Rubin, Gat and Bercovici2015). It should be noted that each equation of (3.4) for pressure and stream function is formally equivalent to Gauss’ law with a position-dependent dielectric function, where the pressure and stream function play the role of electrostatic potential, and the right-hand sides of (3.4) are interpreted as source terms. Thus, similarly to the case of non-uniform permittivity in electrostatics (Landau et al. Reference Landau, Bell, Kearsley, Pitaevskii, Lifshitz and Sykes1984), the matching conditions across the common boundary for the pressure and stream function in the case of a piecewise constant distribution of $\unicode[STIX]{x1D6FC}$ or $\langle \boldsymbol{u}_{w}\rangle$ are derived by line integration along an infinitesimal loop across the boundary. Specifically, these matching conditions are constituted by continuity of the potential and its derivative along the tangential direction (denoted by  $r_{t}$ ),

and discontinuity of its derivative along the normal direction (denoted by $r_{n}$ ), given by where the superscripts $(in)$ and $(out)$ correspond to the two adjacent regions.

The governing equation for the pressure for the case of axially symmetric non-homogeneities in both the viscosity and the Helmholtz–Smoluchowski slip velocity ( $\unicode[STIX]{x1D6FD}=0$ and $\boldsymbol{f}=0$ ) is explicitly written as

where $\unicode[STIX]{x1D711}$ is the polar angle coordinate and $\unicode[STIX]{x1D6FE}=12E\unicode[STIX]{x1D716}/h^{2}$ . The term on the right-hand side suggests a dipole type solution, $p(r,\unicode[STIX]{x1D711})=\unicode[STIX]{x1D6FC}(r)\cos (\unicode[STIX]{x1D711})$ , which upon substitution into (3.7), together with (3.1) and (3.2), leads to an ordinary differential equation for $\unicode[STIX]{x1D6FC}(r)$ , which easily integrates and leads towards the following solution in each of the regions:

where $a^{(in/out)}$ , $b^{(in/out)}$ are yet to be determined coefficients. Regularity at the origin dictates $a^{(in)}=0$ , while $b^{(out)}=\unicode[STIX]{x0394}p/\unicode[STIX]{x0394}\ell$ can be identified as a mean pressure gradient applied along the cell, as inferred from boundary conditions at $x\rightarrow \pm \infty$ . An additional relation between the coefficients can be taken as continuity of pressure or discontinuity of its derivative along the radial direction on the curve $r=R$ , constituted respectively by (3.5a,b ) and (3.6a ), leading to the following expressions for the dipole strength, $a^{(out)}$ , and pressure gradient in the inner region, $b^{(in)}$ :

where $m=\unicode[STIX]{x1D707}^{(in)}/\unicode[STIX]{x1D707}^{(out)}$ . Specific expressions for the pressure and the flow field can be obtained by substituting these coefficients, (3.9b ), into (2.17) and (3.8). Figure 2 presents the pressure and flow fields for several cases of spatial non-uniformities. Figures 2(a) and 2(b) present the cases of increasing ( $m>1$ ) or decreasing ( $m<1$ ) the viscosity in the bulk, with vanishing slip velocity throughout. Under a global pressure gradient, the streamlines are repelled from the central region for $m>1$ and converge for $m<1$ (as could be anticipated due to the formal analogy between position-dependent $1/\unicode[STIX]{x1D707}$ in Stokes flow and position-dependent dielectric constant in an electrostatic problem). In the limit $m\gg 1$ , the radial component of the flow field in the outer region, which is provided in appendix B, vanishes on the surface $r=R$ , and the velocity in the inner region tends to zero. Substituting $m=1$ into (3.9b ) corresponds to a case of uniform viscosity, and reproduces the electro-osmotic dipole obtained in Boyko et al. (Reference Boyko, Rubin, Gat and Bercovici2015) with closed streamlines shown in figure 2(c). In appendix C, we discuss the use of the Green’s function method, which allows us to solve (3.4) for more general cases.

Figure 2.

Analytical results showing the effect of a disk-shaped non-homogeneity on the resulting pressure distribution (normalized grey-scale colourmap) and flow field (white streamlines). (a,b) The ratio of internal to external viscosity, $m$ , affects the pressure distribution in the plane. For $m>1$ , (a), the density of isobars increases within the disk and decreases outside of it. For $m<1$ , (b), the opposite occurs and gradients are smaller within the disk and larger outside of it. In a complementary manner, for $m>1$ , streamlines are repelled from the disk region, whereas $m<1$ shows a focusing effect, with streamlines attracted to the disk. (c) A mismatch in an axial body force, or in EOF slip velocities, between the inner and outer regions of the disk results in a dipole which produces high and low pressures at the edges of the disk.

4 Static elastic deformations

Consider static or quasi-static processes characterized by time scales much larger than the characteristic visco-elastic time $r_{0}^{6}/\unicode[STIX]{x1D6FC}Dh_{0}$ , defined in (2.20). The corresponding governing equation for such processes then allows us to neglect partial derivatives with respect to time in (2.19), and for the particular case of non-uniform Helmholtz–Smoluchowski slip velocity it takes the following form:

Since the spatial derivative and the source terms in (2.19) both depend on $1/\unicode[STIX]{x1D707}$ , the relation (4.1) between the static deformation and the corresponding zeta potential does not depend on viscosity. For an arbitrary zeta-potential distribution, the deformation field obtained from (4.1) is not necessarily localized and may create deformations with the characteristic length scale of the configuration $L$ . Nevertheless, we may require a localized deformation field and calculate the corresponding zeta potential. For an electric field $E$ that is homogeneous and oriented along the $x$ direction, we utilize (4.1) to express the zeta potential as a function of the deformation $\unicode[STIX]{x1D702}$ via

The integration introduces some function of $y$ , which represents the effects of boundary conditions. In our problem, we assume that the boundaries are sufficiently far from the region of interest, and therefore set this function to zero. Hereafter, we consider the set-up described by figure 1, where the non-homogeneous zeta potential, $\unicode[STIX]{x1D701}$ , resides on the bottom rigid plate. Thus, the mean zeta potential, $\langle \unicode[STIX]{x1D701}\rangle$ , relates to the zeta potential on the bottom plate, $\unicode[STIX]{x1D701}$ , via $\langle \unicode[STIX]{x1D701}\rangle =\unicode[STIX]{x1D701}/2$ .

Figure 3.

Analytical results presenting the zeta potential required for generating a desired Gaussian deformation, and the resulting pressure distribution. Panels (a) and (b) present the 1D and 2D cases respectively. In both cases, a central region of positive pressure is bracketed by regions of negative pressure to ensure that the deformation remains localized.

Let us determine the pressure and zeta-potential distribution necessary to generate a local deformation of a Gaussian shape,

Substituting (4.3) into (2.4) yields an expression for $p$ , while utilizing (4.2) yields the following closed form expression for the necessary zeta-potential distribution:

Figures 3(a) and 3(b) respectively present the one- and two-dimensional zeta-potential distributions necessary to generate the corresponding Gaussian deformation. It should be noted that in both cases the central region with positive pressure distribution is surrounded by regions with negative pressure, which ensure that the total deformation decays to zero on a length scale  $\unicode[STIX]{x1D70E}$ .

To determine the zeta-potential distribution necessary to generate an arbitrary deformation, $f(\boldsymbol{r})$ , we take advantage of the linearity of our problem and the finite Gaussian representation introduced in Gabor (Reference Gabor1946), which allows us to approximate any surface as a linear combination of Gaussians centred at different points via $f(\boldsymbol{r})\simeq \sum _{i=1}^{N}c_{i}g(\boldsymbol{r};\boldsymbol{r}_{i})$ , where the index $i$ identifies a Gaussian centred at a point $\boldsymbol{r}_{i}$ . The coefficients are fixed by demanding that the sum of the Gaussian terms is equal to the values of the function $f(\boldsymbol{r})$ at some $N$ auxiliary points, leading to the following system of $N$ algebraic equations for $N$ coefficients $c_{i}$ :

where for convenience we set these auxiliary points as the central points of each Gaussian. The approximation in a given domain becomes more accurate as we increase the number of Gaussian terms in the sum and decrease the distances between their central points. The corresponding zeta potential, $\langle \unicode[STIX]{x1D701}(\boldsymbol{r})\rangle$ , necessary to generate a desired deformation is given by

where each of the $\langle \unicode[STIX]{x1D701}(\boldsymbol{r};\boldsymbol{r}_{i})\rangle$ terms in the sum above corresponds to the zeta potential needed to maintain a single Gaussian deformation, given by (4.4).

For illustration, consider an example for the one- and the two-dimensional cases, where we determine the zeta potential necessary to create desired deformations which have the shape of a human face and the African continent respectively. Figure 4(a) presents the shape of a human face created by a sum of 45 Gaussians, while figure 4(b) presents the set of Gaussians multiplied by the corresponding coefficients $c_{i}$ . Figure 4(c) presents the resultant zeta-potential distribution obtained by summing the contribution needed to create each one of the Gaussians, by utilizing (4.6). Following the same steps, it is possible to obtain the two-dimensional zeta-potential distribution necessary to create arbitrary deformation profiles in 2D. Figure 4(d) presents the result of summing 400 Gaussians, while figure 4(e) presents the corresponding zeta potential obtained by inserting the zeta potential needed to support each one of Gaussians, (4.4), into the sum (4.6).

Figure 4.

Analytical results demonstrating the calculation of zeta-potential distributions required in order to create complex deformations. (a) A smooth curve of a human face profile, which is generated by superposition of 45 Gaussians, presented in (b). (c) The total zeta-potential distribution required to generate this deformation. (d) The sum of 400 Gaussians ( $20\times 20$ array) is used to describe a deformation corresponding to the topography of the African continent. (e) The zeta-potential distribution (colourmap) in the $x{-}y$ plane required to achieve this deformation.

5 Time-dependent elastic deformations

Consider the general elastic deformations equation, (2.19), where the only spatial non-homogeneity is in the Helmholtz–Smoluchowski slip velocity (i.e. $\boldsymbol{f}=0$ and $\unicode[STIX]{x1D6FD}=0$ ), and the source term is then the component of the zeta-potential gradient along the electric field. For convenience, we scale to dimensionless variables via $\boldsymbol{r}\rightarrow r_{0}\boldsymbol{r}$ , $\unicode[STIX]{x1D702}\rightarrow \unicode[STIX]{x1D702}_{0}\unicode[STIX]{x1D702}$ , $t\rightarrow \unicode[STIX]{x1D70F}_{VE}t$ , $E\rightarrow E_{0}E$ and $\unicode[STIX]{x1D701}\rightarrow \unicode[STIX]{x1D701}_{0}\unicode[STIX]{x1D701}$ . The equation then takes the form of a sixth-order diffusion equation,

where $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D716}_{0}u_{0}\unicode[STIX]{x1D70F}_{VE}/\unicode[STIX]{x1D702}_{0}$ , $\unicode[STIX]{x1D716}_{0}=h_{0}/r_{0}$ , $u_{0}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}_{0}E_{0}/\unicode[STIX]{x1D707}$ is the typical slip velocity, and $\unicode[STIX]{x1D70F}_{VE}$ is defined by (2.20). Similarly to (4.2), an explicit expression for the zeta potential required to generate pre-assigned spatio-temporal elastic deformations, $\unicode[STIX]{x1D702}(\boldsymbol{r},t)$ , can be obtained by expressing the corresponding zeta potential, $\unicode[STIX]{x1D701}(\boldsymbol{r},t)$ , from the right-hand side of (5.1) through the following integral:

where we have assumed that the external electric field is uniform and is applied along the $x$ -axis. The resultant expression for $\unicode[STIX]{x1D701}(\boldsymbol{r},t)$ in such cases is a sum of two terms. The first term expresses the zeta potential required to sustain the rate of deformation and the second term represents the zeta potential of the static problem (i.e. $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}t=0$ ), given by (4.2).

6 Concluding remarks

In this work, we studied non-uniform lubrication flows in the gap between a rigid plate and a parallel elastic plate as a mechanism to create time-varying deformation fields of the elastic plate. The approach we presented can be relevant to various applications such as adaptable optics, soft actuators and lab-on-a-chip devices. It is instructive to estimate the resultant deformation for realistic parameters. For actuation due to EOF-driven Helmholtz–Smoluchowski slip velocity, and using typical values of the zeta potential of $\unicode[STIX]{x1D701}=50~[\text{mV}]$ (i.e. nearly twice the thermal voltage $V_{T}\sim 25.6~[\text{mV}]$ ), the plate Young modulus $Y=1~[\text{mPa}]$ corresponding to polydimethylsiloxane (PDMS) (Cheng et al. Reference Cheng, Sun, Meininger and Almasri2010), the magnitude of the applied electric field $E=10^{4}~[\text{V}~\text{m}^{-1}]$ , the gap between the plates $h_{0}=50~[\unicode[STIX]{x03BC}\text{m}]$ , the plate thickness $w=10~[\unicode[STIX]{x03BC}\text{m}]$ and the dielectric constant $\unicode[STIX]{x1D716}=80\unicode[STIX]{x1D700}_{0}$ (where $\unicode[STIX]{x1D700}_{0}$ is the dielectric permittivity of the vacuum), we obtain a deformation $\unicode[STIX]{x1D702}$ of the order of $40~[\unicode[STIX]{x03BC}\text{m}]$ distributed over a region of width $300~[\unicode[STIX]{x03BC}\text{m}]$ . Naively, from the linearity of (5.1), one could expect that an order of magnitude increase of the zeta potential would lead to a similar increase in the deformation. However, as can be directly obtained by the Boltzmann distribution for the monovalent ion concentration in the EDL, $c_{EDL}/c_{0}=\text{e}^{\unicode[STIX]{x1D701}/V_{T}}$ , together with the fact that the ions have a finite size, there exists an upper bound for $\unicode[STIX]{x1D701}$ above which it is not possible to pack additional charges into the EDL. For example, a bulk solution concentration of $c_{0}=1~[\text{mM}]$ would reach the critical packing limit on the surface at approximately $\unicode[STIX]{x1D701}=350~[\text{mV}]$ , and is expected to introduce corrections even at lower values of the zeta potential (Bazant et al. Reference Bazant, Kilic, Storey and Ajdari2009). Indeed, Van Der Wouden et al. (Reference Van Der Wouden, Heuser, Hermes, Oosterbroek, Gardeniers and Van Den Berg2005) reported saturation in EOF velocity for sufficiently high zeta potentials.

In §§ 4 and 5, we examined static and asymptotically static deformations. Such deformations are achieved by requiring a steady deformation field within a defined finite region and then calculating the corresponding zeta-potential distribution. In such configurations, the deformations vanish outside of the deformed region, and thus boundary conditions can be applied near the deformed region without affecting the solution. Thus, the requirement $\ell \ll L$ (stating that the configuration characteristic length $L$ is much greater than the characteristic length $\ell$ of the deformed region) is unnecessary for static and asymptotically static spatially localized solutions.

Our analysis of localized steady zeta-potential actuation yielded a solution that grows indefinitely with time, which indicates that these solutions are no longer valid after a certain time, as they fail to satisfy the underlying assumptions of our model. In practice, there are several stabilizing mechanisms such as physical walls and the liquid’s own weight, which impose restoring forces on the elastic plate and are expected to prevent the diverging behaviour, though investigation of such models lies beyond the scope of this study. Within the small-deformation regime used in this work, the static solutions are valid provided that $\unicode[STIX]{x1D702}_{L}/\unicode[STIX]{x1D702}_{\ell }\sim p_{L}L^{4}p_{\ell }\ell ^{4}\ll 1$ , where the subscript $\ell$ indicates the deformation and pressure associated with the localized solution, whereas $L$ indicates quantities associated with the perturbation field over the entire configuration.

In this work, we focused on deformations that were small compared with the gap between the plates, which allowed us to consider one-way coupling between the fluid pressure and the elastic deformation of the plate. From our order of magnitude analysis, it appears feasible to obtain pressures that would result in deformations of a few tens of microns. In principle, it is possible to achieve larger deformations by utilizing a thinner elastic plate (at least in the Kirchhoff–Love model the deformation depends on the thickness in a power of three) and/or utilizing more compliant materials which may become more available with the advance of materials science (e.g. materials with Young modulus $Y\sim 1~[\text{kPa}]$ , reported in Moraes et al. (Reference Moraes, Labuz, Shao, Fu and Takayama2015)). The modelling of large deformations is an interesting extension of this work where one must consider changes in the flow field and the electric field due to deformations, as well as account for additional elastic forces such as the tension in the plate (which could be captured by the Föppl–von Karman model), and also body forces such as gravity, which are expected to have an additional stabilizing effect on the dynamics of the plate.