Nonlinear EHD Instability of Two-Superposed Walters’ B Fluids Moving through Porous Media

Abstract

The current work examines the application of the viscous potential flow to the Kelvin-Helmholtz instability (KHI) of a planar interface between two visco-elastic Walters’ B fluids. The fluids are fully saturated in porous media in the presence of heat and mass transfer across the interface. Additionally, the structure is pervaded via a uniform, normal electrical field in the absence of superficial charges. The nonlinear scheme basically depends on analyzing the linear principal equation of motion, and then applying the appropriate nonlinear boundary-conditions. The current organization creates a nonlinear characteristic equation describing the amplitude performance of the surface waves. The classical Routh–Hrutwitz theory is employed to judge the linear stability criteria. Once more, the implication of the multiple time scale with the aid of Taylor theory yields a Ginzburg–Landau equation, which controls the nonlinear stability criteria. Furthermore, the Poincaré–Lindstedt technique is implemented to achieve an analytic estimated bounded solution for the surface deflection. Many special cases draw upon appropriate data selections. Finally, all theoretical findings are numerically confirmed in such a way that ensures the effectiveness of various physical parameters.

1. Introduction

Electrohydrodynamics (EHD) can be described by the classical theory of magnetism and electricity and it has drawn a great deal of attention from many authors. Furthermore, EHD is of supreme significance to several problems in practical engineering, especially with the departure of charge particles occurring in colloids, Deoxyribonucleic acid proteins, cells, and numerous additional elements of organic concentration. Devitt and Melcher [1] showed that the field coupled surface wave happens at the interface between two fluids when it is stressed by an electric field. El-Sayed [2] used the multiple-time scales, technique to explore the nonlinear modulation of the interfacial waves of two superposed dielectric fluids with uniform depths and horizontal borders under the influence of a uniform normal electric field. He discovered that a nonlinear Schrödinger equation can describe quasi-monochromatic travelling waves. Papageorgiou and Vanden-Broeck [3] examined the surface waves along with the action of a uniform horizontal field. The fluid was taken as an inviscid, incompressible and nonconducting, and the controlled waves were numerically calculated for arbitrary wavelengths and amplitudes. Furthermore, the numerical technique was grounded on considering a system of integro-differential nonlinear equations. Grandison et al. [4] studied the interfacial capillary waves in the presence of electric fields. Their nonlinear solutions were computed by using the boundary integral methods. The investigation included both the symmetric and anti-symmetric modes. The EHD nonlinear surface wave instability through porous media under the influence of a uniform field was investigated by El-Sayed [5]. It was found that the surface of separation is administrated by a Schrödinger type. Consequently, the nonlinear stability profile was examined along with the coefficients of Schrödinger equation. Aldini and Seyed-Yagoobi [6] investigated how the EHD instability induction pumps can express itself in a rapid drop/jump in pump output. As a result of the instability, alternating/bidirectional flow might occur. A non-dimensional stability analysis of EHD induction pumping of liquid film in a vertical annular configuration in the presence of an external load for the repulsion mode was carried out. Burcham and Saville [7] examined EHD stability of a liquid bridge suspended in a dielectric gas by Taylor–Melcher theory. Amer and Moatimid [8] have recently investigated EHD instability of a streaming dielectric liquid jet. The system was pervaded by a uniform axial electrostatic field. They demonstrated that the ratio between the basic velocities plays a dual role in the instability configuration.

Fluids are all around us. From a practical point of view, the fluid behaviour is indispensable, i.e., aircrafts and ships travel through fluids; lubricants for mechanical maneuvers are fluids; the atmosphere and the weather are administrated by fluid dynamics. One of the most important branches of fluids is non-Newtonian fluids. In recent decades, a great deal of interest has been drawn to the comprehensive examination of the rheological influences happening during the flow of the non-Newtonian fluids which are saturated in porous media. This phenomenon has become of special interest to oil reservoir engineering. Subsequently, it has become the rheological application of non-Newtonian supplant and expatriate fluids is vital to oil extraction. Numerous technical procedures include the analogous flow of fluids with different density, viscosity, and elasticity throughout permeable media. Such flows appear in petroleum engineering, boiling in porous media, packed bed reactors, chemical industry, and in several additional procedures. The unstable interface causes a considerable growth in the resistance to the flow occurs. EHD surface wave stability between two superimposed visco-elastic Walters’ B fluids that are saturated in porous media was investigated by El-Sayed [9]. Under the influence of a tangential electric field, El-Sayed [10] investigated the EHD Kelvin–Helmholtz instability of planner interface between two uniform suspended viscous and flowing dielectric fluids penetrated with suspended particles through porous media. He found that the presence of both streaming and tangential electric field had an effect on the disturbances. A nonlinear stability of surface wave between two conducting fluids was examined by Zakaria et al. [11]. Based on the multiple scale approach, they derived a nonlinear Schrödinger equation with complex coefficients. Consequently, the stability criteria were accomplished. In the same vein, such instabilities lead to emulsion formation in petroleum production engineering. Therefore, the information about the circumstances of the instability will allow us to guess the restrictive process circumstances. As known, the flow in porous media is of utmost importance to oil engineers and geophysical fluid. The stability of a horizontal interface between two fluids of different densities is named the Rayleigh–Taylor instability (RTI). KHI of Oldroydian visco-elastic fluids in a permeable medium finds its main usefulness in biochemical technology and geophysical fluid dynamics. It is also assumed to be more suitable for use in the oil industry. El-Sayed [12] examined EHD atomization instability and Rayleigh regimes for dielectric liquid jet emitting a parabolic velocity profile into a stationary dielectric gas through a porous medium. Moatimid and Zeky [13] examined the EHD nonlinear stability of a cylindrical interface between two visco-elastic fluids of Walters’ B type. The flow was saturated in porous media. The nonlinear instability approach resulted in a Ginzburg–Landau type. Therefore, the instability standards were analytically attained and numerically confirmed. Additionally, along with the nonlinear expanded frequency, an analytic bounded approximate solution of the surface deflection was originated.

Typically, in practical engineering, the understanding of the mechanism of heat and mass transfer has become of great importance. It has bearing on the vehicles, electronic devices, and structures. The inclusion of the energy/concentration equations takes into account the governing momentum equation. The interface stability phenomenon is examined under the presumption that the two-fluids are immiscible. Consequently, no mass transfer at the interface happens. Simultaneously, thermal influence acts as a portion of slight roles, so the impact of heat transmission is canceled. Indeed, mass and heat transfer through the interface is of paramount importance in a wide range of fields of practical engineering and physics. If heat transfer across the interface occurs, mass transfer also happens. Commonly, the latent heat becomes very large. Thomas and Hartnett proposed an adequate approximation of immiscible fluids [14]. Taking these considerations into account, and apart from the treatment of concentration with an energy equation, a shortened formulation of the problem of heat and mass at the interface stability has been modeled by Hsieh [15]. Hsieh [15] applied this simplification to examine the RTI. A complicated dispersion relation was achieved in understanding these problems. Furthermore, despite the validity of the traditional stability criterion, Hsieh [16] formulated the influence of heat and mass transfer, which simplified the stability configuration. It was shown that the stability criterion, in the case of KHI, is substantially modified due to the influence of heat and mass transmission. The nonlinear RTI in the action of heat and mass transfer was examined by Hsieh [17]. It is shown that when this transfer is sufficient enough, the characteristically unstable profile is stabilized in accordance with the nonlinear effect. Many authors followed Hsieh’s simplified formulation to examine the stability problems under the influence of mass and heat transfer; for illustration, see Nayak and Chakraborty [18]. It was shown that KHI has a destabilizing influence on the stability profile. Moreover, KHI is less in a plane than in cylindrical geometry. Lee [19] investigated the nonlinear Rayleigh–Taylor stability of a cylindrical contact between vapor and liquid phases of a fluid. The investigation was conducted using the multiple time scale expansion method. The regions of stability/instability were graphically depicted. Özgen [20] investigated the stability and transition difficulties of two-dimensional boundary-layers with hot malls numerically using linear stability theory. An effective shoot–search strategy was used to solve the equations found.

