2.1. Governing equations

The fluid motion satisfies the continuity and Navier–Stokes equations given by

(2.1a,b)\begin{equation} {\boldsymbol{\nabla}}\boldsymbol{\cdot}{\boldsymbol{v}} = 0 \quad \mbox{and} \quad \rho\left(\frac{\partial{\boldsymbol{v}}}{\partial t} + {\boldsymbol{v}}\boldsymbol{\cdot}{\boldsymbol{\nabla}}{\boldsymbol{v}}\right)= {\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{T}+\rho g \boldsymbol{i}_{z}, \end{equation}

where $\boldsymbol {T}=-p\boldsymbol {I}+\mu ({\boldsymbol {\nabla }}{\boldsymbol {v}}+ {\boldsymbol {\nabla }}{\boldsymbol {v}}^{\textrm {T}})$ is the stress tensor, $\boldsymbol {i}_{z}$ is the unit vector along the positive $z$ direction and $I$ is the identity tensor. The interface speed, $\mathcal {U}$, the unit normal vector, $\boldsymbol {n}$, and the unit tangent vector, $\boldsymbol {t}$, are given by

(2.2a–c)\begin{align} \mathcal{U} = \frac{\dfrac{\partial h}{\partial t}}{\left[1+\left(\dfrac{\partial h}{\partial x}\right)^{2}\right]^{1/2}},\quad \boldsymbol{n}= \frac{-\dfrac{\partial h}{\partial x}i_{x}+i_{z}}{\left[1+\left(\dfrac{\partial h}{\partial x}\right)^{2}\right]^{1/2}}\quad \mbox{and} \quad \boldsymbol{t}= \dfrac{i_{x}+\dfrac{\partial h}{\partial x}i_{z}}{\left[1+ \left(\dfrac{\partial h}{\partial x}\right)^{2}\right]^{1/2}}. \end{align}

At the wavy wall, no slip and no penetration are satisfied. In addition, the free surface which is at $z=h(x,t)$ is a material surface and thus ${\boldsymbol {v}} \boldsymbol {\cdot } {\boldsymbol {n}}- \mathcal {U}= 0$ must hold. The interfacial force balances at the surface are given by

(2.3a,b)\begin{equation} {\boldsymbol{n}} \boldsymbol{\cdot} {\boldsymbol{T}}\boldsymbol{\cdot}{\boldsymbol{t}}= 0 \quad \mbox{and} \quad{\boldsymbol{n}}\boldsymbol{\cdot}{\boldsymbol{T}} \boldsymbol{\cdot}{\boldsymbol{n}}={-}\gamma {\boldsymbol{\nabla}} \boldsymbol{\cdot} {\boldsymbol{n}}. \end{equation}

To analyse the effect of the wavy-wall amplitude, the modelling equations are investigated in the long-wave limit.

2.2. The long-wave model and boundary-layer equations

The instability and subsequent dynamics are investigated by employing a long-wave approximation. Thus we take $H\ll \lambda$, where $H$ is the thickness of the fluid layer and $\lambda$ is a characteristic horizontal length scale such as the wavelength of the disturbance. The long-wave model is obtained via a separation of length scales. To do so the governing equations are made dimensionless by using the following scales, denoted by the subscript ‘$c$’:

(2.4a–f)\begin{equation} x_{c}=\lambda,\quad z_{c}=H,\quad u_{c}=U,\quad w_{c}=\epsilon U,\quad t_{c}=\frac{\lambda}{U},\quad p_{c}=\frac{\mu U}{H}. \end{equation}

Here $U$ is a characteristic velocity scale, and $\epsilon$, the long-wave parameter (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011), is defined by $\epsilon =H/\lambda$. The wavy-wall function is represented by $f(x)=\mathcal {A}\cos (2 k_{w} x)$, where $k_{w}={{\rm \pi} }/{W}$ and $\mathcal {A}$ is the amplitude of the wavy wall. Note that $W$ and $\mathcal {A}$ are scaled quantities, having been scaled with the thickness of the liquid film. The boundary-layer assumption is invoked, i.e. $\epsilon < 1$, and upon retaining terms of ${O}(\epsilon )$, the non-dimensional model, without changing notation for the dependent variables, is given by

(2.5)\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z} = 0, \end{gather}

(2.6)\begin{gather} \epsilon {Re}\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+w \frac{\partial u}{\partial z}\right) ={-} \epsilon \frac{\partial p}{\partial x}+ \frac{\partial^{2} u}{\partial z^{2}},\quad \mbox{where}\ Re=\frac{\rho U H}{\mu}, \end{gather}

