Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Branching dynamics in electrohydrodynamic instabilities of viscoelastic soft gels

Gyandeep Balram and Bhagavatula Dinesh*
Department of Chemical Engineering and Technology, Indian Institute of Technology-BHU, Varanasi, UP 221005, India. E-mail: dinesh.che@iitbhu.ac.in

Received 16th July 2025 , Accepted 6th September 2025

First published on 24th September 2025


Abstract

An electric field imposed on a bilayer of fluids that are stably stratified in the presence of gravity leads to an instability manifested by interfacial deflections. The layer of perfect conductor is simulated using a linear viscoelastic model and the perfect dielectric is considered to be a layer of air. Under neutral conditions, the key dimensionless groups are the dimensionless electric potential, Bond number and the Weissenberg number. The branching behavior upon instability to sinusoidal disturbances is determined by weak nonlinear analysis with the dimensionless potential advanced from its critical value at neutral stability. An analytical expression obtained from weak nonlinear analysis leads to the unintuitive result that sinusoidal deflections can either lead to supercritical saturated waves or lead to subcritical breakup depending on the elasticity of the perfect conductor. The analytical expression also indicates that there is a transition wave number below which supercritical saturation ought to occur, it can be shown that such wave numbers can be geometrically accessed, thus permitting any supercritical saturation to steady waves. In contrast, our results demonstrate that when the perfect conductor is modeled as an Oldroyd-B fluid, the branching remains subcritical in nature, ultimately leading to interface rupture—mirroring the behavior observed in the Newtonian fluid case (as demonstrated by B. Dinesh and R. Narayanan, Phys. Rev. Fluids, 2021, 6, 054001).


1 Introduction and physics

Electrohydrodynamics concerns the interaction between fluids and electric fields. The foundational studies by Taylor and Melcher2,3 established the theoretical basis for analyzing such phenomena. Subsequent research has examined how electric fields influence the morphology and stability of interfaces in systems such as liquid bridges, drops, jets, and thin films. When an electric field is applied perpendicular to an initially flat interface, it typically induces destabilization.

Vertical electric fields have been applied in quiescent multilayer flows for the controlled self-assembly of microscale hierarchical structures in polymer melts and processes like soft lithography.4–8 Other applications include inkjet printing,9 enhancement of mixing and heat transfer, as well as the development of soft devices with adjustable shapes.10,11 For an overview of electrohydrodynamic instability in thin fluid films, see the review by Papageorgiou.12

In this work, we investigate the branching behavior of electrohydrodynamic instability in a system comprising a perfect conductor–dielectric fluid pair. The perfect conductor is modeled in two ways: first, using a linear viscoelastic model, and second, using the Oldroyd-B fluid model. The dielectric fluid is assumed to be a hydrodynamically passive air layer. As shown in Fig. 1, the fluids are confined between two rigid plates, with the bottom plate held at a constant voltage D, while the top plate is grounded. The applied voltage counteracts gravity and drives an instability that leads to the formation of an array of pillars in the conducting fluid.1,12–15 However, when the elastic stresses in the conducting fluid overcome the destabilizing effect of the constant voltage, the interface undergoes saturated deformation without branching.


image file: d5sm00727e-f1.tif
Fig. 1 Schematic of electrostatic instability. The bottom viscoelastic fluid is a perfect conductor and the top fluid is a perfect dielectric. The bottom plate is maintained a constant voltage D and the top plate is grounded.

The instability arises due to competition between the applied potential and stabilizing factors such as gravity, elastic stresses, and surface tension. These competing effects can produce a minimum in the voltage–wave number curve at the onset of instability.1,16

At neutral stability, if the perfect conductor is modeled as a Newtonian fluid, the velocity perturbations vanish. For a viscoelastic perfect conductor, at low wave numbers, the destabilizing voltage is balanced by gravity and elastic stresses, while at high wave numbers, surface tension opposes the destabilizing effect.1,16,17 The high-wave number behavior is similar to that seen in the Rayleigh–Taylor instability, which is typically subcritical in nature. In contrast, the low-wave number balance resembles that of thermocapillary-driven problems, such as the Bénard–Marangoni instability, which generally exhibits supercritical behavior.18

The objective of this study is to determine whether a transition occurs between supercritical and subcritical behavior at a critical wave number, or if the instability consistently exhibits one of these two behaviors across all wave numbers.

We review previous studies relevant to this work. Wu and Chou19 showed that, for Maxwell fluids subjected to DC electric fields, elasticity increases the growth rate of pillar formation. They identified a critical Deborah number above which the growth rate for certain wavenumbers diverges.

Tomar et al.20 examined the surface instability of a confined viscoelastic liquid film under an applied electric field, considering both Maxwell and Jeffreys models. Their analysis showed that incorporating fluid inertia removes the singularity, resulting in large but finite growth rates for all Deborah numbers and identifying the dominant instability wavelength. They also observed that the small inertia limit does not coincide with the zero-inertia limit when describing the instability dynamics and wavelength in polymer melts.

Espín et al.7 studied the effects of viscoelasticity on instabilities under AC and DC electric fields using the Jeffreys model for both perfect and leaky dielectrics. In the DC case, asymptotic methods were applied to resolve a singularity occurring when solvent viscosity is neglected, corresponding to the Maxwell limit. Their results show that elasticity increases both the maximum growth rate and the wavenumber of the instability.

Dinesh et al.1 investigated the electrohydrodynamic instability between a perfect conductor and a dielectric, modeling both as Newtonian fluids. Their study found that the branching behavior is always subcritical, resulting in rupture of the interface. In contrast, the present study demonstrates that both subcritical and supercritical branching behavior can occur when the perfect conductor is described using a viscoelastic model.

Observations from previous studies indicate that pillar formation due to instability spans the gap between the plates, which is characteristic of subcritical interface behavior. However, elasticity in the viscoelastic fluid can also result in supercritical saturation of the interface. The objective of this work is to demonstrate this behavior through mathematical analysis and to identify the physical conditions that lead to either subcritical or supercritical instability.

The analysis presented here follows the methodology outlined in the works of Dinesh and Narayanan.1,21 We consider a system comprising a perfect conductor–dielectric fluid pair and derive an analytical expression to characterize the branching behavior. The perfect conductor is modeled using a linear viscoelastic framework. Additionally, we examine a case where the perfect conductor is represented by an Oldroyd-B fluid.

A perturbation expansion is carried out around the neutral stability state, with the applied voltage slightly exceeding its critical threshold. This approach allows us to investigate the nonlinear evolution of the interface and distinguish between subcritical and supercritical regimes.

2 Mathematical model

The model assumes a hydrodynamically active viscoelastic fluid with constant properties whose free surface is exposed to a passive gas in a stabilizing gravitational field, as depicted in Fig. 1. The gas layer is considered to be air and hydrodynamically inactive. The viscoelastic fluid and the air layers are between two rigid electrically conducting plates located at −h* and h, across which a constant voltage difference, D, is applied (cf. Fig. 1). The potential field in the top fluid is given by
 
2ψ = 0 (1)
The bottom viscoelastic fluid, represented by an asterisk, is assumed to be a perfect conductor while the top fluid is taken to be a perfect dielectric. The bottom fluid is taken to be linearly viscoelastic for algebraic simplicity whilst retaining essential physics in this study. Now, the stress tensor in a viscoelastic fluid may be expressed in terms of the displacement field. The displacement vector, R, is the displacement of the position vector, x, in the current configuration from the position vector, ζ, in the reference configuration. In other words,
 
x = ζ + R(x) (2)