The perturbation methods [21] were initially established for treating the weak nonlinear phenomenon. Therefore, numerous investigators aimed to incorporate them in nonlinear aspects. All of the classical perturbation techniques need the presence of small parameter. In light of perturbation theory, Lindstedt–Poincaré (L–P) method is a technique to obtain a uniform analytic approximate periodic solution to an ordinary differential equation, which cannot be reached by the classical perturbation approaches. The method eliminates secular terms present in the straightforward presentation of perturbation theory that is applied to a weakness nonlinear problem with a finite oscillatory solution. The method is named after Henri Poincaré (1893) and Anders Lindstedt (1882). The L–P technique was the most commonly utilized analytical approach in analyzing some problems of nonlinear oscillations with a small parameter. Conversely, there are several nonlinear problems that do not include small parameters. Consequently, a new analytical technique would be industrialized so that the constraint slight parameters can be disregarded. Burton [22] offered an adapted form of L–P technique aimed to manage strong nonlinear oscillators. Cheung et al. [23] modified the L–P technique. Nevertheless, all the improved varieties remain actual operative for the Duffing oscillator. Alam et al. [24] aimed to use the modified L–P technique for numerous nonlinear oscillators. The beginning of these formulations and their operations for penetrating a proper approximate solution are very easy.

In accordance with the aforementioned aspects, the stability problem of the two visco-elastic fluids is of great importance. Therefore, the current problem is concerned with a nonlinear EHD instability problem. Two fluids are assumed to obey the Walters’ B type. A uniform vertical electric field is acted upon by the two media. Due to the wide practical applications of porous media, the media under consideration represent a limited amount of porous media. Additionally, the upper and lower boundaries remain deliberated to be horizontally rigid porous plates. The presence of constant suction/injection velocities is also included. Therefore, the plates were steaming normal uniform velocities. To facilitate mathematical manipulation, the hypothesis of the viscous potential are considered. To simplify the performance of the problem, the rest of the manuscript is structured as follows: Section 2 presents the methodology of the theoretical model and includes the foremost equations of motion and the consistent suitable nonlinear boundary conditions. Furthermore, the technique of clarification by the resources of the normal mode analysis is given in this Section. The linear dispersion relation and stability analysis are introduced in Section 3. The linear stability analysis of the linear methodology is conducted in this Section. The nonlinear stability yields a Ginzburg–Landau equation, and the hypothetical and numerical calculations are obtained in Section 4. In accordance with the P–L technique, an analytic bounded approximate solution of the surface deflection is derived in Section 5. Furthermore, a numerical approximation of the surface wave elevation is introduced in this Section. The results of the study are summed up as concluding observations in Section 6. This Section includes the physical findings yielded from the analysis of linear/nonlinear stability profile, together with an annotation about a future work.

2. Methodology

A scheme consisting of two standardized, dielectric, incompressible, and streaming visco-elastic fluids of Walters’ B type in infinite and parallel flows is considered. The fluids are considered as uniform, isotropic, and fully saturated in porous media. Dynamic viscosities, dielectric constants, and a single medium permeability are presumed. For simplicity, the porosity is taken as unity in both media. A statically stable situation is supposed in such a way that the upper fluid is assumed to be light (vapor or gas); meanwhile, the lower one is assumed to be heavier (liquid). For more convenience, the Cartesian coordinates will be utilized. Unlike the mainstream in mathematical analyses, only two-dimensional motion is tackled in the current study. The undisturbed interface between the two fluids is presupposed to be well addressed and is originally horizontal to yield the plane $y=0$. Moreover, the two fluids are considered as moving with uniform parallel velocities during the two superimposed permeable media. The upper together with the lower boundaries are considered as horizontally rigid permeable plates. The presence of a constant suction/injection streaming at the lower and upper plates is deliberated. Therefore, the plates permit a steaming velocity in the orthogonal direction of the flow channel. The gravitational force acts perpendicularly downwards. The upper fluid is bounded by plane $y=h_2$, which is raised to a uniform potential. The lower liquid is restricted by an earthed conducting plane at $y=-h_1$. The potential difference between the two rigid plates generates a uniform, vertical electric field $E_{0j}$. Additionally, the rigid plates are reserved at constant temperatures. Meanwhile, the undisturbed interface $y=0$ is kept at a uniform temperature. Therefore, the influence of mass and heat transfer is analyzed, for simplicity, in light of the Hsieh [15,16] simplified formulation. The theoretical physical model is sketched below in Figure 1.

After assuming a slight disturbance, the boundary function develops as follows:

$$\mathsf{\Sigma }(x,y,t)=y-\eta (x;t)$$

(1)

Henceforth, the unit outward normal vector is given as

$$\underset{¯}{n}=\nabla \mathsf{\Sigma }/\left|\nabla \mathsf{\Sigma }\right|=\frac{1}{\sqrt{1+{(\partial \eta /\partial x)}^2}}\left(-\frac{\partial \eta }{\partial x} {\underset{¯}{e}}_x+{\underset{¯}{e}}_y\right)$$

(2)

Characteristically, as argued by Funada and Joseph [25,26,27] and others, it has been demonstrated that the potential flow, in the case of an incompressible fluid, is considered by $\underset{¯}{v}=-\nabla \phi$. It follows that the Navier–Stokes equation exhibits similar behavior to Euler’s. Consequently, the term $\mu {\nabla }^2\underset{¯}{v}$ will disappear from the substance of the fluid. In a few words, along with the potential flow, the pressure and velocity are the same as in the inviscid flow. On the other hand, the viscosity of the fluids enters only at the surface of separation. Consequently, the equation of a viscous incompressible fluid, in the case of the potential flow through the porous media, obeys Brinkman–Darcy equation as given by the following momentum equation:

$${\rho }_j \left(\frac{\partial {\underset{¯}{v}}_j}{\partial t}+({\underset{¯}{v}}_j·\nabla ){\underset{¯}{v}}_j\right)=-\nabla P_j-\frac{{\mu }_j}{\kappa }{\underset{¯}{v}}_j-{\rho }_{j}g {\underset{¯}{e}}_y, \hspace{1em}\hspace{1em}\hspace{1em}j=1, 2.$$

(3)

The frictional forces are a consequence of the interactions between the fluid and the porous media. It is comparable to the flow velocity which is represented by the term $\frac{{\mu }_j}{\kappa }{\underset{¯}{v}}_j$.

The zero-order solution of Equation (3) yields

$$P_{0j}=-\frac{{\mu }_j}{\kappa }(V_{j}y+U_{j}x)-{\rho }_{j}gy+{\lambda }_j,$$

(4)

In accordance with the viscous potential theory and the incompressibility condition of the two fluids, the potentials ${\phi }_1(x,y;t)$ and ${\phi }_2(x,y;t)$ should satisfy the following Laplace’s equations:

$${\nabla }^2{\phi }_1=0, -h_1<y<\eta (x,t)$$

(5)

and

$${\nabla }^2{\phi }_2=0, \eta (x,t) <y<h_2$$

(6)

The integration of the linear equation of motion as given in Equation (3) results in Bernoulli’s formulation, which gives the following distribution function of the pressure:

$$P_j=-{\rho }_j\left(\frac{\partial {\phi }_j}{\partial t}+V_j\frac{\partial {\phi }_j}{\partial y} +U_j\frac{\partial {\phi }_j}{\partial x}\right)-\frac{{\mu }_j}{\kappa }{\phi }_j$$

