Instability and dripping of electrified liquid films flowing down inverted substrates

TL;DR

This study investigates how electric fields influence the stability and dripping behavior of thin liquid films on inverted substrates, demonstrating that electric fields can effectively suppress dripping through linear stability mechanisms.

Contribution

The paper introduces a combined linear stability and numerical simulation approach to predict electric field thresholds for dripping suppression in electrified liquid films.

Findings

Electric fields can significantly delay or prevent dripping.

Linear stability analysis accurately predicts the electric field strength needed for suppression.

Numerical simulations confirm the theoretical predictions.

Abstract

We consider the gravity-driven flow of a perfect dielectric, viscous, thin liquid film, wetting a flat substrate inclined at a non-zero angle to the horizontal. The dynamics of the thin film is influenced by an electric field which is set up parallel to the substrate surface - this nonlocal physical mechanism has a linearly stabilizing effect on the interfacial dynamics. Our particular interest is in fluid films that are hanging from the underside of the substrate; these films may drip depending on physical parameters, and we investigate whether a sufficiently strong electric field can suppress such nonlinear phenomena. For a non-electrified flow, it was observed by Brun et al. (Phys. Fluids 27, 084107, 2015) that the thresholds of linear absolute instability and dripping are reasonably close. In the present study, we incorporate an electric field and analyse the absolute/convective instabilities of a hierarchy of reduced-order models to predict the dripping limit in parameter space. The spatial stability results for the reduced-order models are verified by performing an impulse-response analysis with direct numerical simulations (DNS) of the Navier-Stokes equations coupled to the appropriate electrical equations. Guided by the results of the linear theory, we perform DNS on extended domains with inflow/outflow conditions (mimicking an experimental set-up) to investigate the dripping limit for both non-electrified and electrified liquid films. For the latter, we find that the absolute instability threshold provides an order-of-magnitude estimate for the electric field strength required to suppress dripping; the linear theory may thus be used to determine the feasibility of dripping suppression given a set of geometrical, fluid and electrical parameters.

Tables

Table 1. Table 1 : Physical parameters for Castor Oil (Hänseler AG).

$\rho$	$\nu$	$\mu$	$\sigma$
(kg/m3)	(m2/s)	(kg/m$\cdot$s)	(N/m)
$946$	$9.14\times {10}^{-4}$	$0.865$	$0.0339$

Table 2. Table 2 : Critical E 0 subscript 𝐸 0 E_{0} ( V / m ) for A/C transition with θ = 5 ​ π / 6 𝜃 5 𝜋 6 \theta=5\pi/6 and ℓ = 2.5 ℓ 2.5 \ell=2.5 mm , in the case of Castor Oil (Hänseler AG).

Model	Inertialess Benney	Benney	WIBL1	WIBL2	Stokes flow
$E_0$ (V/m)	$3.07\times {10}^5$	$3.45\times {10}^5$	$3.34\times {10}^5$	$2.48\times {10}^5$	$2.10\times {10}^5$

Table 3. Table 3 : Absolute instability regions and dripping limits for the non-electrified case (lengths are given in mm ).

$\theta$	$\theta -90{(}^{\circ })$	Inertialess Benney	Benney	WIBL1	WIBL2	Stokes flow	${\ell }^{\text{drip}}$ (DNS)
$7\pi /8$	$67.5$	$1.52\le \ell$	$1.52\le \ell$	$1.52\le \ell \le 20.1$	$1.35\le \ell \le 8.76$	$1.36\le \ell$	$1.40$–$1.45$
$5\pi /6$	$60$	$2.19\le \ell$	$2.15\le \ell$	$2.16\le \ell \le 15.8$	$2.18\le \ell \le 5.21$	$2.25\le \ell \le 13.1$	$2.30$–$2.50$
$49\pi /60$	$57$	$2.51\le \ell$	$2.42\le \ell$	$2.44\le \ell \le 14.5$	—	—	$2.40$–$2.60$

Equations

$${\mathbf{u}}_t+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}+\mathbf{g},\nabla \cdot \mathbf{u}=0,$$

$${\nabla }^2{V}^{\iota }=0,\iota =\text{S},\text{I},\text{II}.$$

$$\nabla V^{\text{S}}\to -{\mathbf{E}}_0,\text{as }z\to -\infty ,\nabla V^{\text{II}}\to -{\mathbf{E}}_0,\text{as }z\to +\infty ,$$

$${\left[{ϵ}^{\iota }{V}_z^{\iota }\right]}_{\text{I}}^{\text{S}}=0,{\left[V^{\iota }\right]}_{\text{I}}^{\text{S}}=0,\text{at }z=0,$$

$$w=h_t+u{h}_x,$$

$${\left[{ϵ}^{\iota }\nabla V^{\iota }\cdot \mathbf{n}\right]}_{\text{II}}^{\text{I}}=0,{\left[V^{\iota }\right]}_{\text{II}}^{\text{I}}=0.$$

$${\left[\left({\mathbf{T}}^{\iota }\mathbf{n}\right)\cdot \mathbf{n}\right]}_{\text{II}}^{\text{I}}=\sigma h_{xx}/(1+h_x^2{)}^{3/2},{\left[\left({\mathbf{T}}^{\iota }\mathbf{n}\right)\cdot \mathbf{t}\right]}_{\text{II}}^{\text{I}}=0,$$

$${\mathbf{T}}_{jk}^{\text{I}}$$

$${\mathbf{T}}_{jk}^{\text{II}}$$

$$2\mu \left(u_x(h_x^2-1)-(u_z+w_x)h_x\right)+\frac{{ϵ}_0}{2}{\left[{ϵ}^{\iota }\left((h_x^2-1)((V_x^{\iota }{)}^2-(V_z^{\iota }{)}^2)-4h_x{V}_x^{\iota }{V}_z^{\iota }\right)\right]}_{\text{II}}^{\text{I}}$$

$$+(p_{\text{atm}}-p)(1+h_x^2)$$

$$(1-h_x^2)(u_z+w_x)+4w_z{h}_x$$

$$\left.\begin{array}[]{c}{\displaystyle\overline{h}=\ell,\qquad\overline{u}=\frac{g\sin\theta}{2\nu}(2\ell z-z^{2}),\qquad\overline{w}=0,}\\[8.0pt] {\displaystyle\overline{p}=p_{\textrm{atm}}-\frac{\epsilon_{0}E_{0}^{2}}{2}(\epsilon^{\textrm{I}}-\epsilon^{\textrm{II}})-\rho g(z-\ell)\cos\theta,\quad\overline{V}^{\iota}=-E_{0}x\quad\textrm{for $\iota=$ S, I, II}.}\end{array}\right\}$$

$$\mathrm{Re}=\frac{U_0\ell }{\nu }=\frac{g{\ell }^3\sin \theta }{2{\nu }^2},\mathrm{We}=\frac{{ϵ}_0{E}_0^2\ell }{2\mu U_0}=\frac{{ϵ}_0{E}_0^2}{\rho g\ell \sin \theta },\text{mbox}C=\frac{U_0\mu }{\sigma }=\frac{\rho g{\ell }^2\sin \theta }{2\sigma }.$$

$$\left.\begin{array}[]{c}{\displaystyle x^{*}=\frac{1}{\ell}x,\quad z^{*}=\frac{1}{\ell}z,\quad\bm{u}^{*}=\frac{1}{U_{0}}\bm{u},\quad t^{*}=\frac{U_{0}}{\ell}t,\quad h^{*}=\frac{1}{\ell}h,}\\[8.0pt] {\displaystyle p^{*}=\frac{\ell}{\mu U_{0}}\left(p-p_{\textrm{atm}}+\frac{\epsilon_{0}E_{0}^{2}}{2}(\epsilon^{\textrm{I}}-\epsilon^{\textrm{II}})+\rho gz\cos\theta\right),\quad(V^{\iota})^{*}=\frac{1}{E_{0}\ell}(V^{\iota}+E_{0}x)\quad\textrm{for }\iota=\textrm{S},\textrm{I},\textrm{II},}\end{array}\right\}$$

$$\mathrm{Re}({\mathbf{u}}_t+(\mathbf{u}\cdot \nabla )\mathbf{u})=-\nabla p+{\nabla }^2\mathbf{u}+2{\mathbf{e}}_x,\nabla \cdot \mathbf{u}=0.$$

$$\nabla V^{\text{S}}\to 0,\text{as }z\to -\infty ,\nabla V^{\text{II}}\to 0,\text{as }z\to +\infty ,$$

$${\left[{ϵ}^{\iota }(-h_x(V_x^{\iota }-1)+V_z^{\iota })\right]}_{\text{II}}^{\text{I}}=0.$$

$$u_x(h_x^2-1)-(u_z+w_x)h_x+h(1+h_x^2)\cot \theta -\frac{1}{2}p(1+h_x^2)$$

$$+\frac{\mathrm{We}}{2}{\left[{ϵ}^{\iota }\left((h_x^2-1)((V_x^{\iota }-1{)}^2-(V_z^{\iota }{)}^2)-4h_x(V_x^{\iota }-1)V_z^{\iota }+(1+h_x{)}^2\right)\right]}_{\text{II}}^{\text{I}}=\frac{1}{2\text{mbox}C}\frac{h_{xx}}{(1+h_x^2{)}^{1/2}}.$$

$$\overline{h}=1,\overline{u}=2z-z^2,\overline{w}=0,\overline{p}=0,{\overline{V}}^{\iota }=0\text{for ι= S, I, II}.$$

$$({\partial }_t,{\partial }_x,w)\to \delta ({\partial }_t,{\partial }_x,w),$$

$$H_t+{\left[\frac{2}{3}{H}^3+\delta \left[\frac{8\mathrm{Re}}{15}{H}^6{H}_x-\frac{1}{3}{H}^3{P}_x\right]\right]}_x=0.$$

$$h_t+f_x=0,f(x,t)={\int }_0^{h(x,t)}u(x,z,t)\mathrm{d}z,$$

$$H_t+F_x=$$

$$F+\delta \frac{2\mathrm{Re}}{5}{H}^2{F}_t=$$

$$+{\delta }^2\left[\frac{8}{5}{H}_x^{2}F-\frac{9}{5}H{H}_x{F}_x-\frac{12}{5}H{H}_{xx}F+\frac{9}{5}{H}^2{F}_{xx}\right].$$

$$h\cot \theta -\frac{1}{2}p+\frac{\mathrm{We}}{2}{\left[{ϵ}^{\iota }\left(2\delta V_x^{\iota }+(V_z^{\iota }{)}^2+4\delta h_x{V}_z^{\iota }\right)\right]}_{\text{II}}^{\text{I}}=\frac{{\delta }^2}{2\text{mbox}C}{h}_{xx}+O(\delta ),$$

$$h=h_0+\delta h_1+\dots ,V^{\iota }=V_0^{\iota }+\delta V_1^{\iota }+\dots ,\iota =\mathrm{I},\mathrm{II},\mathrm{S}.$$

$${\delta }^2{V}_{xx}^{\text{I}}+V_{zz}^{\text{I}}=0,V_{xx}^{\iota }+V_{\zeta \zeta }^{\iota }=0,\text{for }\iota =\text{S},\text{II}.$$

$${\hat{V}}^{\text{S}}={\hat{A}}^{\text{S}}(\xi )e^{\mathrm{sign}({\xi }_{\text{r}})\xi \zeta },{\hat{V}}^{\text{II}}={\hat{B}}^{\text{II}}(\xi )e^{-\mathrm{sign}({\xi }_{\text{r}})\xi \zeta },$$

$${\hat{V}}^{\text{I}}={\hat{C}}_0^{\text{I}}(\xi )z+{\hat{D}}_0^{\text{I}}(\xi )+\delta \left({\hat{C}}_1^{\text{I}}(\xi )z+{\hat{D}}_1^{\text{I}}(\xi )\right)+O({\delta }^2).$$

$${\hat{A}}_0^{\text{S}}(\xi )={\hat{D}}_0^{\text{I}}(\xi ),{\hat{A}}_1^{\text{S}}(\xi )={\hat{D}}_1^{\text{I}}(\xi ),{\hat{C}}_0^{\text{I}}(\xi )=0,{ϵ}^{\text{S}}\mathrm{sign}({\xi }_{\text{r}})\xi {\hat{A}}_0^{\text{S}}(\xi )-{ϵ}^{\text{I}}{\hat{C}}_1^{\text{I}}(\xi )=0.$$