For a linear viscoelastic fluid (cf., ref. 22–24) the stress tensor, T, is given by

 
T = −pI + G(∇R + ∇RT) + μg(∇v + ∇vT) (3)
Here G and μg are the shear modulus and viscosity of the fluid material and v is the velocity field. Now the velocity field in the fluid is itself expressed in terms of the displacement field (cf., ref. 25 for a detailed explanation). This expression is given by
 
image file: d5sm00727e-t1.tif(4)

In a linear viscoelastic soft-gel layer, the elastic field is described by the displacement field vector (cf. eqn (2)), which represents the deformation of material points, while the velocity field vector in eqn (4), describes the rate of change of these displacements over time. The relaxation mechanism occurs as the internal elastic stresses, generated by deformation, gradually dissipate through viscous dissipation in the soft-gel. This means the displacement field changes slowly, and the velocity field captures how the material returns toward its undeformed state by relaxing the elastic stresses. Over time, energy stored in elastic deformation is released or converted to viscous dissipation, resulting in a decrease of internal stresses and a gradual return to equilibrium. This relaxation process is characterized by time-dependent decay of stresses and strains governed by the viscoelastic properties of the gel.21,23,24,26–29 The equations of motion in the viscoelastic medium are thus

 
image file: d5sm00727e-t2.tif(5)
where T and v are given by eqn (3) and (4), where in eqn (4) the superscript T on ∇R is the transpose of the tensor ∇R and I is the identity tensor and where iz is the base vector in the positive z-direction in eqn (5). The horizontal and vertical components of the displacement vector, R, are denoted by X and Z.

The viscoelastic fluid is taken to be incompressible and the mass conservation demands (ref. 25 and 30)

 
det(F) = 1 (6)
where F is the deformation tensor given by image file: d5sm00727e-t3.tif. Here ζi represent the components of the position vector in the reference configuration. Upon expansion of eqn (6) and using eqn (2) we get
 
image file: d5sm00727e-t4.tif(7)
eqn (5) and (7) along with the representations for T and v constitute the domain equations. Other models, such as the neo-Hookean25,29,31 and Oldroyd-B models,32 can also describe the soft-gel layer. In this work, we use the linear viscoelastic model to obtain a mathematically tractable solution and focus on the core physics of the electrohydrodynamic instability. A stability analysis using the Oldroyd-B model is included in the Appendix D. The governing equations are complemented by boundary conditions at the rigid wall and at the soft-gel air interface.

At the wall, the displacement fields are taken to be zero. At the interface, z = h(x, t), the normal and the tangential components of the momentum balance hold. They are

 
n·T·t = 0  and [thin space (1/6-em)]n·T·nD(TM·n) = −γ∇·n (8)
where the unit normal vector (n) and the unit tangent vector (t) are given by
 
image file: d5sm00727e-t5.tif(9)
In addition we have an impermeable interface along with its kinematic relation
 
image file: d5sm00727e-t6.tif(10)
In eqn (8), the dimensional form of the Maxwell stress tensor, i.e., TM, (cf. ref. 33), is given by image file: d5sm00727e-t7.tif where E = −∇ψ is the electric field, ε is the relative permittivity of the fluid and ε0 is the permittivity of free space.

The governing equations are made dimensionless by using the following scales denoted by the subscript ‘c’:

 
image file: d5sm00727e-t8.tif(11)
Here U is a characteristic velocity scale. The scaled potential field in the top fluid is given by
 
2ψ = 0 (12)
The potential field is subject to a constant value of unity at the bottom plate, i.e., at z = −1, while the top plate at z = [script letter H] is assumed to be grounded, where [script letter H] = h/h*.

The dimensionless governing equations for the bottom fluid are given by

 
image file: d5sm00727e-t9.tif(13)
where image file: d5sm00727e-t10.tif, image file: d5sm00727e-t11.tif and image file: d5sm00727e-t12.tif and image file: d5sm00727e-t13.tif. The mass conservation eqn (7) remains unchanged and does not induce any dimensionless groups. The interfacial force balance conditions at z = h(x, t) become
 
image file: d5sm00727e-t14.tif(14)
The normal component of the Maxwell stress in dimensionless form can be written as
 
image file: d5sm00727e-t15.tif(15)
Here, [scr D, script letter D] = εε0D2/h*μ*U. For the perfect conductor–dielectric model, the tangential components of the Maxwell stress tensor vanish (cf. ref. 33). From the scaling of the equations and boundary conditions, five key dimensionless groups emerge: Re, Ca, Bo, [scr D, script letter D], and Wi, which are evaluated using representative physical parameters (see Tables 1 and 2). The base state R = v = T = 0 is subject to a stability analysis to explore the onset and nature of the bifurcation. We seek the neutral stability condition and characterize steady-state branching, governed primarily by Wi, Ca, Bo, and [scr D, script letter D]. Notably, under steady-state conditions, Ca appears only as Ca/Wi, allowing U to be chosen such that Ca = 1.

Table 1 Physical properties of polyacrylamide viscoelastic fluid34,35
Physical property Range
Density (ρ) 1000 kg m−3
Viscosity (μ) 1.5 kg m−1 s−1
Shear modulus (G) 10–45 Pa
Surface tension (γ) 73 × 10−3 N m−1
Thickness (h*) 1 × 10−3–20 × 10−3 m


Table 2 Range of key dimensionless quantities for polyacrylamide viscoelastic fluid. The velocity scale U is based on the capillary number such that image file: d5sm00727e-t16.tif
Dimensionless parameter Definition Range
Re ρh*U/μ 1.5
Wi μU/Gh* 0.1–10
Bo ρgh*2/γ 0.1–50


3 Linear stability analysis

A linear stability analysis is carried out by introducing small perturbations ψ1, R1, p1, and ζ1 about the quiescent base state. These are expressed using the expansion
ψ(x, z, t) = ψ0(z) + ψ1(z)cos(kx)eσt,
with similar forms assumed for R, p, and ζ. Here, k denotes the wavenumber of the perturbation and the subscript “0” refers to the base-state quantities. Neutral stability corresponds to σ = 0 in the above expansion.

Assuming constant fluid volume imposes the constraint

image file: d5sm00727e-t17.tif
is the disturbance wavelength.

Since the base-state displacement R0 = 0, the base-state pressure gradient is balanced by the gravitational body force:

image file: d5sm00727e-t18.tif
here [capital G, script] = ρgh*2/μU and δj = ρj/ρ, takes values ρdieltric/ρconductor and 1. The base-state potential field is given by ψ0 = 1 − z/[script letter H]. The perturbation fields, being non-zero, are determined by solving the linearized governing equations for the soft gel. Substituting the assumed forms of the disturbances into the governing equations yields:
 
image file: d5sm00727e-t19.tif(16)
 
image file: d5sm00727e-t20.tif(17)
 
image file: d5sm00727e-t21.tif(18)

Linearizing the boundary conditions at the reference interface, z = 0

 
Z1 = ζ1 (19)
 
image file: d5sm00727e-t22.tif(20)
The no displacement fields X1 = 0 and Z1 = 0, boundary conditions, at the wall i.e., z = −1 are imposed. Upon eliminating p1 and X1 from the governing equations and using the boundary conditions we get an expression for Z1, i.e.,
 
image file: d5sm00727e-t23.tif(21)
The linearized equation for the potential field in the top fluid are given by
 
image file: d5sm00727e-t24.tif(22)
where at z = 0 and z = [script letter H], we have image file: d5sm00727e-t25.tif. Upon solving eqn (22), we get
 
image file: d5sm00727e-t26.tif(23)
To determine the neutral stability criteria, we turn toward the perturbed normal component of the momentum balance along the interface, z = 0. First, we take the derivative of this equation with respect to x to eliminate pressure, p1 using the x-momentum equation. Second, using the continuity equation, we eliminate X1 in terms of Z1, i.e.,
 