(7)

Typically, in EHD problems, it is usually assumed that the quasi-static approximation is effective [28]. Therefore, the electrical equations are

$$\nabla ·\epsilon \underset{¯}{E}=0,$$

(8)

and

$$\nabla \times \underset{¯}{E}=0.$$

(9)

As a result of the potential difference between the planes $y=h_2$ and $y=- h_1$, a normal electric field is created in two regions. Hence, the electric field may be written as:

$${\underset{¯}{E}}_j=-E_{0j}^{} {\underset{¯}{e}}_y-\nabla {\psi }_j .$$

(10)

It follows that the electrostatic potentials satisfy the following Laplace equations:

$${\nabla }^2{\psi }_1=0, -h_1<y<\eta (x,t)$$

(11)

and

$${\nabla }^2{\psi }_2=0, \eta (x,t) <y<h_2.$$

(12)

To complete the considered nonlinear boundary-value problem, the solutions for the above principal equations must satisfy the following boundary conditions.

Nonlinear Boundary Conditions

The overall solutions for the velocity and electric potential distributions must satisfy the applicable nonlinear boundary conditions. These circumstances may be categorized into two distinctive groups as follows:

On the Rigid Restrictions $y=h_2$ and $y=-h_1$.

• As the rigid boundaries are isolated, the tangential components of the electric field must vanish; it follows that

  $$\frac{\partial {\psi }_1}{\partial x}=0,\hspace{1em}\hspace{1em}\hspace{1em}y=-h_1,$$

  (13)

  and

  $$\frac{\partial {\psi }_2}{\partial x}=0,\hspace{1em}\hspace{1em}\hspace{1em}y=h_2.$$

  (14)
• The suction/injection velocities in the normal direction at the permeable boundaries require

  $$\frac{\partial {\phi }_1}{\partial y}=V_1,\hspace{1em}\hspace{1em}\hspace{1em}y=-h_1,$$

  (15)

  and

  $$\frac{\partial {\phi }_2}{\partial y}=V_2,\hspace{1em}\hspace{1em}\hspace{1em}y=h_2.$$

  (16)

At the Perturbed Boundary at $y=\eta (x, t)$.

The tangential components of the electric field are continuous at the interface. This may be written as follows:

$$\underset{¯}{n}\times ‖\underset{¯}{E}‖=\underset{¯}{0},$$

(17)

where $‖\ast ‖$ represents the jump across the interface.

Equation (17) gives

$${\frac{\partial \eta }{\partial x}}_{ }‖E_{0j}+\frac{\partial {\psi }_j}{\partial y}‖+ ‖\frac{\partial {\psi }_j}{\partial x}‖=0$$

(18)

In the absence of the free surface currents, it follows that the normal components of the electric induction vector are continuous, which gives

$$\underset{¯}{n}·‖\epsilon \underset{¯}{E}‖=0$$

(19)

.

Once more, the Equation (19) yields

$$\frac{\partial \eta }{\partial x} ‖ {\epsilon }_j \frac{\partial {\psi }_j}{\partial x}‖-‖ {\epsilon }_j\frac{\partial {\psi }_j}{\partial y}‖=0$$

(20)

One may consider ${\underset{¯}{v}}_j=U_j {\underset{¯}{e}}_x+V_j {\underset{¯}{e}}_y+\nabla {\phi }_j$ to be the velocity vector field of the fluid particles. Following Hsieh [15,16,17], one gets the following conditions:

The conservation of mass across the surface of separation yields

$$‖{\rho }_j \left(\frac{\partial S}{\partial t}+{\underset{¯}{v}}_j·\nabla S\right)‖ =0.$$

(21)

This interfacial condition of energy transfer, in the nonlinear form, was given by Moatimid et al. [6] to yield

$${\rho }_1 \left(\frac{\partial S}{\partial t}+{\underset{¯}{v}}_1·\nabla S\right)={\alpha }_1\left(\eta +{\alpha }_2{\eta }^2+{\alpha }_3{\eta }^3\right),$$

(22)

where heat and transfer coefficients are given as follows:

$${\alpha }_1=\frac{G}{L} \left(\frac{1}{h_1}+\frac{1}{h_2}\right), {\alpha }_2=\frac{1}{h_2}-\frac{1}{h_1}, \hspace{1em}\mathrm{and}\hspace{1em}{\alpha }_3=\frac{1}{h_2^2}-\frac{1}{h_1 h_2}+\frac{1}{h_1^2},$$

(23)

where $G=\frac{K_2(T_0-T_2)}{h_2-y}=\frac{K_1(T_1-T_0)}{h_1+y}$ indicate that, in the steadiness formal, the heat changes are equal through the interface in the two fluids.

In accordance with these conditions, the solutions of the velocities and electric potentials may be written as follows:

$${\phi }_1(x,y,t) =\frac{\cosh \left(k(y+h_1)\right)}{k \left(1-i {\eta }_x\coth k{h}_1\right)\sinh k{h}_1}\left( U_1{\eta }_x+{\eta }_t+\frac{{\alpha }_1}{{\rho }_1}(\eta +{\alpha }_2{\eta }^2+{\alpha }_3{\eta }^3)\right) ,$$

(24)

$${\phi }_2(x,y,t) =-\frac{\cosh \left(k(y-h_2)\right)}{k \left(1+i {\eta }_x\coth k{h}_2 \right)\sinh k{h}_2}\left( U_2{\eta }_x+ {\eta }_t+\frac{{\alpha }_1}{{\rho }_2}(\eta +{\alpha }_2{\eta }^2+{\alpha }_3{\eta }^3)\right),$$

(25)

$${\psi }_1(x,y,t) =\frac{i }{{\epsilon }^{\ast }}{\epsilon }_2({\epsilon }_1-{\epsilon }_2)\sinh [k(y+h_1)] {\eta }_x \left(\cosh k{h}_2+i {\eta }_x\sinh k{h}_2 \right) \mathrm{\Phi }$$

(26)

and

$${\psi }_2(x,y,t) =\frac{i }{{\epsilon }^{\ast }}{\epsilon }_1({\epsilon }_1-{\epsilon }_2)\sinh \left[k(y-h_2)\right]{\eta }_x \left(\cosh k{h}_1-i {\eta }_x \sinh k{h}_1 \right) \mathrm{\Phi }$$

(27)

where

$${\epsilon }^{\ast }=k(h_1{\epsilon }_2+h_2{\epsilon }_1) \left[{\epsilon }_1 \left(1-i {\eta }_x\tanh k{h}_1\right)\left(\tanh k{h}_2+i{\eta }_x\right)+{\epsilon }_2 \left(1+i {\eta }_x\tanh k{h}_2\right)\left(\tanh k{h}_1-i {\eta }_x\right)\right] \cosh k{h}_1 \cosh k{h}_2,$$

and

$$E_{01}-E_{02}=\frac{({\epsilon }_1-{\epsilon }_2)}{({\epsilon }_2{h}_1+{\epsilon }_1{h}_2)}\mathrm{\Phi }$$

It should be noted that the above equations are compatible with the earlier results obtained by Moatimid et al. [6] in the special case where the electric field is taken as a horizontal one. In addition, as the linear terms are ignored, the velocity potential distributions are in agreement with those obtained by Awasthi et al. [29]. It should be noted that the aforementioned distributions of the potential functions ${\psi }_j$ and ${\phi }_j$ contain the nonlinear terms in the restriction $\eta$. This nonlinearity occurs in accordance with the nonlinear boundary conditions that are exemplified overhead.

To examine the stability of the physical model, the remaining boundary condition arises from the normal component of the stress tensor. In agreement with the occurrence of the quantity of surface tension, the normal component must be discontinuous. The visco-elastic force of the Walters’ B type is given by Tonekaboni et al. [30] as