Full text

I Introduction

Thin liquid films are encountered in numerous industrial applications such as coating processes, lab-on-a-chip systems Stone et al. (2004), and liquid film cooling Miyara (1999); Serifi et al. (2004), as well as geological and geomorphological applications Shorts et al. (2005); Camporeale (2017). Generally such flows can become unstable to inertial and gravitational instabilities and can support complicated wave structures including spatiotemporal chaotic dynamics; it is of interest to consider ways to control such phenomena. In this work we investigate the possible flow control of gravity-driven perfect dielectric liquid films by imposing an electric field parallel to the substrate on which the fluid lies (and hence parallel to the undisturbed film surface). With this orientation, the electric field is stabilizing and we show that it can be utilized quite dramatically in films on inverted substrates, for example, to arrest dripping or even damp interfacial waves completely.

The non-electrified falling film problem has been studied extensively since the pioneering experiments of Kapitza and Kapitza Kapitza and Kapitza (1949). For an overlying film flow with a sufficiently shallow inclination angle, a flat interface with semi-parabolic velocity profile is an exact and stable solution – the Nusselt solution Nusselt (1916). Linear stability analyses by Yih Yih (1955, 1963) and Benjamin Benjamin (1957) find a critical Reynolds number, dependent on the inclination angle, above which the Nusselt solution first becomes unstable. Two-dimensional (2D) solitary waves form Liu and Gollub (1994), comprising a humped wave front with leading capillary ripples, which transition to 3D waves due to secondary instabilities Kharlamov et al. (2015). Inlet forcing can be used to induce these 2D waves (as opposed to allowing the instability to develop from noise alone), and can also be employed as a flow control to preserve the 2D phase of the dynamics Alekseenko et al. (1994); Park and Nosoko (2003). For hanging film flows on inverted substrates, the situation is more complicated and there are fewer experimental and theoretical studies. The cross-stream component of gravity now destabilizes the film interface, aiding the inertial (Kapitza) instability in the streamwise direction, and inducing the classical Rayleigh–Taylor instability in the spanwise direction. The latter gives rise to rivulet structures as observed in experiments Rothrock (1968); Charogiannis and Markides (2016); Charogiannis et al. (2018), with waves on the crests which may drip depending on parameters Indeikina et al. (1997); Brun et al. (2015). A recent numerical study of both overlying and hanging liquid films in the 2D setting was performed by Rohlfs et al. Rohlfs et al. (2017), in which the onsets of flow reversal and circulating waves are investigated.

Rothrock Rothrock (1968) was the first to perform a careful experimental study of liquid films on inverted substrates, additionally carrying out a linear stability analysis to obtain onset conditions for waves on pendant rivulets – these conditions agreed with experiments. Film dewetting was observed, as the rivulets which formed near the inlet fed into a wedge-shaped fluid film which ended with a single pendant rivulet. Regimes of dripping and no dripping were observed. Drop pinch-off from pendant rivulets was studied by Indeikina et al. Indeikina et al. (1997), with further experiments and the consideration of an inertialess long-wave lubrication equation coupled with a matched asymptotics procedure. The authors identified two distinct mechanisms for dripping depending on the static contact angle of the rivulet and the substrate inclination angle: a jet mechanism at high flow rates due to failure of axial curvature to counteract gravity, and another due to the failure of azimuthal curvature. Charogiannis and Markides Charogiannis and Markides (2016) and Charogiannis et al. Charogiannis et al. (2018) performed experiments with fully wetting films on the underside of a flat substrate inclined at $15^{\circ},30^{\circ}$, and $45^{\circ}$ beyond vertical, for which no dripping was reported. The interface shapes observed were strongly 3D, with clear rivulet structures in the majority of cases with low to moderate Reynolds numbers. The pulses on the crests of the rivulets increased in amplitude as the Reynolds number increased. For larger Reynolds numbers, wavefronts formed across multiple rivulets, causing loss of fluid mass from the rivulets into the separating troughs. By approximating the transverse wavelength from experimental data, Charogiannis et al. Charogiannis et al. (2018) classified two types of rivulet formation depending on the inclination angle and the Kapitza number. They found that the wavelengths of the transverse rivulets for the more extreme inclinations (i.e. $45^{\circ}$ beyond vertical) and/or lower Kapitza numbers, were as predicted by the linear stability of the flat film solution, matching the wavelength of the most unstable transverse mode arising from the competition between the cross-stream component of gravity and surface tension. However, for much smaller inclinations beyond vertical and/or larger Kapitza numbers, they found that the transverse wavelength matched that of the canonical Rayleigh–Taylor instability for a film hanging from a horizontal substrate (full vertical gravity versus surface tension). They suggested that in the former case, the primary instability was Rayleigh–Taylor, whereas in the latter case, the rivulet formation was due to a secondary instability of suspended 2D wavefronts.

A hierarchy of reduced-order models may be constructed using a long-wave methodology to describe the dynamics of the fluid film; these simplify the problem both analytically and numerically while retaining the relevant physical effects. A so-called Benney equation Benney (1966); Gjevik (1970) for the film thickness arises for Reynolds numbers close to critical. This highly nonlinear equation retains the effects of inertia, gravity, viscosity and surface tension. Analytical and numerical studies of the Benney equation were carried out by a number of authors Pumir et al. (1983); Rosenau et al. (1992); Salamon et al. (1994); Oron and Gottlieb (2002); Gottlieb and Oron (2004); Oron and Gottlieb (2004); Scheid et al. (2005), and finite time blow-up was observed in simulations (regions of parameter space where blow-up occurs coincide closely with parameters for which nonlinear traveling waves cease to exist Pumir et al. (1983); Scheid et al. (2005)). Rosenau et al. Rosenau et al. (1992) considered the Benney equation with a full-curvature regularisation, however this was not effective for all parameters for which finite-time blow-up occurred. Coupled systems of equations for the interface height and fluid flux may be derived using an averaging methodology, giving much better agreement with the true dynamics of liquid films for moderate Reynolds numbers. Such models were first constructed by Kapitza Kapitza (1948a, b) and Shkadov Shkadov (1967), however their systems predicted an incorrect critical Reynolds number. This issue was corrected by the weighted integral boundary layer (WIBL) models developed by Ruyer-Quil and Manneville Ruyer-Quil and Manneville (1998, 2000, 2002); WIBL models show good agreement with DNS and experiment Denner et al. (2018), and better matching with full linear theory (based on the Orr–Sommerfeld problem) than the Benney equation Kalliadasis et al. (2012).

The complete mechanisms for the dripping of a hanging film are not yet fully understood, but their relation to the absolute (or spatial) instability of the Nusselt solution has been the subject of recent research. Overlying and vertical film flows exhibit convective instabilities, whereas films hanging from the underside of a horizontal substrate exhibit an absolute Rayleigh–Taylor instability. Thus, for given fluid parameters, a film flow transitions from convective to absolute instability as the inclination is increased to a critical angle beyond vertical. The connection between absolute instability and dripping of hanging films was investigated by Brun et al. Brun et al. (2015). The authors derived an inertialess Benney equation for the 2D flow, and identified regions in parameter space for which the flat interface solution was convectively or absolutely unstable. They performed experiments of gravity-driven films hanging from inverted substrates, and found good agreement between the region of parameter space for which dripping was observed, and that in which the Benney model exhibited absolute linear instability. The experiments were conducted by pouring castor oil onto a flat substrate, and letting it spread until the fluid layer was roughly uniform with a given thickness. The substrate was then rotated to some angle beyond vertical, and the number of drips falling from a fixed region of the substrate was recorded. It is noticeable from their results that only a small amount of dripping was observed just beyond the absolute–convective (A/C) threshold, with much more intense dripping further into the absolute instability regime. Inertial effects and higher order terms were included by Scheid et al. Scheid et al. (2016) in their study of the A/C transition for WIBL models. The authors reported a large discrepancy with the results of the inertialess Benney equation away from zero Reynolds numbers, and found a fluid-independent critical angle below which only convective instabilities exist for their models. Kofman et al. Kofman et al. (2018) employed DNS for the 2D problem with spatially periodic boundary conditions, finding that the dripping onset did not coincide closely with the A/C transition curve obtained from the WIBL models in Scheid et al. (2016). They attribute dripping to a secondary instability of travelling wave solutions. We note that, the dripping thresholds computed by Kofman et al. Kofman et al. (2018) improve in their agreement with the A/C threshold predictions as the length of the periodic domain increases (in the region where the A/C curve in terms of thickness and inclination angle is monotonic). We emphasize that due to nonlinear effects, dripping and absolute instability should not align exactly, however a predictor for dripping phenomena based on linear theory is useful. The consideration of 2D models to study the inherently 3D process of dripping is not entirely justified, in particular a 2D study cannot capture the second mechanism identified by Indeikina et al. Indeikina et al. (1997). However, for domains which are sufficiently small in the spanwise direction, i.e. below the threshold of the spanwise Rayleigh–Taylor instability, a 2D flow assumption is appropriate. Furthermore, we believe that a 2D model is a reasonable approximation for the flow on the crest of a wide rivulet.

Lin et al. Lin et al. (2012) considered the dynamics of a 3D fluid front on the underside of an inclined plane and performed numerical simulations of a multidimensional Benney equation, regularizing the problem with the addition of a thin precursor film. They found that the fluid fronts were unstable to a transverse fingering instability. Thin rivulets form with approximately equal width in the spanwise direction, and fast moving “drop-like” waves appear on the rivulet crests as in the wetted case. Although the finger formation is not attributed to the Rayleigh–Taylor instability, the dynamics on their crests is of relevance to both the wetted and non-wetted cases. The effect of vertical electric fields and temperature gradients on the linear stability of such fluid fronts was considered by Conroy et al. Conroy et al. (2019).

The use of horizontal electric fields to stabilize the Rayleigh–Taylor instability in stratified systems of dielectric fluids was considered by Cimpeanu et al. Cimpeanu et al. (2014) and Anderson et al. Anderson et al. (2017). The former work investigates an unbounded two-fluid system of viscous dielectric fluids with a less dense lower fluid. The latter study considers the related problem of an upper fluid layer and a hydrodynamically passive lower layer bounded above and below by solid dielectric substrates. Linear theory showed that the growth rates decrease as the field strength is increased, and the band of unstable wavenumbers shrinks in extent. If the domain is finite in the horizontal direction, then complete linear stabilization of the flat interface solution is attained with a sufficiently strong field. DNS of the Navier–Stokes equations coupled with electrostatics was also carried out, and the parameters for which finger formation was suppressed was found to be in good agreement with the linear theory. In addition, Anderson et al. Anderson et al. (2017) derived a nonlinear nonlocal evolution equation for the interfacial dynamics valid for sufficiently thin liquid layers. Numerical solutions of the model accurately capture the primary collar and secondary lobe structures present in the early stages of finger formation (as validated with DNS). In both Cimpeanu et al. (2014) and Anderson et al. (2017), the authors demonstrate numerically the possibility of active control of the underlying Rayleigh–Taylor instability and production of sustained nonlinear interfacial oscillations for arbitrarily long times. Such oscillations can enhance mixing, for example, as in Cimpeanu and Papageorgiou (2015). Importantly, we note that such varying electric fields are required for this particular problem to give bounded nontrivial solutions at large times; for constant field strengths, either the flat state is stable or finger formation accompanied by film rupture occurs. Of interest, therefore, is the recent study of Kord and Capecelatro Kord and Capecelatro (2019), on optimal perturbations for controlling the Rayleigh–Taylor instability and the induced mixing.