image file: d5sm00727e-t27.tif(24)
Upon substituting the solution for ψ1 and Z1, in eqn (24), we get
 
image file: d5sm00727e-t28.tif(25)
 
image file: d5sm00727e-t29.tif(26)
Here, [capital G, script]Ca (also known as the Bond number, Bo), Wi/Ca and [scr D, script letter D]Ca, hereafter termed [scr D, script letter D]0Ca at the neutral point, can be grouped as pairs and are observed to be independent of the characteristic velocity scale, U, and the viscosities. In addition, here δ is the ratio of the density of the perfect dielectric to the density of the perfect conductor.
 
image file: d5sm00727e-t30.tif(27)
Observe that for low wavenumbers and shallow film thicknesses, term I is of image file: d5sm00727e-t31.tif and term II is of [scr O, script letter O](1). This is obtained by expanding cosh(2k), sinh(k), and cosh(k) in term I in series of the wavenumber k and neglecting higher order terms of wavenumbers. For example, coth(k[script letter H]) when expressed in a series of wavenumbers can be written as 1/k[script letter H]. Upon multiplying eqn (27) by k2, the Bo term is of [scr O, script letter O](k2), the term consisting of image file: d5sm00727e-t32.tif is of [scr O, script letter O](1), the curvature term is of [scr O, script letter O](k4), and the term containing [scr D, script letter D]Ca is of [scr O, script letter O](k2). This indicates that for low wavenumber, the stabilizing effect of elastic stresses is dominant over the stabilizing effects of gravity, surface tension, and the destabilizing effect of the potential.

It is noteworthy that the term associated with the elasticity (Wi) of the perfect conductor has an opposite sign compared to the term associated with the critical potential ([scr D, script letter D]). This indicates that the elastic stresses, like gravity, cause a stabilizing effect, whereas the potential leads to a destabilizing effect on the interfacial deformations. It is noteworthy that when a similar linear stability analysis is performed by employing the Oldroyd-B model to capture the behavior of the soft-gel layer, the following expression relating Bo and [scr D, script letter D]Ca is obtained (cf. Appendix D for the derivation), i.e.,

 
image file: d5sm00727e-t33.tif(28)
This expression matches the critical potential for a Newtonian fluid (perfect conductor) and a hydrodynamically passive gas (perfect dielectric).1 The elasticity of the Oldroyd-B fluid does not change the critical potential for the onset of instability. As a result, the Oldroyd-B model does not affect the outcome of the linear stability analysis under neutral conditions. This result contrasts with our current findings, where decreasing the Weissenberg number, Wi, raises the critical potential required for instability. This is the main reason for using the linear viscoelastic model to represent the soft-gel layer in both the linear stability and weakly nonlinear analyses.

Typical neutral stability curves obtained from eqn (27) are shown in Fig. 2 drawn for the example of air on top of polyacrylamide gel. Observe that on the falling branch of the neutral stability curve, i.e., for low wave numbers the electrostatic potential driving the instability is balanced by gravity and elastic stresses, while on the rising branch, i.e., for high wave numbers the voltage applied across the horizontal conducting walls is stabilized by interfacial tension. Interfacial tension is a force acting tangentially along the interface but creates a net restoring force normal (perpendicular) to the interface through curvature effects. The force pulls the interface inward to minimize its surface area and restore stability. Interfacial tension stabilizes the interface by acting as a restoring force that resists deformation, minimizing the surface area between two phases and suppressing disturbances that would increase irregularities. This helps maintain a smooth and stable boundary. Interfacial tension enters the mathematical formulation through the normal stress balance equation in eqn (8) (dimensional form) and via the Capillary number Ca in eqn (14). This stabilizing effect of interfacial tension for high wave numbers is analogous to Rayleigh–Taylor instability of a fluid, wherein the instability is always subcritical in nature. The balance between potential, elastic stresses and gravity for low wave numbers is reminiscent of the Bénard instability which is supercritical in nature. Note that for high Weissenberg numbers i.e., Wi = 1 and Wi = 10, the neutral curve tends towards the critical curve corresponding to that of a Newtonian perfect conductor, i.e., the curve with the solid dots (cf. Fig. 2(a) of Dinesh et al.1). For low Weissenberg numbers, i.e., Wi = 0.1, the critical potential [scr D, script letter D]0Ca for the onset of an instability increases. This increment in the critical potential is attributed to the stabilizing elastic stresses in the soft-gel layer.


image file: d5sm00727e-f2.tif
Fig. 2 [scr D, script letter D]0Ca vs. k2 obtained from eqn (26). The case of air on top of water is assumed. Figure (a) is drawn for [script letter H] = 1, Bo = 5 and Wi = 0.1, 1 and 10. Note that the solid dots in the figure correspond to the case where the perfect conductor is a Newtonian fluid. Figure (b) is drawn for [script letter H] = 1, Bo = 5, and Wi = 0.1 taking image file: d5sm00727e-t199.tif, with n being the horizontal mode index.

Not all wave numbers in Fig. 2a are accessible. When a minimum appears in the curves, the accessible wave numbers, k are related to the horizontal modal index, n, as shown in Fig. 2b for a one-dimensional system of width, w, where image file: d5sm00727e-t34.tif. where n takes integer values 1, 2, 3,…, with n = 1 representing one wave at the interface, n = 2 representing two waves, and so on. Here, w is the width of the container confining the perfect conductor and perfect dielectric between flat plates. To obtain the curve of the critical dimensionless potential versus the wavenumber, the waveform (i.e., the value of n) is fixed while the width “w” of the container is varied. This variation results in a continuous and smooth curve for the dimensionless potential as a function of the wavenumber. There are specific values of “w” where two consecutive modes can coexist, and around these points, the curves are non-monotonic. These points and the wave numbers between them are crucial to our study's conclusions. Our aim is to determine whether the instability is supercritical or subcritical for each accessible wave number by performing a weak nonlinear analysis around the critical state.

4 Weak nonlinear analysis and discussion: linear viscoelastic model

A weak nonlinear analysis is used to study how an interface behaves just beyond the onset of an instability, whether it stabilizes at a finite amplitude (supercritical saturation) or grows toward rupture (subcritical rupture). For the case of the perfect conductor behaving like a linear viscoelastic fluid, the ensuing analysis follows the procedure outlined in the work of Dinesh and Narayanan.1,21 We start with a linear stability analysis to find the critical conditions, the governing equations are then expanded in a small parameter (ε, given below in eqn (29)) that measures deviation from this threshold. Next, we solve the equations order by order, i.e., [scr O, script letter O](ε2) and [scr O, script letter O](ε3). This leads to an equation in terms of the amplitude of the interfacial deformation (which we refer as ζ1 or [scr A, script letter A]). To determine the nature of the bifurcation, the dimensionless potential, [scr D, script letter D], is advanced from its critical value, [scr D, script letter D]0 by an amount ε, such that ε is defined by
 
image file: d5sm00727e-t35.tif(29)
The response of the potential to an increase of the control variable [scr D, script letter D] is given by
 
image file: d5sm00727e-t36.tif(30)
The displacement fields, pressure field and the interface deflection change from their base state by
 
image file: d5sm00727e-t37.tif(31)
 
image file: d5sm00727e-t38.tif(32)
 
image file: d5sm00727e-t39.tif(33)
and
 
image file: d5sm00727e-t40.tif(34)
At the interface, the interior displacement field, X, is expressed as ref. 36
 
image file: d5sm00727e-t41.tif(35)
and likewise for Z(x, z) and p(x, z). The equations at various orders can thus be obtained. The equations at [scr O, script letter O](ε) are given by
 
image file: d5sm00727e-t42.tif(36)
 