$${\sigma }_{ij}^{viscoelastic}=-P{\delta }_{ij}+2\left({\mu }_j-{{\mu }^{\prime }}_j \frac{\partial }{\partial t}\right) {\underset{¯}{e}}_{ij},$$

(28)

$${\underset{¯}{e}}_{ij}=\frac{1}{2} \left(\nabla v+\nabla v^T\right) ,$$

(29)

The interfacial condition for the preservation of momentum is referred to by Moatimid et al. [6] and Kumar [31]. Therefore, one gets

$$\begin{array}{l}{\rho }_1\left[v_1·\nabla S\right] \left(\frac{\partial S}{\partial t}+v_1\cdot \nabla S\right)={\rho }_2\left[v_2·\nabla S\right] \left(\frac{\partial S}{\partial t}+v_2\cdot \nabla S\right)+\\ \left( ‖P_j‖-2 ‖\left({\mu }_j-{{\mu }^{\prime }}_j \frac{\partial }{\partial t}\right) \underset{¯}{n}· \left[\underset{¯}{n} \cdot {\underset{¯}{e}}_{ij}\right] ‖- \frac{1}{2}‖\epsilon (E_n^2-E_t^2)‖ +\sigma \nabla \cdot \underset{¯}{n}\right) {\left|\nabla S\right|}^2\end{array} \mathrm{at}\hspace{1em}y=\eta ,$$

(30)

The pressure will be eliminated by making use of Bernoulli’s equation. Considering the Darcy’s model of the flow throughout porous media, Equation (30) may be written as

$$a_0{\eta }_{tt}+a_1{\eta }_{xt}+a_2 {\eta }_{xx}+(a_3+i b_1) {\eta }_x+(a_4+i b_2) {\eta }_t+(a_5+i b_3) \eta =N_1(\eta )+N_2(\eta ),$$

(31)

where the coefficients $a_i$ and $b_i$ are real constants; they are given in the Appendix A. In addition, the nonlinear terms $N_1(\eta )$ and $N_2(\eta )$ represent the quadratic and cubic nonlinear terms in the variable $\eta$.

The zero-order interfacial condition of the conservation of momentum gives

$$\frac{x}{\kappa }\left({\zeta }_1{U}_1-{\zeta }_2{U}_2\right)+\frac{1}{2} \left({\epsilon }_2 E_{02}^2-{\epsilon }_1 E_{01}^2\right)+{\lambda }_2-{\lambda }_1=0$$

(32)

In what follows, ignoring the right hand side of Equation (31), the linear stability approach will be introduced throughout the next Section.

3. Linear Stability Analysis

In accordance with of this approach, the linearized analysis of the nonlinear equation given in Equation (31) is obtained when the nonlinear terms of the surface deflection are disregarded. Consequently, the linearized dispersion equation may be formulated as follows:

$${\eta }_{tt}+a_1{\eta }_{xt}+a_2 {\eta }_{xx}+(a_3+i b_1) {\eta }_x+(a_4+i b_2) {\eta }_t+(a_5+i b_3) \eta = 0 .$$

(33)

Therefore, in view of the normal modes analysis, one may accept a constant wave train solution for Equation (33) in the following form:

$$\eta (x;t)=\gamma e^{i(k x-\omega t)}+c·c ,$$

(34)

For a nontrivial solution of Equation (33), the dispersion relation then becomes

$${\omega }^2+(F_1+i G_1) \omega +(F_2+i G_2)=0,$$

(35)

where the coefficients $F_1,G_1, F_2$ and $G_2$ are known of the background of the problem.

To shorten the numerical evaluations, it is preferable to work out the distinguishing equation in terms of the following appropriate non-dimensional quantities. This can be completed in different ways depending mainly on the selectivity of the physiognomies of time and length. For this purpose, consider the parameters $\sqrt{{\rho }_1{h}_1^3/\sigma }$, $h_1$ and ${\rho }_1{h}_1^3$ to indicate the characteristics of time, length and mass, respectively. The other non-dimensional quantities might be assumed as:

$${\rho }_2=\rho {\rho }_1, {\mu }_2=\mu {\mu }_1, {{\mu }^{\prime }}_2={\mu }^{\prime } {{\mu }^{\prime }}_1, U_2=\tilde{U}{U}_1,V_2=V{V}_1, {\epsilon }_2=\epsilon {\epsilon }_2,h_2=h h_1, \mathrm{and} {\mathrm{\Phi }}^{\ast 2}={\mathrm{\Phi }}^2 {\epsilon }_1{h}_1/\sigma .$$

For straightforwardness, the “$\ast$” mark may be ignored in the subsequent analysis. Subsequently utilizing this selection, the succeeding non-dimensional numbers will seem in the dispersion relation as follows:

• Weber numbers $W{e}_{\tilde{U}}={\rho }_1 U_1^2 h_1/\sigma$ and $W{e}_V={\rho }_1 V_1^2 h_1/\sigma$, for the horizontal and vertical velocities, respectively, distinguish the proportion of the disturbing hydrodynamic forces from the soothing surface tension;
• Ohnesorge number $Oh={\mu }_1/\sqrt{\sigma {\rho }_1{h}_1}$ relates the viscous forces to an inertial and surface tension force, where, ${\mu }_j={\rho }_j{\zeta }_j$ represents the dynamic viscosity and ${\zeta }_j$ represents the kinematic viscosity;
• Darcy number $Da=\kappa /h_1^2$ signifies the comparative conclusion of the permeability of the medium against its cross-sectional area;
• Elasticity numbe $El={\mu }^{\prime }/{\rho }_1{h}_1^2$ characterizes the capability of a frame to attack a miss representing effect and to return to its original magnitude and form, when the impact of force is disregarded;
• Bond number $Bo=h_1^2 {\rho }_1 g/\sigma$ represents the ratio of gravitational force to surface tension force;
• Potential Bond number ${\mathrm{\Phi }}^2 {\epsilon }_1{h}_1/\sigma$ represents the ratio between the electric potential and the surface tension force.

The stability standards of the dispersion relation (35) are referred to by the Routh–Hrutwitz theory; for illustration, see Zahreddin and El-Shehawey [32]. It follows that the stability criteria may be written as:

$$F_1>0$$

(36)

and

$$F_1^2{F}_2+F_1{G}_1{G}_2-G_2^2>0.$$

(37)

Although the calculation showed that $F_1$ is independent of ${\mathrm{\Phi }}^2$, it will be taken into consideration in the forgoing numerical analysis. According to the non-dimensional chosen data, the inequality (36) is satisfied in domains $k<0.0119$ and $k>0.0604$. The distributions are discussed numerically and illustrated graphically in the region $k>0.0605$. Instantaneously, the additional inequality of Equation (37) might be correlated to ${\mathrm{\Phi }}^2$ in the subsequent equation:

$$A {\mathrm{\Phi }}^2+B>0,$$

(38)

where $A$ and $B$ are known from the context.

As specified beforehand, the repercussion of the inequality in Equation (38) must be encompassed. Consequently, all the succeeding figures are planned in a definite domain, where the condition of Equation (37) is automatically verified. In addition, the evaluations designated that the parameter $A$ is continuously of negative implication. This displays a destabilizing effect of the normal electric field, which is an early result. It is confirmed by numerous authors; for instance, see Dvitt and Melcher [1], and many references cited therein.

To this end, our interest focuses on the inequality of Equation (38). For this objective, the potential Bond number $Log {\mathrm{\Phi }}^2$ will be designed against the wavenumber $k$ of the surface waves. In the following figures, the stable region is referred to by letter $S$. Meanwhile, letter $U$ stands for the unstable region. It is suitable to designate the impact of several physical parameters on the stability configuration. The following figures are designed for a system, taking the subsequent specifics:

$$Bo=1, \rho =0.25, \mu =0.5, \epsilon =0.1, {\mu }^{\prime }=10, \tilde{U}=-15, h=2.5, {\alpha }_1=2.4, Oh=0.8, W{e}_{\tilde{U}}=2,W{e}_V=4,Da=4,El=30, V=-0.5.$$

Figure 2 illustrates the stability profile drawn over the range of the wavenumber. It is found that the domain of the wave number is portioned into three parts as shown in the figure: in the left region, where $0<k<0.0119$, it is found that $F_1>0$ and $A>0$, which shows that the potential Bound number plays a stability influence. In the middle part, where $0.0119<k<0.0604$, the term $F_1<0$ shows that the stability criterion is not satisfied. On the other hand, where $k>0.0604$, the term $F_1>0$ shows a destabilizing influence of the normal electric field. As a conclusion, contrary to the previous studies of the inviscid fluids, the normal electric field plays a dual role in the stability portrait.

Figure 3 indicates the inspection of the ratio between the two thicknesses on the stability profile. As seen from this figure, as the thickness of the upper layer is greater than the lower one, the stable region increases. This mechanism is enhanced for large values of the wavenumber. This displays a stabilizing influence of this parameter of the selected input parameters. Similar results were earlier gained by Moatimid and Zekry [13].

As given throughout the non-dimensional procedure, the ratio between the viscosities $\mu$ represents the ratio of the viscosity of the second fluid to the first one. The behavior of this parameter is sketched in Figure 4. As seen in this figure, as the parameter $\mu$ increases, the stable region also increases. This indicates a stabilizing effect of this limitation on the designated input parameters. As in the previous figure, the stability mechanism is enhanced for large values of the wavenumber of the surface waves. Similar results were earlier obtained.

From the definition of $W{e}_{\tilde{U}}$, it is noted that its increase may be produced by either a growth of the liquid stream velocity, thicknesses, and density, or by reducing surface tension. From the examination of Figure 5, the influence of $W{e}_{\tilde{U}}$ on the stability profile is revealed. All physical parameters are held fixed except $W{e}_{\tilde{U}}$. This figure shows that $W{e}_{\tilde{U}}$ plays a destabilizing influence. Actually, the increase of this parameter causes the initial streaming to increase. The destabilizing influence of KHI is an early phenomenon verified by many authors cited herein. More specifically, the consequence is in conformity with those obtained earlier.

Figure 6 designates the impacts of $W{e}_V$ on the stability representation. The physical restrictions are entirely fixed, excluding $W{e}_V$. As realized, the figure displays that $W{e}_V$ exerts a stabilizing influence. Indeed, the increase of this parameter means an increase in the suction/injection in the lower layer. This consequence agrees with the finding that has been previously attained. The comparison between Figure 4; Figure 5 displays that the Weber numbers have a dual role in the stability profile.

Figure 7 indicates the influence of Darcy number $Da$ on the stability profile. The inspection of the figure shows that this parameter exerts a destabilizing influence. This mechanism is enhanced great values of the wavenumber. This result is in conformity with the findings that were previously found by Moatimid and Amer [33], El-Sayed et al. [34], and many other researchers.

Figure 8 indicates the presence of the linear mass and heat transfer parameters ${\alpha }_1$ on the stability picture. As illustrated from this figure, the stable region grows as mass and heat transfer increases. This mechanism shows a stabilizing influence ${\alpha }_1$, especially at large values of wave number $k$. This mechanism corresponds to the findings that were reported earlier by Hsieh [16] and Elfenawy [35].

Figure 9 illustrates the effect of Bond number on the interfacial stability. It is found that when the values of Bond number increase, the stability region increases too. As seen from the mathematical formula of this parameter, it has the same properties of gravity. Physically, this is an acceptable result which has been in agreement with that which was earlier obtained by Chen and Chen [36].

The velocity $V$ represents the ratio of the velocities of the upper plate to the lower one of the suction/injection. The positive standards of $V$ denote an injection at the permeable upper layer $y=h_2$ with a consistent suction at the wall $y=-h_1$; meanwhile, the negative values mean an injection at the permeable plate $y=-h_1$ with a corresponding suction on the other plate. To indicate the suction/injection phenomena, Figure 10; Figure 11 are drawn. Therefore, Figure 10 illustrates the result of suction/injection velocities on the porous boundaries of the flow channel. The inspection of this figure shows that the suction/injection at both boundaries has a stabilizing effect. On the other hand, if $V_1$ and $V_2$ have the same sign, the value of the ratio $V$ has a positive value. The effect of the vertical velocity $V$ ratio has been displayed in Figure 11. As the parameter $V$ increases, the disturbance will grow faster and the system becomes a destabilizing situation. Therefore, this restriction has a destabilizing effect on the stability behaviour. This shows that the suction/injection plays a dual role in the stability picture.

Figure 12 designates the influence of the Ohnesorge number $Oh$ on the stability examination. All physical restrictions remain except the parameter $Oh$. This figure shows that $Oh$ has a stabilizing influence. A similar effect has been recently reported by Amer and Moatimid [8] and many other researchers; see Kourmatzis et al. [37].

Figure 13 designates the influence of the elasticity quantity $El$ on the stability outline. All physical structures are maintained except for $El$. As represented in this figure, as the values of elasticity number $El$ increase, the stability region also increases. Therefore, the elasticity number $El$ has a stabilizing influence on the considered system in the occurrence of a vertical electric field. This result is in conformity with that which was gained earlier by El-Sayed et al. [38].

4. Nonlinear Ginzburg–Landau Equation

As specified throughout the linear stability analysis, the interface surface elevation $\eta =\eta (x,t)$ has a special form by substituting Equation (34) into Equation (31). Therefore, the nonlinear characteristic equation in terms of the interfacial distance $\eta$ may be formulated as follows:

$$D(k,\omega )\eta =\gamma (k,\omega ){\eta }^2+\beta (k,\omega ){\eta }^3$$

(39)

The terms $D\left(\omega , k\right)$,$\gamma \left(\omega , k\right)$ and $\beta \left(\omega , k\right)$ are recognized from the background of the problem. To avoid the excessive length of the article, they are all excluded.

The following exploration depends mainly on the multiple time scale approach [39]. This performance is contingent principally on a slight parameter. It processes the ratio of a characteristic wavelength and/or periodic time of variation. Consequently, one may accept that $\delta$ is a trivial limitation that describes the deliberate modulation. In the view of this method, the independent variables $x$ and $t$, which are addressed on the scale of the characteristic wavelength and period time, can be represented as alternate and independent variables,

$$X_n={\delta }^{n}x\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}\hspace{1em}{T}_n={\delta }^{n}t\hspace{1em}\hspace{1em}\hspace{1em}n=0,1,2,\dots$$

(40)

where

$$L\left(\frac{\partial }{\partial x},\frac{\partial }{\partial t}\right)\eta =0.$$

(41)

Therefore, the independent variables $X_0, T_0$ characterize the fast distinctions. Meanwile, $X_1, T_1, X_2, T_2$ refer to the deliberate ones. The differential operators can be currently connected as the following derivative expansions:

$$\frac{\partial }{\partial x}\equiv k\frac{\partial }{\partial {\theta }_0}+\delta \frac{\partial }{\partial X_1}+{\delta }^2\frac{\partial }{\partial X_2}+\dots \hspace{1em}\hspace{1em}\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}\hspace{1em}\frac{\partial }{\partial t}\equiv -\omega \frac{\partial }{\partial {\theta }_0}+\delta \frac{\partial }{\partial T_1}+{\delta }^2\frac{\partial }{\partial T_2}+\dots ,$$