The studies described above were conducted in the absence of a mean flow, which adds to the complexity of the physical system and further enriches the interplay between the competing mechanisms. In the present work, we consider the application of parallel electric fields to gravity-driven flows on inclined flat substrates. We work with a 2D formulation of the problem, and allow the fluid to be either overlying or hanging. A long-wave approximation allows us to construct a fully nonlinear Benney equation with nonlocal electric field effects. We also provide electrified versions of two WIBL models derived by Ruyer-Quil and Manneville Ruyer-Quil and Manneville (1998, 2000, 2002). Temporal linear stability analyses of all three models show the stabilizing effect of the electric field, and complete stabilization of finite length systems (even in the absence of surface tension) as found in Cimpeanu et al. (2014); Anderson et al. (2017). The major effort focuses on hanging films and the effect of the electric field on the A/C transition. Unsurprisingly, we find that an increased field strength restricts the range of parameters that yield absolutely unstable systems, i.e. promotes convective instability (spatial stability) of the flat film solution. Using a WIBL model, we obtain a minimum critical angle depending on the electric field strength, below which only convective instabilities occur, extending the findings in Scheid et al. (2016). Comparisons with results for the full stability problem at zero Reynolds number are also provided. A similar stability study of a nonlocal problem was performed for the Benney equation with (linearly destabilizing) normal electric field effects by Blyth et al. Blyth et al. (2018). They investigated the stability of solitary wave pulses by analysing absolute/convective instabilities in a reference frame travelling with the pulses. They found that increased field strengths promoted convective instability; their finding appears to be linked to the increased pulse speed as the electric field strength is raised, so that the pulses escape the expanding wave packets generated by localized disturbances. For large field strengths, the pulse solutions are hence much more stable and attracting structures, in agreement with full time-dependent computations. Although the results are not presented here, we found that normal electric fields promoted absolute instability of the flat film solution (in a fixed reference frame).

Following Brun et al. (2015); Kofman et al. (2018), we investigate the relationship between absolute instability and dripping of hanging films, both with and without electric fields. The working fluid for the DNS is the same as used in experiments by Brun et al. Brun et al. (2015), and we took care to faithfully reproduce the experiments in silico, consequently investigating electric field effects through extensive computations that could be hard to sample experimentally. To this end, and in contrast to Kofman et al. (2018) where spatially periodic boundary conditions are used, we initialize with a flat interface and excite the most unstable waves with time-periodic forcing at the inflow (as used in Denner et al. (2018) for the computation of solitary waves on overlying films). Hence, we both mimic experimental conditions and are able to produce solutions with a regular spacing of waves/drops. In the non-electrified case, we find a reasonably close agreement between the onset of dripping and the lower threshold of absolute instability. We believe that this is partly due to the low Reynolds numbers we consider, but also due to the use of long domains in our DNS coupled with appropriate inflow/outflow boundary conditions. In borderline cases, it can take many wavelengths from the inlet for a drip to develop. Our simulations indicate that (non-transient) dripping takes place as an instability of individual solitary pulses, agreeing with the discussions in Kofman et al. (2018). For the electrified problem the agreement is less striking, but the A/C threshold still remains a good order-of-magnitude estimate for complete dripping suppression, something that cannot be obtained from a temporal stability analysis. Increased electric field strengths delay the first dripping event, and induce temporally and spatially irregular dripping that does not appear to be transient (as in the non-electrified case). We believe that the nonlocality of the problem is an important factor here. We emphasize that, unlike the electric field stabilization of the classical Rayleigh–Taylor instability Cimpeanu et al. (2014); Anderson et al. (2017), the mean flow gives rise to nontrivial bounded solutions (wave trains that do not drip) below the threshold of stabilization of the flat interface solution. For this reason, we do not investigate active control strategies.

The structure of the rest of the paper is as follows. In Section II, we present the physical model and full set of governing equations for the 2D flow. In Section III, we give the long-wave models for the electrified gravity-driven film flow, and perform a temporal linear stability analysis. In Section IV we specialize to hanging film flows, and perform a spatial stability analysis of the long-wave models to determine parameter regimes of absolute or convective instabilities. We obtain a minimum critical angle, depending on the electric field strength, below which all instabilities are convective for all flow parameters. We also compare the linear results of the long-wave models at zero Reynolds number to those of the full Stokes flow problem. In Section V, we carry out DNS of the full problem with small-amplitude pulse initial data to verify the regimes of absolute and convective instability predicted by the reduced-order models. Section VI provides DNS of the full problem to determine the dripping onset for both the non-electrified and electrified case, and investigates its relation to absolute instability. A discussion and conclusions are given in Section VII.

II Physical Model and Governing Equations

We consider a Newtonian fluid with constant density $\rho$ and dynamic viscosity $\mu$ (the kinematic viscosity is $\nu=\mu/\rho$), flowing under gravity on a flat substrate inclined at some angle $\theta$ to the horizontal. For the 2D problem considered here, we take coordinates $(x,z)$ with $x$ in the streamwise direction and $z$ perpendicular to the substrate surface – see the schematic in figure 1. We have overlying film flows for $\theta\in(0,\pi/2)$, a vertical film flow for $\theta=\pi/2$, and hanging flows for $\theta\in(\pi/2,\pi)$ as shown in figure 1. The surface tension coefficient between the fluid and the surrounding medium is $\sigma$, and the gravitational force in our chosen coordinate system is $\bm{g}=g\sin\theta\,\bm{e}_{x}-g\cos\theta\,\bm{e}_{z}$, with $\bm{e}_{x}$ and $\bm{e}_{z}$ denoting unit vectors in the $x$ and $z$ directions, respectively. The local film thickness is denoted by $h(x,t)$, a function of space and time, with unperturbed thickness $\ell$. The liquid layer is denoted by Region I, with Region II being the hydrodynamically passive medium defined by $z>h(x,t)$ (here, the pressure is constant and denoted by $p_{\text{atm}}$). The fluid in Region I is governed by the usual Navier–Stokes equations

[TABLE]

where $\bm{u}=(u,w)$ is the velocity field and $p$ is the pressure. At the substrate surface, we have the no-slip and impermeability conditions, $\bm{u}|_{z=0}=\bm{0}$.

We assume that the substrate (Region S), fluid (Region I), and gas phase (Region II) are perfect dielectrics with dielectric permittivities $\epsilon_{0}\epsilon^{\textrm{S}}$, $\epsilon_{0}\epsilon^{\textrm{I}}$, and $\epsilon_{0}\epsilon^{\textrm{II}}$, respectively, where $\epsilon_{0}$ is the permittivity of free space. We denote the voltage potentials in each region by $V^{\textrm{S}}$, $V^{\textrm{I}}$, and $V^{\textrm{II}}$, each with corresponding electric fields $\bm{E}^{\iota}=-\bm{\nabla}V^{\iota}$. Gauss’ law provides harmonic problems for the voltages, i.e.

[TABLE]

These arguments follow from the electrostatic limit of Maxwell’s equations appropriate to this study – see the review by Papageorgiou Papageorgiou (2019) for details. The imposed field is parallel to the solid substrate surface as shown in figure 1, and hence the boundary conditions far from the liquid phase in the normal direction are

[TABLE]

where $\bm{E}_{0}=E_{0}\,\bm{e}_{x}$ and $E_{0}=V_{0}/L_{0}$ measures the strength of the imposed field – here $V_{0}$ is the voltage drop across a system of length $L_{0}$. Making the usual assumption of zero impressed charges at interfaces, we have the following electrical boundary conditions at the substrate

[TABLE]

subscripts denote partial derivatives and the jump notation $\left[\;\cdot\;\right]_{\textrm{I}}^{\textrm{S}}=(\;\cdot\;)_{\textrm{S}}-(\;\cdot\;)_{\textrm{I}}$ has been introduced. The first condition corresponds to continuity of the displacement field and the second to continuity of the voltage potentials – see Papageorgiou (2019).

To calculate the boundary conditions at the gas–liquid interface, we first define the unit tangent and outward-pointing normal vectors there, $\bm{t}=(1,h_{x})^{T}/\sqrt{1+h_{x}^{2}}$ and $\bm{n}=(-h_{x},1)^{T}/\sqrt{1+h_{x}^{2}}$. For the remainder of the section, all $z$-dependent expressions are evaluated at the interface $z=h(x,t)$. The kinematic condition reads

[TABLE]

and the conditions on the voltage potentials, analogous to (4), are

[TABLE]

The final two conditions are the continuity of normal and tangential stresses at the interface

[TABLE]

where the stress tensors include hydrodynamic and/or Maxwell stresses as appropriate, and (employing the usual subscript notation) are given by

[TABLE]

The stress balances (7) written out in full become

[TABLE]

where (6) is used to simplify the tangential stress balance (10). Electrohydrodynamic coupling is present through the normal stress balance alone as expected for interfaces between perfect dielectrics (and also between a perfect dielectric and a perfect conductor Papageorgiou and Petropoulos (2004)). For finite conductivities, consideration of the Taylor–Melcher leaky dielectric model is appropriate Melcher and Taylor (1969); Saville (1997); Pease and Russel (2002); Papageorgiou (2019).

II.1 Exact solution and non-dimensional equations

An exact solution to the full formulation (extending the Nusselt solution to the electrified problem) is

[TABLE]

The velocity profile is semi-parabolic in $z$, and the voltage potential is linear in $x$. We non-dimensionalize velocities with the base velocity at the free surface, $U_{0}=\overline{u}|_{z=\ell}=g\ell^{2}\sin\theta/2\nu$, and make use of the non-dimensional parameters

[TABLE]

Here, ${\operatorname{\mathit{R\kern-1.49994pte}}}$ is the Reynolds number measuring the ratio of inertial to viscous forces, ${\operatorname{\mathit{W\kern-3.00003pte}}}$ is the electric Weber number measuring the ratio of electrical to fluid pressures, and C is the capillary number measuring the ratio of viscous to surface tension forces. To non-dimensionalize we write

[TABLE]

substitute into the governing equations and boundary conditions, and drop the stars. In Region I, the Navier–Stokes equations transform to

[TABLE]

The no-slip and impermeability conditions, and Laplace’s equations for the voltage potentials are unchanged under the change of variables (13), the far-field conditions (3) become

[TABLE]

while the conditions (4) for the voltage potentials at the substrate surface are unchanged. For the interfacial conditions, the kinematic condition (5), continuity of voltage (6b), and the tangential stress balance (10) are also unchanged. Equation (6a) transforms to

[TABLE]

Finally, the normal stress balance becomes

[TABLE]

The exact Nusselt solution to the Navier–Stokes equations (14) and the above boundary conditions in non-dimensional form is

[TABLE]

In Appendix A, we provide the Orr–Sommerfeld system for the linearization of the dynamics about (18); the electrostatics component is analytically tractable and the remaining problem is a modification of the usual Orr–Sommerfeld system for thin films Kalliadasis et al. (2012). We also consider the Stokes flow limit in Appendix B, for which we obtain an exact dispersion relation and compare this with the linear behavior of the long-wave models given next.

III Hierarchy of nonlinear long-wave models

We utilize a hierarchy of long-wave models to analyse the nature of the linear instabilities present in the electrified flow – the spatial stability of these models is investigated in Section IV, and detailed comparisons with DNS are undertaken in Section V. The models considered are the Benney and two WIBL models. We omit details of their derivations since they appear in Tseluiko and Papageorgiou (2006); Tomlin et al. (2017) for the Benney equation and Ruyer-Quil and Manneville (1998, 2000, 2002) for the WIBL models (albeit using a different non-dimensionalisation), but give the calculations used to obtain the electric field contribution that enters through the normal stress balance (17).

To derive the models, we assume that (in the original dimensional variables) the typical wavelength of interfacial deformations $\lambda$ is long in comparison to the undisturbed liquid height $\ell$, so that $\delta=\ell/\lambda\ll 1$. We introduce the rescalings

[TABLE]

and apply a systematic asymptotic procedure. For the derivations, it is also assumed that ${\operatorname{\mathit{R\kern-1.49994pte}}}=O(1)$, although this has been found to be an unnecessary restriction for the WIBL models. The Benney equation arises from an asymptotically correct elimination of the flow field variables in the kinematic equation (5). For an approximation of the interface thickness $H$ (the first two terms of an asymptotic expansion of $h$), the Benney model with errors of $O(\delta^{2})$ is

[TABLE]