image file: d5sm00727e-t43.tif(37)
 
image file: d5sm00727e-t44.tif(38)
The boundary conditions at [scr O, script letter O](ε) at the rigid wall, i.e., z = −1 are
 
Z1 = 0 and X1 = 0 (39)
The interfacial conditions, i.e., at z = 0 at [scr O, script letter O](ε) become
 
Z1ζ1 = 0 (40)
 
image file: d5sm00727e-t45.tif(41)
and
 
image file: d5sm00727e-t46.tif(42)
Observe that the [scr O, script letter O](ε) equations are precisely the neutral stability equations. Hence ζ1(x) becomes ζ1(x) = [scr A, script letter A][thin space (1/6-em)]cos(kx). Our job now is to determine the sign of [scr A, script letter A]2, noting that a positive value of [scr A, script letter A]2 implies a supercritical bifurcation and a negative value implies a subcritical bifurcation. When the system consisting of the perfect conductor and perfect dielectric becomes unstable, interfacial deformations start to grow. If the square of the amplitude, i.e., [scr A, script letter A]2, is positive, the interfacial deformations grow gently and the amplitude of these deformations increases smoothly as we go past the critical potential—this is a supercritical bifurcation. This finally leads to the formation of saturated interfacial deformation. If the square of the amplitude, i.e., [scr A, script letter A]2, is negative, the interface resists disturbances but then suddenly jumps to a large amplitude deformation once pushed away from the critical potential, often showing hysteresis—this is a subcritical bifurcation. In short, supercritical bifurcation means a smooth transition to a saturated interfacial deformation, while subcritical bifurcation means an abrupt jump leading to the rupture of the interface, which ultimately leads to pillars of the perfect conductor spanning the gap between the plates.12 To determine [scr A, script letter A]2, we proceed to the next order.

The governing equations at image file: d5sm00727e-t47.tif are given by

 
image file: d5sm00727e-t48.tif(43)
 
image file: d5sm00727e-t49.tif(44)
and
 
image file: d5sm00727e-t50.tif(45)

The boundary conditions at z = −1 are

 
Z2 = 0[thin space (1/6-em)] and [thin space (1/6-em)]X2 = 0 (46)
and the interfacial conditions, i.e., at z = 0, are now
 
image file: d5sm00727e-t51.tif(47)
 
image file: d5sm00727e-t52.tif(48)
and
 
image file: d5sm00727e-t53.tif(49)
The potential field at image file: d5sm00727e-t54.tif is given by
 
2ψ2 = 0 (50)
subject to the following conditions
 
image file: d5sm00727e-t55.tif(51)
Observe the pattern of eqn in (36)–(42) and compare them with eqn (43)–(51). The right-hand sides are comprised of forcing terms that are quadratic combinations of [scr O, script letter O](ε) terms and directly proportional to [scr A, script letter A]2. We call these quadratic combinations (1, 1) terms due to their bi-linear combinations of first order variables. The solution to the second order problem must therefore be the sum of several [scr A, script letter A]2 dependent terms. In addition, the boxed term in eqn (51) evolves due to the expansion image file: d5sm00727e-t56.tif at image file: d5sm00727e-t57.tif. This term leads to the determination of [scr A, script letter A]2 at image file: d5sm00727e-t58.tif. Upon observing that eqn (36)–(42) are homogeneous and employing solvability conditions on eqn (36)–(42), we see that solvability of eqn (43)–(51) is automatically satisfied. Thus we cannot determine [scr A, script letter A]2 at this order and must advance to the next order. In weak nonlinear analysis, the solvability condition ensures that the perturbation expansions remain bounded by preventing unphysical growth of perturbations. Applied to the normal stress balance equation at the interface, this condition requires that higher-order nonlinear terms satisfy a compatibility criterion, correlating factors like disturbance amplitude, surface tension, elasticity, and electrostatic potential. This balance leads to an amplitude equation that governs the evolution and saturation of interfacial deformations beyond the linear stability threshold, ensuring a physically consistent description of the interface evolution.

The governing equations at image file: d5sm00727e-t59.tif are given by

 
image file: d5sm00727e-t60.tif(52)
 
image file: d5sm00727e-t61.tif(53)
and
 
image file: d5sm00727e-t62.tif(54)
The boundary conditions at z = −1 are
 
Z3 = 0 and [thin space (1/6-em)]X3 = 0 (55)

The interfacial conditions at z = 0 are

 
image file: d5sm00727e-t63.tif(56)
 
image file: d5sm00727e-t64.tif(57)
and
 
image file: d5sm00727e-t65.tif(58)

At the third order the potential field is governed by

 
2ψ3 = 0 (59)

subject to ψ3 = 0 at z = [script letter H] and the following condition at z = 0

 
image file: d5sm00727e-t66.tif(60)
In the absence of electrostatic potential and upon reversing the direction of gravity, the governing equations for the weakly nonlinear regime reduce to the classical Rayleigh–Taylor instability problem. In this limiting case, the Bond number advances beyond its critical value, i.e., Bo = Bo0 + ε2/2. In Rayleigh–Taylor instability, gravity destabilizes the system when a denser fluid overlies a lighter fluid (cf. ref. 21 for weak nonlinear analysis of Rayleigh–Taylor instability of a soft-gel layer). In our current study, the destabilization arises from Maxwell stresses at the soft-gel fluid interface, which causes interfacial deformations. However, the stabilization of the interfacial deformations arises from elastic stresses, gravity and surface tension. It is noteworthy that from eqn (27) term II arises from the Maxwell stresses in the normal force balance equation. This term is solely responsible for driving an instability when a critical potential is applied across the perfect conductor and perfect dielectric fluid layers. However, this term is balanced by the Bond number term, term I, and surface tension. The Bond number term indicates that gravity opposes the potential, term I shows that elastic stresses act against the potential, and the surface tension term also acts in opposition to the potential.

If [scr A, script letter A]2 is positive, the branch is supercritical, but if it is negative, [scr D, script letter D] would not be advanced as indicated by eqn (29) but reduced instead and the branch would become a subcritical pitchfork. The calculations that involve solvability require the use of symbolic manipulation, which was carried out in Mathematica®. The algebraic complications can be reduced if the nonlinear terms of eqn (7) are dropped.21 It has been observed by the present authors that the results do not change qualitatively.

4.1 Solution of the nonlinear equations and tracing the cause of the transition from super to subcritical branching

The solution to the first order governing equations is obtained from the linear stability analysis discussed in the Section 3. We now proceed to seek the solution for the governing equation at image file: d5sm00727e-t67.tif. At this second order, the displacement fields, pressure field and interface deformation must be expressed as
 
image file: d5sm00727e-t68.tif(61)
Here the displacement fields and the pressure field consist of an x-dependent part and an x-independent part. The domain equations for the x-dependent part of the problem are
 
image file: d5sm00727e-t69.tif(62)
 
image file: d5sm00727e-t70.tif(63)
and
 
image file: d5sm00727e-t71.tif(64)
As before, upon eliminating image file: d5sm00727e-t72.tif from the governing eqn (63) and (64), by using the continuity eqn (62), we get
 
image file: d5sm00727e-t73.tif(65)
At z = 0, we have the tangential stress and kinematic conditions. Eliminating image file: d5sm00727e-t74.tif from the x-dependent part of the tangential stress condition, eqn (48), by taking the horizontal derivative and using the continuity equation, (62), we get
 
image file: d5sm00727e-t75.tif(66)
while the x-dependent part of kinematic condition, (47), at z = 0 gives
 
image file: d5sm00727e-t76.tif(67)
Solving eqn (65) and applying the tangential stress condition, (66), along with the kinematic condition, (67), we get
 