(42)

where $\theta =k{X}_0-\omega T_0$ denotes the lowermost order.

It is satisfactory to enlarge the operator $L$ in the following form:

$$L\left((ik,-i\omega )+i\delta (\frac{\partial }{\partial X_1},\frac{\partial }{\partial T_1})+i{\delta }^2(\frac{\partial }{\partial X_2},\frac{\partial }{\partial T_2})+\dots \right)$$

(43)

The operator $L$ can be extended powers of $\delta$. By means of Taylor’s explanation around $(k,-\omega )$, we recall merely the terms up to $O({\delta }^2)$. Based on this explanation, one finds

$$L\to L_0+\delta L_1+{\delta }^2{L}_2+\dots ,$$

(44)

where

$$L_0\equiv (-\omega \frac{\partial }{\partial {\theta }_0},k\frac{\partial }{\partial {\theta }_0}),$$

(45)

$$L_1\equiv i\left(\frac{\partial L_0}{\partial \omega }\right)\frac{\partial }{\partial T_1}-i\left(\frac{\partial L_0}{\partial k}\right)\frac{\partial }{\partial X_1} ,$$

(46)

and

$$L_2\equiv i\left(\frac{\partial L_0}{\partial \omega }\right)\frac{\partial }{\partial T_2}-i\left(\frac{\partial L_0}{\partial k}\right)\frac{\partial }{\partial X_2}-\frac{1}{2}\left(\frac{{\partial }^2{L}_0}{\partial {\omega }^2}\right)\frac{{\partial }^2}{\partial T_1^2}-\frac{1}{2}\left(\frac{{\partial }^2{L}_0}{\partial k^2}\right)\frac{{\partial }^2}{\partial X_1^2}+\frac{1}{2}\left(\frac{{\partial }^2{L}_0}{\partial k \partial \omega }\right)\frac{{\partial }^2}{\partial X_1\partial T_1}.$$

(47)

Representing the development of the operator as given in Equation (44) into Equation (34), one finds

$$\left(L_0+\delta L_1+{\delta }^2{L}_2\right) \eta =0 .$$

(48)

The abovementioned analysis, reviews the perturbation technique to attain a uniform effective solution. Certainly, this concept needs invalidation of the secular terms.

The stability process of Equation (39) examined in detail in our preceding paper given by Moatimid et al. [13]. Based on the previous findings in this reference, one obtains the following Ginzburg–Landau equation:

$$i\frac{\partial \gamma }{\partial \tau }+(P_r+i P_i)\frac{{\partial }^2\gamma }{\partial {\zeta }^2}=(Q_r+i Q_i) {\gamma }^2\overline{\gamma },$$

(49)

where $\overline{\gamma }$ is the complex conjugate of $\gamma$,

$$\begin{array}{cc}{P}_r+i P_i=& -\frac{1}{2}{\left(\frac{\partial D}{\partial \omega }\right)}^{-1}\left(V_g^2\frac{{\partial }^{2}D}{\partial {\omega }^2}+2V_g\frac{{\partial }^{2}D}{\partial \omega \partial k}+\frac{{\partial }^{2}D}{\partial k^2}\right),\\ & Q_r+i Q_i={\left(\frac{\partial D}{\partial \omega }\right)}^{-1}\left(\frac{2\gamma }{\mathrm{\Omega }}+3\beta \right),\\ & \zeta =\delta (x-V_{g}t), \hspace{1em} \hspace{1em} \tau ={\delta }^{2}t,\end{array}$$

(50)

and the group velocity is given by $V_g=-\frac{\partial D}{\partial k}{\left(\frac{\partial D}{\partial \omega }\right)}^{-1}.$

The non-zero denominator $\mathrm{\Omega }$ results from the linear dispersion function $D(\omega , k)$ by replacing both $\omega$ and $k$ by $2\omega$ and $2k$, correspondingly. The stability standards of the Ginzburg–Landau as given in Equation (49) have been previously established by Lange and Newell [40]. By putting the assumptions of linear perturbation into practice for this equation, they postulated that the perturbation is stable throughout the subsequent conditions:

$$Q_i<0,\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{P}_r Q_r+P_i Q_i>0.$$

(51)

On the other hand, the system develops unstable. Accordingly, the evolution curves that distinguish the stable from the unstable areas are equivalent to

$$Q_i=0,\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{P}_r Q_r+P_i Q_i=0.$$

(52)

Before distributing the numerical calculations, the transition curves given by Equation (52) must be inscribed in a convenient non-dimensional procedure. This can be prepared in a number of methods fundamentally based on the special appearances of time, length and mass. Consider that the parameters $1/\omega$, $h_1$ and $\sigma /{\omega }^2$ denote the appearances of time, length and mass, respectively.

The calculation of the transition curve given by $Q_i=0$ may be represented in a third-degree polynomial of ${\mathrm{\Phi }}^2$ as

$$L_3 {({\mathrm{\Phi }}^2)}^3+L_2 {({\mathrm{\Phi }}^2)}^2+L_1 {\mathrm{\Phi }}^2+L_0=0 .$$

(53)

Meanwhile, the second transition curve given $P_r Q_r+P_i Q_i=0$ can be arranged in a fifth-degree polynomial on ${\mathrm{\Phi }}^2$ as:

$$G_5 {({\mathrm{\Phi }}^2)}^5+G_4 {({\mathrm{\Phi }}^2)}^4+G_3 {({\mathrm{\Phi }}^2)}^3+G_2 {({\mathrm{\Phi }}^2)}^2+G_1 {\mathrm{\Phi }}^2+G_0=0 ,$$

(54)

where the coefficients given in Equations (53) and (54) are well recognized from the background. To decrease the length of the paper, they will be crossed out.

In order to exemplify the stability standards in this nonlinear stability methodology, the characteristic curves that are arranged in Equations (53) and (54) will be numerically analyzed. In what follows, the numerical calculations consider two different data choices as follows:

Data A consider the following particular system:

$$\rho =3, \mu =0.5, \epsilon =0.1, {\mu }^{\prime }=0.5, \tilde{U}=-15, h=2.5, {\alpha }_1=0.4, \widehat{\kappa }=4, V=-0.5, {\alpha }_2=10, \mathrm{and} {\alpha }_3=25$$

Based on the choice of Data A, Equation (53) yields only one positive real root; meanwhile, the additional double roots are complex conjugates. Essentially, this is an algebraic concept. Therefore, only one transition curve can be graphed. On the other hand, Equation (54) also results in one, real root. For this objective, one seeks a common domain of the wave number at which both of the two roots occur. The other complex conjugate roots have no influence on the stability profile. Consequently, $Log {\mathrm{\Phi }}^2$ is planned versus the wave number of the surface waves $k$ throughout Figure 14. The transition curve that is given by Equation (53) is signified by a purple curve. Meanwhile, the red curve refers to the transition curves that are given in Equation (54). In contrast with the linear stability approach, these transition curves divide the stability diagram into three regions. To address the nature of each region, we check the occurrence of the inequalities that are given in (51). As seen from Figure 14, the plane is divided into several parts of stability/instability regions.

Data B considers the following particular system:

$$\rho =3, \mu =0.25, \epsilon =4.5, {\mu }^{\prime }=15, \tilde{U}=7, h=2.5, {\alpha }_1=2.5, \widehat{\kappa }=0.3, V=5, {\alpha }_2=30, \mathrm{and} {\alpha }_3=1.6$$

Along with this choice, the numerical calculation showed, as before, that there is only one real positive root of Equation (53). On the other hand, Equation (54) yields three positive roots as follows: two of them are integrated to yield only one curve. The third one gives a separate curve. The other roots of Equations (53) and (54) are complex conjugates. Therefore, they have no implication in the stability diagram. Figure 15 shows these different transition curves. As shown, the stability diagram is separated by numerous portions of stability/instability constituencies. In order to illustrate the influences of the different physical parameters along with the considered nonlinear approach, and for the sake of simplicity, the chosen sample system will be compatible with data A.

