1. Introduction
Porous materials, owing to their high surface-to-volume ratio, serve as vital platforms in diverse applications, including Micro-Electro-Mechanical Systems/Nano-Electro-Mechanical Systems (Wang et al. Reference Wang, Zheng, Li, Huck and Sirringhaus2004), biomimetics (Vincent Reference Vincent2009), separation and sensing technologies (Tasangtong et al. Reference Tasangtong, Pholsaptanakorn, Tapsawut, Wiwekwin, Mettakoonpitak, Na Nongkhai and Sameenoi2024) and energy storage devices (Pourzolfaghar et al. Reference Pourzolfaghar, Wang, Jiang, Kositsarakhom, Jirasupcharoen, Suwantri, Jyothi, Prabhakaran and Li2024). Interfacial instabilities in such media critically influence functional performance, particularly in lab-on-chip diagnostics (Kapoor et al. Reference Kapoor, Ramamoorthy, Sundaramurthy, Vaishampayan, Sridhar, Balasubramanian and Ponnuchamy2024), organ-on-chip systems (Tiwari & Chauhan Reference Tiwari and Chauhan2019) and solid-state batteries (Hüter et al. Reference Hüter, Fu, Finsterbusch, Figgemeier, Wells and Spatschek2017). These instabilities, often triggered by thermomechanical or electrostatic gradients, yield self-organised patterns (Chou & Zhuang Reference Chou and Zhuang1999; Mukherjee et al. Reference Mukherjee, Sharma, Patil, Faruqui and Pattader2008), central to miniaturised fabrication strategies. Among these mechanisms, electrohydrodynamic (EHD) instabilities stand out for enabling remote, tuneable pattern formation on soft, confined substrates. In solid-state electrochemical interfaces, such instabilities can deteriorate wettability (Chen et al. Reference Chen, Che, Feng, Liao, Wang, Yin and Ma2019) and thereby undermine performance and reliability.
However, most prior studies treat the substrate either as rigid or purely viscoelastic, neglecting the interplay between porous substrate and solid viscoelasticity. Such coupling is critical in soft, confined environments – ranging from solid-state electrochemical systems to hydrated biological tissues – where both electrostatic and mechanical responses co-regulate interface evolution. This study develops a unified EHD framework that explicitly accounts for deformable solid–porous (SP) interface coupled to viscoelastic solids, addressing this gap.
The foundational works by Herminghaus (Reference Herminghaus1999) and Chou and co-researchers (Chou, Zhuang & Guo Reference Chou, Zhuang and Guo1999; Deshpande, Sun & Chou Reference Deshpande, Sun and Chou2001), theoretically and experimentally established that thin dielectric polymer films, when confined between conductive plates and subjected to an external electric field, can undergo spontaneous electrostatic destabilisation. Building on this, Schäffer et al. (Reference Schäffer, Albrecht, Russell and Steiner2000) demonstrated that polymer–air interfaces could be patterned into periodic columnar structures under a direct current (DC) field, indicating the robustness of EHD instabilities as a non-contact patterning technique. Later studies expanded the patterning repertoire by incorporating chemically or topographically structured electrodes, enabling a wide array of tuneable micro- and nanostructures with a large-range ordering (Verma et al. Reference Verma, Sharma, Kargupta and Bhaumik2005; Wu & Russel Reference Wu and Russel2005; Dickey et al. Reference Dickey, Gupta, Leach, Collister, Willson and Russell2006; Pattader et al. Reference Pattader, Banerjee, Sharma and Bandyopadhyay2011).
The physical mechanism underlying these instabilities involves Maxwell stresses generated due to charge accumulation at the deformable dielectric interface. These stresses scale nonlinearly with the applied voltage and are further modulated by surface tension, leading to either stable or unstable configurations depending on system parameters such as film thickness, dielectric contrast and confinement geometry (Saville Reference Saville1997; Pease & Russel Reference Pease and Russel2003; Verma et al. Reference Verma, Sharma, Kargupta and Bhaumik2005; Wu & Chou Reference Wu and Chou2005). The viscoelastic liquid films with frequency-dependent elasticity exhibit accelerated deformation and longer relaxation times (Tomar et al. Reference Tomar, Shankar, Sharma and Biswas2007; Bandyopadhyay, Reddy & Sharma Reference Bandyopadhyay, Reddy and Sharma2012). However, viscoelastic solids display fundamentally different instability modes. Such systems exhibit two distinct regimes: a viscous-dominated regime with unconditionally unstable long-wave modes, and an elasticity-dominated regime characterised by finite-wavenumber instabilities that require a threshold voltage for onset (Sarkar, Sharma & Shenoy Reference Sarkar, Sharma and Shenoy2008; Arun et al. Reference Arun, Sharma, Pattader, Banerjee, Dixit and Narayan2009). The transition between these regimes is captured by the dimensionless ratio
$\gamma / \mu h$
, where
$\gamma$
is the surface tension,
$\mu$
the shear modulus and
$h$
the film thickness (Arun et al. Reference Arun, Sharma, Pattader, Banerjee, Dixit and Narayan2009). Additionally, studies on leaky dielectrics (Pease & Russel Reference Pease and Russel2003; Mondal, Kumar & Bandyopadhyay Reference Mondal, Kumar and Bandyopadhyay2013), anisotropic systems (Mondal et al. Reference Mondal, Ghosh, Chaudhuri and Bandyopadhyay2018), and alternating current (AC)-driven or thermocapillary-EHD modes (Roberts & Kumar Reference Roberts and Kumar2009; Nazaripoor et al. Reference Nazaripoor, Koch, Sadrzadeh and Bhattacharjee2016) have expanded the design space for nanoscale patterning.
Multilayer polymer systems further extend this potential by enabling embedded, encapsulated and phase-inverted miniaturised micro-structures due to the presence of multiple unstable interfaces and reduced interfacial tension (Lin et al. Reference Lin, Kerle, Baker, Hoagland, Schärfer, Steiner and Russell2001; Morariu et al. Reference Morariu, Voicu, Schäffer, Lin, Russell and Steiner2003; Dickey et al. Reference Dickey, Gupta, Leach, Collister, Willson and Russell2006; Bandyopadhyay et al. Reference Bandyopadhyay, Sharma, Thiele and Reddy2009; Reddy, Bandyopadhyay & Sharma Reference Reddy, Bandyopadhyay and Sharma2012). These configurations allow for precise control over instability wavelength and morphology through variations in loss and storage moduli, layer thickness and dielectric contrast (Bandyopadhyay et al. Reference Bandyopadhyay, Sharma, Thiele and Reddy2009). Additionally, patterned electrodes enable lateral tuning of structure dimensions (Pattader et al. Reference Pattader, Banerjee, Sharma and Bandyopadhyay2011; Mahajan et al. Reference Mahajan, Hutter, Steiner and Oppenheimer2013; Lee et al. Reference Lee, Jung, Kang and Lee2016; Liu et al. Reference Liu, Li, Yu, Wang, Liu and Desmulliez2016).
In confined bilayers, the replacement of free surfaces with polymer–polymer interfaces further reduces the instability length scale (
$\lambda$
). Since polymer–polymer interfacial tension is typically an order of magnitude lower than the surface tension at polymer–air interfaces (Bandyopadhyay, Sharma & Shankar Reference Bandyopadhyay, Sharma and Shankar2008). Application of external AC or DC fields to such confined surfaces – particularly those incorporating leaky dielectric fluids – substantially reduces pattern wavelengths (Thaokar & Kumaran Reference Thaokar and Kumaran2005; Ozen et al. Reference Ozen, Aubry, Papageorgiou and Petropoulos2006; Gambhire & Thaokar Reference Gambhire and Thaokar2011; Mondal & Bandyopadhyay Reference Mondal and Bandyopadhyay2014). Experiments (Wang, Zhang & Zhao Reference Wang, Zhang and Zhao2011) reveal that confined elastic films produce creases and craters with
$\lambda \sim h$
or
$1.5h$
, in contrast to
$\lambda \sim 3h$
in polymer–air systems. These morphological transitions – wrinkling to creasing to cratering – are governed by the dimensionless parameter
$\gamma /\mu h$
(Wang & Zhao Reference Wang and Zhao2013), reflecting the viscoelastic properties of the film and aligning with similar transitions in other deformable systems. Moreover, leaky dielectric interfaces enable pattern formation at lower field strengths due to enhanced interfacial charge accumulation (Mondal & Bandyopadhyay Reference Mondal and Bandyopadhyay2014). Further, electrokinetic instabilities in confined bilayer flows – particularly in microchannels – are modulated by fluid rheology (Newtonian/non-Newtonian/viscoelastic) and field type (DC or AC), impacting mixing, phase separation and droplet dynamics (Ozen et al. Reference Ozen, Aubry, Papageorgiou and Petropoulos2006; Zahn & Reddy Reference Zahn and Reddy2006; Mandal et al. Reference Mandal, Ghosh, Bandopadhyay and Chakraborty2015; Navarkar et al. Reference Navarkar, Amiroudine, Mayur and Demekhin2015; Dutta et al. Reference Dutta, Ghosh, Pattader and Bandyopadhyay2019; Kaykanat & Uguz Reference Kaykanat and Uguz2023).
Recent studies have extended the scope of EHD instability analyses to confined porous systems, revealing the significant influence of substrate permeability on flow destabilisation. Linear stability analyses have shown that the critical Reynolds number for instability in porous-supported flows can drop substantially, from
${\textit{Re}}_{cr} = 5772$
in classical non-porous Poiseuille flow to as low as
${\textit{Re}}_{cr} = 2500$
– underlining the destabilising role of porous substrates (Chang, Chen & Straughan Reference Chang, Chen and Straughan2006). The stability of such systems is governed by several parameters, including viscosity stratification (El-Sayed, Amer & Alfayzi Reference El-Sayed, Amer and Alfayzi2023), the Darcy number (
$Da$
, the ratio of permeability to the square of the porous-layer thickness) (Sharma, Aleria & Bera Reference Sharma, Aleria and Bera2024) and stress-jump conditions, which intensify shear perturbations at fluid–porous boundaries (Barman, Aleria & Bera Reference Barman, Aleria and Bera2024). Incorporation of porous media into microscale flows enhances mixing, emulsification and interfacial transport – as shown by Kumar and co-researchers (Goyal et al. Reference Goyal, Kumar, Bandyopadhyay, Usha and Banerjee2013; Kumar et al. Reference Kumar, Goyal, Banerjee and Bandyopadhyay2013a
; Kumar et al. Reference Kumar, Usha, Banerjee and Bandyopadhyay2013b
) – through interactions between porosity and the multiphase dynamics. In anisotropic porous substrates, depth-dependent permeability gradients further modulate flow stability, leading to distinct instability regimes as a function of
$Da$
(Deepu, Anand & Basu Reference Deepu, Anand and Basu2015). Importantly, EHD forcing amplifies this destabilisation, making porous supports highly effective in tuning interfacial responses in liquid films (Sadiq & Usha Reference Sadiq and Usha2010).
Further, the rheology of confined fluids – particularly non-Newtonian and viscoelastic characteristics – enhances interfacial instabilities in porous systems by facilitating momentum transport and altering the balance between elasticity and permeability, thereby significantly modifying the critical Reynolds number (Sengupta & De Reference Sengupta and De2019; Bharathi & Kudenatti Reference Bharathi and Kudenatti2022). In practical scenarios like enhanced oil recovery, viscoelastic polymer displacements trigger ‘Saffman–Taylor instabilities’ even at low Reynolds number (
${\textit{Re}}$
), aiding in the mobilisation of trapped fluids (Jung, Shakeri & Seemann Reference Jung, Shakeri and Seemann2023). Interestingly, when confined within anisotropic porous media, these instabilities become increasingly stable with rising Deborah number (
$De$
), underscoring the critical role of elasticity in tuning the interfacial dynamics (Shokri, Kayhani & Norouzi Reference Shokri, Kayhani and Norouzi2018).
While much of the literature considers rigid porous substrates, deformable porous layers more accurately represent real-world soft systems, such as hydrated gels, biofilms and compressible polymer electrolytes. Unlike rigid layers, deformable substrates exhibit non-uniform stress distributions and coupling between solid deformation and fluid flow (Barry, Parker & Aldis Reference Barry, Parker and Aldis1991; Alexiou & Kapellos Reference Alexiou and Kapellos2013), which are key to understanding cardiovascular dynamics (Sreenadh et al. Reference Sreenadh, Prasad, Vaidya, Sudhakara, Krishna and Krishnamurthy2017), blood perfusion in tissues and synovial lubrication in joints (Salmi et al. Reference Salmi, Madkhali, Arif, Alharbi and Malik2022). Recent elastohydrodynamic models show that even moderate pressure gradients can significantly deform such porous layers, altering flow fields (Kumar & Sekhar Reference Kumar and Sekhar2020). A striking example is the work by Vanarse et al. (Reference Vanarse, Thakur, Ghosh, Parmar and Bandyopadhyay2024), which explored how ‘Rayleigh–Taylor’ and ‘Saffman–Taylor’ instabilities coexist and interact during viscous fluid invasion into deformable porous media, demonstrating the rich multi-physics of these systems.
Schematic representation of an external electric field (voltage,
${\psi }_{0}$
) mediated thin confined bilayer consisting of a viscoelastic solid gel atop a saturated porous medium undergoing small interfacial deformation. The local and mean thicknesses of the porous layer are
${h}_{0}$
and
$h( x,t )$
, respectively. The electrodes are separated by distance, d.
Despite extensive studies on EHD instabilities in soft confined systems, the coupled influence of porous medium-moderated solid viscoelasticity remains largely unexplored, particularly in the context of interfacial wettability (Huo et al. Reference Huo, Luo, Thangadurai, Guo, Nan and Sun2020; Wang et al. Reference Wang, Li, Li, Wang, Zhu, Sun and Yi2024), miniaturised self-organised pattern formation and its implications for emerging technologies. Soft-solid-saturated porous materials such as hydrated polymer gels (Mori, Shimizu & Hara Reference Mori, Shimizu and Hara2012; Xu, Li & Zhang Reference Xu, Li and Zhang2013), endothelial surface layers in biological tissues (Khaled & Vafai Reference Khaled and Vafai2003; Liu, Patil & Narusawa Reference Liu, Patil and Narusawa2007), engineered soft materials (Kashani & Cho Reference Kashani and Cho2024) and porous solid electrolytes used in next-generation batteries (Liu et al. Reference Liu2018; Thomas-Alyea Reference Thomas-Alyea2018), can all be described as a viscoelastic solid-saturated porous framework. Understanding the interfacial wettability in such framework is crucial for designing stable interfaces in hydrogel-based biomedical systems (Zahin et al. Reference Zahin, Al Barghouthi and Dickerson2024), thin-film coatings (Al Akhrass et al. Reference Al Akhrass, Reiter, Hou, Yang, Chang, Chang, Wang and Yang2008) and flexible electronics (Zhang, Gu & Chen Reference Zhang, Gu and Chen2023), thereby improving interfacial stability and overall device performance.
To address this gap, the present study emulates a deformable interface between a polymer gel and saturated porous substrate confined under an externally applied electric field across the gel layer, as schematically illustrated in figure 1. These confined configurations, particularly saturated porous medium-supported systems, inherently support the emergence of encapsulated porous patterns – whose formation would be governed by interfacial stresses and substrate deformation.. In our formulation, the polymeric gel is modelled as a zero-frequency Kelvin–Voigt viscoelastic solid characterised by constant storage and loss moduli (Kumar & Matar Reference Kumar and Matar2004; Matar, Gkanis & Kumar Reference Matar, Gkanis and Kumar2005; Ghosh, Bandyopadhyay & Sharma Reference Ghosh, Bandyopadhyay and Sharma2018; Grazioli et al. Reference Grazioli, Verners, Zadin, Brandell and Simone2019; Qin et al. Reference Qin, Yang, He, Zhang and Hu2023), saturating a Brinkman porous medium (Goyal et al. Reference Goyal, Kumar, Bandyopadhyay, Usha and Banerjee2013; Straughan Reference Straughan2021, Reference Straughan2025), while the saturated porous layer is modelled using a Darcy–Brinkman–Kelvin–Voigt similar framework (Straughan Reference Straughan2023, Reference Straughan2025). Both layers are considered purely dielectric; the term ‘dielectric’ here refers to their electrical permittivity rather than their mechanical state. Several previous studies (Kumar & Matar Reference Kumar and Matar2004; Matar et al. Reference Matar, Gkanis and Kumar2005; Bandyopadhyay et al. Reference Bandyopadhyay, Reddy and Sharma2012; Ghosh et al. Reference Ghosh, Bandyopadhyay and Sharma2018) have demonstrated that cross-linked elastomeric dielectrics exhibiting viscoelastic solid-like behaviour can be effectively modelled within a continuum EHD framework using Kelvin–Voigt-type constitutive relations under small-deformation conditions.
The present formulation extends classical interfacial EHD theory by incorporating the coupled effects of viscoelasticity and porous medium, thereby generalising the instability framework to the deformable interfaces of the soft-condensed systems. We theoretically investigate the salient features of electric-field-induced instabilities at such an interface between the soft-solid and the soft-solid-saturated porous substrates. By systematically exploring the role of the SP interfacial deformation in such bilayer systems, we aim to elucidate how the interplay between confinement, viscoelastic and Maxwell stresses and permeability dictates the spatio-temporal evolution of instability modes. The hydrodynamic and electrostatic governing equations are linearised and coupled with appropriate boundary conditions to characterise the characteristic length and time scales of emerging patterns. Numerical results derived from a generalised linear stability analysis (GLSA) are validated in the long-wave asymptotic limit for several representative regimes. The theoretical framework developed in this study reveals that the electric-field-mediated bilayer system integrates three synergistic mechanisms that govern interfacial instability and pattern formation. First, the interplay between substrate permeability and the rheological response of the saturated porous layer – combined with the applied electric field – enables tuneable control over instability growth rates and dominant wavelengths. Second, the viscoelastic film exhibits two distinct regimes: a viscous-dominated mode that is unconditionally unstable and an elastic-dominated mode characterised by finite-wavenumber bifurcations. The onset of the finite-wavenumber mode dictates interfacial wettability, while the emerging dominant wavelengths determine the degree of pattern miniaturisation. Third, geometric confinement and dielectric contrast drive the system toward shorter-wavelength instabilities compared with free-surface analogues, enabling refined morphological control. These insights not only elucidate the mechanisms underlying wettability transitions but also establish the operational limits for maintaining interfacial stability in flexible and solid-state electrochemical systems. Most notably, the porous medium itself emerges as a key design lever – especially in the short-wavelength regime – facilitating the precise tuning of interfacial patterns down to the sub-micron and nanoscale, which remains challenging with conventional fabrication approaches.
The remainder of the paper is structured as follows. Section 2 presents the theoretical formulation, including governing equations and boundary conditions. Non-dimensionalisation, linear stability analysis and the numerical methodology are detailed in §§ 3–5. Simulation results and physical interpretations are discussed in § 6, and concluding insights are summarised in § 7.
2. Theoretical formulation
Figure 1 schematically shows an interface between a viscoelastic solid gel and the same solid gel-saturated porous medium of infinite lateral extent, in the presence of an external electric field applied perpendicularly with a constant potential of
${\psi }_{0}$
. The permeability of the porous medium is denoted as
$\kappa$
. The loss (or viscosity) and storage (or elasticity) moduli of the upper viscoelastic layer (l) are
${\eta }_{l}$
and
${\mu }_{l}$
, respectively. The effective loss and storage moduli of the saturated porous medium (m) are
${\eta }_{e}$
and
${\mu }_{e}$
, respectively. The subscripts (l) and (m) are used to identify the viscoelastic and porous-layer variables. The dielectric constant of the ith layer is
${\varepsilon }_{i}$
(i = l and m). The interfacial tension of the deformable interface is denoted as
$\gamma$
. The geometrical variables of the system are shown in figure 1.
A two-dimensional (2-D) Cartesian coordinate system is employed for the formulation, with the origin fixed at the porous-layer bottom electrode interface (
$z=0$
), while x- and z-coordinates are parallel and normal to the bottom electrode, respectively. The symbol t represents time, and the over dots represent the total time derivatives. The bold variables denote vectors and tensors, whereas the variables within the curly brackets, followed by the bold variables, represent the vector or tensor components of the preceding variable. The subscripts or the symbol
$\partial$
that are followed by a comma represent the partial differentiations with respect to the subscripted variable. The displacement and velocity vectors, stress tensors and hydrostatic pressures of the ith layer are denoted as
$\boldsymbol{{u}}_{i} \left \{u_{i}, w_{i}\right \}$
,
$\boldsymbol{{\dot u}}_{i}$
,
${\boldsymbol{\tau }}_i$
and
$p_i$
, respectively.
2.1. Governing equations
The dynamics of the thin bilayer considered in the present study is assumed to be thick enough for the use of continuum approximation, and thin enough to neglect the inertial and gravitational forces. Also, the magnitude of the electric field is considered much stronger than that of gravity. The layers are assumed to be incompressible and isothermal. The dynamics of the upper layer (l) is governed by the following continuity equation and equation of motion, respectively:
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{{\boldsymbol{u}}_{l}}=\boldsymbol{\nabla }\boldsymbol{\cdot }\,{{\boldsymbol{\dot {u}}}_{l}}=0, \\[-28pt] \nonumber \end{align}
\begin{align} {-}\boldsymbol{\nabla }{{p}_{l}}+\left ( {{\mu }_{l}}+{{\eta }_{l}}{{\partial }_{,t}} \right ){{{\nabla} }^{2}}{{\boldsymbol{u}}_{l}}=0. \\[2pt] \nonumber \end{align}
Equation (2.2) is obtained by considering the layer (l) as a zero-frequency linear viscoelastic solid following the Kelvin–Voigt model,
${{\boldsymbol{\tau }}_{l}}={{\mu }_{l}} ( \boldsymbol{\nabla }{{\boldsymbol{u}}_{l}}+\boldsymbol{\nabla }{{\boldsymbol{u}}_{l}}^{\mathbf{T}} )+{{\eta }_{l}} ( \boldsymbol{\nabla }{{{\boldsymbol{\dot {u}}}}_{l}}+\boldsymbol{\nabla }{{{\boldsymbol{\dot {u}}}}_{l}}^{\mathbf{T}} )$
(Kumar & Matar Reference Kumar and Matar2004; Matar et al. Reference Matar, Gkanis and Kumar2005; Patra et al. Reference Patra, Bandyopadhyay, Tomar, Sharma and Biswas2011; Ghosh et al. Reference Ghosh, Bandyopadhyay and Sharma2018). The constitutive relation considers the lth layer material’s response to be equivalent to a parallel linear connection of a dashpot (purely viscous) and a spring (purely elastic). The constitutive relation describes the rheological behaviour of any soft polymeric solids when the applied stress and subsequent deformations are small (Yoon & Huang Reference Yoon and Huang2011; Grazioli et al. Reference Grazioli, Verners, Zadin, Brandell and Simone2019; Qin et al. Reference Qin, Yang, He, Zhang and Hu2023).
The dynamics of the lower layer (m) of the viscoelastic solid saturated porous medium is modelled using mass conservation and a modified Darcy–Brinkman–Kelvin–Voigt (mDBKV) framework as
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{{\boldsymbol{u}}_{m}}=\boldsymbol{\nabla }\boldsymbol{\cdot }\,{{\boldsymbol{\dot {u}}}_{m}}=0, \\[-28pt] \nonumber \end{align}
\begin{align} {-}\boldsymbol{\nabla }{{p}_{m}}+\left ( {{\mu }_{e}}+{{\eta }_{e}}{{\partial }_{,t}} \right ){{{\nabla} }^{2}}{{\boldsymbol{u}}_{m}}-{{\kappa }^{-1}}\left ( {{\mu }_{l}}+{{\eta }_{l}}{{\partial }_{,t}} \right ){{\boldsymbol{u}}_{m}}=0. \\[2pt] \nonumber \end{align}
Equation (2.4) represents momentum conservation, where
$\boldsymbol{{u}}_{m}$
is considered the effective displacement field of the lower layer which captures the pore-averaged displacement of the soft viscoelastic phase relative to the porous skeleton (Tilton & Cortelezzi Reference Tilton and Cortelezzi2008; Carrillo & Bourg Reference Carrillo and Bourg2019; Barman et al. Reference Barman, Aleria and Bera2024). The time derivative,
$\boldsymbol{{\dot u}}_{m}$
, represents the corresponding seepage velocity. It is important to emphasise that the present formulation is restricted to the ‘small-strain, small-deformation regime’ and considers no base-state deformation or flow in either layer, i.e.
$\bar {\boldsymbol{u}}_i = 0$
. In this context, the Darcy resistance term
$\kappa ^{-1}(\mu _l+\eta _l\partial _t)\boldsymbol{u}_m$
in (2.4) correctly captures the viscoelastic drag associated with relative motion between the permeating soft solid and the porous skeleton. This is fully consistent with Darcy–Brinkman–Biot-type formulations, where the drag force depends on the relative velocity (Barry et al. Reference Barry, Parker and Aldis1991; Carrillo & Bourg Reference Carrillo and Bourg2019). In the present study, however, the porous skeleton with much higher storage modulus is assumed to be undergoing an infinitesimal strain compared with the soft gel and does not undergo finite base-state deformation (Kumar & Sekhar Reference Kumar and Sekhar2020). Therefore, even when the interface at
$z=h(x,t)$
varies in
$O(1)$
amplitude, skeleton strain remains
$O(\delta )$
with
$\delta \ll 1$
, allowing interface motion while neglecting internal skeleton strain. Hence, a single effective displacement field,
$\boldsymbol{{u}}_{m}$
, should serve the purpose at leading order. The above-mentioned consideration also enables us to implement the porous medium to be homogeneous and mesoporous, which allows the use of averaged physical continuum quantities such as
$\kappa$
,
$\xi$
,
$\varepsilon _i$
,
$\mu _i$
,
$\eta _i$
and
$\xi$
in the macroscopic equations.
Equation (2.4) is formulated by analogy to the Kelvin–Voigt constitutive relation, which has been widely used to model thin cross-linked polymer films under small deformations (Kumar & Matar Reference Kumar and Matar2004; Matar et al. Reference Matar, Gkanis and Kumar2005). In porous systems, similar analogical extensions of classical constitutive laws have been employed, such as the Darcy–Brinkman–Jeffrey model (Rudraiah, Radhadevi & Kaloni Reference Rudraiah, Radhadevi and Kaloni1990), Darcy–Oldroyd models (Tan & Masuoka Reference Tan and Masuoka2005; Zhang, Fu & Tan Reference Zhang, Fu and Tan2008) and more recently the Darcy–Brinkman–Kelvin–Voigt framework for viscoelastic convection (Straughan Reference Straughan2023, Reference Straughan2025). Our current formulation follows this lineage, replacing velocity with pore-averaged displacement to represent the effective displacement of the lower layer. The terms in (2.4) have direct physical meaning. The terms
$(\kappa ^{-1}\eta _l\partial _{,t}\boldsymbol{u}_m)$
and
$(\kappa ^{-1}\mu _l\boldsymbol{u}_m)$
, as stated above, denote Darcy resistance due to viscous dissipation and elastic restoring stress, respectively, and the terms (
$\eta _e\partial _{,t}{\nabla} ^2\boldsymbol{u}_m$
) and (
$\mu _e{\nabla} ^2\boldsymbol{u}_m$
) represent Brinkman viscous dissipation and elastic resistance, respectively. The effective stress tensor of the lower layer is represented as,
${{\boldsymbol{\tau }}_{m}}={{\mu }_{e}} ( \boldsymbol{\nabla }{{\boldsymbol{u}}_{m}}+\boldsymbol{\nabla }{{\boldsymbol{u}}_{m}}^{\mathbf{T}} )+{{\eta }_{e}} ( \boldsymbol{\nabla }{{{\boldsymbol{\dot {u}}}}_{m}}+\boldsymbol{\nabla }{{{\boldsymbol{\dot {u}}}}_{m}}^{\mathbf{T}} )$
. The inclusion of the Brinkman term in the present formulation is essential to account for interfacial shear dissipation and elastic resistance, ensuring smooth stress and velocity continuity across the deformable interface (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995; Straughan Reference Straughan2023).
Importantly, (2.4) reduces to known limits: (i) the Darcy–Brinkman model for Newtonian fluids when
$(\mu _e, \mu _l) \to 0$
; (ii) Darcy’s law when
$(\mu _e, \mu _l, \eta _e ) \to 0$
; and (iii) the Navier–Stokes–Voigt model for viscoelastic solids as
$\kappa \to \infty$
. We further validated this formulation by considering limiting scenarios (viscous, viscoelastic and porous asymptotics) where analytical dispersion relations exist (figure 2), obtaining consistent results. Thus, (2.3)–(2.4) provide a flexible yet realistic description of soft-solid-saturated porous layers.
In the case of a soft solid with a low storage modulus, such as hydrated polymer gels (Mori et al. Reference Mori, Shimizu and Hara2012; Xu et al. Reference Xu, Li and Zhang2013) or porous solid electrolytes (Liu et al. Reference Liu2018; Thomas-Alyea Reference Thomas-Alyea2018), the limited elastic resistance permits viscous deformation and partial permeation into the porous medium. This is precisely the regime considered in the present study. Within the mDBKV framework, permeability acts as a resistance term regulating the evolution of the viscoelastic displacement field inside the porous skeleton, analogous to its role in classical Darcy–Brinkman formulations. This description captures the essential features of porous networks where both the interstitial soft solid and elastic matrix contribute to transport, including polymer gels such as alginate, polyacrylamide and collagen, as well as endothelial surface layers in blood vessels (Khaled & Vafai Reference Khaled and Vafai2003; Liu et al. Reference Liu, Patil and Narusawa2007). These examples demonstrate that the mDBKV framework provides a physically grounded and versatile constitutive relation for real soft-solid-saturated porous systems. The subsequent stability analysis in § 6, therefore, directly addresses the parameter ranges relevant to such materials.
We assumed that both layers are perfectly dielectric in nature, with constant dielectric permittivity that ensures electroneutrality in the bulk and there is no free charge accumulation at the interface (Saville Reference Saville1997). In (2.2) and (2.4), charge neutrality (
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{M}_i = 0$
) in the bulk is considered. Here,
$\boldsymbol{M}_i$
represents the Maxwell stresses originating from the external electric field for the ith layer
\begin{equation} {\boldsymbol{{M}}}_i={\varepsilon }_i{\varepsilon }_0\left ({\boldsymbol{{E}}}_i{\boldsymbol{{E}}}_i-\mathrm{0.5}\left ({\boldsymbol{{E}}}_i\boldsymbol{\cdot }{\boldsymbol{{E}}}_i\right )\boldsymbol{{I}}\right ) . \end{equation}
The symbol
$\varepsilon _{0}$
represents the dielectric permittivity of free space. The notation
$\textrm {I}$
denotes the identity matrix. The variable
${\boldsymbol{{E}}}_i$
is the irrotational electric field
$ (\boldsymbol{\nabla }\times {\boldsymbol{{E}}}_i =0 )$
in the ith layer which can be expressed in terms of electric-field potentials,
$\psi _i$
, as
${\boldsymbol{{E}}}_i =-\boldsymbol{\nabla }\psi _{i}$
. Thus, both the perfectly dielectric (
$\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{{E}}}_i =0$
) layers (i = l and m) are governed by the following Laplace equation (Shankar & Sharma Reference Shankar and Sharma2004):
\begin{equation} {\nabla} ^{2} \psi _{i} =0. \end{equation}
2.2. Boundary conditions
The governing equations (2.1)–(2.6) are supplemented by the following boundary conditions. No-slip and impermeability conditions are employed at both the bottom electrode-porous (EP),
$z=0$
, and top electrode-solid (ES),
$z=d$
, interfaces
\begin{equation} {\boldsymbol{{u}}}_i={\dot {\boldsymbol{{u}}}}_i=0. \end{equation}
Continuity of displacements, normal and tangential stress balances together with the kinematic condition are implemented at the deformable SP,
$z=h (x,t )$
, interface
\begin{align} {\left (\boldsymbol{{u}}\boldsymbol{\cdot }\boldsymbol{{n}}\right )}_m&={\left (\boldsymbol{{u}}\boldsymbol{\cdot }\boldsymbol{{n}}\right )}_l, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{{n}}\boldsymbol{\cdot }\left (-p_l\boldsymbol{{I}}+{\boldsymbol{{\tau }}}_l+{\boldsymbol{{M}}}_l\right )\boldsymbol{\cdot }\boldsymbol{{n}}-\boldsymbol{{n}}\boldsymbol{\cdot }\left (-p_m\boldsymbol{{I}}+{\boldsymbol{{\tau }}}_m+{\boldsymbol{{M}}}_m\right )\boldsymbol{\cdot }\boldsymbol{{n}}&=\gamma (\mathrm{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{{n}}), \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{{n}}\boldsymbol{\cdot }\left (-p_l\boldsymbol{{I}}+{\boldsymbol{{\tau }}}_l+{\boldsymbol{{M}}}_l\right )\boldsymbol{\cdot }\boldsymbol{{t}}-\boldsymbol{{n}}\boldsymbol{\cdot }\left (-p_m\boldsymbol{{I}}+{\boldsymbol{{\tau }}}_m+{\boldsymbol{{M}}}_m\right )\boldsymbol{\cdot }\boldsymbol{{t}}&=-\frac {\xi }{\sqrt {\kappa }} ({\dot {\boldsymbol{{u}}}}_m\boldsymbol{\cdot }\boldsymbol{{t}}) , \\[-28pt] \nonumber \end{align}
\begin{align} \dot {h}+{ \dot u}_i\boldsymbol{\cdot }{\mathrm{\boldsymbol{\nabla }}}_sh &= \dot w_i. \\[2pt] \nonumber \end{align}
Here, the stress-jump coefficient,
$\xi$
, is employed to account for the spatial heterogeneity across the SP interface, (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995). The symbol,
${\mathrm{\boldsymbol{\nabla }}}_s$
is the surface gradient operator which is defined as,
$ ({\mathrm{\boldsymbol{\nabla }}}_s=\boldsymbol{{I}}-\boldsymbol{{nn}} )$
. The unit outward normal (
$\boldsymbol{{n}}$
) and tangential (
$\boldsymbol{{t}}$
) vectors are defined as,
$\boldsymbol{{n}}={\mathrm{\boldsymbol{\nabla }} (z-h )}/{\left |\mathrm{\boldsymbol{\nabla }} (z-h )\right |}$
and
$\boldsymbol{{n}}\boldsymbol{\cdot }\boldsymbol{{t}}\boldsymbol{=}0$
, respectively.
The interfacial boundary conditions employed above follow from the general stress balance relations for a deformable EHD interface. The normal and tangential stress balances, respectively, equate the total (hydrostatic + shear (viscous + elastic) + electrostatic (Maxwell)) stress jumps to the capillary pressure and shear resistance at the interface. For the deformable SP interface, the tangential stress continuity is modified through a stress-jump term
$-(\xi /{\sqrt {\kappa }}) ({\dot {\boldsymbol{{u}}}}_m\boldsymbol{\cdot }\boldsymbol{{t}} )$
, accounting for finite shear transmission across the porous layer (Beavers & Joseph Reference Beavers and Joseph1967; Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995; Vanarse et al. Reference Vanarse, Thakur, Ghosh, Parmar and Bandyopadhyay2024). The stress tensors in both layers
${\boldsymbol{\tau }}_{l}$
and
${\boldsymbol{\tau }}_{m}$
incorporate the respective storage and loss moduli, ensuring that both viscous and elastic responses are represented consistently. While the form of these interfacial boundary conditions (2.8)–(2.11) remains consistent with those used in classical fluid–porous systems. In the limiting case when the storage moduli
${(\mu }_l,{\mu }_m\to 0)$
in
$({\boldsymbol{{\tau }}}_l,\ {\boldsymbol{{\tau }}}_m)$
and Maxwell stress terms
${\boldsymbol{{(}}\boldsymbol{{M}}}_l,{\boldsymbol{{M}}}_m\boldsymbol{\to }0)$
, our boundary conditions reduce exactly to the well-established fluid–porous interfacial boundary conditions employed in earlier studies (Hirata, Goyeau & Gobin Reference Hirata, Goyeau and Gobin2007; Kumar et al. Reference Kumar, Goyal, Banerjee and Bandyopadhyay2013a
). The continuity of displacement (2.8), the kinematic condition (2.11) and the stress balances ((2.9) –(2.10)) collectively capture the characteristics of the deformable SP interface. A detailed derivation of (2.11) is provided in Appendix B.
Constant potential electrostatics boundary conditions are enforced at the bottom EP,
$z=0$
, and top ES,
$z=d$
, electrode interfaces
\begin{align} \psi _m &= 0, \\[-10pt] \nonumber \end{align}
\begin{align} \psi _l &= \psi _0. \end{align}
The normal and the tangential electric-field components are balanced at the deformable SP,
$z=h (x,t )$
, interface
\begin{align} {\varepsilon }_0{\varepsilon }_l ({\boldsymbol{{E}}}_l\boldsymbol{\cdot }\boldsymbol{{n}})-{\varepsilon }_0{\varepsilon }_m ({\boldsymbol{{E}}}_m\boldsymbol{\cdot }\boldsymbol{{n}})&=0, \\[-28pt] \nonumber \end{align}
\begin{align} {\boldsymbol{{E}}}_l\boldsymbol{\cdot }\boldsymbol{{t}}&={\boldsymbol{{E}}}_m\boldsymbol{\cdot }\boldsymbol{{t}}. \\[2pt] \nonumber \end{align}
3. Non-dimensional equations
3.1. Equations for general analysis
The governing equations (2.1)–(2.6) and the boundary conditions (2.7)–(2.14) are non-dimensionalised using the following scalings:
$(X,Z,H)=(x,z,h) / h_{0} $
,
$T= (\varepsilon _{0} \psi _{0}^{2} / \eta _e h_{0}^{2} )t$
,
$ (U, W ) = (u, w ) / h_0$
,
$ P= (h_{0}^{2}/ \varepsilon _{0} \psi _{0}^{2} )p $
,
$\varPsi = \psi / \psi _{0} $
. The consequent non-dimensional continuity equations for both the layers (i = l and m) can be written as
\begin{equation} U_{i,X}+W_{i,Z}=0. \end{equation}
The dimensionless X- and Z-components of equations of motion for the lth layer are
\begin{align} -P_{l,X}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )\left (U_{l,\textit{XX}}+U_{l,\textit{ZZ}}\right )&=0, \\[-8pt] \nonumber \end{align}
\begin{align} -P_{l,Z}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )\left (W_{l,\textit{XX}}+W_{l,\textit{ZZ}}\right )&=0. \\[10pt] \nonumber \end{align}
The dimensionless parameter,
${\textit{Ra}} = {\varepsilon _0 {\psi _0}^2} / {\mu _l h_0}$
, is the modified electric Rayleigh number signifies the ratio of the electrostatic force (
$\varepsilon _0\psi _0^2 /{h_0^2}$
) to the elastic stiffness (
$\mu _l/{h_0}$
) of the lth layer.
The dimensionless X- and Z-components of equations of motion for the mth layer are
\begin{align} -P_{m,X}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+{\partial }_{,T}\big )\left (U_{m,\textit{XX}}+U_{m,\textit{ZZ}}\right )-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )U_m&=0, \\[-8pt] \nonumber \end{align}
\begin{align} -P_{m,Z}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+{\partial }_{,T}\big )\left (W_{m,\textit{XX}}+W_{m,\textit{ZZ}}\right )-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )W_m&=0. \\[10pt] \nonumber \end{align}
The dimensionless number,
$Da=\kappa /{{h_0}^{\mathrm{2}}}$
, signifies the relative effect of the permeability of the medium to its thickness. The condition,
$Da\to 0$
, indicates an impermeable porous medium layer, whereas
$Da\to \infty$
indicates that the porous layer behaves as a viscoelastic solid layer. The other dimensionless parameters obtained in (3.2)–(3.5) are storage and loss modulus ratios,
${\mu }_r={\mu }_l/{{\mu }_e}$
and
${\eta }_r={\eta }_l/{{\eta }_e}$
, respectively.
The dimensionless electrostatics governing equation for the ith layer (i = l and m) is obtained as
\begin{equation} {{\varPsi }}_{i,\textit{XX}}+{{\varPsi }}_{i,\textit{ZZ}}=0. \end{equation}
The non-dimensional no-slip and impermeability boundary conditions at the EP,
$Z=0$
, and ES,
$Z=D$
, interfaces are
\begin{equation} U_i=W_i=0. \end{equation}
The dimensionless continuity of displacements, normal and tangential stress balances, and kinematic condition at the SP interface,
$Z=H(X, T)$
, are obtained as
\begin{align} U_m = U_l, \\[-28pt] \nonumber \end{align}
\begin{align} W_m = W_l, \\[-28pt] \nonumber \end{align}
\begin{align} & -\left ( {{P}_{l}}-{{P}_{m}} \right )+2\big ( R{{a}^{-1}}+{{\eta }_{r}}{{\partial }_{,T}} \big ){{W}_{l,Z}}-\big ( {{\left ( {{\mu }_{r}}Ra \right )}^{-1}}+{{\partial }_{,T}} \big ){{W}_{m,Z}} \nonumber \\&\qquad +\left ( {{{\varepsilon }_{m}}}/{2}\; \right )\big ( {{\varepsilon }_{r}}{{\varPsi }^{2}_{l,Z}}-{{\varPsi }^{2}_{m,Z}} \big )+\varGamma {{H}_{,\textit{XX}}}=0, \\[-28pt] \nonumber \end{align}
\begin{align} & H_{,X}{\varepsilon }_m\big ({\varepsilon }_r{{{\varPsi }}^{2}_{l,Z}}-{{{\varPsi }}^{2}_{m,Z}}\big )+{\varepsilon }_m\left ({\varepsilon }_r{{\varPsi }}_{l,X}{{\varPsi }}_{l,Z}-{{\varPsi }}_{m,X}{{\varPsi }}_{m,Z}\right ) +\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )\nonumber \\& \quad \times\left (U_{l,Z}+W_{l,X}\right ) -\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+{\partial }_{,T}\big )\left (U_{m,Z}+W_{m,X}\right )=-\chi Da^{{-\mathrm{1}}/{\mathrm{2}}}U_{m,T} , \\[-28pt] \nonumber \end{align}
\begin{align} H_{,T}+U_{i,T}H_{,X}=W_{i,T}. \\[2pt] \nonumber \end{align}
The dimensionless number,
$\mathrm{\varGamma }=\gamma h_0 \big /{{\varepsilon }_0{{\psi }_0}^{\mathrm{2}}}$
, signifies the ratio of the capillary to electrostatic stresses, whereas,
$\chi =\xi /{{\eta }_e}$
, is the dimensionless stress-jump coefficient. The other non-dimensional parameter,
${\varepsilon }_r={\varepsilon }_l/{{\varepsilon }_m}$
, is defined as the ratio of the dielectric permittivities. The dimensionless electrostatic boundary conditions at both the ES,
$Z=0$
, and EP, interfaces,
$Z=D$
, are obtained as
\begin{align} {{\varPsi }}_m &= 0 , \\[-10pt] \nonumber \end{align}
\begin{align} {{\varPsi }}_l &= \mathrm{1}. \end{align}
The non-dimensional normal and tangential electric-field component balances at the SP interface,
$Z=H(X,T)$
, are
\begin{align} {\varepsilon }_r{{\varPsi }}_{l,Z}-{{\varPsi }}_{m,Z} =0, \\[-28pt] \nonumber \end{align}
\begin{align} \left ({{\varPsi }}_{m,X}-{{\varPsi }}_{l,X}\right )+H_{,X}\left ({{\varPsi }}_{m,Z}-{{\varPsi }}_{l,Z}\right )=0. \\[2pt] \nonumber \end{align}
The governing equations (2.1)–(2.6) together with the boundary conditions (2.7)–(2.14) presented in § 2 are expressed in the non-dimensionalised form (3.1)–(3.6) and (3.7)–(3.14) in § 3. The above non-dimensional governing equations (3.1)–(3.6) and boundary conditions (3.7)–(3.14) are employed to perform linear stability analysis (LSA) in § 4. Two different approaches are followed to carry out the LSA: (i) general linear stability analysis (GLSA) and (ii) long-wave linear stability analysis (LWLSA-Analytical). Equations (3.1)–(3.14) are reduced in the long-wave limit to carry out LWLSA-Analytical, which provides an asymptotic framework to obtain the length and time scales of the instabilities. Equations (3.1)–(3.14) are employed to carry out GLSA numerically. Long-wave equations are also numerically solved to validate the eigenvalues obtained by performing GLSA in the long-wave limit. Here, we have termed the numerically solved long-wave equations as LWLSA-Numerical. The reduced long-wave equations are shown in § 3.2.
However, the reduced long-wave equations can be further simplified to a few asymptotic cases, such as electric-field-induced viscous and viscoelastic confined bilayers and electric-field-induced porous medium-moderated viscous confined bilayers implementing a certain parametric variation (refer to § 6.1), where comparatively simpler analytical solutions are available. These analytical solutions (LWLSA-Analytical) are used to check the validity of the implemented numerical methodology by comparing the eigenvalues obtained from LWLSA-Numerical.
3.2. Equations for long-wave analysis
The reduced long-wave governing equations corresponding to (3.1)–(3.5) are as follows:
\begin{align} U_{i,X}+W_{i,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -P_{l,X}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )U_{l,\textit{ZZ}}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} P_{l,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -P_{m,X}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+{\partial }_{,T}\big )U_{m,\textit{ZZ}}-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )U_m&=0, \\[-28pt] \nonumber \end{align}
\begin{align} P_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} {{\varPsi }}_{i,\textit{ZZ}}&=0. \\[2pt] \nonumber \end{align}
The reduced solid dynamics and electrostatics boundary conditions at the EP and ES interfaces,
$Z=0$
and
$Z=D$
, are
\begin{align}& U_i=W_i=0, \\[-9pt] \nonumber \end{align}
\begin{align}&\quad {{\varPsi }}_m = 0, \\[-9pt] \nonumber \end{align}
\begin{align}&\quad\, {{\varPsi }}_l =\mathrm{1}. \end{align}
The other reduced boundary conditions at
$Z=H(X,T)$
are
\begin{align} U_m &= U_l, \\[-28pt] \nonumber \end{align}
\begin{align} W_m &= W_l , \\[-28pt] \nonumber \end{align}
\begin{align} -\left (P_l-P_m\right )+\left ({{\varepsilon }_m}/{\mathrm{2}}\right )\big ({\varepsilon }_r{{{\varPsi }}^{2}_{l,Z}}-{{{\varPsi }}^{2}_{m,Z}}\big )+\mathrm{\varGamma }H_{,\textit{XX}}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+{\partial }_{,T}\big )U_{m,Z}-\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r{\partial }_{,T}\big )U_{l,Z}&=\chi Da^{{-\mathrm{1}}/{\mathrm{2}}}U_{m,T}, \\[-28pt] \nonumber \end{align}
\begin{align} H_{,T}+U_{i,T}H_{,X}&=W_{i,T}, \\[-28pt] \nonumber \end{align}
\begin{align} {\varepsilon }_r{{\varPsi }}_{l,Z}-{{\varPsi }}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \left ({{\varPsi }}_{m,X}-{{\varPsi }}_{l,X}\right )+H_{,X}\left ({{\varPsi }}_{m,Z}-{{\varPsi }}_{l,Z}\right )&=0. \\[-8pt] \nonumber \end{align}
Equations (3.15)–(3.28) are used to perform LWLSA-Numerical.
4. Linear stability analysis
The general ((3.1)–(3.14)) and reduced ((3.15)–(3.28)) non-dimensional equations are linearised by employing the following normal linear modes,
$ [U_i,W_i,H,P_i,{{\varPsi }}_i ]= [0, 0, \overline {H},{\overline {P}}_i,{\overline {{\varPsi }}}_i ]+[{\tilde {U}}_i,{\widetilde {W}}_i,\tilde {H},{\tilde {P}}_i,{\widetilde {{\varPsi }}}_i](Z)e^{{\varOmega }T+iKX}$
. We have considered that the deformations of the SP interface are small in the presence of small-amplitude perturbations to the thin films’ static and space-invariant base state to perform GLSA and LWLSA-Numerical. Consequently, the X- and Z-directional displacements are set to
${\overline {U}}_i={\overline {W}}_i=0$
, and the non-dimensional base-state film thickness is employed as
$\overline {H}=1$
. Here, the symbols
$ {\varOmega }$
and K represent the dimensionless linear growth coefficient and wavenumber of the perturbations, respectively. The over-bar (over-tilde) variables represent the base-state (perturbed-state) physical quantities. The system is unstable (stable) when
${\textit{Re}} [{{\varOmega }} ]\gt 0$
$ (Re [{{\varOmega }} ]\lt 0 )$
.
4.1. Base-state analysis
The aforementioned normal linear modes are substituted into (3.1)–(3.28) to linearise the equations about the base-state quantities by neglecting the higher-order nonlinear terms in the perturbations. In this section, we lay out the base-state equations, whereas the linearised perturbed-state equations are laid out in the following section. The base-state equations for both GLSA and LWLSA-Numerical remain the same. There are no significant effects of the solid dynamics base states on the stability of the systems as the solid displacements are set to zero. However, base-sate electrostatic potentials significantly affect the stability of the system by modifying the base-sate pressure difference between the films, obtained from the normal stress balance ((3.9) and (3.24)) as
\begin{equation} \left ({\overline {P}}_l-{\overline {P}}_m\right )+\mathrm{0.5}{\varepsilon }_m\big ({\varepsilon }_r{{\overline {{\varPsi }}}^{2}_{l,Z}}-{{\overline {{\varPsi }}}^{2}_{m,Z}}\big )=0. \end{equation}
The base-state electrostatics governing equations for the ith layer is obtained as
\begin{equation} {\overline {{\varPsi }}}_{i,\textit{ZZ}}=0. \end{equation}
The corresponding boundary conditions at the EP and ES interfaces,
$Z=0$
and
$Z=D$
, are derived as
\begin{align} {\overline {{\varPsi }}}_m &= 0, \end{align}
\begin{align} {\overline {{\varPsi }}}_l &= \mathrm{1}, \end{align}
\begin{align} {\overline {{\varPsi }}}_l-{\overline {{\varPsi }}}_m &= 0, \end{align}
\begin{align} {\varepsilon }_r{\overline {{\varPsi }}}_{l,Z}-{\overline {{\varPsi }}}_{m,Z} &=0. \end{align}
The set of ordinary differential electrostatic governing equations (4.2) is solved analytically with the boundary conditions (4.3) and (4.4) to obtain the following base-state solutions:
\begin{align} {\overline {{\varPsi }}}_l&=\frac {Z+\overline {H}\left ({\varepsilon }_r-\mathrm{1}\right )}{\overline {H}\left ({\varepsilon }_r-\mathrm{1}\right )+D}, \\[-8pt] \nonumber \end{align}
\begin{align} {\overline {{\varPsi }}}_m&=\frac {{\varepsilon }_rZ}{\overline {H}\left ({\varepsilon }_r-\mathrm{1}\right )+D}. \end{align}
4.2. General perturbed-state analysis
The linearised perturbed-state mass and momentum balance equations (3.1)–(3.5) lead to the following ordinary differential equations (ODEs):
\begin{align} iK{\tilde {U}}_i+{\widetilde {W}}_{i,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -{\tilde {P}}_lK^{\mathrm{2}}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big )\big ({\widetilde {W}}_{l,\textit{ZZZ}}-K^{\mathrm{2}}{\widetilde {W}}_{l,Z}\big )&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -{\tilde {P}}_{l,Z}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big )\big ({\widetilde {W}}_{l,\textit{ZZ}}-K^{\mathrm{2}}{\widetilde {W}}_l\big )&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -{\tilde {P}}_mK^{\mathrm{2}}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big )\big ({\widetilde {W}}_{m,\textit{ZZZ}}-K^{\mathrm{2}}{\widetilde {W}}_{m,Z}\big )-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big ){\widetilde {W}}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} & -{\tilde {P}}_{m,Z}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big )\big ({\widetilde {W}}_{m,\textit{ZZ}}-K^{\mathrm{2}}{\widetilde {W}}_m\big )-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big ){\widetilde {W}}_m =0. \\[2pt] \nonumber \end{align}
Eliminating the perturbed pressures from (4.8)–(4.11) and using (4.7), the following fourth-order bi-harmonic ODEs are obtained for the lth and mth layers
\begin{align} {\widetilde {W}}_{l,\textit{ZZZZ}}-\mathrm{2}K^{\mathrm{2}}{\widetilde {W}}_{l,\textit{ZZ}}+K^{\mathrm{4}}{\widetilde {W}}_l=0, \\[-28pt] \nonumber \end{align}
\begin{align}& \big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big ) \big ({\widetilde {W}}_{m,\textit{ZZZZ}}-2K^{2}{\widetilde {W}}_{m,\textit{ZZ}}+K^{4}{\widetilde {W}}_m\big ) \nonumber\\&\quad - Da^{-1}\big ({\textit{Ra}}^{-1}+{\eta }_r {\varOmega }\big ) \big ({\widetilde {W}}_{m,\textit{ZZ}}-K^{2}{\widetilde {W}}_m\big ) = 0. \\[2pt] \nonumber \end{align}
The linearised electrostatic governing equations for the potential of the ith layer is
\begin{equation} {\widetilde {{\varPsi }}}_{i,\textit{ZZ}}-K^{\mathrm{2}}{\widetilde {{\varPsi }}}_i=0. \end{equation}
The boundary conditions are also linearised following similar pathways. The linearised no-slip and non-permeable boundary conditions at
$Z=0$
and
$Z=D$
are
\begin{equation} {\widetilde {W}}_{i,Z}={\widetilde {W}}_i=0. \end{equation}
The linearised normal and tangential stress balance, continuity of displacements and kinematic condition at the deformed interface
$Z=H$
projecting onto the base interface
$Z=\overline H$
are
\begin{align} & -\big ( \tilde {P}_l - \tilde {P}_m \big ) + 2 \big ( {\textit{Ra}}^{-1} + \eta _r {\varOmega } \big ) \widetilde {W}_{l,Z} - 2 \big ( \left ( \mu _r Ra \right )^{-1} + {\varOmega } \big ) \widetilde {W}_{m,Z} \nonumber\\&\quad + \varepsilon _m \left ( \varepsilon _r \overline {{\varPsi }}_{l,Z} \widetilde {{\varPsi }}_{l,Z} - \overline {{\varPsi }}_{m,Z} \widetilde {{\varPsi }}_{m,Z} \right ) - \varGamma \widetilde {W}_m K^2 = 0, \\[-28pt] \nonumber \end{align}
\begin{align} & \tilde {H} K^2 \varepsilon _m \big ( \varepsilon _r \overline {{\varPsi }}_{l,Z}^2 - \overline {{\varPsi }}_{m,Z}^2 \big ) + K^2 \varepsilon _m \big ( \varepsilon _r \widetilde {{\varPsi }}_l \overline {{\varPsi }}_{l,Z} - \widetilde {{\varPsi }}_m \overline {{\varPsi }}_{m,Z} \big ) + \big({\textit{Ra}}^{-1} + \eta _r {\varOmega }\big)\nonumber\\&\quad \times \big(\widetilde {W}_{l,\textit{ZZ}} + K^2 \widetilde {W}_l\big) - \big ( (\mu _r Ra)^{-1} + {\varOmega } \big ) \big(\widetilde {W}_{m,\textit{ZZ}} + K^2 \widetilde {W}_m\big) - {\varOmega } \chi Da^{-1/2} \widetilde {W}_{m,Z} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {W}}_{m,Z} = {\widetilde {W}}_{l,Z}, \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {W}}_m = {\widetilde {W}}_l, \\[-28pt] \nonumber \end{align}
\begin{align} \tilde {H}={\widetilde {W}}_i. \\[2pt] \nonumber \end{align}
The linearised electrostatic boundary conditions at
$Z=0$
and
$Z=D$
are
\begin{align} {\widetilde {{\varPsi }}}_m &= 0, \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {{\varPsi }}}_l &= 0. \\[2pt] \nonumber \end{align}
The linearised normal and tangential electric-field component balances at
$Z=H$
projecting onto the base interface
$Z=\overline H$
are obtained as
\begin{align} {\varepsilon }_r{\widetilde {{\varPsi }}}_{l,Z}-{\widetilde {{\varPsi }}}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \left ({\overline {{\varPsi }}}_{l,Z}-{\overline {{\varPsi }}}_{m,Z}\right ){\widetilde {W}}_m+\left ({\widetilde {{\varPsi }}}_l-{\widetilde {{\varPsi }}}_m\right )&=0. \\[2pt] \nonumber \end{align}
The analytical solution obtained for the variable
${\widetilde {W}}_m$
from (4.13) is cumbersome. Also, the variables,
${\widetilde {W}}_l$
,
${\widetilde {W}}_m$
,
${\widetilde {{\varPsi }}}_l$
and
${\widetilde {{\varPsi }}}_m$
are coupled through boundary conditions (4.15)–(4.22). Thus, the whole set of governing equations (4.12)–(4.14), together with the boundary conditions (4.15)–(4.22), is solved numerically to obtain the GLSA results to envisage the time and length scales for the porous medium-moderated electric-field instabilities in viscoelastic solid confined bilayers. The details of the numerical methodology are provided in § 5.
4.3. Long-wave perturbed-state analysis
The reduced long-wave equations (3.15)–(3.28) are also linearised, employing the same normal linear modes. The reduced linearised governing equations corresponding to (3.15)–(3.19) are
\begin{align} iK{\tilde {U}}_i+{\widetilde {W}}_{i,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -{\tilde {P}}_lK^{\mathrm{2}}+\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big )\left ({\widetilde {W}}_{l,\textit{ZZZ}}\right )&=0, \\[-28pt] \nonumber \end{align}
\begin{align} {\tilde {P}}_{l,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} -{\tilde {P}}_mK^{\mathrm{2}}+\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big )\left ({\widetilde {W}}_{m,\textit{ZZZ}}\right )-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big ){\widetilde {W}}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} {\tilde {P}}_{m,Z}&=0. \\[-9pt] \nonumber \end{align}
The following reduced fourth-order bi-harmonic equations are obtained using (4.23)–(4.27)
\begin{align} {\widetilde {W}}_{l,\textit{ZZZZ}}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big ){\widetilde {W}}_{m,\textit{ZZZZ}}-Da^{-\mathrm{1}}\big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big ){\widetilde {W}}_{m,\textit{ZZ}}&=0. \\[2pt] \nonumber \end{align}
The reduced linearised electric-field potential governing equations are obtained as
\begin{equation} {\widetilde {{\varPsi }}}_{i,\textit{ZZ}}=0. \end{equation}
The linearised boundary conditions corresponding to (3.21)–(3.22) at
$Z=0$
and
$Z=D$
are
\begin{align}& {\widetilde {W}}_{i,Z}={\widetilde {W}}_i=0, \end{align}
\begin{align}&\quad\, {\widetilde {{\varPsi }}}_m = 0, \end{align}
\begin{align}&\quad\,\, {\widetilde {{\varPsi }}}_l = 0. \end{align}
The linearised boundary conditions at the deformable interface corresponding to (3.23)–(3.28) at
$Z=\overline H$
are
\begin{align} -\big ({\tilde {P}}_l-{\tilde {P}}_m\big )+{{\varepsilon }_{m}}\left ({\varepsilon }_r{\overline {{\varPsi }}}_{l,Z}{\widetilde {{\varPsi }}}_{l,Z}-{\overline {{\varPsi }}}_{m,Z}{\widetilde {{\varPsi }}}_{m,Z}\right )-\mathrm{\varGamma }{\widetilde {W}}_mK^{\mathrm{2}}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \big ({\textit{Ra}}^{-\mathrm{1}}+{\eta }_r {\varOmega }\big ){\widetilde {W}}_{l,\textit{ZZ}}-\big ({\left ({\mu }_r\textit{Ra}\right )}^{-\mathrm{1}}+ {\varOmega }\big ){\widetilde {W}}_{m,\textit{ZZ}}+ {\varOmega }\chi Da^{{-\mathrm{1}}/{\mathrm{2}}}{\widetilde {W}}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {W}}_{m,Z} &= {\widetilde {W}}_{l,Z} , \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {W}}_m &= {\widetilde {W}}_l , \\[-28pt] \nonumber \end{align}
\begin{align} \tilde {H}&={\widetilde {W}}_i, \\[-28pt] \nonumber \end{align}
\begin{align} {\varepsilon }_r{\widetilde {{\varPsi }}}_{l,Z}-{\widetilde {{\varPsi }}}_{m,Z}&=0, \\[-28pt] \nonumber \end{align}
\begin{align} \left ({\overline {{\varPsi }}}_{l,Z}-{\overline {{\varPsi }}}_{m,Z}\right ){\widetilde {W}}_m+\left ({\widetilde {{\varPsi }}}_l-{\widetilde {{\varPsi }}}_m\right )&=0. \\[-15pt] \nonumber \end{align}
Equations (4.28)–(4.30), together with the boundary conditions (4.31)–(4.38) are also solved using the same numerical methodology laid down in § 5 to obtain LWLSA-Numerical results. In certain asymptotic limits, however, the general analytical solutions of (4.28)–(4.30) become significantly more straightforward. Section 6.1 presents three such asymptotic cases used to validate the implemented numerical schemes LWLSA-Numerical and GLSA). The analytical solutions of the bi-harmonic equations for these cases are provided in Appendix A. Each solution involves twelve integration constants (
$Q_1$
–
$Q_{12}$
), which are uniquely determined by enforcing the boundary conditions in (4.31)–(4.38).
For any individual asymptotic case, expressions of
${\tilde {P}}_l$
and
${\tilde {P}}_m$
((4.24)–(4.26)), and the solutions of the variables,
${\widetilde {W}}_l$
,
${\widetilde {W}}_m$
,
${\widetilde {{\varPsi }}}_l$
and
${\widetilde {{\varPsi }}}_m$
((A1)–(A3c
)), are replaced into (4.31)–(4.35), (4.37) and (4.38) to obtain twelve homogeneous linear algebraic equations involving twelve unknowns,
$Q_1-Q_{12}$
. The determinant of the coefficient matrix of the set of these algebraic equations is equated to zero to obtain the general dispersion relation for the corresponding system. The determinant is solved analytically (LWLSA-Analytical). The derivation, subsequent calculations, and solution of the determinant are done using the commercial package Mathematica
$\mathrm{{}^{TM}}$
. The expression of the general dispersion relation for each case is too inconvenient to be provided within the text. The eigenvalue solutions yielded from this approach are expected to reproduce similar eigenvalues obtained from the LWLSA-Numerical of the corresponding asymptotic cases.
5. Numerical methodology
In this section, GLSA (LWLSA-Numerical) is performed by solving the eigenvalue problem resulting from the set of coupled Orr–Sommerfeld governing equations (4.12)–(4.14) ((4.27)–(4.29)) and boundary conditions (4.15)–(4.18) and (4.20)–(4.22) ((4.31)–(4.35) and (4.37)–(4.38)) to obtain the linear growth rate (
$ {\varOmega }$
) and the corresponding wavenumber (
$K$
) for the unstable modes (Orszag Reference Orszag1971; Weideman & Reddy Reference Weideman and Reddy2000). Subsequently, perturbation growth rate is written as,
$ {\varOmega } = -\textit{iKC}$
, with
$C$
denoting the complex wave speed, yielding
$ {\varOmega } = -\textit{iKC}_{{\textit{Re}}} + \textit{KC}_{\textit{Im}}$
. Here, the real part,
${\textit{Re}}[ {\varOmega }] = \textit{KC}_{\textit{Im}}$
, represents the exponential growth (
${\textit{Re}}[ {\varOmega }] \gt 0$
) or decay (
${\textit{Re}}[ {\varOmega }] \lt 0$
) rate of the perturbation amplitude. In contrast, the imaginary part,
$\mathrm{Im}[ {\varOmega }] = -\textit{KC}_{{\textit{Re}}}$
, corresponds to the oscillatory component or angular frequency of the disturbance. The results obtained in this study indicate that
$\mathrm{Im}[ {\varOmega }]$
is negligible (
$O(10^{-13})$
) compared with
${\textit{Re}}[ {\varOmega }]$
(
$O(10^{-3})$
), confirming that the disturbances correspond to standing-wave modes rather than travelling waves. Consequently, only the real part of
$ {\varOmega }$
, which governs the onset and strength of interfacial instability, is analysed and presented hereafter.
The spectral collocation method is used to obtain solutions numerically using the help of the commercial software package MATLAB. The method is implemented by mapping the computational domain into
$(-1, 1)$
for each layer by employing the linear transformations,
$\alpha =-\mathrm{2}Z+\mathrm{1}$
and
$\beta ={ (-\mathrm{2}Z+D+\mathrm{1} )}/{ (D-\mathrm{1} )}$
for the porous and solid layers, respectively. The transformed governing equations of the eigenvalue problem to perform GLSA are evaluated as
\begin{align} \big ({\mathrm{16}}/{{(D-\mathrm{1} )}^{\mathrm{4}}}\big ){\widetilde {W}}_{l,\beta \beta \beta \beta }-\big ({\mathrm{8}K^{\mathrm{2}}}/{{(D-\mathrm{1})}^{\mathrm{2}}}\big ){\widetilde {W}}_{l,\beta \beta }+K^{\mathrm{4}}{\widetilde {W}}_l=0, \\[-28pt] \nonumber \end{align}
$\begin{align} & \big ({ ({\mu }_r\textit{Ra} )}^{-1}-{iKC}\big )\big (16{\widetilde {W}}_{m,\alpha \alpha \alpha \alpha }-8K^{2}{\widetilde {W}}_{m,\alpha \alpha }+K^{4}{\widetilde {W}}_m\big )\nonumber \\&\qquad -Da^{-1}\big ({\textit{Ra}}^{-1}-\textit{iKC}{\eta }_r\big )\big (4{\widetilde {W}}_{m,\alpha \alpha }-K^{2}{\widetilde {W}}_m\big ), \\[-28pt] \nonumber \end{align}$
\begin{align} \big ({\mathrm{4}}/{{(D-\mathrm{1})}^{\mathrm{2}}}\big ){\widetilde {{\varPsi }}}_{l,\beta \beta }-K^{\mathrm{2}}{\widetilde {{\varPsi }}}_l=0, \\[-28pt] \nonumber \end{align}
\begin{align} \mathrm{4}{\widetilde {{\varPsi }}}_{m,\alpha \alpha }-K^{\mathrm{2}}{\widetilde {{\varPsi }}}_m=0. \\[2pt] \nonumber \end{align}
The transformed eigenvalue problem boundary conditions at the EP interface,
$\alpha = 1$
, are obtained as
\begin{align} {\widetilde {W}}_m &= 0, \end{align}
\begin{align} {\widetilde {W}}_{m,\alpha } &= 0, \end{align}
\begin{align} {\widetilde {{\varPsi }}}_m &= 0. \end{align}
The transformed eigenvalue problem boundary conditions at the SP interface,
$\alpha = -1$
and
$\beta =1$
, are derived as
\begin{align} & \big ( {{ ( {{\mu }_{r}}Ra )}^{-1}}-\textit{iKC} \big )8{{{\widetilde {W}}}_{m,\alpha \alpha \alpha }}-6{{K}^{2}}\big ( {{( {{\mu }_{r}}Ra)}^{-1}}-\textit{iKC} \big ){{{\widetilde {W}}}_{m,\alpha }}\nonumber \\&\quad -2D{{a}^{-1}}\big ( R{{a}^{-1}}-\textit{iKC}{{\eta }_{r}} \big ){{{\widetilde {W}}}_{m,\alpha }} -\big ( R{{a}^{-1}}-\textit{iKC}{{\eta }_{r}} \big )\big ( {8}/{{{ ( 1-D )}^{3}}}\; \big ){{{\widetilde {W}}}_{l,\beta \beta \beta }}\nonumber \\ &\quad +\big ( R{{a}^{-1}}-\textit{iKC}{{\eta }_{r}} \big )\big ( {6{{K}^{2}}}/{( 1-D )}\; \big ){{{\widetilde {W}}}_{l,\beta }} \nonumber \\ &\quad +2{{K}^{2}}{{\varepsilon }_{m}}\big ( {{{\bar {{\varPsi } }}}_{l,Z}} ( {{{\varepsilon }_{r}}}/{ ( 1-D )}\; ){{{\tilde {{\varPsi } }}}_{l,\beta }}-{{{\bar {{\varPsi } }}}_{m,Z}}{{{\tilde {{\varPsi } }}}_{m,\alpha }} \big )-\varGamma {{{\widetilde {W}}}_{m}}{{K}^{4}}=0, \\[-28pt] \nonumber \end{align}
$\begin{align} & {{{\widetilde {W}}}_{m}}{{K}^{2}}{{\varepsilon }_{m}}\big ( {{\varepsilon }_{r}}{{{\bar {{\varPsi } }}}_{l,Z}}^{2}-{{{\bar {{\varPsi } }}}_{m,Z}}^{2} \big )+{{K}^{2}}{{\varepsilon }_{m}}\big ( {{\varepsilon }_{r}}{{{\tilde {{\varPsi } }}}_{l}}{{{\bar {{\varPsi } }}}_{l,Z}}-{{{\tilde {{\varPsi } }}}_{m}}{{{\bar {{\varPsi } }}}_{m,Z}} \big ) \nonumber \\&\quad -\big ( {{ ( {{\mu }_{r}}Ra)}^{-1}}-\textit{iKC} \big )\big ( 4{{{\widetilde {W}}}_{m,\alpha \alpha }}+{{K}^{2}}{{{\widetilde {W}}}_{m}} \big ) +2D{{a}^{-1/2}}\chi ({{iKC}} ){{{\widetilde {W}}}_{m,\alpha }} \nonumber \\&\quad +\big ( R{{a}^{-1}}-\textit{iKC}{{\eta }_{r}} \big )\big ( \big ( {4}/{{{( 1-D)}^{2}}} \big ){{{\widetilde {W}}}_{l,\beta \beta }}+{{K}^{2}}{{{\widetilde {W}}}_{l}} \big )=0, \\[-28pt] \nonumber \end{align}$
\begin{align} \left ({\mathrm{2}}/{\left (\mathrm{1}-D\right )}\right ){\widetilde {W}}_{l,Z}-\mathrm{2}{\widetilde {W}}_{m,Z} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} {\widetilde {W}}_l-{\widetilde {W}}_m = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \left ({{\varepsilon }_r}/{\left (\mathrm{1}-D\right )}\right ){\widetilde {{\varPsi }}}_{l,\beta }-{\widetilde {{\varPsi }}}_{m,\alpha }=0 , \\[-28pt] \nonumber \end{align}
\begin{align} \big({\overline {{\varPsi }}}_{l,Z}-{\overline {{\varPsi }}}_{m,Z}\big ){\widetilde {W}}_m+\big ({\widetilde {{\varPsi }}}_l-{\widetilde {{\varPsi }}}_m\big )=0. \\[0pt] \nonumber \end{align}
The transformed eigenvalue problem boundary conditions implemented at the ES interface,
$\beta =-1$
, are
\begin{align} \widetilde {W}_l &= 0, \end{align}
\begin{align} \widetilde {W}_{l,\beta } &= 0, \end{align}
\begin{align} \widetilde {{\varPsi }}_l &= 0. \end{align}
Similarly, (4.28)–(4.35), (4.37) and (4.38) are also mapped using the same linear transformation to perform LWLSA-Numerical.
The subscripts
$\alpha$
and
$\beta$
preceded with a comma in the transformed eigenvalue problem equations denote ordinary differentiation with respect to
$\alpha$
and
$\beta$
, respectively. The base-state solutions (4.5) and (4.6) are also mapped into the transformed (−1,1) domain and replaced into (5.6), (5.7) and (5.10). Subsequently, the Chebyshev polynomials
$T_n (Z )$
are used to expand the set of transformed eigenvalue problem governing ODEs (5.1)–(5.4) and boundary conditions (5.5)–(5.11). Four ODEs with twelve boundary conditions form a
$ (\mathrm{4}N )\times (\mathrm{4}N )$
matrix for N Chebyshev polynomials. The derived matrix is used to obtain the eigenvalues. The accuracy and the simulation time trade-off of the code is achieved by using
$\mathrm{\sim }$
30 Chebyshev polynomials and efficaciously removing any spurious eigenvalues. The number of Chebyshev polynomials is varied to confirm that the solutions remain the same and are independent of N. The corresponding convergence data and comparative plots are provided in the Electronic Supplementary Information (available at https://doi.org/10.1017/jfm.2026.11475). All results corresponding to the LWLSA-Numerical and GLSA presented in this study are obtained using the spectral collocation-based numerical scheme outlined in this section.
The time and length scales of the instability are determined by evaluating the dominant growth coefficient,
${{{\varOmega }}}_{{M}}$
, by finding the maximum value of
$ {\varOmega }$
, and the corresponding wavenumber,
$K_{\kern-1pt {M}}$
, associated with the fastest growing mode. The dominant wavelength (
${{\varLambda }}_{{M}}={2\pi }/{K_{\kern-1pt {M}}}$
) and inverse of the dominant growth coefficient (
${{\varOmega }}^{-1}_{{M}}$
) dictate the length and time scales of such porous-moderated electric-field-induced confined viscoelastic bilayer instabilities.
6. Results and discussion
6.1. Validations
We begin by validating the theoretical framework established in the previous section. As mentioned in § 3.1, three asymptotic scenarios (cases 1–3) are emulated, as illustrated in figure 2, where simpler analytical solutions exist. Cases 1–3 correspond to electric-field-induced viscous confined bilayers (
$Da\to \infty$
,
${\textit{Ra}}\to \infty$
,
$\chi =0$
), electric-field-induced viscoelastic confined bilayers (
$Da\to \infty$
,
$\chi =0$
) and electric-field-induced porous medium-moderated viscous confined bilayers (
${\textit{Ra}}\to \infty$
,
$\chi \not =0$
), respectively. Case 4 represents the current scenario.
Results obtained from LWLSA-Analytical, LWLSA-Numerical and GLSA showing the variation of
$ {\varOmega }$
with
$K$
. The circular symbols are obtained from the analytical solutions (LWLSA-Analytical), whereas the solid lines are obtained from the numerical solutions (LWLSA-Numerical and GLSA). (ai–ci) Correspond to the schematic representation of asymptotic configurations, namely, viscous–viscous (case 1) (Thaokar & Kumaran Reference Thaokar and Kumaran2005), viscoelastic–viscoelastic (case 2) and viscous–porous (case 3) bilayers, respectively, and (di) schematically shows the configuration of the present study (case 4). (a–c) Show LWLSA-Numerical and LWLSA-Analytical results for cases 1–3, respectively. Cases 1–3 are emulated as follows. The present system corresponds to case 1 when
$Da\to \infty$
(
$Da=\mathrm{1}0^{\mathrm{4}}$
),
${\textit{Ra}}\to \infty$
(
${\mu }_l,{\mu }_e\to 0$
),
$\chi =0$
(
$\xi =0$
); case 2 when
$Da\to \infty$
(
$Da=\mathrm{1}0^{\mathrm{4}}$
),
${\mu }_l=\mathrm{25}$
KPa,
${\mu }_e=\mathrm{30}$
KPa,
$\chi =0$
(
$\xi =0$
); and case 3 when
$Da=\mathrm{0.1}$
,
${\mu }_l,{\mu }_e\to 0$
(
${\textit{Ra}}\to \infty$
),
$\chi =\mathrm{1}0^{\mathrm{-4}}$
(
$\xi \neq 0$
). (d) Shows the LWLSA-Numerical and GLSA results for case 4. The other necessary dimensional parameters required for the simulations are provided in table 1.
Figure 2 compares and contrasts the variations in the linear growth coefficient (
$ {\varOmega }$
) with wavenumber (
$K$
) obtained from the LWLSA-Numerical, LWLSA-Analytical and GLSA for cases 1–4. The steps followed to carry out LWLSA-Analytical, LWLSA-Numerical and GLSA for all the asymptotic cases are provided at the ends of §§ 4.3 and 5. The solid curves with darker (lighter) shades in panels (a)–(d) correspond to the LWLSA-Numerical (GLSA) results, while the symbols represent analytical dispersion relations in the long-wave limit (LWLSA-Analytical). In panel (a), the LWLSA-Numerical results match the analytically obtained LWLSA-Analytical dispersion relation by Thaokar & Kumaran (Reference Thaokar and Kumaran2005), emulating case 1. The dispersion relation curve (
$ {\varOmega }$
vs
$K$
) exhibits unconditionally unstable behaviour, analogous to that of single-layer viscous films undergoing EHD instabilities, primarily due to the presence of a single deformable interface (Pease & Russel Reference Pease and Russel2002). For case 2, the anticipated finite-wavenumber characteristics of confined viscoelastic bilayers are obtained by executing both LWLSA-Numerical and LWLSA-Analytical, as shown in panel (b). These results in panels (a) and (b) confirm that the proposed LWLSA-Numerical numerical scheme accurately predicts the confined bilayer characteristics in the long-wave limit, similar to the analytical results (LWLSA-Analytical).
Dimensional values used in the simulations, if not stated otherwise. See ESI for detailed justification of the chosen parameter ranges.
Panel (c) shows the agreement between the LWLSA-Analytical and LWLSA-Numerical results for a porous medium-moderated viscous confined bilayer, emulating case 3, confirming the accuracy of the numerical scheme for predicting dispersion curves even in the presence of porous media. The results from panels (a)–(c) establish the capability of the numerical scheme to quantitatively predict eigenvalues for the current scenario (case 4), as shown in panel (d). Here, the LWLSA-Numerical and GLSA results agree well in the long-wave limit but deviate significantly at shorter wavelengths. This agreement between LWLSA-Numerical and GLSA once again confirms the accuracy of the numerical implementation. However, the deviation at shorter wavelengths highlights the necessity of using GLSA over LWLSA for more accurate results in these cases. Previous studies also suggest that electric-field-induced instabilities in viscoelastic films are dominated by short-wave behaviour at higher field intensities and film stiffness (Sarkar et al. Reference Sarkar, Sharma and Shenoy2008; Bandyopadhyay et al. Reference Bandyopadhyay, Reddy and Sharma2012). Additionally, the presence of porous media and elasto-capillarity induces short-wave hydrodynamic instabilities (Thiele, Goyeau & Velarde Reference Thiele, Goyeau and Velarde2009; Goyal et al. Reference Goyal, Kumar, Bandyopadhyay, Usha and Banerjee2013; Shargatov et al. Reference Shargatov, Tsypkin, Gorkunov, Kozhurina and Bogdanova2022). Hence, GLSA implementation is indispensable for accurately predicting the eigenvalues in the present work. Moving forward, all results presented have been obtained through GLSA computations. The range of values used in the simulations was selected to reflect realistic physical conditions and to maintain generality across a wide parameter space (Pelrine et al. Reference Pelrine, Kornbluh, Pei and Joseph2000; Mora et al. Reference Mora, Phou, Fromental, Pismen and Pomeau2010; Xu et al. Reference Xu, Li and Zhang2013). Permeability space used in porous–viscoelastic models of soft materials such as composite biological membrane gels (Mattern, Nakornchai & Deen Reference Mattern, Nakornchai and Deen2008) and hydrated collagen gel scaffolds (Serpooshan et al. Reference Serpooshan, Quinn, Muja and Nazhat2013) often spans
$10^{-13}-10^{-15}$
${\mathrm{m}}^{\mathrm{2}}$
along with elastic moduli in the range relevant (1–1000 Pa) (Cacheux et al. Reference Cacheux, Ordonez-Miranda, Bancaud, Jalabert, Alcaide, Nomura and Matsunaga2023).
6.2. Porous media modulated EHD instabilities
Figure 3 illustrates the influence of the applied electrostatic force
$(\varepsilon _0 \psi _0^2 / h_0^2)$
and the elastic stiffness of the solid layer
$(\mu _l / h_0)$
– represented through variations in
${\textit{Ra}}$
– on the deformable SP interface. The characteristic of the porous medium is held constant at
$Da = 0.1$
throughout. Panel (a) presents the dispersion relation, showing the variation of
$ {\varOmega }$
with
$K$
for different values of
${\textit{Ra}}$
, corresponding to varying electric-field potentials
$\psi _0$
. As
$\psi _0$
increases, the dispersion curves shift towards higher wavenumbers, indicating that the instability modes progressively favour shorter wavelengths. This trend is consistent with the EHD instabilities in thin confined bilayers previously reported by Mondal & Bandyopadhyay (Reference Mondal and Bandyopadhyay2014). Panel (b) further quantifies this behaviour, showing that
$ {\varOmega }_{{M}}$
increases and
${\varLambda }_{{M}}$
decreases with increasing
${\textit{Ra}}$
, in agreement with the trends observed in panel (a). Panel (c) examines the effect of increasing the elastic modulus
$\mu _l$
of the solid layer. The corresponding curves exhibit a shift toward finite-wavenumber instability, indicating the emergence of elasticity-driven mode selection. In panel (d),
$ {\varOmega }_{{M}}$
rises sharply for stiffer layers (lower
${\textit{Ra}}$
), but the growth rate saturates for softer layers once
${\textit{Ra}} \gtrsim 10^4$
. Similar behaviour is observed when the porosity of the medium is varied by adjusting the elasticity ratio between the porous and solid layers. Interestingly, panel (d) also shows that
${\varLambda }_{{M}}$
remains nearly invariant with changes in
$\mu _l$
. However, the onset conditions – i.e. the critical strength of the destabilising stresses required to initiate instability – are significantly influenced by
$\mu _l$
, thereby altering the critical wavelength
${\varLambda }_{{C}}$
. This indicates that, although elasticity has minimal influence on the dominant mode’s spatial scale, it significantly alters the critical threshold for the onset of instability.
The GLSA results to unveil the effect of electric Rayleigh number (
${\textit{Ra}}$
). Panels (a) and (c) show the variation of
$ {\varOmega }$
with
$K$
for different
${\textit{Ra}}$
by varying
$\psi _0$
and
$\mu _l$
, respectively, while other parameters are kept constant. The curves 1–3 in panel (a) (c) represent
$\psi _0 =$
70, 100 and 130 V (
$\mu _l = 4.0\times 10^3$
,
$2.4\times 10^4$
and
$4.7\times 10^4$
Pa), respectively. Panels (b) and (d) show the variation of the dominant growth coefficient,
$ {\varOmega }_{ {M}}$
(solid colour lines), and length scale,
${\varLambda }_{{M}}$
(blue dashed line and symbols), with
${\textit{Ra}}$
corresponding to panels (a) and (c), respectively. In panel (d), three different curves correspond to
$\mu _r =$
0.3, 0.5 and 0.8. The other essential parameters used for these simulations are shown in table 1.
Results obtained from GLSA to show the effect of
$Da$
. Panels (a) and (c) show the variation of
$ {\varOmega }$
with
$K$
for different
$Da$
by varying
$\kappa$
(
$h_0(=100$
nm) is kept constant) and
$h_0$
(
$\kappa (=1.00$
$\mathrm{\times }$
10
${}^{-15}{}^{\ }{\mathrm{m}}^{\mathrm{2}}$
is kept constant), respectively. The curves 1–5 in panel (a)
$\kappa =\mathrm{1.00}\times \mathrm{1}0^{-\mathrm{12}}\,{\mathrm{m}}^{\mathrm{2}}$
,
$\mathrm{1.38}\times \mathrm{1}0^{-\mathrm{14}}\,{\mathrm{m}}^{\mathrm{2}}$
,
$\mathrm{3.91}\times \mathrm{1}0^{-\mathrm{16}}\,{\mathrm{m}}^{\mathrm{2}}$
,
$\mathrm{1.58}\times \mathrm{1}0^{-\mathrm{17}}\,{\mathrm{m}}^{\mathrm{2}}$
and
$\mathrm{7.53}\times \mathrm{1}0^{-\mathrm{18}}\,{\mathrm{m}}^{\mathrm{2}}$
, respectively. Whereas curves 1–3 in panel (c) represent
$h_0=\mathrm{115}$
,
$\mathrm{125}$
and
$\mathrm{135}$
nm, respectively. Panels (b) and (d) show the variation of
${{\varOmega }}_{{M}}$
(black solid line) and
${{\varLambda }}_{{M}}$
(blue dashed line) with
$Da$
corresponding to panels (a) and (c), respectively. The other essential parameters used for these simulations are shown in table 1.
Figure 4 illustrates the effect of the porous medium, characterised by the Darcy number
$ Da = \kappa / h_0^2$
, on the interfacial instability of the SP system. The applied electrostatic forcing is held constant at
$ Ra = 1.77 \times 10^4$
. Panels (a) and (b) present the variation in instability behaviour with permeability, while panels (c) and (d) explore the impact of varying the porous medium thickness. Panel (b) reveals a sigmoidal response in both
${\varLambda }_{{M}}$
and
$ {\varOmega }_{{M}}$
, corresponding to the dispersion curves in panel (a). This sigmoidal behaviour divides the system response into three distinct regions: (I) nearly impermeable, (II) intermediately permeable and (III) highly permeable media. In region III
$(Da \rightarrow \infty )$
, the porous medium mimics the behaviour of an additional viscoelastic layer, resulting in
$ {\varLambda }_{{M}} \sim 2.5$
, with a marginal increase in
$ {\varOmega }_{{M}}$
, as seen in curves 1 and 2 of panel (a). As permeability decreases (region II), both
$ \varLambda _{{M}}$
and
$ {\varOmega }_{{M}}$
gradually decline, demonstrating transitionary behaviour (curves 2–4). Further reduction in permeability leads to stabilisation of the interface (region I), as observed in curves 4–5 of panel (a). Panel (b) confirms that introducing a porous substrate with moderate permeability
$(Da \approx 1 - 0.001)$
can reduce
$ {\varLambda }_{{M}}$
to as low as
$\sim \kern-1pt 1.6$
, thus enabling pattern miniaturisation via porous-medium-regulated EHD instabilities. This complements the earlier observation in figure 3 that elastic softening enables shorter-wavelength instabilities; here, the same is achieved through permeability modulation. Panels (c) and (d) show that increasing the thickness of the porous medium (reducing
$ Da$
) leads to a progressive decline in
$ {\varLambda }_{{M}}$
, while the growth rate
$ {\varOmega }_{{M}}$
increases with decreasing layer thickness – behaviour consistent with previous observations in saturated rigid porous media (Goyal et al. Reference Goyal, Kumar, Bandyopadhyay, Usha and Banerjee2013; Kumar et al. Reference Kumar, Goyal, Banerjee and Bandyopadhyay2013a
; Kumar et al. Reference Kumar, Usha, Banerjee and Bandyopadhyay2013b
). These findings collectively emphasise that the Darcy number
$ Da$
functions as a critical control parameter for controlling the length and time scales of interfacial instabilities in confined viscoelastic–porous systems.
Neutral stability curves from GLSA showing the variation of
$K$
with
${\textit{Ra}}$
by varying
$\psi _{0}$
(
$h_{0}(=100$
nm) and
$\mu _{l}(=500$
Pa) are kept constant),
$\mu _{l}$
(
$h_{0}(=100$
nm) and
$\psi _{0}(=100$
V) are kept constant) for
$Da = 10^{-1}$
in panels (a) and (b), respectively. Panels (d) and (e) show the variation of
$K$
with
$Da$
by varying
$\kappa$
(
$h_{0}(=100$
nm) is kept constant) and
$h_{0}$
(
$\kappa =1.00\times 10^{-15}$
$\textrm {m}^{\textrm {2}}$
) for
${\textit{Ra}} = 1.77 \times 10^{4}$
, respectively. The circle indicates the critical point
$(K_{{C}}, Ra_{{C}} )$
of the curves shown in (a) and (b) and the critical point
$(K_{{C}}, Da_{{C}} )$
of the curves shown in (d) and (e). Panels (c) and (f) show the variation of critical wavelength
${{\varLambda }}_{{C}}={\mathrm{2}\pi }/{K_{{C}}}$
with critical Rayleigh number (
${\textit{Ra}}_{{C}}$
) and critical Darcy number
$(Da_{{C}} )$
, respectively.
Figure 5 further elucidates the onset conditions for instability through neutral stability curves obtained using GLSA. The neutral curve, defined by
$ ( {\varOmega } = 0 )$
, delineates the boundary between stable
$ ( {\varOmega } \lt 0 )$
and unstable
$ ( {\varOmega } \gt 0 )$
regions. The shaded region under the curve corresponds to unstable configurations, where interfacial instabilities can emerge. Along this curve, the destabilising electrostatic stresses are exactly counterbalanced by the stabilising storage modulus and interfacial tension – whose relative contributions are modulated by the underlying porous substrate, as formalised in §§ 3 and 4. Panels (a), (b), (d) and (e) display the characteristic finite-wavenumber bifurcation associated with confined SP interfaces. The neutral wavenumbers
$K_{{C1}}$
(black/darker shade) and
$K_{{C2}}$
(brown/lighter shade) bound the range of unstable wavenumbers for which
$ {\varOmega } \gt 0$
, indicating a finite window of unstable zone. Wavenumbers
$K \lt K_{{C1}}$
are suppressed by elastic stresses, whereas those with
$K \gt K_{ {C2}}$
are dampened by surface tension. As the destabilising electrostatic force increases, the area enclosed by the neutral stability curve broadens, indicating enhanced instability. The intersection points
$ (K_{{C}}, Ra_{{C}} )$
and
$ (K_{{C}}, Da_{{C}} )$
denote critical thresholds beyond which instabilities emerge (Balram & Dinesh Reference Balram and Dinesh2025).
Importantly, the results show that the system remains entirely stable for
${\textit{Ra}} \lt Ra_{{C}}$
or
$Da \lt Da_{{C}}$
, identifying these as bifurcation parameters. Panels (a, b) indicate that increasing the viscoelastic stiffness of the layer raises
${\textit{Ra}}_{{C}}$
, requiring stronger electric fields for destabilisation (Bandyopadhyay et al. Reference Bandyopadhyay, Reddy and Sharma2012). Likewise, panels (d, e) emphasise that both permeability (
$\kappa$
) and porous-layer thickness (
$h_0$
) influence the onset of instability through their effect on
$Da_{{C}}$
. These observations are reinforced by panels (c) and (f). In panel (c), the critical wavelength (
${{\varLambda }}_{{C}}$
) decreases sharply with increasing
${\textit{Ra}}_C$
because the growing influence of
$\psi _0$
which enhances electrostatic destabilisation, promoting shorter-wavelength instabilities. Beyond a threshold, however, the increase in
$h_0$
induces a stabilising effect, leading to a subsequent rise in
${{\varLambda }}_{{C}}$
. This non-monotonic variation reflects a competition between the destabilising EHD stresses and the moderating influence of the porous-layer confinement. Panel (f) recovers the sigmoidal-like dependence of
${{\varLambda }}_{{C}}$
on
$Da_{ C}$
, as observed earlier in figure 4(b). Collectively, these findings reveal that, while
${{\varLambda }}_{{M}}$
determines the morphology of self-organised patterns under strongly unstable conditions,
${{\varLambda }}_{{C}}$
sets the threshold for the emergence of interfacial wettability – both length scales being critically governed by porous medium characteristics and elastic confinement. Additionally, these results can be extended for
$h_{0} =10\,\mu \mathrm{m}$
(added to the supplementary information is available at https://doi.org/10.1017/jfm.2026.11475).
The GLSA results showing the effect of the filling ratio
$\upsilon = h_0/d$
. Panels (a) and (b) show the variation of
$ {\varOmega }_{{M}}$
and
${\varLambda }_{{M}}$
with
$\upsilon$
(by varying
$d$
) for a constant
$Da=10^{-2}$
. Panels (c)–(h) show the contours of
${\varLambda }_{{M}}$
in a
${\textit{Ra}}-Da$
plane for
$\upsilon =0.4$
(c–e) and
$\upsilon =0.6$
(f–h). The other essential parameters used for these simulations are shown in table 1.
Figures 3–5 collectively underscore the efficacy of a porous medium as a modular element for controlling both the critical and dominant length scales of deformable SP interfaces. The ability of porous substrates to modulate EHD instability characteristics offers a compelling alternative to traditional multi-layer configurations, enabling substantial pattern miniaturisation in confined micro- and nanoscale systems. Building on the mechanistic insights from these figures, we now focus on strategies to push the dominant pattern length scale below unity (
${{\varLambda }}_{{M}} \lt 1.0$
), thereby extending the design space for soft-material interfaces driven by porous-medium-modulated instabilities.
6.3. Miniaturisation pathways – influence of geometric confinement
Figure 6 illustrates the influence of the filling ratio of the porous medium
$ (\upsilon = h_0/d )$
on interfacial instabilities across system sizes ranging from
$h_0 = 10$
nm to
$1\,\mu$
m. Panel (a) shows that increasing both system size and filling ratio enhances the instability of the SP interface, particularly when the dielectric permittivity of the porous medium exceeds that of the solid layer. This indicates that the porous layer can amplify EHD instabilities under such a dielectric contrast. Panel (b) demonstrates a clear transition from longer to shorter dominant wavelengths (
${\varLambda }_{{M}}$
) as both
$h_0$
and
$\upsilon$
increase. This transition is further explored in panels (c)–(h), which map
${\varLambda }_{{M}}$
over the
${\textit{Ra}}$
–
$Da$
parameter space. For smaller systems (e.g.
$h_0 = 10$
nm), shown in panels (c) and (f), only long-wave EHD instabilities are present, with negligible influence from the porous medium – indicating that the electric field dominates the instability mechanism in this regime. As the system size increases to
$h_0 = 500$
nm, panels (d) and (g) show that the porous layer begins to moderate the instability, leading to a sigmoid-like response and an apparent reduction in
${\varLambda }_{{M}}$
. This effect becomes more pronounced at
$h_0 = 1000$
nm, as illustrated in panels (e) and (h), where a strong sigmoid trend emerges and
${\varLambda }_{{M}}$
decreases substantially. Furthermore, increasing the filling ratio from
$\upsilon = 0.4$
to
$\upsilon = 0.6$
results in a notable reduction of
${\varLambda }_{{M}}$
to approximately 0.2. These results confirm that the porous medium plays a critical role in facilitating short-wave instabilities, thereby offering a promising pathway for achieving extreme pattern miniaturisation.
Three-dimensional contours of
${\varLambda }_{{M}}$
(colour bar) in
${\textit{Ra}}$
–
$Da$
–
$\upsilon$
space. Panels (b)–(e) show the contours of
${\varLambda }_{{M}}$
in a
${\textit{Ra}}-Da$
plane for
$\upsilon =0.34$
,
$\upsilon =0.43$
,
$\upsilon =0.53$
and
$\upsilon =0.67$
, respectively, for
$h_0 = 100$
nm (varying
$d$
). The other essential parameters used for the presented simulation results are shown in table 1.
Figure 7 succinctly synthesises the insights from figure 6 by mapping the combined influence of
${\textit{Ra}}$
,
$Da$
and
$\upsilon$
on
${\varLambda }_{{M}}$
. Panel (a) presents a 3-D contour volume that vividly illustrates the multivariate dependency of
${\varLambda }_{{M}}$
on the electrostatic forces, permeability and confinement geometry. These high-dimensional visualisations offer a powerful framework for interpreting the coupled effects of system parameters and for identifying regimes conducive to desired pattern scales. The contours confirm that the dominant wavelength is tuneable over a broad range – from
${\varLambda }_{{M}} \sim 0.2$
to
${\varLambda }_{{M}} \sim 9$
– by systematically navigating the
${\textit{Ra}}$
–
$Da$
–
$\upsilon$
space. The signature sigmoidal transition is observed across intermediate values of
$\upsilon$
, with long-wave instabilities emerging at lower
$\upsilon$
and short-wave instabilities dominating at higher
$\upsilon$
, consistent with the underlying modulation by the porous layer. To elucidate these trends further, panels (b)–(e) provide 2-D contour slices corresponding to fixed values of
$\upsilon$
spanning 0.34–0.67. Panel (b), representative of lower filling ratios, shows that the system favours long-wave modes due to minimal influence from the porous layer. As
$\upsilon$
increases, panels (c)–(e) reveal a marked enhancement in destabilising EHD stresses, driven by increasing porous-layer participation. This results in a progressive shift toward shorter wavelengths, confirming the porous medium’s role in fuelling short-wave instabilities through increased field coupling and stress localisation.
Figures 6 and 7 demonstrate that precise control over instability length scales – from long to short – is achievable via judicious tuning of porous medium characteristics, electrostatic forcing and confinement geometry. These findings reinforce the pivotal role of the porous substrate as an active control element, enabling design strategies for functional interfaces in soft-matter systems with application-tailored micro/nanopatterning in porous substrate.
The GLSA results to exhibit the influence of dielectric contrast
$\varepsilon _r=\varepsilon _l/\varepsilon _m$
. Results presented in panels (a–c) (d–f) are obtained by varying
$\varepsilon _l$
(
$\varepsilon _m$
) and keeping
$\varepsilon _m (= 4.5)$
(
$\varepsilon _l (=2.5)$
) constant. Panels (a) and (d) show the variation of
$ {\varOmega }$
with
$K$
for different
$\varepsilon _{r}$
. Curves 1–5 in panel (a) correspond to
$\varepsilon _r=0.02$
, 0.20, 0.71, 1.31 and 1.80, respectively. Curves 1–5 in the panel (b) correspond to
$\varepsilon _r=0.35$
, 0.80, 2.08, 5.00 and 50.00, respectively. Panels (b) and (e) show the variation of
$ {\varOmega }_{{M}}$
and
${\varLambda }_{{M}}$
with
$\varepsilon _r$
. Panels (c) and (f) present the neutral stability curves showing the variation of
$K$
with
$\varepsilon _{r}$
. The other essential parameters used for these simulations are shown in table 1.
6.4. Miniaturisation pathways – influence of dielectric contrast
Figure 8 presents the influence of
$\varepsilon _r (= \varepsilon _l/\varepsilon _m )$
, on the length and time scales of instability. Before analysing the plotted results, we examine the expression for the electrostatic stress
$ ( \phi _{\textit{EL}} )$
acting on the deformable SP interface, which provides critical insight into the effect of
$\varepsilon _r$
. Substituting the base-state electric-field expressions from (4.5) and (4.6) into the normal stress balance (4.16) yields
\begin{equation} \phi _{\textit{EL}} = -\frac {\varepsilon _m \varepsilon _r (\varepsilon _r - 1)}{2 (D + (\varepsilon _r - 1))^{2}}. \end{equation}
As evident from (6.1), the electrostatic stress vanishes when the
$\varepsilon _r$
approaches unity
$ ( \varepsilon _l \to \varepsilon _m )$
, leading to
$\phi _{\textit{EL}} \to 0$
and thereby stabilising the SP interface
$ ( {\varOmega }_{{M}} \lt 0 )$
(Kumar & Matar Reference Kumar and Matar2004; Bandyopadhyay et al. Reference Bandyopadhyay, Reddy and Sharma2012). In contrast, a non-unity
$ \varepsilon _r (\neq 1 )$
induces interfacial instabilities
$ ( {\varOmega }_{{M}} \gt 0 )$
by enhancing the interfacial Maxwell stress. That signifies that increasing the magnitude of the dielectric contrast,
$ |\Delta \varepsilon | = |\varepsilon _l-\varepsilon _m|$
achieved by ensuring either
$\varepsilon _r \lt 1$
or
$\varepsilon _r \gt 1$
, amplifies the growth rate while reducing the instability wavelength, as observed in panels (a), (b), (d) and (e). Panels (a) and (d) depict dispersion curves for various values of
$\varepsilon _r$
, confirming that instabilities are promoted on both sides of
$\varepsilon _r = 1$
. Panels (b) and (e) further illustrate the variation of
$ {\varOmega }_{{M}}$
and
${\varLambda }_{{M}}$
with
$\varepsilon _r$
, revealing a singularity at
$\varepsilon _r = 1$
, where both the growth rate and dominant wavelength exhibit extremal behaviour. The trends in
$ {\varOmega }_{{M}}$
(increasing) and
${\varLambda }_{{M}}$
(decreasing) closely follow the variation in
$\phi _{\textit{EL}}$
, with extrema in
$\phi _{\textit{EL}}$
analytically verified at
$\varepsilon _r = 0.3$
and
$\varepsilon _r = 3.0$
via
${\rm d}\phi _{\textit{EL}}/{\rm d}\varepsilon _r$
. However, the peak in
$ {\varOmega }_{{M}}$
and the corresponding trough in
${\varLambda }_{{M}}$
are shifted relative to these extrema, due to the modulation of EHD instabilities by the underlying porous medium – an effect captured through the stress balances ((4.16) and (4.17)). Curves 2 and 4 in panels (a) and (d), corresponding to
$\varepsilon _r = 0.2$
and
$\varepsilon _r = 5$
, respectively, exemplify this deviation. Finally, panels (c) and (f) present the neutral stability boundaries as a function of
$\varepsilon _r$
, confirming trends similar to those in panels (b) and (e). These results underscore the critical role of dielectric contrast
$ ( |\Delta \varepsilon | )$
in governing both the onset and morphology of EHD instabilities at the deformable SP interface.
Three-dimensional contours of
${{\varLambda }}_{{M}}$
(colour bar) in
${\textit{Ra}}$
–
$Da$
–
$\varepsilon _r$
space. Panels (a)–(d) show the contours of
${\varLambda }_{{M}}$
in a
${\textit{Ra}}$
–
$Da$
plane for different
$\varepsilon _r$
while
$\varepsilon _m (=4.5)$
is kept constant. The other essential parameters used in the simulations are shown in table 1.
Figure 9 further elucidates how
$\varepsilon _r$
modulates the dominant instability length scale
$ ( {\varLambda }_{{M}} )$
in conjunction with
${\textit{Ra}}$
and
$Da$
. As shown in panel (a), near-unity permittivity ratio (
$\varepsilon _r \to 1$
) yields long-wave instabilities with
${\varLambda }_{{M}} \sim 9$
, corresponding to vanishing electrostatic stresses and a near-stable dynamics. In contrast, panels (b)–(e) demonstrate that a significant departure from
$\varepsilon _r = 1$
, i.e. the increase in
$|\Delta \varepsilon |$
– drastically shortens
${\varLambda }_{{M}}$
, leading to tighter, more localised instability modes (
${\varLambda }_{{M}} \sim 3$
). When
$|\Delta \varepsilon |$
becomes extreme, i.e.
$\varepsilon _r \gg 1$
or
$\varepsilon _r \ll 1$
, miniaturised patterns with
${\varLambda }_{{M}} \lt 1$
emerge, showcasing the system’s sensitivity to interfacial electrostatics. This whole range of variation in
${\varLambda }_{{M}}$
is captured in panel (a). Such dielectric tuning not only reduces the critical electric potential (
$\psi _0$
) required to trigger instability (Bandyopadhyay et al. Reference Bandyopadhyay, Sharma, Thiele and Reddy2009), but also allows morphological control – columns versus holes – based on permittivity asymmetry across the interface (Mondal et al. Reference Mondal, Kumar and Bandyopadhyay2013).
Figures 8 and 9 underscore the pivotal role of dielectric contrast
$ ( |\Delta \varepsilon | )$
in shaping both the stability threshold and dominant length scales of EHD instabilities in confined porous–viscoelastic systems. While figure 8 establishes the theoretical basis through electrostatic stress analysis, figure 9 demonstrates the practical utility of tuning
$\varepsilon _r$
to access a wide range of miniaturised pattern morphologies. The coupling of dielectric mismatch with porous medium and elasticity parameters thus offers a powerful design strategy for engineering controllable, energy-efficient encapsulated porous patterning in soft micro/nano-scale systems.
7. Conclusions
This study presents a comprehensive GLSA to investigate EHD instabilities at a deformable interface between a viscoelastic solid and a solid-saturated porous medium confined under an externally applied electric field. By coupling the mDBKV framework for the solid-saturated porous substrate with Maxwell stress formulations for the electric field, we derive a complete set of governing equations and interfacial boundary conditions – including stress-jump constraints – to describe the dynamics of the deformable interface.
To the best of our knowledge, this is the first study to formulate and apply a GLSA framework that systematically captures the interplay between SP interfacial deformation, electrostatic forcing and viscoelastic confinement, uncovering controllable instability pathways in SP bilayer systems. Key conclusions emerging from this work include:
-
(i) The GLSA framework offers a full-spectrum characterisation of interface stability in resolving finite-wavenumber bifurcations. Non-dimensional parameters – Darcy number (
$Da$
), electric Rayleigh number (
${\textit{Ra}}$
), dielectric permittivity ratio (
$\varepsilon _r$
) and filling ratio (
$\upsilon$
) – are identified as key regulators of pattern length scales and growth kinetics. -
(ii) External electrostatic force destabilises the SP interface, leading to finite-wavenumber bifurcations in the presence of solid layers with higher elastic modulus. While the dominant wavelength
$ ({\varLambda }_{{M}} )$
decreases with increasing
${\textit{Ra}}$
, the elastic modulus mainly influences the critical onset of instability rather than the dominant wavelength itself. This critical mode governs interfacial wettability and thus impacts the performance and reliability of solid-state electrochemical devices. -
(iii) Neutral stability curves reveal critical thresholds (
${\textit{Ra}}_C$
,
$Da_C$
) separating stable and unstable regimes. These bifurcations are sensitive to permeability and thickness of the porous medium (
$\kappa$
,
$h_0$
), applied potential (
$\psi _0$
) and dielectric asymmetry (
$\varepsilon _r$
), enabling predictive control of instability onset. -
(iv) The porous medium acts as an effective modulator of instability, promoting miniaturisation through its higher dielectric permittivity (
$\varepsilon _m$
). A sigmoidal dependence of both dominant wavelength
$ ({\varLambda }_{{M}} )$
and time scale
$ ( {\varOmega }_{{M}} )$
on
$Da$
delineates three functional regimes – impermeable, transitional and highly permeable – mirroring a crossover from solid-like to fluid-like substrate behaviour. -
(v) The filling ratio (
$\upsilon$
) and dielectric permittivity ratio (
$\varepsilon _r$
) emerge as key tuneable parameters for tailoring the instability dynamics. Increased
$\upsilon$
intensifies EHD stresses at the interface, while high dielectric contrast
$ ( |\Delta \varepsilon | )$
(either
$\varepsilon _r \gt 1$
or
$\varepsilon _r \lt 1$
) accelerates growth rates and compresses instability wavelengths. In contrast,
$\varepsilon _r = 1$
leads to complete suppression of instability due to vanishing electrostatic stresses. -
(vi) The system size plays a pivotal role in accessing different instability regimes. As the system transits from nano- to micro-scale, the interfacial instabilities shift from long-wave to short-wave and eventually
${\varLambda }_{{M}} \lt 1$
is achieved, enabling high-aspect-ratio pattern miniaturisation across large-area surfaces. -
(vii) Multivariate 3-D contour analyses across
${\textit{Ra}}$
–
$Da$
–
$\upsilon$
–
$\varepsilon _r$
reveal that a wide spectrum of interfacial patterns (
${\varLambda }_{{M}} \sim 0.2-10$
) can be achieved through synergistic parameter tuning. This establishes a quantitative design framework for engineering self-organised, encapsulated porous patterns with embedded morphological confinement, offering a robust handle for hierarchical structuring within the bilayer system. A subsequent future work will focus on extending the present linear framework toward fully nonlinear simulations to visualise the morphological evolution of perturbation modes and validate the predicted instability pathways.
These results yield actionable design principles:
-
(i) amplify instability by introducing dielectric contrast
$ ( |\Delta \varepsilon | )$
across the interface; -
(ii) exploit substrate permeability (
$Da$
) to access short-wave morphologies; -
(iii) tailor geometric filling ratio (
$\upsilon$
) to intensify EHD destabilisation; -
(iv) tune
$\upsilon$
and
$\varepsilon _r$
in the
${\textit{Ra}}$
–
$Da$
space to transition across longer (
${\varLambda }_{{M}} \sim 10$
) to ultra-short (
${\varLambda }_{{M}} \sim 0.2$
) wavelength regimes; -
(v) use GLSA-predicted thresholds to determine safe and unstable operating windows.
Together, these findings provide the first theoretical foundation for porous medium-modulated EHD pattern formation in confined viscoelastic bilayers. Rather than relying on complex multilayer architectures, this work shows how electrostatic forcing, saturated porous-layer compliance and confinement geometry can be harnessed to design functional soft interfaces.
Outlook. This work marks a significant step toward the rational design of electrohydrodynamically active soft-matter platforms. In particular, the mDBKV framework provides a tractable way to capture the wettability dynamics in real soft porous systems such as hydrated gels and porous solid electrolytes, where EHD instabilities strongly impact interfacial stability and performance. Thus, the present stability analysis directly addresses parameter regimes of practical interest in both biological and electrochemical soft-matter systems.
The insights and design strategies proposed here offer a foundation for advancing next-generation solid-state energy technologies, nanostructured EHD devices and influencing biological systems through informed control of interfacial instabilities.
Future work will extend the present analysis to finite base-state deformation using a Darcy–Brinkman–Biot framework (Carrillo & Bourg Reference Carrillo and Bourg2019) in order to investigate nonlinear bulk–interface coupling in soft porous systems.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.11475.
Acknowledgements
The insightful discussions with Dr V.B. Vanarse are gratefully acknowledged. The authors thank the Computer Centre at IIT (BHU) Varanasi for providing access to the ParamShivay Supercomputing Facility, supported under the National Supercomputing Mission (NSM), Government of India.
Funding
A.G. acknowledges financial support from the Anusandhan National Research Foundation (ANRF), formerly known as the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India, through the Start-up Research Grant (SRG/2022/000089).
Declaration of interests
The authors declare no conflicts of interest.
Data availability statement
Data are available within the article or its supplementary materials.
Appendix A. Analytical solutions to asymptotic cases
The general analytical solutions to (4.28) and (4.30) corresponding to the three asymptotic cases (cases 1–3) are as follows:
\begin{align} {\widetilde {W}}_l&=Q_{\mathrm{1}}+Q_{\mathrm{2}}Z+Q_{\mathrm{3}}Z^{\mathrm{2}}+Q_{\mathrm{4}}Z^{\mathrm{3}} ,\\[-9pt] \nonumber \end{align}
\begin{align} {\widetilde {{\varPsi }}}_l &= Q_{\mathrm{5}}Z+Q_{\mathrm{6}}, \\[-9pt] \nonumber \end{align}
\begin{align} {\widetilde {{\varPsi }}}_m &= Q_{\mathrm{7}}Z+Q_{\mathrm{8}}. \end{align}
However, the general analytical solutions to (4.29) for cases 1–3 are as follows:
\begin{align} \text{Case 1:} \; \; {\widetilde {W}}_m &= Q_{\mathrm{9}}+Q_{\mathrm{10}}Z+Q_{\mathrm{11}}Z^{\mathrm{2}}+Q_{\mathrm{12}}Z^{\mathrm{3}}, \\[-9pt] \nonumber \end{align}
\begin{align} \text{Case 2:} \; \; {\widetilde {W}}_m &= Q_{\mathrm{9}}+Q_{\mathrm{10}}Z+Q_{\mathrm{11}}Z^{\mathrm{2}}+Q_{\mathrm{12}}Z^{\mathrm{3}}, \\[-9pt] \nonumber \end{align}
\begin{align} \text{Case 3:} \; \; {\widetilde {W}}_m &= {\left (\eta _r/Da\right )}^{\mathrm{-1}} \big (Q_{\mathrm{9}}e^{\beta Z}+Q_{\mathrm{10}}e^{-\beta Z}\big ) + Q_{\mathrm{11}}Z + Q_{\mathrm{12}}. \\[10pt] \nonumber \end{align}
Appendix B. Derivation of kinematic condition (2.11)
Any interface can be represented as
\begin{equation} F\left (x,z,t\right )=z-h\left (x,t\right ) .\end{equation}
In the case of a deformable interface, a material point on the interface remains on the interface for all time, which means its material derivative
$(D/Dt)$
is zero, i.e.
\begin{align} \frac {\textit{DF}}{Dt} &=0 , \\[-9pt] \nonumber \end{align}
\begin{align} \frac {\textit{DF}}{Dt} &=\frac {\partial F}{\partial t} +\frac {\partial u}{\partial t} \frac {\partial F}{\partial x} +\frac {\partial w}{\partial t} \frac {\partial F}{\partial z} =0 . \end{align}
Here,
$({\partial u}/{\partial t})$
and
$({\partial w}/{\partial t})$
are the velocity components as time derivatives of the displacements in x and z-directions, respectively.
Now, substituting (B1) in (B3), we get
\begin{equation} \frac {\partial h}{\partial t} +\frac {\partial u}{\partial t} \left (\frac {\partial h}{\partial x} \right )=\frac {\partial w}{\partial t} .\end{equation}
This can be written for any solid layer (l) as
\begin{equation} \dot {h}+\dot {u}_{l} \boldsymbol{\cdot }\boldsymbol{\nabla} _{s} h=\dot {w}_{l} .\end{equation}