image file: d5sm00727e-t77.tif(68)
Employing eqn (68) in the continuity eqn (62), and the x-momentum eqn (63), yields
 
image file: d5sm00727e-t78.tif(69)
and
 
image file: d5sm00727e-t79.tif(70)

At the second order we are left with obtaining the solutions for the x-independent parts of X2, Z2 and p2. To do this, we eliminate pressure from the momentum eqn (44) and (45), and further eliminate X2 by using the continuity eqn (43), resulting in an equation for Z2. The x-independent part of this equation gives

 
image file: d5sm00727e-t80.tif(71)
At z = 0, we have the tangential stress condition and the kinematic conditions. The tangential stress condition for Z20 is obtained by eliminating X2 from the eqn (48) by taking the horizontal derivative and using the continuity equation, (43). The x-independent part of the resulting equation gives
 
image file: d5sm00727e-t81.tif(72)
and at z = 0 the x-independent part of the kinematic condition, (47), yields
 
image file: d5sm00727e-t82.tif(73)
The solution to eqn (71) using the conditions (72) and (73) along with the no displacements conditions at z = −1, gives
 
Z20 = 0 (74)
Applying this solution to the z-momentum eqn (45), yields
 
image file: d5sm00727e-t83.tif(75)
Noting this, we turn toward the remaining equations at image file: d5sm00727e-t84.tif. The potential field is given by
 
2ψ2 = 0 (76)
subject to the following conditions
 
image file: d5sm00727e-t85.tif(77)
The normal component of the momentum balance at the interface, z = 0, is
 
image file: d5sm00727e-t86.tif(78)
where TMaxwell2f consists of the forcing terms which are bilinear, also called (1, 1) terms, because they are products of terms with subscript 1 (cf. Appendix A for the expression of TMaxwell2f). The terms associated with [scr D, script letter D]0 arrive from the expansion of the Maxwell stresses at the second order. We note that all of the forcing terms at this order are bilinear combinations of first order terms, i.e., (1,1) terms with the exception of the boxed term in eqn (77). This implies that the forcing terms are superposition of second harmonics, i.e., cos(2kx) terms and x-independent terms. Their projection on to the eigenspace, i.e., cos(kx) is zero and thus solvability at second order is automatically satisfied. The boxed term in eqn (77) is so identified as it is the sole reason for us to determine [scr A, script letter A]2 at the third order. At the second order, ψ2 is expressed as,
 
image file: d5sm00727e-t87.tif(79)
where image file: d5sm00727e-t88.tif is given by
 
image file: d5sm00727e-t89.tif(80)
and
 
image file: d5sm00727e-t90.tif(81)
and where ψ20(z) in eqn (79) is given by
 
image file: d5sm00727e-t91.tif(82)
The reason for us to be able to calculate [scr A, script letter A]2 at the third order is due to the term ψf20. This term arises from the correction of the base state potential at the second order.

Having obtained the solutions to the first and second order equations in terms of [scr A, script letter A], we now turn to the equations at the third order, i.e. image file: d5sm00727e-t92.tif to determine [scr A, script letter A]2. As noted earlier, the forcing terms in eqn (52)–(60) are bi-linear combinations of second order and first order terms i.e., (2, 1) terms and trilinear combinations of first order terms i.e., (1, 1, 1) terms. From these combinations we can infer that the displacement fields, pressure field and the interface deformation at this order can be expressed as

 
image file: d5sm00727e-t93.tif(83)
and
 
image file: d5sm00727e-t94.tif(84)
At this order only the cos(kx) part of the displacement, pressure and the interface deformation fields (the terms underlined in eqn (83) and (84)) play a role in determining the amplitude, [scr A, script letter A]. Therefore, we see that the domain equations for the cos(kx) part of the variables are
 
image file: d5sm00727e-t95.tif(85)
 
image file: d5sm00727e-t96.tif(86)
and
 
image file: d5sm00727e-t97.tif(87)
Again as before, eliminating [X with combining circumflex]3 from (86) and (87) and using the continuity equation, (85), we get
 
image file: d5sm00727e-t98.tif(88)
At z = 0 we have the tangential stress and the kinematic conditions. The tangential stress is modified by eliminating [X with combining circumflex]3 from the cos(kx) part of the eqn (57). This is done by taking the horizontal derivative of the cos(kx) part of the eqn (57) and by using the continuity equation, (52). This gives
 
image file: d5sm00727e-t99.tif(89)
and at z = 0, the cos(kx) part of the kinematic condition, (56), yields
 
image file: d5sm00727e-t100.tif(90)
Solving eqn (85)–(90) gives
 
image file: d5sm00727e-t101.tif(91)
We split image file: d5sm00727e-t102.tif into two parts, one that is free of ĥ3 and the other that is homogeneous in ĥ3. Thus,
 
image file: d5sm00727e-t103.tif(92)
where image file: d5sm00727e-t104.tif and image file: d5sm00727e-t105.tif.

At the third order the potential field is governed by

 
2ψ3 = 0 (93)
subject to ψ3 = 0 at z = [script letter H] and the following condition at z = 0
 
image file: d5sm00727e-t106.tif(94)
The normal component of the momentum balance along the interface at z = 0 gives
 
image file: d5sm00727e-t107.tif(95)
where TMaxwell3f consists of the (1, 2) and (1, 1, 1) forcing terms (cf. Appendix A for the expression of TMaxwell3f). The boxed term in eqn (95) gives rise to a term in the expression for [scr A, script letter A] which is always negative indicating the subcritical nature of the bifurcation and plays a key role in the high wave number regime, leading to subcritical branching upon instability. It is analogous to a similar term in the Rayleigh–Taylor problem.

At this order, ψ3 is expressed as

 
image file: d5sm00727e-t108.tif(96)
Here, the cos(kx) part of the third order problem alone plays a role in the determination of [scr A, script letter A]2 due to the requirement of the solvability condition at the third order (cf. Appendix B for the solution of image file: d5sm00727e-t109.tif). At the third order, solvability requires that the inhomogeneous terms reside in the null space of the homogeneous problem. This yields an expression for [scr A, script letter A]2, the details of which may be found in the Appendix C. In the limit of Wi → ∞, [scr A, script letter A]2 is given by
 
image file: d5sm00727e-t110.tif(97)
with β = [k(cosh(2k[script letter H]) + 2)csch2(k[script letter H])]/[(2[thin space (1/6-em)]tanh(k[script letter H])(Bo + k2) − 6k2coth(k[script letter H]))/(Bo + k2)]. The above expression for 1/[scr A, script letter A]2 is exactly the same that was reported in Dinesh et al.1 Note that in the limit of Wi → ∞, the bottom fluid mimics the behavior of a Newtonian perfect conductor. We now focus in the expression for 1/[scr A, script letter A]2 when the bottom fluid is considered to be a viscoelastic perfect conductor, i.e.,
 
image file: d5sm00727e-t111.tif(98)
 
image file: d5sm00727e-t112.tif(99)
It is noteworthy that the above expression in eqn (98) for 1/[scr A, script letter A]2 helps us glean the physics of electrostatic potential acting on a viscoelastic perfect conductor-dielectric fluid pair. A positive sign for 1/[scr A, script letter A]2 indicates that the branching behavior for the interface deformation will lead to a supercritical saturation and a negative sign signifies that the branching behavior leads to a subcritical rupture of the interface. Also note that this expression consists of terms containing Wi, these terms play a key role in stabilization of interfacial deformations driven by the electrostatic potential. This stabilization is primarily attributed to the elastic stresses in the perfect conductor that counteract the destabilizing potential. This can be seen from the expression for [scr A, script letter A]2 in eqn (99), which is obtained by expanding [scr A, script letter A]2 from eqn (98) as a series in the wavenumber k. At low wavenumbers, the [scr O, script letter O](k2) term is dominant. This term is positive and includes the Weissenberg number (Wi), indicating that the interface undergoes supercritical saturation. The presence of Wi shows that elastic stresses contribute to the saturation of interfacial deformations. This term arises from the tangential and normal elastic stresses in eqn (41), (42), (48), (49), (57) and (58).