Figure 16, Figure 17 and Figure 18 designate the influences of mass and heat transfer parameters ${\alpha }_1$, ${\alpha }_2$ and ${\alpha }_3$, respectively, on the stability picture. As shown from the theoretical procedure, the parameter ${\alpha }_1$ appears along with the linear stability approach. Meanwhile, the parameters ${\alpha }_2$ and ${\alpha }_3$ emerge due to the nonlinear sense. The inspection of Figure 16 shows that ${\alpha }_1$ plays a stabilizing role. This influence remains compatible with the linear stability approach given in Figure 8. Depending on the nonlinear approach, Figure 17 shows that ${\alpha }_2$ also has a stabilizing effect. On the other hand, Figure 18 demonstrates that ${\alpha }_3$ having a destabilizing influence. Therefore, one can say that mass and heat transfer parameters have a dual role in the stability criterion. These results are in agreement with the results that were earlier obtained by Hsieh [16].

Figure 19 indicates the influences of the ratio of Walters’ B visco-elasticity ${\mu }^{\prime }$ on the stability picture. As understood from this figure, the increase of this stricture enhances the instability regions. Therefore, it consumes a destabilizing consequence. This consequence is in arrangement with the findings that were previously reported by El-Sayed et al. [41].

Figure 20 and Figure 21 illustrate the suction/injection influence throughout the nonlinear stability approach. As previously seen, when $V_1$ and $V_2$ have the same sign, the ratio $V$ has a positive value. It is seen in Figure 20 that this parameter has a destabilizing effect on the stability behaviour as revealed by our discussion in linear stability analysis.

On the other hand, Figure 21 indicates the consequence of suction/injection velocities on the porous borders of the flow channel. We know that if $V_1$ takes a negative value at the lower plate and $V_2$ takes a positive value at the upper plate; this condition is called the suction velocity, while the inverse occurs for the injection velocity. The inspection of this figure shows that the suction/injection at both boundaries has a stabilizing effect. This displays that the suction/injection plays a dual role through the non-linear stability profile. Similar results were obtained earlier.

5. Profile Interface by Means of Poincaré–Lindstedt Technique

The aim of this Section is to achieve an approximate bounded solution for the surface elevation. As shown in the aforementioned section, the nonlinear methodology results in the nonlinear distinguishing relation that is given in Equation (31). It characterizes a nonlinear second-order differential equation through complex coefficients of the interface displacement $\eta (x,t)$. Really, the analysis of this equation, in its current form, is slightly difficult. Physically, the behavior of the amplitude deflection $\eta (x,t)$ is essentially a real function. To simplify the subsequent calculations, one may consider only the time dependent, i.e., $\eta =\eta (0,t)=\gamma (t)$. Consequently, the previous distinguishing Equation (31) may be divided into its real and imaginary parts as follows:

$${\gamma }^{″}+l_1 {\gamma }^{\prime }+l_2 \gamma +l_{3 }{\gamma }^2+l_4 \gamma {\gamma }^{\prime }+l_5 \gamma {\gamma }^{″}+l_6 {\gamma }^3+l_7 {\gamma }^2 {\gamma }^{\prime }+l_8 {\gamma }^2 {\gamma }^{″}=0,$$

(55)

and

$$m_1 {\gamma }^{\prime }+m_2 \gamma +m_{3 }{\gamma }^2+m_4 \gamma {\gamma }^{\prime }+m_5 {\gamma }^3+m_6 {\gamma }^2 {\gamma }^{\prime }=0,$$

(56)

where $l_i, m_i (i=1,2,....)$ are recognized from the context. To restrict the length of the paper, they will be crossed out from the text.

To be more understandable, the mixture of Equations (55) and (56) may be attained by terminating the term ${\gamma }^{\prime }$, and the subsequent equation can be written as follows:

$${\gamma }^{″}+{\mathrm{\Omega }}^2 \gamma +s_{1 }{\gamma }^2+s_2 {\gamma }^3=0,$$

(57)

where ${{\Omega }}^2,s_i, (i=1,2)$ is recognized from the context. This amount affects the ordinary frequency of the problem. Characteristically, it is necessarily positive. To decrease the length of the paper, they will be omitted.

At this stage, the nonlinear amplitude as given in Equation (57) has real coefficients. It signifies a comprehensive cubic nonlinear differential equation. It is sometimes called Rayleigh-Duffing equation. The solution for this equation needs initial conditions. For this purpose, the following original conditions may be introduced:

$$\gamma \left(0\right)=1; \mathrm{and} {\gamma }^{\prime }(0)=0.$$

(58)

For this objective, the homotopy formulation of the considered parametric equation becomes

$${\gamma }^{″}+{\mathrm{\Omega }}^2 \gamma +\chi \left(s_{1 }{\gamma }^2+s_2 {\gamma }^3\right)=0; \chi \in \left[0, 1\right].$$

(59)

Typically, a new variable $\tau =\varpi t$ is presented to substitute $t$, and then Equation (59) develops

$${\varpi }^2\ddot{\gamma }+{\mathrm{\Omega }}^2 \gamma +\chi \left(s_{1 }{\gamma }^2+s_2 {\gamma }^3\right)=0.$$

(60)

Herein $\varpi$ is recognized as the frequency of the oscillator, and the prime denotes the differentiation with respect to $\tau$. Allowing the L–P technique, see [42,43,44,45,46,47,48], and the parameters $\gamma$ and $\varpi$ may be expended as follows:

$$\gamma (\tau )={\displaystyle \sum _{j=0}^{\infty }{\chi }^j{\gamma }_j(\tau )=}{\gamma }_0(\tau )+\chi {\gamma }_1(\tau )+{\chi }^2{\gamma }_2(\tau )+.... \mathrm{and} \varpi =1+\chi {\varpi }_1+{\chi }^2{\varpi }_2+........$$

(61)

Substituting from Equations (61) into Equation (60) and identifying the coefficients of the same powers of $\chi$ on both sides, one finds the subsequent grading equations:

$${\chi }^0: {\gamma }_0(\tau )=\cos (\mathrm{\Omega } \tau ) ,$$

(62)

$$\chi : {\gamma }_1(\tau )=L_T^{-1}\left\{\frac{S}{S^2+{\mathrm{\Omega }}^2}\right\}-L_T^{-1}\left\{\frac{S}{S^2+{\mathrm{\Omega }}^2}{L}_T\left\{s_1{\gamma }_0^2+s_2{\gamma }_0^3+2{\varpi }_1{\ddot{\gamma }}_0\right\}\right\},$$

(63)

and

$${\chi }^2: {\gamma }_2(\tau )=L_T^{-1}\left\{\frac{S}{S^2+{\mathrm{\Omega }}^2}\right\}-L_T^{-1}\left\{\frac{1}{s^2+{\mathrm{\Omega }}^2}{L}_T\left\{2s_1{\gamma }_0{\gamma }_1+3s_2{\gamma }_0^2{\gamma }_1+{\varpi }_1^2{\ddot{\gamma }}_0+2{\varpi }_2{\ddot{\gamma }}_0+2{\varpi }_1{\ddot{\gamma }}_1\right\}\right\}.$$

(64)

On replacing Equation (60) into Equation (61), we find