Here, $P$ is the leading-order pressure at the interface, taken to be $O(1)$ so that it enters the dynamics. Equation (20) is effective at modelling flows for Reynolds numbers just beyond critical; increasing $\operatorname{\mathit{R\kern-1.49994pte}}$ further, singular phenomena such as finite time blow-up are observed in numerical simulations Pumir et al. (1983); Rosenau et al. (1992).

Even though WIBL models rely on closure assumptions rather than rational asymptotic approximations, they are far more accurate than Benney equations at describing thin films for Reynolds numbers away from critical; they correctly capture the dynamics beyond the drag–gravity regime to which Benney equations are restricted Kalliadasis et al. (2012). In order to obtain the WIBL models, we rewrite the kinematic condition (5) as

[TABLE]

where $f$ is the fluid flux through a slice of the film in the $z$-plane. Following a weighted residuals strategy, in which the flow field is expanded in polynomials of $z$, the WIBL2 (simplified second-order) model, which comprises an approximation of the kinematic condition (21) coupled to an equation for the time evolution of $F$ (an approximation of $f$), is given by

[TABLE]

The WIBL1 (first-order) model is obtained by omitting the $O(\delta^{2})$ term in (23). The dependence on the flow field variables is not completely eliminated as in the Benney equation (20); the latter may be obtained from WIBL1 by using the leading order relation $F=2H^{3}/3$ in the $O(\delta)$ terms of (23), and substituting the result into (22). Furthermore, Benney and WIBL1 are identical at zero Reynolds numbers, and all models reduce to the same nonlinear hyperbolic equation for $\delta=0$. Denner et al. Denner et al. (2018) reported that simulations of WIBL2 correctly captured the main humps in solitary wave trains even up to ${\operatorname{\mathit{R\kern-1.49994pte}}}=100$, agreeing with DNS and experiments, but overestimated the amplitude of the leading capillary ripples for ${\operatorname{\mathit{R\kern-1.49994pte}}}\gtrsim 10$ (note that our ${\operatorname{\mathit{R\kern-1.49994pte}}}$ is $1.5$ times greater than the Reynolds number used in Denner et al. (2018)). We consider Reynolds numbers up to these values in the linear theory to follow – it is known that there is good agreement with the Orr–Sommerfeld linear theory for these values Kalliadasis et al. (2012). For all models it remains to compute $P$ in terms of $H$, and this is undertaken next to clarify the additional effect due to the electric field.

The normal stress balance (17) under the lubrication scalings (19) becomes

[TABLE]

where we have retained leading and next-order terms inside the bracket corresponding to the electric field effect. We expand the solutions as

[TABLE]

We also introduce a stretched normal variable $z=\zeta/\delta$, $\zeta=O(1)$, in the non-slender Regions S and II, so that the Laplace equations in the three regions become

[TABLE]

The solutions for $V^{\textrm{S}}$, $V^{\textrm{II}}$ can be found by taking a Fourier transform in $x$ and applying the far-field conditions (15); these are

[TABLE]

where the complex wavenumber $\xi=\xi_{\textrm{r}}+i\xi_{\textrm{i}}$ has been introduced in anticipation of the spatial linear stability analysis to follow. The multipliers $\widehat{A}^{\textrm{S}}(\xi)$ and $\widehat{B}^{\textrm{II}}(\xi)$ have asymptotic expansions similar to (25). The solution in Region I has the asymptotic expansion

[TABLE]

The electric field boundary conditions at the substrate (4) yield the relations

[TABLE]

From the interfacial boundary conditions for the voltage potentials (6b,16), we also obtain

[TABLE]

These give the Fourier transform of the voltage potential at the interface to leading order as (c.f. equation (40) in Anderson et al. (2017) with $\epsilon_{B}=1$, corresponding to the case of no lower bounding solid)

[TABLE]

for $\iota=$ I, II, where $\mathcal{H}$ denotes the Hilbert transform with Fourier symbol $\widehat{\mathcal{H}}(\xi)=-i\operatorname{sign}(\xi)$ for $\xi\in\mathbb{R}$. The expression for $\widehat{V}^{\iota}$ for complex wavenumbers given in (31a) will be important for the spatial stability analysis in the next section. It follows from this result and (24) that to retain electrical and capillary effects in the leading order evolution, we require ${\operatorname{\mathit{W\kern-3.00003pte}}}=O(\delta^{-1})$ and ${\mbox{{C}}}=O(\delta^{2})$. The normal stress balance (24) gives

[TABLE]

where

[TABLE]

It follows that the field has no effect on the interfacial dynamics if the electrical permittivities in Regions I and II are equal as ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ is always zero. The impact of the electric field is maximized for $\epsilon^{\textrm{I}}$ large, with $\epsilon^{\textrm{S}}$ and $\epsilon^{\textrm{II}}$ both small.

With the rescaling $\delta(\partial_{t},\partial_{x})\rightarrow(\partial_{t},\partial_{x})$, we may return to the original time and space variables (as in the non-dimensional Navier–Stokes setting), formally setting $\delta=1$ in the long-wave models (assumed from this point). As we will see in the following linear stability analysis, both surface tension and the electric field have a stabilizing effect on the interface dynamics; for $\mbox{{C}}^{-1/2}\sim\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$, both effects are relevant. For the results of the models to be valid, parameters should be chosen so that at least one of these effects is retained; if both surface tension and the electric field are negligible in the leading-order dynamics, the instability moves towards the short waves, invalidating the initial long-wave assumption. Although the small parameter is now scaled out of the equations, it remains in the problem implicitly, indicating the subset of parameter space which is compatible with the long-wave assumption. For non-zero field strengths, the models are well-posed in the limit of weak surface tension, i.e. $\mbox{{C}}^{-1/2}=o(\operatorname{\mathit{W\kern-3.00003pte}}^{\prime})$, in the sense that high wavenumbers remain damped. The corresponding weakly nonlinear evolution for electrified flows with weak surface tension is described by the nonlocal Kuramoto–Sivashinsky type equation

[TABLE]

which is also well-posed. A fourth-order term is included if surface tension effects are not negligible. We do not consider (34) in the current work.

III.1 Temporal linear stability analysis

The dispersion relation for the Benney model is obtained by substituting $H=1+\tilde{H}e^{i\xi x+\omega t}$, with $\xi$ real and complex $\omega=\omega_{\textrm{r}}+i\omega_{\textrm{i}}$, into (20) and linearizing for small $\tilde{H}$. This gives

[TABLE]

The phase velocity is $-\omega_{\textrm{i}}/\xi=2$, inertial forces are linearly destabilizing, with both electric and surface tension forces being linearly stabilizing. Gravitational forces are stabilizing for overlying films ($\theta<\pi/2$), and destabilizing for hanging films. The critical Reynolds number at which the flat film state destabilizes is recovered, ${\operatorname{\mathit{R\kern-1.49994pte}}}_{\textrm{c}}\equiv 5\cot\theta/4$, a quantity which is zero for vertical films. It follows from (35), and also holds for the WIBL models and the full system, that if a flow is above critical, ${\operatorname{\mathit{R\kern-1.49994pte}}}>{\operatorname{\mathit{R\kern-1.49994pte}}}_{\textrm{c}}$ (which is always verified for hanging flows since $\operatorname{\mathit{R\kern-1.49994pte}}_{\textrm{c}}<0\leq\operatorname{\mathit{R\kern-1.49994pte}}$), then increasing surface tension and/or the strength of the parallel electric field cannot prevent linear instability in all wavenumbers; we have $\omega_{\textrm{r}}\sim\xi^{2}$ locally near $\xi=0$, and thus, for any $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$ and C, the flat state will be linearly unstable on a sufficiently large spatial domain. For the temporal linear stability of the WIBL models, we substitute $H=1+\tilde{H}e^{i\xi x+\omega t}$ and $F=2/3+\tilde{F}e^{i\xi x+\omega t}$ into (22,23), and linearize to obtain a quadratic in $\omega$:

[TABLE]

The two double-underlined terms in the above expression correspond to $O(\delta^{2})$ terms in the WIBL2 model, and are dropped for the consideration of WIBL1. Both WIBL models have the same critical Reynolds number ${\operatorname{\mathit{R\kern-1.49994pte}}}_{\textrm{c}}$ as the Benney model and Orr–Sommerfeld theory.

Figure 2 plots growth rates for the three models in the vertical substrate case, $\theta=\pi/2$, with ${\mbox{{C}}}=0.01$ and $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}=0,10,20$, for ${\operatorname{\mathit{R\kern-1.49994pte}}}=1$ in panel (a) and ${\operatorname{\mathit{R\kern-1.49994pte}}}=4$ in panel (b). The growth rates for the Benney model, WIBL1 and WIBL2 are plotted with dotted, dashed and solid curves, respectively. The results from the three models are in closer agreement for $\operatorname{\mathit{R\kern-1.49994pte}}=1$, which is nearer to the critical Reynolds number ${\operatorname{\mathit{R\kern-1.49994pte}}}_{\textrm{c}}=0$, with more visible discrepancies for $\operatorname{\mathit{R\kern-1.49994pte}}=4$ in panel (b). The stabilizing effect of the electric field is evident – in panel (a), the band of modes with wavenumbers between approximately $0.06$ and $0.12$ become linearly stable as $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$ is increased from [math] to $10$, and the maximum growth rate decreases.

The band of unstable modes with positive wavenumbers $(0,\xi_{\textrm{c}})$ is the same for the Benney and WIBL1 models, where $\xi_{\textrm{c}}$ is the critical wavenumber satisfying $\omega_{\textrm{r}}=0$, with linear stability obtained if the non-dimensional system length $L=L_{0}/\ell$ satisfies $L<2\pi/\xi_{\textrm{c}}$.

In figure 3 we plot the variation of the critical wavenumber $\xi_{\textrm{c}}$ against the Reynolds number for a vertical film flow. For a given $\operatorname{\mathit{R\kern-1.49994pte}}$, instability is found for wavenumbers below the curves. The long-wave models are compared to the full stability results based on the Orr–Sommerfeld equation discussed in Appendix A, taking ${\mbox{{C}}}=0.01$ and ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0,\,40$. We note that the Benney and WIBL1 models give the same $\xi_{\textrm{c}}$, and thus are represented in the figure with the same (dotted) line. Unlike the long-wave models, dependence of the Orr–Sommerfeld system on the dielectric constants is not completely absorbed into the non-dimensional parameter ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ – see expression (61) in Appendix A. Thus, we must also prescribe electrical permittivities; we took $\epsilon^{\textrm{S}}=1.5$, $\epsilon^{\textrm{I}}=2$, and $\epsilon^{\textrm{II}}=1$ to obtain the solid curve in figure 3. It is clear that WIBL2 agrees with the Orr–Sommerfeld theory much more closely than Benney/WIBL1 in both the non-electrified (as reported in Kalliadasis et al. (2012); see their Fig. 4.2, Ch. 4, p. 72) and electrified case.

The critical electric field strength beyond which the flat film solution is linearly stabilized may also be computed from the above dispersion relations; of particular interest is its dependence on the system length. For example, taking the Benney model on a spatial domain of non-dimensional length $L=L_{0}/\ell$, we require

[TABLE]

for linear stability. This implies that for $L_{0}\gg 1$, the critical $E_{0}$ required to stabilize the flat film solution scales with $L_{0}^{1/2}$. Nonlinearly, we expect that linear stability will imply that the flat film solution is reached after large times, as found by Anderson et al. Anderson et al. (2017) in the related problem without a mean flow. Thus, a sufficiently strong electric field can fully stabilize an overlying or hanging film on a finite spatial domain (at least in 2D). We do not provide nonlinear simulations of the long-wave models here as this is not the focus of the present study. We expect that the field strength required for dripping suppression should be independent of the system length (for $L_{0}\gg 1$), and thus the temporal stability results cannot be used as a predictor of this nonlinear phenomena.

IV Absolute and convective instabilities

In this section we investigate the absolute or convective nature of the linear instabilities for the various long-wave models. We compare these results with (analytical) full linear computations in the case of zero Reynolds number. We also perform DNS with pulse initial conditions in section V to validate the linear results, and to investigate the nonlinear development of the instabilities.