To further glean the physics of electrostatic instability, we calculate the sign of 1/[scr A, script letter A]2 from eqn (98), along the neutral stability curve for Wi = 0.1. In addition, we also calculate the wavenumber beyond which the sign of 1/[scr A, script letter A]2 transits from a positive to negative sign. This transition wavenumber is obtained from the expression for 1/[scr A, script letter A]2 by equating [scr A, script letter A] 2 = 0. The neutral curve and the curve for the transition wavenumber (straight line) are shown in Fig. 3. Note that the region on the neutral stability curve above the straight line which is the transition wavenumber always leads to supercritical saturation of the interface deformation and the region below the straight line leads to subcritical rupture of the interface. The intersection of the neutral curve and the straight line corresponding to the transition wavenumber is a co-dimension 2 point. Observe that in the region for the supercritical behavior, the elastic stresses along with gravity counter act the Maxwell stresses due to the electrostatic potential. These stresses balance out leading to a saturated interface deformation. However, if the fluid were to be Newtonian the Maxwell stresses dominate over gravity and always lead to a subcritical behavior. The region on the neutral curve below the straight line depicts the subcritical branching behavior of the interface. This behavior is attributed to both the Maxwell stresses and the surface tension which ultimately lead to the rupture of the interface. The surface tension term arises from the diminished curvature, i.e., the boxed term in eqn (95) which is of [scr O, script letter O](k3). In short, if the wavenumber k, is chosen such a way that the critical potential lies to the left of the intersection point of the neutral curve and the transition wavenumber, this will lead to a supercritical saturation of the interface. On the contrary, if the critical potential lies to the right of the intersection point, the interfacial deformations lead to a subcritical rupture.


image file: d5sm00727e-f3.tif
Fig. 3 The critical potential versus wavenumber given by the neutral stability curve and the transition wavenumber is represented by the straight line. The critical potential versus wavenumber to the left of the intersection point lead to a supercritical saturation of the interface and the region to the right of the intersection point represents the subcritical rupture of the interface. The results are presented for Wi = 0.1.

4.2 Comments on the use of the Oldroyd-B fluid model for the perfect conductor: linear stability and weakly nonlinear analysis

A similar stability analysis is performed for the case where the perfect conductor is simulated using the Oldroyd-B fluid model (the details of the linear stability analysis are relegated to Appendix D). For an Oldroyd-B fluid modeled as a perfect conductor, the neutral stability relationship involving the critical dimensionless potential, Bond number, and curvature terms remains identical to the relationship observed when the perfect conductor is modeled as a Newtonian fluid. In addition, the expression for [scr A, script letter A]2 in the case of an Oldroyd-B fluid is identical to that obtained when the perfect conductor behaves as a Newtonian fluid (the derivation of the weak nonlinear analysis is provided in Appendix D). This equivalence arises because, in both cases, the perturbed velocities and pressure vanish at orders [scr O, script letter O](ε), image file: d5sm00727e-t113.tif, and image file: d5sm00727e-t114.tif. As a result, the governing equations remain the same for both Newtonian and Oldroyd-B fluids. Consequently, the bifurcation behavior is also identical, with both fluids exhibiting subcritical rupture of the interface.1

5 Conclusions

Using weakly nonlinear analysis, it is shown that an electric field imposed on a linear viscoelastic fluid, modeled as a perfect conductor, adjacent to a dielectric air layer in the presence of gravity can lead to either supercritical or subcritical instability of the interface. The analysis is based on a regular perturbation expansion, where the perturbation parameter is defined in terms of the deviation of the applied potential from its critical value.

The results indicate that for low wavenumbers (i.e., wide containers), the interface deforms smoothly and exhibits supercritical saturation. In contrast, at high wavenumbers (i.e., narrow containers), subcritical branching occurs. An analytical expression shows that supercritical behavior arises due to the stabilizing influence of elastic stresses and gravity, which dominate over the destabilizing Maxwell stresses. On the other hand, the subcritical nature of the instability at high wavenumbers is primarily driven by electrostatic forcing and the reduced influence of curvature, which scales as [scr O, script letter O](k3). These effects are sufficient to overcome the stabilizing contributions from gravity and elasticity.

Additionally, when the perfect conductor is modeled as an Oldroyd-B fluid, the analysis reveals that the branching remains subcritical, resulting in interface rupture. This behavior is consistent with that observed in Newtonian fluids.

Author contributions

G. B.: conceptualization, formal analysis, investigation, methodology, software, visualization, writing; D. B.: conceptualization, methodology, investigation, resources, supervision, writing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available within the article.

Appendices

A Maxwell stresses at image file: d5sm00727e-t115.tif and at image file: d5sm00727e-t116.tif

TMaxwell2f is given by the following expression
 
image file: d5sm00727e-t117.tif(A1)
TMaxwell3f is given by the following expression
 
image file: d5sm00727e-t118.tif(A2)

B Solution to image file: d5sm00727e-t119.tif at image file: d5sm00727e-t120.tif

We only need to determine the solution of the cos(kx) part of the third order problem, because this part alone plays a role in the determination of [scr A, script letter A]2 at this order. The governing equation for image file: d5sm00727e-t121.tif is given by
 
image file: d5sm00727e-t122.tif(B1)
The corresponding boundary conditions for image file: d5sm00727e-t123.tif at z = 0 and at z = [script letter H], are
 
image file: d5sm00727e-t124.tif(B2)
and
 
image file: d5sm00727e-t125.tif(B3)
The solution to the governing equation for image file: d5sm00727e-t126.tif along with its boundary conditions is,
 
image file: d5sm00727e-t127.tif(B4)

C Solvability condition at image file: d5sm00727e-t128.tif

The normal force balance at the third order, yields
 
image file: d5sm00727e-t129.tif(C1)
Recall that the cos(kx) part of the normal force balance at the first order is given by
 
image file: d5sm00727e-t130.tif(C2)
Multiplying eqn (C1) with [small zeta, Greek, circumflex]1[thin space (1/6-em)]cos(kx) and eqn (2) with [small zeta, Greek, circumflex]3[thin space (1/6-em)]cos(kx) and upon subtracting the resulting equations, we get an expression for [scr A, script letter A]2.

D Linear stability and weak nonlinear analysis using the Oldroyd-B model for the perfect conductor

We now perform a weakly nonlinear analysis of the electrohydrodynamic instability when the perfect conductor is represented by an Oldroyd-B model. The governing equations are expanded and solved at successive orders: [scr O, script letter O](ε), image file: d5sm00727e-t131.tif, and image file: d5sm00727e-t132.tif. The ensuing analysis follows the procedure outlined in the work of Dinesh and Narayanan.1,21 In addition, the solution methodology is outlined in Section 4. For a perfect conductor, the dimensionless governing equations using the Oldroyd-B model are given as follows, beginning with the continuity equation:
 
image file: d5sm00727e-t133.tif(C3)
the x-momentum equation
 
image file: d5sm00727e-t134.tif(C4)
and the z-momentum equation
 
image file: d5sm00727e-t135.tif(C5)
where τxx, τxz, τzx, τzz are components of the polymeric stress tensor. The xx-component of shear stress is given by
 
image file: d5sm00727e-t136.tif(C6)
The zz-component of shear stress yields
 
image file: d5sm00727e-t137.tif(C7)
The xz-component of shear stress is represented by
 