$${\gamma }_1(\tau )=L_T^{-1}\left[\frac{1}{S^2+{\mathrm{\Omega }}^2}\right]-L_T^{-1}\left[\frac{1}{S^2+{\mathrm{\Omega }}^2}{L}_T\left\{\frac{1}{4}\left(2s_1+3s_2\cos ( \tau )-8 {\mathrm{\Omega }}^2{\varpi }_1\cos ( \tau )+2s_1\cos (2\mathrm{\Omega } \tau )+s_2\cos (3\mathrm{\Omega } \tau )\right)\right\}\right].$$

(65)

The development necessitates the termination of the secular terms. Consequently, the coefficient of the function $\cos (\mathrm{\Omega } \tau )$ must be omitted. Hence, the parameter ${\varpi }_1$ is determined as

$${\varpi }_1=\frac{3 s_2}{8{\mathrm{\Omega }}^2}.$$

(66)

It follows that the periodic solution ${\gamma }_1(\tau )$ becomes

$${\gamma }_1(\tau )=\frac{s_1}{2{\mathrm{\Omega }}^2}+\left(1+\frac{s_2}{32 {\mathrm{\Omega }}^2}-\frac{s_1}{3 {\sigma }^2}\right)\cos (\mathrm{\Omega }\tau )+\frac{s_1}{6 {\mathrm{\Omega }}^2}\cos (2\mathrm{\Omega }\tau )-\frac{s_2}{32 {\mathrm{\Omega }}^2}\cos (3\mathrm{\Omega }\tau ).$$

(67)

Once more, replacing Equations (62) and (67) into Equation (64), the invalidation of the secular term needs

$${\varpi }_2=\frac{320s_1^2-192s_1{s}_2-45s_2^2+576s_2{\mathrm{\Omega }}^2}{768 {\mathrm{\Omega }}^4}.$$

(68)

By similar arguments as assumed before, after a long, but straightforward calculations, one finds the solution ${\gamma }_2(\tau )$ is

$$\begin{array}{l}{\gamma }_2(\tau )=\frac{s_1\left(32s_1+63s_2+96{\mathrm{\Omega }}^2\right)}{96{\mathrm{\Omega }}^4}+\left(1+\frac{29s_1^2}{144{\mathrm{\Omega }}^4}-\frac{35s_1{s}_2}{96{\mathrm{\Omega }}^4}+\frac{23s_2^2}{1024{\mathrm{\Omega }}^4}-\frac{2s_1}{3{\mathrm{\Omega }}^2}+\frac{3s_2}{32{\mathrm{\Omega }}^2}\right)\cos (\mathrm{\Omega }\tau )+\frac{s_1\left(s_1-3s_2-3{\mathrm{\Omega }}^2\right)}{9{\mathrm{\Omega }}^4}\cos (2\mathrm{\Omega }\tau )+\\ \hspace{1em}\hspace{1em}\frac{8s_1^2+12s_1{s}_2-9s_2^2-36s_2{\mathrm{\Omega }}^2}{384{\mathrm{\Omega }}^4}\cos (3\mathrm{\Omega }\tau )+\frac{s_1{s}_2}{96 {\mathrm{\Omega }}^4}\cos (4\mathrm{\Omega }\tau )+\frac{s_2^2}{1024 {\mathrm{\Omega }}^4}\cos (5\mathrm{\Omega }\tau ).\end{array}$$

(69)

By the same token, the estimated circumscribed solution of the equation of motion that is given in Equation (57) may be formulated as follows:

$$\gamma (t)=\underset{\chi \to 1}{\lim }\left({\gamma }_0(t)+\chi {\gamma }_1(t)+{\chi }^2{\gamma }_2(t)\right),$$

(70)

Though there are alternative approaches to Equation (57) numerically or analytically (see [43,44,45,46,47,48]), the present method provides a simple but effective tool for the present analysis.

In what follows, a numerical calculation of the theoretic consequences that are reported in this section is completed. As stated before, the non-dimensional amounts are actually useful in shortening the obtained results. For this purpose, consider the same non-dimensional quantities given in the case presented in Section 4.

Therefore, Figure 22 graphs the perturbed solution as given in Equation (70). This figure characterizes the perturbed solution for a system, taking the subsequent specifics:

$$\rho =9, \mu =0.1, \epsilon =0.4, {\mu }^{\prime }=5, \tilde{U}=2, h=0.9, {\alpha }_1=0.5, \widehat{\kappa }=0.5, V=8, {\alpha }_2=2, {\alpha }_3=0.6, k=0.02 \mathrm{and} E_0=5.$$

For more convenience, the analytical perturbed solution of the surface wave is graphed along with the three-dimensions as shown in Figure 23.

6. Concluding Remarks

In accordance with the theory of the viscous potential, the behavior of nonlinear interfacial waves between two horizontal immiscible liquids in the occurrence of a uniform normal electric field is analyzed. A better quantitative and quantitative understanding of heat and mass transfer mechanisms behind various food production, processing, preservation, and storage. Therefore, the current paper includes the presence of mass and heat transfer. For straightforwardness, the shortened model of Hsieh [15,16,17] has been employed. Recently, Zanutto et al. [49] involved the interfacial mass and heat transfer in their work. Additionally, the coupling of the mass and heat transfer with viscoelastic Walters B fluid was addressed by Ghasemi et al. [50]. Their governing equations of motion were scrutinized by means of the homotopy analysis method. Many models explore fluids of both elastic and viscous characteristics. These visco-elastic fluids are becoming increasingly important in current technologies and industries. Therefore, the present study considered visco-elastic Walters’ B fluid. Pandey et al. [51] investigated the initiation of convection in a horizontal layer of Walter’ B visco-elastic nanofluid. Equation (57) can also be solved by other analytical methods [52,53,54,55,56], especially He’s frequency formulation [55,56]. Due to the practical implications of the porous media, the media are presumed to be fully saturated in porous structure. Furthermore, the two rigid barriers are considered as permeable, where the suction and injection velocities are taken into account. The nonlinear approach is concerned with the linear equations of motion with the appropriate nonlinear boundary conditions. This technique resulted in a nonlinear distinguishing differential equation, which judges the surface elevation. By means of the Routh–Hrutwitz technique, the stability criteria are obtained. Typically, a non-dimensional practice is utilized for a good presentation of the stability benchmarks. Several non-dimensional numbers, such as Weber, Ohnesorge, Darcy, Elasticity, Bond, and potential Bond numbers, are investigated. The multiple scale method with the aid of the Taylor expansion are employed to provide a Landau–Ginzburg equation. Therefore, the stability criteria are theoretically achieved and numerically confirmed. Additionally, the profile of the surface waves is obtained like a Rayleigh–Duffing equation. By resources from the L–P technique, an analytic approximate bounded solution is derived.

Finally, the influences of various non-dimensional numbers on the different cases can be summarized in the following table:

• In light of the linear approach

Physical Parameters	Behavior
Weber numbers $W{e}_{\tilde{U}}$	U
Ohnesorge number $Oh$	S
Darcy number $Da$	U
Elasticity number $El$	S
Bond number $Bo$	S
The ratio of velocity $V$	Dual role
Mass and heat transfer parameters ${\alpha }_1$	S
The ratio of viscosities $\mu$	S
The ratio between the two thicknesses $(h=h_2/h_1)$	S

• In light of the nonlinear approach

Physical Parameters	Behavior
Mass and heat transfer parameters ${\alpha }_1$	S
Mass and heat transfer parameters ${\alpha }_2$	S
Mass and heat transfer parameters ${\alpha }_3$	U
The ratio of Walters’ B visco-elasticity ${\mu }^{\prime }$	U
The ratio of velocity $V$	Dual role

Due to the impact of suspended particles on the interface stability of superposed visco-elastic fluids in porous media have importance in industrial and chemical engineering, as a future work, it is it is worthy to examine the nonlinear stability analysis of well-known models. The fact that information about the fluid–particle mixture is not commensurate with their industrial and scientific significance is another motivator for future work.