and

(2.7)\begin{equation} - \epsilon \frac{\partial p}{\partial z}+\epsilon G=0,\quad \mbox{where}\quad G=\frac{\rho g H^{2}}{\mu U}. \end{equation}

The above equations result when considering $Re$ and $G$ to be at least ${O}(1)$. At the free surface $z=h(x,t)$, kinematic and tangential force balance conditions give

(2.8a,b)\begin{equation} \epsilon \frac{\partial h}{\partial t}={-}\epsilon u\frac{\partial h}{\partial x}+\epsilon w\quad \mbox{and} \quad \frac{\partial u}{\partial z}=0, \end{equation}

and the normal force balance at $z=h(x,t)$ results in

(2.9)\begin{equation} -\epsilon p=\frac{\epsilon^{3}}{Ca}\frac{\partial^{2} h}{\partial x^{2}},\quad \mbox{where}\ Ca=\frac{\mu U}{\gamma}. \end{equation}

The equations (2.5)–(2.9) are called the $1+\epsilon$ model (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011). The capillary number, $Ca$, is assumed to be no higher than ${O}(\epsilon ^{2})$, so that the effect of curvature at the order of interest can be seen. In the current problem, $Ca$ is set equal to $\epsilon ^{3}$, thereby also setting the velocity scale. Other choices for the velocity scale are possible; in this work the choice is motivated by computations for the fluids of interest, e.g. water/air.

Integration of the vertical component of the momentum balance, i.e. (2.7), from the bulk of the fluid to the free surface, $h(x,t)$, is performed, and the resulting equation is substituted into the horizontal component of the momentum balance, i.e. (2.6), to express pressure in the bulk of the fluid in terms of pressure at the free surface. This yields

(2.10)\begin{equation} \epsilon Re\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+w \frac{\partial u}{\partial z}\right) ={-} \epsilon\left.\frac{\partial p}{\partial x}\right\vert_{h}+\frac{\partial^{2} u}{\partial z^{2}}+\epsilon G\frac{\partial h}{\partial x}. \end{equation}

The normal force balance equation (2.9) is differentiated in the horizontal direction and combined with (2.10) to give

(2.11)\begin{equation} \epsilon Re\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+w \frac{\partial u}{\partial z}\right) = \frac{\epsilon^{3}}{Ca} \frac{\partial^{3} h}{\partial x^{3}}+\frac{\partial^{2} u}{\partial z^{2}}+\epsilon G\frac{\partial h}{\partial x}. \end{equation}

To arrive at the evolution equations, the weighted-residual integral boundary-layer (WRIBL) technique is employed (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015; Dietze et al. Reference Dietze, Picardo and Narayanan2018; Pillai & Narayanan Reference Pillai and Narayanan2018) by first decomposing $u$ into an ${O}(1)$ and ${O}(\epsilon )$ part, i.e.

(2.12)\begin{equation} u(x,z,t)=\underbrace{\hat{u}(x,z,t)}_{{O}(1)}+\underbrace{\tilde{u}(x,z,t)}_{{O}(\epsilon)}, \end{equation}

and then introducing this into (2.11). A Galerkin projection of the resulting equation at ${O}(\epsilon )$ with a suitable weight function, $F$, yields the dynamic evolution equations that only depend on the surface position $h(x,t)$. Thus we first get

(2.13)\begin{align} \int_{f}^{h} \epsilon Re\left(\frac{\partial \hat{u}}{\partial t}+\hat{u} \frac{\partial \hat{u}}{\partial x}+\hat{w} \frac{\partial \hat{u}}{\partial z}\right) F\,\textrm{d} z= \int_{f}^{h} \frac{\partial^{2} \hat{u}}{\partial z^{2}} F\,\textrm{d} z+ \left(\epsilon G\frac{\partial h}{\partial x}+\frac{\epsilon^{3}}{Ca} \frac{\partial^{3} h}{\partial x^{3}} \right) \int_{f}^{h} F\,\textrm{d} z. \end{align}

We then choose the leading-order velocity $\hat {u}(x,z,t)$ to be parabolic along the horizontal direction (cf. Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011). It is determined so that the following conditions are satisfied:

(2.14a–d)\begin{equation} \frac{\partial^{2}\hat{u}}{\partial z^{2}}=K,\quad \int_{f(x)}^{h(x,t)} \hat{u} \,\textrm{d} z=q,\quad \hat{u}\vert_{f(x)}=0 \quad \mbox{and}\quad \left.\frac{\partial \hat{u}}{\partial z}\right\vert_{h(x,t)}=0, \end{equation}