image file: d5sm00727e-t138.tif(C8)
Here image file: d5sm00727e-t139.tif, image file: d5sm00727e-t140.tif and η = μp/μ* are the Reynolds number, Deborah number and the viscosity ratio, i.e., the viscosity of the polymer divided by the viscosity of the solvent.

We now perform a linear stability analysis of the governing equations following a similar procedure presented in Section 3. The linearized form of the above equations under neutral conditions, yield,

 
image file: d5sm00727e-t141.tif(C9)
the x-momentum equation
 
image file: d5sm00727e-t142.tif(C10)
and the z-momentum equation
 
image file: d5sm00727e-t143.tif(C11)
The xx-component of shear stress is
 
image file: d5sm00727e-t144.tif(C12)
The zz-component of shear stress is given by
 
image file: d5sm00727e-t145.tif(C13)
The xz-component of shear stress is represented by
 
image file: d5sm00727e-t146.tif(C14)
It is noteworthy that the Deborah number, De, is absent in the governing equations corresponding to the shear stress components of the Oldroyd-B model. The Deborah number which takes into the elastic nature of the fluid an important characteristic of the viscoelastic fluid is absent under the neutral conditions, indicating that the elastic nature of the fluid does not alter the critical potential required for an instability. The Deborah number measures how fast stresses relax compared to the flow time scale. In neutral stability, perturbations are marginal with no net growth or decay, so relaxation does not influence the neutral boundary but only affects the rate at which instabilities grow or decay away from it. Consequently, De appears in transient growth rate problems but not in neutral stability curves. It can also be shown that under neutral stability conditions, the perturbed velocity and pressure are zero. Toward this end, we calculate the perturbed velocity in the z-direction by eliminating pressure and the x-component velocity in the governing equations, eqn (C9)–(C14). This yields,
 
image file: d5sm00727e-t147.tif(C15)
The solution for ŵ can be written as
 
ŵ = c1ekz + c2ekz + c3zekz + c4zekz (C16)
The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,
 
image file: d5sm00727e-t148.tif(C17)
Upon solving the governing equations along with the boundary conditions, we get ŵ = 0. Since the perturbed velocities ŵ and û are zero, the perturbed pressure [p with combining circumflex] is also zero. Therefore, upon substituting the perturbed quantities in the linearized normal stress balance yields an expression for the critical potential, i.e.,
 
image file: d5sm00727e-t149.tif(C18)
Note that this expression is similar to the critical potential for a Newtonian fluid (recall eqn (28) in Section 3). This is attributed to the fact that the perturbed velocities and pressure are zero at [scr O, script letter O](ε). Therefore, the Oldroyd-B model for the soft-gel layer results in a critical potential that is independent of the Deborah number. We now focus on the solution of the governing equations at image file: d5sm00727e-t150.tif.

D.1 image file: d5sm00727e-t151.tif Governing equations and solution. The governing equations at image file: d5sm00727e-t152.tif are given by the continuity equation, i.e.,
 
image file: d5sm00727e-t153.tif(C19)
the x-momentum equation
 
image file: d5sm00727e-t154.tif(C20)
and the z-momentum equation
 
image file: d5sm00727e-t155.tif(C21)
where (τxx)2, (τxz)2, (τzx)2, (τzz)2 are components of the polymeric stress tensor at image file: d5sm00727e-t156.tif. The xx-component of shear stress is given by
 
image file: d5sm00727e-t157.tif(C22)
The zz-component of shear stress is
 
image file: d5sm00727e-t158.tif(C23)
and the xz-component of shear stress can be written as
 
image file: d5sm00727e-t159.tif(C24)
It is noteworthy that in the above equations [scr O, script letter O](ε) terms are zero, i.e., the first order velocity components in the viscoelastic fluid layer are zero. In addition, all the time derivative terms are also zero since we are conducting a weak nonlinear analysis around the neutral stability point. Therefore, [scr O, script letter O](ε) terms and the time derivative terms vanish from the momentum and in the polymeric stress tensor equations.

We now calculate the perturbed velocity in the z-direction by eliminating the pressure and the x-component velocity in the momentum equations. This yields,

 
image file: d5sm00727e-t160.tif(C25)
At this order, velocities, pressure and interface formation are written as
 
image file: d5sm00727e-t161.tif(C26)
 
image file: d5sm00727e-t162.tif(C27)
The solution for image file: d5sm00727e-t163.tif can be written as
 
image file: d5sm00727e-t164.tif(C28)
The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,
 
image file: d5sm00727e-t165.tif(C29)
Upon solving the governing equations along with the boundary conditions, we get image file: d5sm00727e-t166.tif. Note that image file: d5sm00727e-t167.tif will play a key role in the calculations at [scr O, script letter O](ε3/6).
 
image file: d5sm00727e-t168.tif(C30)
The solution for image file: d5sm00727e-t169.tif can be written as
 
w20 = c1z3 + c2z2 + c3z + c4 (C31)
The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of w20 are given by,
 
image file: d5sm00727e-t170.tif(C32)
Upon solving the governing equations along with the boundary conditions, we get w20 = 0. It is noteworthy that w20 = 0, will play a key role in the calculations at [scr O, script letter O](ε3/6). The potential field is given by
 
2ψ2 = 0 (C33)
subject to the following conditions
 
image file: d5sm00727e-t171.tif(C34)
The normal component of the momentum balance at the interface, z = 0, is
 
image file: d5sm00727e-t172.tif(C35)
(cf. Appendix A for the expression of TMaxwell2f) where image file: d5sm00727e-t173.tif is given by
 
image file: d5sm00727e-t174.tif(C36)
and
 
image file: d5sm00727e-t175.tif(C37)
and where ψ20(z) is given by
 
image file: d5sm00727e-t176.tif(C38)
A detailed solution for the potential is provided in Dinesh et al.1 We now proceed to the solution of the governing equations at image file: d5sm00727e-t177.tif.

D.2 image file: d5sm00727e-t178.tif Governing equations and solution. The governing equations at image file: d5sm00727e-t179.tif are given by
 
image file: d5sm00727e-t180.tif(C39)
the x-momentum equation
 
image file: d5sm00727e-t181.tif(C40)
and the z-momentum equation
 
image file: d5sm00727e-t182.tif(C41)
where (τxx)3, (τxz)3, (τzx)3, (τzz)3 are components of the polymeric stress tensor at image file: d5sm00727e-t183.tif. The xx-component of shear stress is expressed as
 
image file: d5sm00727e-t184.tif(C42)
The zz-component of shear stress can be written as
 
image file: d5sm00727e-t185.tif(C43)
and the xz-component of shear stress is represented by
 
image file: d5sm00727e-t186.tif(C44)
It is noteworthy that in the above equations [scr O, script letter O](ε) and image file: d5sm00727e-t187.tif terms are zero, i.e., the first and second order velocity components and pressure in the viscoelastic fluid layer are zero. Therefore, we do not see the product of these [scr O, script letter O](ε) terms and product of [scr O, script letter O](ε) term with image file: d5sm00727e-t188.tif in the momentum equations and the polymeric stress tensor equations. Toward this end, we calculate the perturbed velocity in the z-direction by eliminating pressure and the x-component velocity in the momentum equations. This yields,
 
image file: d5sm00727e-t189.tif(C45)
At this order, velocities, pressure and interface deformation are written as
 
image file: d5sm00727e-t190.tif(C46)
 
image file: d5sm00727e-t191.tif(C47)
The solution for ŵ3 can be written as
 
ŵ3 = c1ekz + c2ekz + c3zekz + c4zekz (C48)
The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,
 
image file: d5sm00727e-t192.tif(C49)
Upon solving the governing equations along with the boundary conditions, we get ŵ3 = 0. At the third order the potential field is governed by
 