Briefly, a system is classified as convectively unstable if its response to a small localized perturbation is convected away from the location of the initial disturbance, whereas absolute instabilities lead to growth both upstream and downstream of the disturbance location – see Huerre and Monkewitz (1990) for details. Such instabilities are of interest in thin film flows and it is reasonable to expect that an absolutely unstable hanging film is more susceptible to dripping than a convectively unstable one. However, the evolution of liquid films, and in particular the dripping process, is very nonlinear and our approach in addressing such questions is a combination of linear theory and allied DNS. Moreover, linear theory is not in a position to predict whether a convectively unstable film will drip far downstream from the inlet (assuming inlet forcing is driving the system for example), or converge to a bounded wave-train of nonlinearly saturated pulses. Brun et al. Brun et al. (2015) considered this scenario in the case of a 2D flow with small Reynolds numbers and used an inertialess Benney equation given in terms of our non-dimensional parameters as (20) with ${\operatorname{\mathit{R\kern-1.49994pte}}}={\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$. They determined regions of absolute and convective instability analytically in terms of the inclination angle and the ratio of mean film thickness to capillary length (this can be rewritten as a relationship between $\theta$ and C). They found reasonable agreement with experiments, where films in the convectively unstable regime (of their model) were found not to drip at all, or yielded relatively few drips compared to the absolutely unstable regime. An increased number of drips emerged deeper into the absolute instability regime. The linear theory was extended to more complex WIBL models by Scheid et al. Scheid et al. (2016) without the low Reynolds number restriction. They calculated a fluid-dependent critical angle in terms of the Kapitza number (fixed for a given fluid) above which only convective instabilities may exist for any Reynolds number; the minimum of these is a fluid-independent critical angle (corresponding to the limit of zero Kapitza number). In what follows, we show that an electric field may be used to increase these critical angles towards the horizontal arrangement, i.e. increase the parameter space of convectively unstable systems.

We apply the theory of Fokas and Papageorgiou Fokas and Papageorgiou (2005) for the linearized model equations at hand. It is useful to extend the dispersion relation $\omega(\xi)$ to complex wavenumbers $\xi=\xi_{\textrm{r}}+i\xi_{\textrm{i}}\in\mathbb{C}$. For the long-wave models considered here, $\omega(\xi)$ is given by either (35) or (36) with $|\xi|^{3}$ replaced by its correct generalisation to complex wavenumbers, $\operatorname{sign}(\xi_{\textrm{r}})\xi^{3}$. Note also the conjugate symmetry $\omega(-\xi_{\textrm{r}})=\overline{\omega(\xi_{\textrm{r}})}$, where the bar denotes the complex conjugate of a quantity, since the solution of the PDE is real-valued. Given an initial condition $H(x,0)=H_{0}(x)$, the solution is given by the Fourier transform pair

[TABLE]

where $\widehat{H_{0}}(\xi)$ is also conjugate symmetric since $H_{0}(x)$ is real. The first integral in (38) is split into two parts, one integral over the stable frequencies, $\mathbb{R}\backslash(-\xi_{\textrm{c}},\xi_{\textrm{c}})$, where $\xi_{\textrm{c}}\in\mathbb{R}$ is the critical wavenumber, and another integral over the interval $(-\xi_{\textrm{c}},\xi_{\textrm{c}})$ which is denoted $J(x,t)$. The former integral decays to zero as $t\rightarrow\infty$ while keeping $x=O(t)$, thus $J(x,t)$ alone contains information regarding the absolute or convective nature of the instability. Using the conjugate symmetry of $\omega$ and $\widehat{H_{0}}(\xi)$, we have

[TABLE]

As shown in Fokas and Papageorgiou (2005), the type of instability is determined by the topology of the curves $\omega_{\textrm{r}}(\xi)$ in the right half-plane. Cauchy’s theorem is applicable since the dispersion relation is analytic and bounded in this region, where $\operatorname{sign}(\xi_{\textrm{r}})\xi^{3}=\xi^{3}$. If the integral over $(0,\xi_{\textrm{c}})$ can be deformed to part of the curve $\omega_{\textrm{r}}(\xi)=0$ in the right half-plane, then $J(x,t)$ decays for $t$ fixed with $|x|\rightarrow\infty$, and for $x$ fixed with $t\rightarrow\infty$, thus the instability is convective. In the case that the curves of $\omega_{\textrm{r}}(\xi)=0$ do not connect the points [math] and $\xi_{\textrm{c}}$ on $\mathbb{R}_{+}$, then the instability is absolute. The Briggs–Bers pinching criterion can then be applied to our problem – see Brun et al. (2015); Scheid et al. (2016) for its application to the non-electrified flow. Solving $\omega_{\textrm{r}}=0$ and $\partial_{\xi}\omega=0$ for complex wavenumbers in the right half-plane is an algebraic problem, the solutions of which separate the regions between absolute and convective instability. In order to find the transition curve in parameter space, we proceed with numerical continuation of solutions to the algebraic problem using the continuation software AUTO-07P Doedel et al. (2007).

The Benney and WIBL1 models are identical for ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$, and without electric fields they yield the equation studied by Brun et al. Brun et al. (2015). We have a unique non-trivial exact solution to the Briggs–Bers criterion for this system with $\xi_{\textrm{r}}>0$. The critical angle $\theta\in(\pi/2,\pi)$ for this exact solution is defined by $\cot^{3}\theta=-243/4{\mbox{{C}}}(17+7^{3/2}),$ with wavenumber $\xi$ given by

[TABLE]

and $\omega_{\textrm{i}}=-3\xi_{\textrm{r}}/2-2\cot(\theta)\xi_{\textrm{r}}\xi_{\textrm{i}}/3$. In the case of ${\mbox{{C}}}=0.01$, this exact solution is $\theta\approx 2.963$ (about $80^{\circ}$ beyond vertical), $\xi_{\textrm{i}}\approx-0.0873$, $\xi_{\textrm{r}}\approx 0.280$ and $\omega_{\textrm{i}}\approx-0.510$.

Given a capillary number C, the curve in parameter space separating convective from absolute instability for either the Benney, WIBL1 or WIBL2 models can be continued from this exact solution. These transition curves can be validated from the topology of $\omega_{\textrm{r}}=0$ for the given parameters. In figure 4(a,b), we give examples of the $\omega_{\textrm{r}}=0$ curves in the right half-plane (the diagram is symmetric about the imaginary axis) for the Benney and WIBL1 models, respectively, and ${\mbox{{C}}}=0.01$, ${\operatorname{\mathit{R\kern-1.49994pte}}}=4$, $\cot\theta=-8$, with ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0,20,40$. For the Benney model in panel (a), the non-trivial real root is connected to the origin when ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=40$ implying convective instability, but for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0,20$ we predict absolute instability. The “pinching” of the roots corresponding to the Briggs–Bers criterion occurs between ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=20$ and ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=40$ for the Benney model. The results for WIBL1 are different as shown in panel (b), with absolute instability only for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$ out of the three cases depicted (WIBL2 with these parameters gives very similar results). In figure 4(c) we collect the A/C transition curves for the different models as a function of ${\operatorname{\mathit{R\kern-1.49994pte}}}$ and $\cot\theta$, for fixed ${\mbox{{C}}}=0.01$ and ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0,20,40$. The dotted, dashed, and solid curves correspond to the Benney, WIBL1, and WIBL2 models, respectively. Absolute instability is found to the left of the curves with convective instability to their right. The results show that an increase in the applied electric field strength promotes convective instability, with the absolute instability threshold found for larger negative values of $\cot\theta$, i.e. for substrate inclinations tending to the horizontal hanging film configuration – note that $\cot\theta\to-\infty$ as $\theta\to\pi_{-}$ when the inclination becomes horizontal and the film is hanging vertically. The Benney model has a linear transition curve due to its particular dispersion relation, and incorrectly predicts that a vertical film flow ($\cot\theta=0$) undergoes the transition from convective to absolute instability at ${\operatorname{\mathit{R\kern-1.49994pte}}}=6.939$, and that overlying film flows may become absolutely unstable for higher Reynolds numbers.

In contrast, the more accurate WIBL1 and WIBL2 models yield fold points, marked with triangles in the figure, and implying the existence of a fluid-dependent critical angle above which all instabilities are convective irrespective of Reynolds number. The predicted critical angle corresponds to a hanging film, so that WIBL models for overlying films can only be convectively unstable. The star in panel (c) corresponds to $\operatorname{\mathit{R\kern-1.49994pte}}=4$, $\cot\theta=-8$ used in panels (a,b), and verifies the A/C conclusions of the root plots. Figure 4(d) plots the zero Reynolds number transition lines and fold points as the electric Weber number increases; it is evident that there is good agreement between the long-wave models at ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$, and that WIBL1 and WIBL2 agree closely on the location of the critical angle fold point for these parameters; increasing $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$ from [math] to $40$ moves the ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$ transition and fold points by over $5^{\circ}$. Our parameters are related to those in Scheid et al. Scheid et al. (2016) (denoted with a superscript Sch) by $\delta^{\textrm{Sch}}=2^{4/3}{\mbox{{C}}}^{1/3}{\operatorname{\mathit{R\kern-1.49994pte}}}$, $\zeta^{\textrm{Sch}}=2^{1/3}{\mbox{{C}}}^{1/3}\cot\theta$, $\eta^{\textrm{Sch}}=2^{2/3}{\mbox{{C}}}^{2/3}$. These are the reduced Reynolds number, the reduced inclination number, and the viscous extensional number, respectively. For the WIBL1 (and Benney) model, viscous extensional stresses are ignored, i.e. $\eta^{\textrm{Sch}}=0$ or dropping the $O(\delta^{2}$) terms in (23), but for the computations for WIBL2 shown in figure 4(c,d) we have $\eta^{\textrm{Sch}}=0.0737$ – see figure 3 in Scheid et al. (2016). From the results shown in figure 4, we recover their critical values of $\zeta_{A}^{\textrm{Sch}}\approx-1.507$ corresponding to the transition point for zero Reynolds numbers, and $\zeta_{C}^{\textrm{Sch}}\approx-1.221$, corresponding to the fold point for the WIBL1 transition curve; the latter of these is the maximum value of $\zeta^{\textrm{Sch}}$ even for non-zero extensional stresses.

The results in figure 4 are computed for a fixed capillary number C, and varying it gives a 2D surface in $(\theta,\operatorname{\mathit{R\kern-1.49994pte}},\mbox{{C}})$-space, separating the regions of absolute and convective instability; this surface changes with the electric field strength through variations in $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$. In the non-electrified case, Scheid et al. Scheid et al. (2016) showed the existence of a fluid-independent critical angle $\theta_{\textrm{c}}=147.4^{\circ}$ ($\theta_{\textrm{c}}-90^{\circ}=57.4^{\circ}$ beyond vertical) for the full WIBL2 model; below $\theta_{c}$ absolute instabilities cannot exist (see figures 4, 5 and 6 in Scheid et al. (2016)), providing a lower bound on the A/C transition surface in terms of the inclination angle. For ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$, this lower bound is attained at $\operatorname{\mathit{R\kern-1.49994pte}}=0$. We recover the minimum critical angle with the simplified WIBL2 system, but not WIBL1 – inclusion of viscous extensional stresses is vital as absolute instability can be found for WIBL1 at any angle beyond vertical.

Figure 5(a,b) plots the A/C transition surfaces of the WIBL2 model for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$ and ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=5$, respectively, with absolute instability obtained above the surfaces (we give the critical inclination angle from the vertical, $\theta_{\textrm{c}}-90^{\circ}$, expressed in degrees for comparison with Scheid et al. (2016)). The superimposed dotted curves depict the transition curve for fixed ${\operatorname{\mathit{R\kern-1.49994pte}}}$ as C varies, and the solid lines track the minima of $\theta_{\textrm{c}}-90^{\circ}$ over all C as ${\operatorname{\mathit{R\kern-1.49994pte}}}$ is varied. As can be seen from the results in panel 5(a), there is a global minimum value of $\theta_{\textrm{c}}-90^{\circ}=57.4^{\circ}$ at ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$ as found by Scheid et al. Scheid et al. (2016) in their study of the non-electrified flow. Interestingly, the minima are found at $O(1)$ values of C; this non-monotonic A/C threshold is not an artefact of long-wave modelling, as confirmed by DNS (see figure 7 in Scheid et al. (2016)). Unsurprisingly, the strong surface tension limit $\mbox{{C}}\to 0$ yields critical angles which are very close to horizontal for any ${\operatorname{\mathit{R\kern-1.49994pte}}}$. We turn next to the inclusion of a field with ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=5$ shown in figure 5(b). The solid curve tracing the minima diverges to ${\mbox{{C}}}=\infty$ as the Reynolds number is decreased from $100$, corresponding to zero surface tension effects. The dotted curves for small ${\operatorname{\mathit{R\kern-1.49994pte}}}$ shown in panel (b) are monotonically decreasing in C, and their minimum value is obtained in the limit ${\mbox{{C}}}\rightarrow\infty$, which is qualitatively different to the non-electrified case. This is verified by performing continuation of the algebraic system corresponding to the Briggs–Bers criterion in the parameter $1/{\mbox{{C}}}$, where it is found that the fold point passes to prohibited negative values of $1/{\mbox{{C}}}$ as the Reynolds number is decreased.

Panel (c) gives the projection of the minimum curves in panels (a) and (b) with other choices of ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ against ${\operatorname{\mathit{R\kern-1.49994pte}}}$. The region above a given line corresponds to either absolute or convective instability, whereas all parameters below the line give convectively unstable systems. We find three regimes depending on the value of ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$: For ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}<1.85$, the critical angle for each ${\operatorname{\mathit{R\kern-1.49994pte}}}$ is attained at finite C, and the minimum critical angle is found in the zero Reynolds number limit. For ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}\in(1.85,3.62)$, the critical angle for small Reynolds numbers is found in the limit of zero surface tension (indicated with a dotted line section), yet the minimum is still found at ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$. The former is also true for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}>3.62$, except the minimum critical angle is found for ${\operatorname{\mathit{R\kern-1.49994pte}}}>0$, shown with diamonds in panel (c). The minimum critical angle (the minimum of the A/C transition surface over both C and ${\operatorname{\mathit{R\kern-1.49994pte}}}$) is plotted as a function of ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}\sin\theta$ (a $\theta$-independent quantity) in figure 5(d), where the three sections of the solid curve (separated with stars) correspond to the three cases in panel (c). The dotted lines are continuations of the middle section of the solid curve corresponding to ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$ and ${\mbox{{C}}}\rightarrow\infty$.