where $K$ is obtained in terms of the flow rate, $q$, using the integral condition. Solving (2.14a–d) for $\hat {u}$, we get

(2.15)\begin{equation} \hat{u}=\frac{3 z^2 q}{2 (\,f-h)^3}+\frac{3 z h q}{(h-f)^3}-\frac{3 q (\,f^2-2 f h)}{2 (\,f-h)^3}. \end{equation}

The weight function, $F$, used in the Galerkin projection in (2.13) has the same functional form as the leading-order velocity, $\hat {u}$. It is determined by solving

(2.16a–c)\begin{equation} \frac{\partial^{2}F}{\partial z^{2}}=1,\quad F\vert_{f(x)}=0 \quad \mbox{and}\quad \left.\frac{\partial F}{\partial z}\right\vert_{h(x,t)}=0, \end{equation}

and is given by

(2.17)\begin{equation} F=\frac{1}{2} \left(2 f h-f^2\right)-z h+\frac{z^2}{2}. \end{equation}

The vertical component of velocity, $w$, required in (2.13) is obtained from the continuity equation (2.5). This, along with expressions for $\hat {u}$ and $F$, is then substituted into (2.13) to obtain the evolution equations in terms of $h(x,t)$ and $q(x,t)$ alone, i.e.

(2.18)\begin{align} & \epsilon Re\left[\left(\frac{-2}{5}(\,f-h)^{2} \frac{\partial q}{\partial t} +\frac{23}{40}(\,f-h)q \frac{\partial q}{\partial x}\right)+\frac{3q}{280} \left(q \left(67\frac{\partial h}{\partial x}-32\frac{\partial f}{\partial x}\right)+48 (\,f-h) \frac{\partial q}{\partial x}\right) \vphantom{\frac{3 q\dfrac{\partial q}{\partial x} (\,f-h) \left({-}39 f^2 h-66 f h^2+22 h^3+48 f^3\right)}{560(\,f-h)^3}}\right.\nonumber\\ &\qquad +\frac{3 q^{2} \left(\dfrac{\partial f}{\partial x} ({-}9 f^2 h+114 f h^2+32 h^3-32 f^3)+\dfrac{\partial h}{\partial x} ({-}96 f^2 h-114 f h^2+38 h^3+67 f^3)\right)}{560(\,f-h)^3} \nonumber\\ &\left.\qquad +\frac{3 q\dfrac{\partial q}{\partial x} (\,f-h) \left({-}39 f^2 h-66 f h^2+22 h^3+48 f^3\right)}{560(\,f-h)^3}\right] \nonumber\\ &\quad =\frac{-1}{3}(h-f)^{3} \left(\epsilon G\frac{\partial h}{\partial x}+\frac{\epsilon^{3}}{Ca}\frac{\partial^{3} h}{\partial x^{3}}\right)+q. \end{align}

Equation (2.18) is used in conjunction with the following equation, obtained by integrating the continuity equation across the fluid layer:

(2.19)\begin{equation} \frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=0, \end{equation}

where Leibniz's integration rule and the kinematic condition are employed. Note that in (2.18), $f$ may be represented by $f=\mathcal {A}\cos (2k_{w}x)$, $k_{w}={{\rm \pi} }/{W}$, for the case of two full waves over a wall length of $2W$. If the wall is given by $\mathcal {A}\cos (4k_{w}x)$, $k_{w}={{\rm \pi} }/{W}$, we have four full waves over a length of $4W$, and so forth. It should be observed that in (2.18), the left-hand side consists of the inertial terms. On the right-hand side the first term arises due to shape of the wavy wall, surface curvature and gravitational effects, while the second term is attributed to transverse viscous diffusion. In the creeping flow limit, the evolution equations (2.18) and (2.19) reduce to a simple equation in terms of $h(x,t)$ alone:

(2.20)\begin{equation} \frac{\partial h}{\partial t}={-}\frac{1}{3}\frac{\partial }{\partial x}\left[\left(h-f\right)^{3} \left(\epsilon G\frac{\partial h}{\partial x}+\frac{\epsilon^{3}}{Ca}\frac{\partial^{3} h}{\partial x^{3}}\right)\right]. \end{equation}

Observe that in the creeping flow limit, the shape of the wavy wall, surface curvature and gravitational effects alone play a key role in the dynamics of the free surface. An analogous procedure for two active fluids with a common interface system may also be employed. This is done in detail in the supplementary material. We now turn to the results obtained from calculations using the reduced-order model, i.e. (2.18) and (2.19), and (2.20), as well as the analogous equations for the two-phase system.