2ψ3 = 0 (C50)
subject to ψ3 = 0 at z = [script letter H] and the following condition at z = 0
 
image file: d5sm00727e-t193.tif(C51)
The normal component of the momentum balance along the interface at z = 0 gives
 
image file: d5sm00727e-t194.tif(C52)
where TMaxwell3f consists of the (1, 2) and (1, 1, 1) forcing terms (cf. Appendix A for the expression of TMaxwell3f). Therefore, at this order, ψ3 and ζ3 are expressed as
 
image file: d5sm00727e-t195.tif(C53)
Here, the cos(kx) part of the third order problem alone plays a role in the determination of [scr A, script letter A]2 due to the requirement of the solvability condition at the third order. The solvability condition yields an expression for [scr A, script letter A]2, i.e.,
 
image file: d5sm00727e-t196.tif(C54)
with β = [k(cosh(2k[script letter H]) + 2)csch2(k[script letter H])]/[(2[thin space (1/6-em)]tanh(k[script letter H])([scr B, script letter B] + k2) − 6k2[thin space (1/6-em)]coth(k[script letter H]))/([scr B, script letter B] + k2)]. The expression for [scr A, script letter A]2 in the case of an Oldroyd-B fluid is identical to that obtained when the perfect conductor behaves as a Newtonian fluid. This equivalence arises because, in both cases, the perturbed velocities and pressure vanish at orders [scr O, script letter O](ε), image file: d5sm00727e-t197.tif, and image file: d5sm00727e-t198.tif. As a result, the governing equations remain the same for both Newtonian and Oldroyd-B fluids. Consequently, the bifurcation behavior is also identical, with both fluids exhibiting subcritical rupture of the interface.1 The Deborah number characterizes the relative importance of elastic relaxation compared to the characteristic time scale of the flow, indicating whether the material responds more like a solid (large De) or more like a fluid (small De). In the case of neutral stability, perturbations are marginal and neither grow nor decay. At this threshold, stress relaxation does not contribute to shifting the onset of instability, since the system is balanced precisely between stable and unstable behavior. Instead, the Deborah number influences the dynamics away from neutrality by determining how rapidly instabilities amplify or decay. Consequently, De plays a central role in governing transient growth rates of perturbations, but it does not explicitly appear in the neutral stability curves that define the onset of instability.

Acknowledgements

The authors gratefully acknowledge funding from Ministry of Higher Education (MHRD), Anusandhan National Research Foundation (ANRF) via the grant numbers ITS/2024/005066 and ANRF/ECRG/2024/001265/ENS.

References

  1. B. Dinesh and R. Narayanan, Phys. Rev. Fluids, 2021, 6, 054001 CrossRef.
  2. G. I. Taylor, Proc. R. Soc. London, Ser. A, 1966, 291, 159–166 Search PubMed.
  3. J. R. Melcher and G. I. Taylor, Annu. Rev. Fluid Mech., 1969, 1, 111–146 CrossRef.
  4. E. SchaÈffer, T. Thurn-Albrecht, T. P. Russell and U. Steiner, Nature, 2000, 403, 874–877 CrossRef PubMed.
  5. V. Shankar and A. Sharma, J. Colloid Interface Sci., 2004, 274, 294–308 CrossRef CAS PubMed.
  6. R. M. Thaokar and V. Kumaran, Phys. Fluids, 2005, 17, 084104 CrossRef.
  7. L. Espn, A. Corbett and S. Kumar, J. Non-Newtonian Fluid Mech., 2013, 196, 102–111 CrossRef.
  8. K. Mondal, P. Kumar and D. Bandyopadhyay, J. Chem. Phys., 2013, 138, 1–15 CrossRef.
  9. J.-U. Park, M. Hardy, S. J. Kang, K. Barton, K. Adair, D. K. Mukhopadhyay, C. Y. Lee, M. S. Strano, A. G. Alleyne and J. G. Georgiadis, et al., Nat. Mater., 2007, 6, 782–789 CrossRef CAS PubMed.
  10. D. Riccobelli and P. Ciarletta, Philos. Trans. R. Soc., A, 2017, 375, 20160421 CrossRef PubMed.
  11. J. Marthelot, E. F. Strong, P. M. Reis and P. T. Brun, Nat. Commun., 2018, 9, 1–7 CrossRef CAS PubMed.
  12. D. T. Papageorgiou, Annu. Rev. Fluid Mech., 2019, 51, 155–187 CrossRef.
  13. S. Y. Chou and L. Zhuang, J. Vac. Sci. Technol., B: Microelectron. Process. Phenom., 1999, 17, 3197–3202 CrossRef CAS.
  14. S. A. Roberts and S. Kumar, J. Fluid Mech., 2009, 631, 255 CrossRef.
  15. P. Gambhire and R. M. Thaokar, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 036301 CrossRef CAS.
  16. K. Ward, S. Matsumoto and R. Narayanan, J. Fluid Mech., 2019, 862, 696–731 CrossRef CAS.
  17. D. Koulova-Nenova, J. Electrost., 1997, 40, 185–190 CrossRef.
  18. P. Colinet, J. C. Legros and M. G. Velarde, Nonlinear dynamics of surface-tension-driven instabilities, Wiley-VCH, Berlin, 2001 Search PubMed.
  19. L. Wu and S. Y. Chou, J. Non-Newtonian Fluid Mech., 2005, 125, 91–99 CrossRef CAS.
  20. G. Tomar, V. Shankar, A. Sharma and G. Biswas, J. Non-Newtonian Fluid Mech., 2007, 143, 120–130 CrossRef CAS.
  21. B. Dinesh and R. Narayanan, J. Fluid Mech., 2021, 915, 1–27 Search PubMed.
  22. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Pergamon Press, 1989 Search PubMed.
  23. V. Shankar and V. Kumaran, J. Fluid Mech., 2000, 407, 291–314 CrossRef.
  24. B. Dinesh and S. Pushpavanam, Phys. Rev. E, 2017, 96, 013119 CrossRef CAS PubMed.
  25. R. Patne, D. Giribabu and V. Shankar, J. Fluid Mech., 2017, 827, 31–66 CrossRef CAS.
  26. V. Kumaran, J. Fluid Mech., 1995, 294, 259–281 CrossRef CAS.
  27. V. Kumaran, J. Fluid Mech., 1995, 302, 117–139 CrossRef.
  28. V. Kumaran, J. Fluid Mech., 1998, 357, 123–140 CrossRef CAS.
  29. B. Dinesh and S. Pushpavanam, Phys. Fluids, 2018, 30, 104103 CrossRef.
  30. P. Howell, G. Kozyreff and J. Ockendon, Applied solid mechanics, Cambridge University Press, 2009 Search PubMed.
  31. P. Ramkarn and V. Shankar, J. Fluid Mech., 2018, 860, 837–885 Search PubMed.
  32. J. Snoeijer, A. Pandey, M. Herrada and J. Eggers, Proc. R. Soc. A, 2020, 476, 20200419 CrossRef CAS.
  33. D. A. Saville, Annu. Rev. Fluid Mech., 1997, 29, 27–64 CrossRef.
  34. S. Mora, T. Phou, J.-M. Fromental and Y. Pomeau, Phys. Rev. Lett., 2014, 113, 178301 CrossRef PubMed.
  35. H. W. Müller and W. Zimmermann, EPL, 1999, 45, 169 CrossRef.
  36. L. E. Johns and R. Narayanan, Interfacial instability, Springer Science & Business Media, 2002, pp. 255–304 Search PubMed.

This journal is © The Royal Society of Chemistry 2025
Click here to see how this site uses Cookies. View our privacy policy here.