The A/C instability regions of WIBL2 in the weak surface tension limit are shown in figure 6. Panel (a) plots the transition surface in terms of ${\operatorname{\mathit{R\kern-1.49994pte}}}$ and ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$. The dotted lines are the transition curves for fixed ${\operatorname{\mathit{R\kern-1.49994pte}}}$ as ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ varies, with the solid curve tracking their minima. We note that the minima are attained at ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$ for $O(1)$ values of ${\operatorname{\mathit{R\kern-1.49994pte}}}$, but at non-zero ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ for larger Reynolds numbers. It is interesting to note that figures 5(a) and 6(a) where the electric field or surface tension is absent, respectively, are qualitatively different since either physical mechanism is stabilizing albeit with a different spectral content. It may be of theoretical interest to consider generalizations of the stabilizing terms in the interfacial pressure (32), but this is beyond the scope of the current study. The projection of the solid curve in panel (a) onto the ${\operatorname{\mathit{R\kern-1.49994pte}}}$-direction is given in panel (b), where the dotted line indicates that the minimum is attained at ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0$. As before, parameters above the curve yield absolute or convective instability, with only convective instability below the curve. The minimum critical angle of $59.2^{\circ}$ below vertical is obtained for ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$ – as shown above, inclusion of surface tension can increase or decrease this value. The minimum of the surface in figure 6(a) over ${\operatorname{\mathit{R\kern-1.49994pte}}}$ (this is not the projection of the solid curve onto the ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$-direction) can be obtained from figure 5(d); it comprises of the initial dotted curve for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}\sin\theta<0.675$, and the two other solid curve segments – these three sections correspond to minima attained in the zero surface tension limit.

IV.1 Comparison of long-wave models with Stokes flow results

We now provide comparisons between the A/C regions for WIBL2 at ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$ with those of the corresponding Stokes flow problem. The regions of validity in parameter space of the long-wave models are restricted due to their asymptotic derivation, whereas for the full Stokes flow problem, there are no restrictions on C, ${\operatorname{\mathit{W\kern-3.00003pte}}}$, or the electrical permittivities. From the full Orr–Sommerfeld system for the electrified problem given in Appendix A, the exact dispersion relation at $\operatorname{\mathit{R\kern-1.49994pte}}=0$, computed in Appendix B, reads

[TABLE]

where the electric field term is given by (61). As in the case of the long-wave models, we are able to take advantage of symmetries of the dispersion relation in order to apply the Briggs–Bers criterion in the right-half $\xi$-plane. We again utilize AUTO-07P to solve this algebraic problem and perform solution continuation. The expressions for the real and imaginary parts of $\omega$ and $\partial_{\xi}\omega$ are exceedingly long – we used Maple to expand these expressions and convert them into a format which was appropriate for AUTO-07P.

Figure 7 gives a comparison of the Stokes results with those of WIBL2 restricted to zero Reynolds number. Panel (a) plots the critical inclination angle for the Stokes problem (solid curve) and WIBL2 model (dotted curve) against C, for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=0,5$. The curves agree closely in the case of strong surface tension, $\mbox{{C}}\rightarrow 0$, with deviation for larger capillary numbers. This is expected as the system becomes unstable to shorter waves for these parameters, thus invalidating the long-wave model. In the non-electrified case, the minimum is again found at an $O(1)$ value of C, although at a different value of $\theta_{\textrm{c}}$, further confirming that the non-monotonic A/C behavior is not a fabrication of the long-wave methodology. A minimum critical angle of $58.0^{\circ}$ below vertical is obtained, differing from the result of $57.4^{\circ}$ for WIBL2 (and the full second-order WIBL model). It may be the case that a higher-order long-wave model captures this value more accurately (WIBL3), but this is beyond the scope of the current study. This reveals a slight downfall in the use of long-wave models to predict the A/C behavior of the full problem. The results for ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}=5$ also agree qualitatively, although with a discrepancy for negligible surface tension – ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}$ is not large enough to maintain the long-wave assumption in this limit. Panel (b) shows a continuation of the minimum critical angle against ${\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}\sin\theta$. The curves intersect and do not converge together as the field strength increases; this is expected as the instability moves towards shorter waves, and the capillary number where the minima are attained increases (to the right of the stars, ${\mbox{{C}}}\rightarrow\infty$). To compute the results for the Stokes flow problem, we used permittivities $\epsilon^{\textrm{S}}=1.5$, $\epsilon^{\textrm{I}}=2$, and $\epsilon^{\textrm{II}}=1$ (as in figure 3).

V Direct Numerical Simulations: Pulse initialization

In this section, we perform DNS of the full Navier–Stokes problem coupled with electrostatics, adding a small-amplitude pulse disturbance to the Nusselt solution as the initial condition. This is utilized to check the accuracy of the predictions of A/C instability for the full problem with that of the low-dimensional models. A similar numerical experiment was performed for the non-electrified case by Scheid et al. Scheid et al. (2016), who showed that the full second-order WIBL model captured the correct A/C behavior of the Navier–Stokes problem. This section provides a similar validation with the addition of the stabilizing electric field.

The hydrodynamical component of the results presented in the previous section can be recast in terms of the mean interface height $\ell$, inclination angle $\theta$, and a geometry-independent Kapitza number $\sigma/\rho g^{1/3}\nu^{4/3}$, which is fixed for a given fluid. Fixing the dielectric constants also means that $E_{0}$ alone parameterizes the electrical component of the problem. Thus, we may obtain a surface in $(\ell,\theta,E_{0})$-space separating the regions of absolute and convective instability for each model. This surface is not monotonic (for a model including viscous extensional stresses) – this is visible from the results for the non-electrified problem by taking a constant-Kapitza number slice of figure 5(a), as well as figure 7 in Scheid et al. (2016) and figure 3 in Kofman et al. (2018). For brevity, we do not re-plot the results in terms of these parameters, but provide A/C thresholds in tables when comparing to DNS.

We fix the working fluid to be the Castor Oil (Hänseler AG) used in experiments by Brun et al. Brun et al. (2015). They measured the capillary and viscous lengths to be $1.91$ mm and $4.4$ mm, respectively, and also give the dynamic viscosity to be $865\pm 5$ cP. The fluid properties corresponding to these measurements are given in SI units in Table 1, and we also fix $g=9.81$ m/s2 for all calculations.

We assign to the fluid the dielectric constant $\epsilon^{\textrm{I}}=5$, which is roughly the value found for castor oils. We note that, although castor oils typically have a very low electrical conductivity (less than $10^{-10}$ Siemens/m), they have been observed in certain experiments to behave as leaky rather than perfect dielectrics Burcham and Saville (2000). Setting $\epsilon^{\textrm{II}}=1$, the value for air, we have $\epsilon^{\textrm{I}}>\epsilon^{\textrm{II}}$, and so the parallel field still has a linearly stabilizing effect on the interface even if the conductivity of the castor oil is accounted for (Papageorgiou and Petropoulos, 2004; Uguz et al., 2008). We ignore conductivities for the purposes of the present study; if conductivities are included, our models become coupled to an equation for the interfacial charge, as in Craster and Matar (2005) – the study of such models is warranted for comparison with experiments. By choosing the electrical permittivity of the solid substrate to match that of the liquid phase, $\epsilon^{\textrm{S}}=5$ – a realistic value for Pyrex glass, we do not need to solve for $V^{\textrm{S}}$ (the numerical implementation is discussed later). For our calculations, we take the permittivity of free space to be $\epsilon_{0}=8.85\times 10^{-12}$ A${}^{2}\cdot$s4/kg$\cdot$m**3. We note that the dielectric breakdown of air occurs for electric field strengths beyond approximately $E_{0}=3\times 10^{6}$ V/m Tipler (1987), and thus we do not consider values of $E_{0}$ beyond this.

For our numerical simulation of the full system with pulse initial conditions, we choose $\theta=5\pi/6$ and $\ell=2.5$ mm, in which case all models (long-wave and Stokes) predict the non-electrified case to be absolutely unstable. However, this is sufficiently close to a region of convective instability as the field strength is increased, so that the required $E_{0}$ is feasible. The critical values of $E_{0}$ above which the models exhibit convective instablility are summarized in Table 2.

The critical value of $E_{0}$ decreases with increasing accuracy of the inertial long wave models, with the inertialess Benney and Stokes flow dispersion relations also giving reasonable values. The field strengths in Table 2 were produced by continuation of the dominant saddle point, and additionally verified with root plots as in figure 4(a,b).

For the DNS, we employ the open-source volume-of-fluid solver $\mathpzc{Gerris}$ Popinet (2003, 2009); López-Herrera et al. (2011). We use viscous–gravity length and time scales as described by Samanta et al. Samanta et al. (2011) in order to non-dimensionalize the system, as it provides a convenient unit effective Reynolds number, embedding the parameter variation inside the dimensionless film thickness and surface tension coefficient instead. We consider a finite computational domain of length $L_{0}=400\ell$ in the streamwise variable, with an extent of $50\ell$ in the perpendicular $z$-direction. For an undisturbed interface, the region $\ell<z<50\ell$ of the computational domain is occupied by the second fluid, in our case we take air (at room temperature). For DNS, Region II is hydrodynamically active; its fluid properties are such that this has a negligible effect, and a comparison with long-wave models derived assuming that Region II is hydrodynamically passive remains valid. We set the density in Region II to $1.17$ kg/m3 and take a dynamic viscosity of $1.83\times 10^{-3}$ kg/m$\cdot$s, both of which are significantly lower than their liquid film counterparts given in Table 1. A spatial filtering technique smoothing out the interface is hence employed to improve the performance of the underlying projection solver. The no-slip and impermeability conditions are applied at the solid substrate boundary ($z=0$), and a no-stress condition is imposed on the opposite side of the computational domain at $z=50\ell$. Dirichlet boundary conditions for the voltage potential are prescribed at the inlet and outlet, chosen to give the desired field strength $E_{0}$. Neumann conditions for the voltage potential are imposed at the upper and lower boundaries of the computation domain – these are found to be an effective equivalent of the far-field conditions since we have chosen $\epsilon^{\textrm{S}}=\epsilon^{\textrm{I}}$. The Nusselt base profile, given by (11) in dimensional variables, is prescribed as the initial condition, on top of which we add a Gaussian pulse whose amplitude is $1\%$ of the liquid film thickness, with center located $150\ell$ from the inlet (far enough to prevent unwanted boundary noise/effects). Guided by the linear theory of the long-wave models, we choose the pulse width to be approximately $4\ell$ to induce the desired instability and minimize transient effects. The Nusselt solution is also used as the inlet condition, while standard outflow conditions are considered at the downstream boundary. We prescribe the maximum level of refinement around the interfacial shape, as well as close to the substrate in order to capture the details of the underlying boundary layer. Changes in the velocity field otherwise dictate the adaptive mesh refinement strategy. Considering the above, the discretized domain is described by approximately $10^{5}$ computational grid cells, with the target evolution to the first dripping events requiring several days of runtime on local high performance computing facilities. The code is executed in parallel on up to $16$ CPUs, with load balancing capabilities resulting in reasonably good scalability properties.

The results of the DNS are shown in figure 8. Panels (a,b) show the evolution of the pulse initial data for the cases of $E_{0}=0,6\times 10^{5}$ V/m, respectively, with the flow from left to right. The profiles are shifted so that the initial pulse is centred at $x=0$. Note also the vector showing the direction in which gravity is acting for these hanging films. From the height scales on the left hand sides of each panel, there is a clear difference in growth rate across these two cases, and the non-electrified dynamics reaches a nonlinear phase of evolution at around $t=4$ s. In order to determine the absolute or convective nature of the instability, we monitor the upstream edge of the wavepacket – plotted in figure 8(c) for a range of $E_{0}$. These curves were computed using the region of $(x,t)$ space for which $|h(x,t)-\ell|$ was above a given tolerance. The boundary of this region is wavy, and so the curves were computed by fitting through the farthest upstream points (similar to figure 8 in Avitabile et al. (2017)). There is an initial transient phase in which the wavepacket is convected downstream in all cases – we found this is difficult to avoid by choosing better initial data due to DNS constraints. For the non-electrified case, the wavepacket edge begins to move back upstream at $t\approx 2$ s, and similarly for $E_{0}=1.5\times 10^{5}$ V/m at $t\approx 3$ s. We are confident that these cases show linear (and nonlinear) absolute instability after a transient phase, as predicted by the low-dimensional models (see Table 2), rather than a linear convective instability followed by a nonlinear absolute instability as discussed by Delbende and Chomaz Delbende and Chomaz (1998). The choices of $E_{0}\geq 3\times 10^{5}$ V/m all showed convective instabilities (in both the linear and nonlinear phase), in agreement with the predictions of WIBL2 and the Stokes dispersion relation.

VI Direct Numerical Simulations: Dripping

In this section, we perform DNS of the full Navier–Stokes equations coupled with electrostatics on a spatial domain of length $L_{0}$, and determine parameters for which dripping occurs. We consider inflow and outflow boundary conditions to mimic an experimental set-up, unlike the spatially periodic study of Kofman et al. Kofman et al. (2018). The majority of the details of the DNS set-up are exactly as in the previous section. We initialize our simulations with the Nusselt solution, and apply a time-periodic forcing at the inlet (located at $x=0$) of the form

[TABLE]

where $\overline{u}(z)$ is the base Nusselt velocity defined in (11), and $\mathfrak{f}$ is a frequency chosen to excite the most unstable interfacial waves using predictions from the long-wave models. The strength of the time-periodic inlet perturbation, $\alpha$, was varied to reduce the transient phase of the dynamics and ensure the development of a clean wave-train close to the inlet. This form of inlet forcing was used effectively by Denner et al. Denner et al. (2018) in their computation of solitary wave profiles on overlying falling films. Furthermore, we find that the resulting solutions are robust to variations in the forcing frequency $\mathfrak{f}$. For all numerical experiments, we fix the aspect ratio by taking the system length (between inlet and outlet) to be $L_{0}=300\ell$. While the domain imposed to capture the dripping dynamics is marginally smaller than in the case of the pulse evolution, the timescale required to comprehensively investigate the instability is at least one order of magnitude larger due to a significant transient period, thus resulting in a considerably increased computational effort. The computational framework described in the previous section naturally allows for interfacial break-up, however we have opted to remove resulting drops as soon as they detach from the main body of fluid to avoid interactions with the opposite boundary, as well as to concentrate computational resources at the level of the investigated liquid film.

VI.1 Non-electrified case

For non-electrified flows we considered three choices of the inclination angle $\theta$, and varied the film thickness $\ell$. These are given in Table 3 along with absolute instability regions in $\ell$ for all of the models discussed in the present paper.

In the case of the WIBL2 model, for example, the values are obtained by taking a constant-Kapitza number slice of the surface shown in figure 5(a); the resulting curve is then considered for the three values of $\theta$. Both WIBL2 and the Stokes flow models exhibit a non-monotonic A/C threshold, in agreement with DNS of the Navier–Stokes equations. The dashes in the last row indicate that the models are convectively unstable for all choices of $\ell$. The other models do not include viscous extensional stresses.

For each case of $\theta$, we performed DNS for a range of values of $\ell$, guided by the results of the spatial stability analysis. We do not present simulation results for large values of $\ell$ near the upper bounds of the absolute instability regions; we found immediate and extensive dripping at such parameters. In agreement with the results of Kofman et al. Kofman et al. (2018) (see their figure 3), we are confident that only the lower A/C threshold can be used as a predictor of dripping; in this regime, viscous extensional stresses are not important and the critical angle for absolute instability depends monotonically on $\ell$. It is expected that the dripping angle should also be monotonic in $\ell$. It is worth noting that in our DNS we varied non-dimensional parameters by an order of magnitude (although ${\operatorname{\mathit{R\kern-1.49994pte}}}$ still remains small); they range from ${\operatorname{\mathit{R\kern-1.49994pte}}}=0.005$ and ${\mbox{{C}}}=0.0885$ for $\ell=1.3$ mm and $\theta=7\pi/8$, to ${\operatorname{\mathit{R\kern-1.49994pte}}}=0.05$ and ${\mbox{{C}}}=0.466$ for $\ell=2.5$ mm and $\theta=49\pi/60$, thus affording a fairly complete view of the physical system.

For $\theta=7\pi/8$, we find a very clear separation between three cases which are plotted in figure 9. From simulations with $\ell=1.25,1.30,1.35$ mm, we find a linear relationship between the mean film thickness and maximum wave height, approximately $h_{\textrm{max}}=1.649\ell$. This linear relationship is expected and in agreement with the results in Denner et al. (2016, 2018) for overlying films. No drips are observed for these cases, even after the coalescence of neighbouring pulses. In figure 9(a) we plot a snapshot of such a bounded wave train for $\ell=1.35$ mm; some of the pulses are wider as they result from the coalescence of two smaller pulses, but all waves very clearly attain the maximum bound sufficiently far from the inlet. Increasing the inflow film thickness to $\ell=1.40$ mm, we find similar bounded wave trains (the dotted line in figure 9(b) indicates the estimate for $h_{\textrm{max}}$), but occasionally in the transient phase of the dynamics, dripping occurs due to the coalescence of two adjacent pulses. Interactions and coalescence between pulses is common for smaller values of $\ell$, however not leading to dripping events. Further increasing $\ell$ to $1.45$ mm, we see dripping due to a different mechanism as displayed in the snapshot in figure 9(c). The individual solitary waves destabilize and drip, rather than dripping due to coalescence or pulse interaction (however, the latter is still possible in this case). The pulses appear to attain the predicted $h_{\textrm{max}}$ briefly, before further increasing in amplitude until pendant drips form. We thus define $\ell^{\textrm{drip}}$ to be the length beyond which dripping occurs due to instability of individual pulses far from the inlet (not transient effects or pulse interactions/coalescence). We believe this to be a sensible definition for this phenomenon. We give predicted ranges for $\ell^{\textrm{drip}}$ in Table 3 based on our DNS results – much longer computations would be required to refine these further.

In all cases we found that the transient dynamics were significantly long, often requiring prohibitively large timescales to observe convergence to a regular state. For $\theta=7\pi/8$, a shallow angle close to the horizontal set-up, the dynamics are relatively well-behaved and the prevailing regime is clearly apparent, with wave-trains developing not too far downstream as shown in figure 9. On the other hand, the inclinations further from horizontal exhibited more complex pulse interactions with wave-trains developing far from the inlet, and hence a more challenging detection of features at the level of the individual pulses is required to inform our previously defined metric for dripping. The predicted intervals of $\ell^{\textrm{drip}}$ for $\theta=5\pi/6$ and $\theta=49\pi/60$ in Table 3 are correspondingly much wider. We note that the regularity of dripping is a useful indicator; for $\ell>\ell^{\textrm{drip}}$, drips often occur at regular time intervals, at a fixed spatial location. For the transient dripping with $\ell<\ell^{\textrm{drip}}$, there is very little regularity in the temporal or spatial location of the dripping events. We see from Table 3 that the predicted $\ell^{\textrm{drip}}$ appears to be monotonic in $\theta$, as expected.

VI.2 Electrified case

In order to study the effect of the electric field on the dripping dynamics, we focus on the case of $\theta=5\pi/6$ and $\ell=2.5$ mm as considered in the previous section – see Table 2 for the critical values of $E_{0}$ for which the various models transition from convective to absolute instability. Since we fix the aspect ratio to $L_{0}=300\ell$, we may express the required voltage in terms of the electric field strength as $V_{0}=300\ell E_{0}$. We note from Table 3 that the film drips without an electric field. For the chosen parameters, it can be calculated that the electric field strength required to linearly stabilize the flat film solution is approximately $6\times 10^{6}$ $\emph{V}/\emph{m}$, which is double the dielectric breakdown of air and thus physically unfeasible. However, we show that dripping suppression may be attained for much smaller and physically feasible field strengths.

In figure 10, we present time series for the maximum deviation of the interface from the mean thickness, $h_{\textrm{max}}-\ell$, for a range of $E_{0}$. We show the first 50 s of longer simulations in the figure; note that the axes are reversed so that downward spikes indicate dripping events. For $E_{0}=0,1.5\times 10^{5},3\times 10^{5}$ V/m, we see very regular dripping events, both temporally and spatially (from inspection of the interfacial evolution). Delay of the first dripping event, and temporal and spatial irregularity of dripping events is clear for $E_{0}=4.5\times 10^{5},6\times 10^{5}$ V/m. This appears to be due to time-periodic disturbances that travel from the inlet downstream, disrupting the regular structure of the pulses and causing considerable pulse interaction. This behavior does not appear to be transient and thus our previously defined metric of dripping is invalid for the electrified problem. However, this transition between dripping due to loss of stability of pulses and dripping due to pulse interaction does agree well with the A/C threshold. For $E_{0}=7.5\times 10^{5},9\times 10^{5}$ V/m, bounded wave trains emerge (with only a few transient drips). Dripping suppression is obtained for $E_{0}$ at least double the convective instability threshold – we expect that the increased disparity between the linear and nonlinear phenomena is due to the nonlocality of the electric field effect. For nonlocal systems, long-range interactions between pulses that are not immediate neighbours are important, with interactions between neighbouring pulses strengthened (Lin et al., 2015; Blyth et al., 2018). Because of this, it may be possible for inlet dynamics to propagate downstream, impeding the formation of a bounded and stable pulse-train.

Finally, in figure 11 we present snapshots of the flow with $E_{0}=3\times 10^{5}$ V/m in panel (a), and $E_{0}=6\times 10^{5}$ V/m in panel (b). These were obtained from the same DNS simulations used to produce figure 10, with the snapshots taken at $t\approx 67$ s (beyond the limits of the time series presented in figure 10). The plots depict the film position, along with a pressure colour map and voltage equipotential curves in the fluid and air regions. The lower plots of each panel (these include the equipotential lines) are enlargements of the smaller rectangular regions in the upper plots, and contain five and seven waves for panels (a) and (b), respectively. For the upper panel, the field is not strong enough to suppress pinching and the leftmost finger depicted is about to break off due to a Rayleigh–Taylor instability rather than the wave merging mechanism. As expected, we see large values of the pressure near the tips of the fingers due to capillary effects. The enlargement in figure 11(b) shows a developed wave train of pulses with amplitudes much smaller than those in panel (a). Dripping is still found in this case – see the time series in figure 10(e) – but is due to pulse coalescence as is visible in the upper zoomed out plot to the left of the enlarged region.

We note that, although the A/C threshold does not agree as closely with the dripping limit as in the non-electrified case, it provides a good order-of-magnitude approximation of the required $E_{0}$ for dripping stabilization. Furthermore, a temporal linear stability does not provide a good estimate of the dripping stabilization threshold since the required $E_{0}$ for linear stability of the system is dependent on $L_{0}$.

VII Conclusions and discussion

Absolute linear instability of the flat interface solution and the highly nonlinear process of dripping are not equivalent, however, it has been found that their thresholds are close in parameter space for non-electrified film flows, but not identical Brun et al. (2015); Kofman et al. (2018). In the present work, we considered the addition of a stabilizing electric field to the problem, and performed a linear stability analysis for different long-wave models and in the Stokes flow limit. In agreement with Cimpeanu et al. (2014); Anderson et al. (2017), we found that a sufficiently strong electric field may be employed to ensure linear and nonlinear stability of the fluid interface. We also found that an absolutely unstable system will become convectively unstable if the imposed electric field is sufficiently strong. Employing DNS to study the full physical system (Navier–Stokes coupled with electrostatics) with pulse initial conditions, we showed that the A/C predictions of the long-wave models are closely aligned with the A/C behavior of the full Navier-Stokes system.

Following an extensive DNS study on long domains using inflow/outflow boundary conditions that mimic experiments, we observed a separation of the fully nonlinear dynamics of hanging film flows into three types: no dripping, dripping due to coalescence and pulse interaction, and dripping due to instability of individual pulses. For both non-electrified and electrified flows, we found that the threshold of absolute instability was a good predictor for the onset of the latter of these three phenomena. For the non-electrified problem, dripping due to coalescence appears to be transient, and thus it is reasonable to define the dripping threshold as the loss of stability of pulses. However, for electrified flows, the regime of dripping due to coalescence and pulse interaction is large in parameter space, and the behavior appears to be persistent – accordingly, we cannot define the dripping limit as in the non-electrified case. The A/C threshold remains a good order-of-magnitude estimate for the critical electric field strength required to prevent dripping (the temporal linear theory provides no such prediction). We expect this to be true for hanging film flows with a range of external effects imposed, e.g. magnetic fields, thermal effects, surfactants. This study shows that long-wave modelling approaches and spatial linear theory can provide valuable predictions for very nonlinear processes, so that computational and experimental resources can be more appropriately targeted. There are no details in the experimental work of Brun et al. Brun et al. (2015) whether the dripping events observed are due to coalescence of waves, but the increased number of drips deeper into the absolute instability regime is congruent with our findings. The inclination angles considered in the experiments of Charogiannis and Markides Charogiannis and Markides (2016) and Charogiannis et al. Charogiannis et al. (2018) are well below the values for which we predict absolute instability or dripping from our 2D study – their results are in agreement since no dripping is observed. Further experimental work at more extreme inclinations, such as those considered in Brun et al. (2015), would be of great interest. For the electrified flow, extension of the models to allow for weakly conducting fluid phases is warranted for comparison with future experiments and also for relevance to applications – Papageorgiou and Petropoulos Papageorgiou and Petropoulos (2004) observed that a parallel field can even become destabilizing for some choices of conductivities and permittivities.

The threshold of absolute instability has been found to successfully predict the onset of nonlinear phenomena in other interfacial flow problems. Rietz et al. Rietz et al. (2017) performed an experimental study of a liquid film on the exterior of a vertical rotating cylinder. The problem is parameterized by a Reynolds number and a ratio of body forces. Reduced-order modelling approaches are viable, and the authors obtained regions of absolute and convective instability for a WIBL model. Similarly to the present study, they found three regimes of the nonlinear dynamics. No dripping was observed for convectively unstable systems. For flows just beyond the A/C threshold, they found that 2D wave-fronts destabilized into rivulets which emitted drips – this is due to wave coalescence on the rivulets. The onsets of the linear and nonlinear phenomena agree most closely for $O(1)$ Reynolds numbers. Going further into the absolute instability regime decreases the inception length of the dripping process, until a regime is reached in which drips form immediately at the inlet. Vellingiri et al. Vellingiri et al. (2015) studied the problem of a gravity-driven thin liquid film sheared by a counter-current turbulent gas flow. For sufficiently large gas flow rates, the usually downward falling film begins to flow back upstream towards the inlet, a phenomenon known as “flooding”. They find that the “flooding point”, the critical value of the gas flow rate at which standing waves form on the interface, is close to the upper limit of absolute instability.

It would also be of interest to compare the A/C predictions and DNS results with the matched asymptotics theory of Indeikina et al. Indeikina et al. (1997) for individual static (pinned) rivulets. The fully wetting rivulet structures observed in the experiments of Charogiannis and Markides Charogiannis and Markides (2016) and Charogiannis et al. Charogiannis et al. (2018) at moderate inclination angles (not close to vertical) have wavelength given by the most unstable spanwise mode, resulting from the balance of destabilizing cross-stream gravity and stabilizing surface tension; this is $\sqrt{2}$ larger than the wavelength of the static rivulets considered in Indeikina et al. (1997) which correspond to the neutrally stable spanwise mode. Lin et al. Lin et al. (2012) considered both the dynamics on a static pinned (neutral-mode) rivulet as in Indeikina et al. (1997), and the development of a front on a thin precursor film (no dewetting or contact lines), finding in the latter case the most unstable mode dominates. The collection of experimental and numerical studies indicate a shift to shorter (stable) wavelengths when the film dewets. However, we also note the hysteresis effect observed in Indeikina et al. (1997), whereby the pinning points of a rivulet with an imposed fluid flux contract as the inclination is increased, but remain fixed as the angle is decreased. This indicates that the rivulets may exist for a range of base widths and static contact angles for a given flow rate. Can fully wetting rivulet structures survive into the dripping regime, or do they give way to one or more static rivulets as observed in Rothrock (1968); Indeikina et al. (1997)? The edges of the fluid film are observed to contract towards the centre in experiments in Charogiannis and Markides (2016); Charogiannis et al. (2018). If dripping occurs in the wedge shaped transition between the inlet and the single rivulet, the mass and fluid flux of the latter are unknown. Thus, it is unclear how to extract a useful prediction from the work of Indeikina et al. Indeikina et al. (1997) for the present 2D setting, since it is not obvious what the dimensions of an effective underlying static rivulet should be, in either the wetted or dewetted regime. We thus leave such comparisons to future work on the 3D problem.

The electrified aspect of the fully 3D problem is also an interesting extension. For this, we allow the electric field to be skewed at an angle $\phi$ to the streamwise flow direction, yet still parallel to the substrate surface. The far-field boundary condition for the voltage potentials in this problem is then

[TABLE]

If $\phi=0$, we recover the arrangement considered in the present work, and if $\phi=\pi/2$ the field is directed purely in the transverse direction. Following through with the long wave analysis, we find that the electric field contribution, which has Fourier symbol $|\xi|^{3}$ in the 2D formulation, now appears as the operator with symbol

[TABLE]

for wavenumber vector $(\xi_{1},\xi_{2})$. From this expression, it is apparent that the electric field has the strongest stabilizing effect on waves which are in the same direction as the undisturbed field lines, and has no influence on interfacial waves which are perpendicular, with a smooth variation in-between. Thus, a sufficiently strong electric field ensures that an interface can only develop instabilities of a select few modes. In particular, a field set up transverse to the flow direction may be employed in experiments to preserve the 2D phase of the dynamics for overlying films, i.e. impede secondary 3D instabilities, or to prevent the formation of rivulets for hanging films without affecting the streamwise dynamics. A detailed study of the fully 3D problem will be presented in future work. Additionally, work is underway on the full 3D problem with an electric field set up normal to the substrate, investigating the electrostatically-induced rivulet formation predicted in Tomlin et al. (2017) for overlying liquid films.

Appendix A Orr–Sommerfeld system

Here we provide details of the Orr–Sommerfeld system we solved to make comparisons between the full and long-wave linear theories. Linearizing (14) about the exact Nusselt solution (18), with perturbations denoted by tildes, we find

[TABLE]

with the no-slip and impermeability conditions at the substrate surface unchanged. At $z=1$, the kinematic condition, tangential stress balance, and normal stress balance become

[TABLE]

where we have used (45c) in deriving (46c). The Laplace equations (2), far-field conditions (15), and the solid–fluid boundary conditions (4) remain unchanged, whereas the electrostatic boundary conditions at the free surface (6) linearize to

[TABLE]

The problem is rearranged to obtain a system for $\tilde{h}$, $\tilde{w}$, and $\tilde{V}^{\iota}$ alone. This is achieved by eliminating the velocity $\tilde{u}$ and pressure $\tilde{p}$ from (45) to find

[TABLE]

and in turn eliminating $\tilde{p}$ in the normal stress balance (46c) gives

[TABLE]

Lastly, a derivative of the tangential stress boundary condition with respect to $x$ yields

[TABLE]

We consider perturbations of the form

[TABLE]

where c.c. denotes the complex conjugate; the wavenumber $\xi$ is allowed to be complex. Defining $D\equiv\mathrm{d}/\mathrm{d}z$ casts equation (48) and the Laplace equations for the voltages into

[TABLE]

The boundary conditions at the substrate $z=0$ are

[TABLE]

The far-field conditions become

[TABLE]

and the interfacial conditions at $z=1$ are

[TABLE]

The electrostatics problem comprising of (53,54c,d,55,58) can be solved analytically,

[TABLE]

where

[TABLE]

Thus, we can analytically give the electric field contribution in (57) at $z=1$,

[TABLE]

(recalling the relationship (33) between $\operatorname{\mathit{W\kern-3.00003pte}}$ and $\operatorname{\mathit{W\kern-3.00003pte}}^{\prime}$), so that only the fluid dynamics problem in the region $0\leq z\leq 1$ remains. The leading-order term of (61) for $\xi\ll 1$ is $2{\operatorname{\mathit{W\kern-3.00003pte}}}^{\prime}\xi^{3}\operatorname{sign}(\xi_{\textrm{r}})$ in agreement with the long-wave calculations. This Orr–Sommerfeld system was solved using the bifurcation and continuation software AUTO-07P, as discussed in Kalliadasis et al. (2012). With this, we produced the solid line in figure 3 showing the dependence of the critical wavenumber $\xi_{\textrm{c}}$ on the Reynolds number and electric field strength.

Appendix B Stokes flow

At zero Reynolds number, the Orr–Sommerfeld equation (52) reduces to

[TABLE]

with boundary conditions (54a,b,56b) and (57) at ${\operatorname{\mathit{R\kern-1.49994pte}}}=0$. Equation (62) has the general solution

[TABLE]

The boundary conditions (54a,b) give $B_{3}(\xi)=0$ and $B_{1}(\xi)+\xi B_{4}(\xi)=0$, respectively, and conditions (56b,57) at $z=1$ imply that

[TABLE]

With these, the dispersion relation in the Stokes flow limit can be obtained from (56a) as

[TABLE]

where the electric field contribution is given by (61). With the small wavenumber approximations

[TABLE]

and that of the electric field term given above, we recover, to leading order, the Benney dispersion relation (35) at zero Reynolds number.

Acknowledgements.

RJT acknowledges a PhD studentship from the Engineering and Physical Sciences Research Council (EPSRC). DTP was partly supported by EPSRC grant EP/L020564/1. RC acknowledges the support of the Mathematical Institute of the University of Oxford, as well as the staff actively maintaining the Imperial College Research Computing Service (DOI:10.14469/hpc/2232).