Stability and Bifurcation Analysis of Two-Immiscible Liquids Film Down an Inclined Slippery Solid Substrate

In this work, the dynamic behavior of linear and nonlinear waves propagating at the separating surface between two thin layers of viscous Newtonian fluids is studied in the presence of the effect of insoluble surface surfactant. The two liquids are confined between two infinite rigid parallel plates and assumed to have different densities and viscosities. The equations of evolution for surface-wave elevation and concentration of surfactant are derived using the lubrication approximation. In the linear stage, by utilizing the normal mode approach, we have derived the dispersion relation that relates the wave angular frequency to the wave number and other parameters that is solved numerically to inspect the influences of some selected parameters on the stability criteria of the fluid flow. Also, analytical expressions for the growth rate as well as its maximum value with corresponding wave number are obtained in the special case of long-wave limiting. It is concluded that the Marangoni number Ma\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ma}$$\end{document} has acquired a significant stabilizing influence on the fluid flow, whereas the inverse of the slippery length of substrate plate β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}, resorts to the destabilize the motion of the interfacial waves. Consequently, both of the Marangoni number and the substrate slippy coefficient can be utilized to control the film flow regime, where they preserve the film laminar flow and tend to prevent the film breakdown. These can be useful in many industrial applications such as coating processes, heat exchangers, cooling microelectronic devices, chemical reactors, food processing, thermal protection design of combustion chambers in rocket engines and operation of Laser cutting and heavy casting production processes.


Introduction
The phenomenon of falling films in a wide range of industrial and applications has received much attention from many authors, for example, Chang (1994), Oron et al. (1997), Weinstein and Ruschak (2004), Craster and Matar (2009), Zakaria (2012), Kalliadasis et al. (2013), Batchvarov et al. (2021) and others. The open-flow systems under the influence of various factors are encountered in many industrial fields, like coating processes as reported by Alekseenko et al. (1994), heat exchangers by Salvagnini and Taqueda (2004), cooling microelectronic devices by Squires and Quake (2005), chemical reactors by Bender et al. (2017), food processing and thermal protection design of combustion chambers in rocket engines by Shine and Nidhi (2018). Many practical experiences have been made continually by researchers to validate and complement the theoretical studies on falling liquid films, for examples, Liu et al. (1994), Alekseenko et al. (2005), Kharlamov et al. (2015), Adebayo et al. (2017), Charogiannis et al. (2017) and Charogiannis and Markides (2019).
Insoluble surfactants are surface active agents that have been shown to influence the surface tension of fluids Chang and Franses (1995). Some experimental measurements of Song et al. Song et al. (2006) for equilibrium surface tension have evidenced that the surface tension is reduced in response to the surfactant concentration in the fluid bulk. The influence of insoluble surfactants on a falling film was first studied theoretically by Benjamin (1964) and Whitaker (1967). They remarked that the increase of the critical Reynolds number when surface elasticity is included, concluding 1 3 30 Page 2 of 32 that there was a stabilizing effect of surfactants, which was qualitatively in agreement with many experimental observations (Stirba and Hurt 1955;Tailby and Portalski 1961;Blyth and Pozrikidis 2004 and others) identified a new Marangoni mode when surfactants are present where a stabilizing effect of surfactants is remarked. According the perturbation theory, Oron and Edwards (1993) derived a nonlinear evolution equation for the free-surface displacement of the liquid film in the presence of interfacial viscous stress. Pereira and Kalliadasis (2008) studied the problem of moved liquid film in both linear and nonlinear regimes. In linear instability, they solved the Orr-Sommerfeld eigenvalue problem and obtained three different modes, where they stated that the Marangoni stress caused by surfactants reduces the domain of film instability. Using the weighted residual model (WRM), they found that the amplitude and velocity of travelling waves on the film are decreased by the Marangoni effect. The linear stability Couette-Poiseuille flow with insoluble surfactants is studied by Frenkel and Halpern (2002) in zero gravity environment and obtained asymptotic long-wave expressions for the growth rates. They found a caused destabilization by the interfacial surfactant. Halpern and Frenkel (2003) studied linear instability of creeping Couette-Poiseuille flow for a two-layer system with a monolayer of an insoluble surfactant in the Stokes approximation. They demonstrated that the unstable waves are not necessarily much longer than the smaller of the two layer thicknesses, there are parametric regimes for which the instability has a mid-wave character and the flow being stable at both sufficiently large and small wavelengths and unstable in between. Frenkel and Halpern (2017) investigated the Couette flow of a horizontal channel of two-layer system subject to vertical gravitational forces and with insoluble interfacial surfactants. They introduced linear and nonlinear instability of the inertialess long-wave with surfactant and gravity dependent using the lubrication approximation which leads to coupled nonlinear evolution equations. In their study, the effect of density stratification included in the system to examine the interacting effects of gravity and Marangoni forces. Linear stability regime of an immiscible bounded two-layer plane Couette flow film in the presence of an insoluble surfactant monolayer along the interface and the boundary parallel plates moving at a constant relative velocity to each other is produced by Frenkel et al. (2019). The impact of insoluble surfactant on the instability of interfacial waves between a thin film and a much thicker fluid has been analyzed in the studies of Bassom et al. (2010), Kalogirou et al. (2012), Kalogirou and Papageorgiou (2016) and Kalogirou (2018) through analysis and numerical computations of a system of evolution equations that includes a non-local contribution from the thicker fluid. The effect of slippery boundary condition of substrate on the dynamics of overlying liquid film flow is discussed by many anthers, for example, Pascal (1999), Samanta et al. (2011), Zakaria (2012, Bhat and Samanta (2018), Sani et al. (2020) and Samanta (2020). Pascal investigated the linear stability of a fluid flowing down an inclined permeable plane by studying the evolution in time of infinitesimal disturbances of long wavelength. He assumed a flow through the porous medium which is governed by Darcy's law, with producing critical conditions for the onset of instability. Also, Pascal remarked that increasing the permeability of the inclined plane destabilizes the flow of the fluid layer flowing above. Samanta et al. (2011) considered the gravity-driven film flow on a slippery inclined plane within the framework of long-wave and boundary layer approximations. They observed the non-trivial effect for Navier slip condition at the wall on the primary instability. Also, the observed role of Navier slip condition on the behavior of backflow phenomenon in the capillary region of the solitary waves. Zakaria studied the linear and nonlinear instability for the evolution of superposed layers of fluid flowing down inside an inclined permeable channel using the Kármán-Pohlhausen method through a long-wave approximation. Bhat & Samanta studied The linear stability analysis of a fluid flow down a slippery inclined plane when the free surface of the fluid is contaminated by a monolayer of insoluble surfactant. According to their analysis, wall slip exhibits a stabilizing effect on the surface mode as opposed to the result in the long-wave regime, while the insoluble surfactant exhibits a stabilizing effect on the surface mode as in the result of the long-wave regime at moderate Reynolds number. At high Reynolds number, both wall slip and insoluble surfactant exhibit a stabilizing effect. The both surface and shear modes compete with each other to dominate the primary instability in case of sufficiently small inclination angle. Sani et al. analyzed the linear instability of a surfactant-laden two-layer falling film over an inclined slippery wall under the influence of external shear stress the top free surface. They found that the overall stabilization of the surface mode by wall velocity slip for the stratified two-fluid flow is contrary to that of the single fluid case in the presence of strong imposed shear. Samanta studied the temporal and spatial disturbance growths for three-dimensional viscous incompressible fluid flows with slippery walls. He investigated that for a free surface flow over a slippery inclined substrate, the maximum temporal energy amplification increases in severity with the effect of wall slip for the spanwise perturbation, but it attenuates with the wall slip for both streamwise and spanwise wavenumbers. The aim of present manuscript is to study the linear and nonlinear instability with the bifurcation view for an overlying two-layer plane Couette flow with different densities, viscosities, thicknesses and insoluble interfacial surfactants on inclined fixed substrate which satisfies the slippery condition and bounded by an upper moved solid wall. The rest of this manuscript is organized as follows: Mathematical equations are posed in section "Mathematical Formulation". Lubrication approximation for the considered problem is discussed in section "Lubrication Approximation". The evolution equations are explored in section "Evolution Equations". Section "Linear Instability Analysis" is devoted for linear instability thorough the growth rate and neutral curves is in investigated. In section "Weakly Nonlinear Analysis of Surface Resonated Waves", Weakly nonlinear analysis of surface resonated waves is browsed through the bifurcation view. The concluding remarks are presented in section "Conclusion".

Mathematical Formulation
Consider a two-dimensional system (x * , z * ) film of twoliquids with different densities, viscosities past a slippery inclined substrate. The film is bounded by a parallel plate which is moved at a constant velocity U 2 with U 2 = V * at a distance z * = Z * as shown in Fig. 1. The z * -axis is the spanwise, coordinate perpendicular to the two plates with the substrate located at z * = 0 where z * = d defines the location of the unperturbed liquid-liquid interface. The symbol * indicates a dimensional quantity.
At the interface, the surface tension, * depends on the concentration of the insoluble surfactant monolayer, * . The basic flow is driven by the combination of the steady motion of the the upper plate and the inclination of substrate.
In the base state, the interface is flat, and the surfactant concentration is uniform.
Once disturbed, the surface tension and surfactant concentration are no longer uniform and the deflection of the liquidliquid interface is represented by the function H * (x * , t * ) where t * represents the time.
The continuity equation and the Navier-Stokes equations govern the fluid film motion in the two layers (with j = 1 for the lower liquid layer and j = 2 for the upper liquid layer). They are where j is the density, v * j = {u * j , w * j } is the fluid velocity vector; p * j is the pressure, j is the viscosity; g is the gravity and is the inclination angle of the substrate.
The interfacial boundary conditions will be as follows. The velocity must be continuous: where [f ] 2 1 = f 2 − f 1 denotes the jump in f across the interface, z * = H * (x * , t * ).
We assume the dependence of surface tension on the surfactant concentration given by the Langmuir isotherm relation (Edwards et al. 1991), which becomes the linear Gibbs isotherm when the surfactant concentration * where 0 and 0 are the basic surface tension and surfactant and Ee is the elasticity parameter. Here we have considered dilute surfactant concentrations only, and therefore a linear equation of state is expected to be valid.
According to the gradient of surface tension, the interfacial balances of the tangential and normal stresses are given by, respectively, and where H = 1 + H * 2 x * 1∕2 . The kinematic interfacial condition is Surface concentration value of the insoluble surfactant on the interface, * [x * , t * ] , obeys the following transport equation Pozrikidis 2004 andHalpern 2017) where D s is the diffusivity constant for the interfacial surfactant. The surface tension gradient in the tangential stress balance Eq. (5) emerges due to local changes in the surfactant concentration and give rise to Marangoni forces.
At the slippery substrate, z * = 0 , the lower layer fluid velocity satisfies Navier-slip and no penetration boundary conditions (Pascal 1999;Samanta 2020) where * is the inverse dimensional sliplength. At the upper plate, we have The deriving of the evolution equation will be made wherein the sum of the flow rates of two-layers Q, is held constant where The governing equations with their boundary conditions are normalized by introducing the following dimensionless variables: (11) where U is the characteristic velocity scale and , d is the thickness of lower layer and is the typical wavelength.

Lubrication Approximation
The interest of this work is to follow the temporal evolution of a given perturbation at the interface, and in particular to explore the effects of surfactants, the slippery of the substrate and film thickness on the emerging linear and nonlinear developments on the interface. In what follows we assume that the wavelength of the disturbance is much larger than the film thickness, which suggests a parameter (= d∕ ) ≪ 1 . The evolution of long waves can be described efficiently working under the lubrication approximation. Under these conditions, the Navier-Stokes equation is greatly simplified by assuming that the fluids are in a state of nearly unidirectional motion (Blyth and Pozrikidis 2004;Bassom et al. 2010;Frenkel and Halpern 2017;Zakaria and Sirwah 2018;Frenkel et al. 2019 and others). The lubrication approximation was used by Halpern and Frenkel (2008), Kalogirou (2018) and Kalogirou and Blyth (2020) to build a model valid for interfacial deformations of the same order as the film thickness. Considering a weak inclination and Reynolds number Re, the inertial effects are negligible and thereby we can apply the concept of lubrication approximation theory. For apply this, we consider cot = S∕ and Reynolds number is assumed to be O( ) Mavromoustaki et al. (2010); Zakaria and Sirwah (2018). Substitution of the normalizations in (12) into the governing equations and boundary conditions, taking into acount the concept of lubrication theory, yields the following set of leading order equations: (12) (x, z, H) = (x * ∕ , z * ∕d, H * ∕d); t = (U∕ )t * ; (u j , w j ) = (u * j , (d∕ )w * j )∕U; = * ∕ 0 ; where Ma = Γ 0 Ee 1 U = −1M a , Ca = 1 U 0 = 3C a , Pe = Ud Ds =Pe are the Marangoni, capillary and Péclet numbers. The parameter scalings considered here are in order to retain gravity, surface tension and Marangoni contributions in the leadingorder dynamics and are similar to those applied in Kalogirou and Blyth (2020).

Evolution Equations
Using the leading-order normal stress balance, the solution of the pressure perturbation namely Using the solutions of the pressure and the boundary conditions, the momentum equations in (13) can be integrated in z twice to give where C(x, t) is the constant of the integration. Finally, the leading-order kinematic equation and leading-order transport equation for interfacial surfactant are given by where p 1 = −Sz + C(x, t), The system of evolution Eqs. (14), (15) and (16), where the tildes are dropped, are derived under the assumptions of lubrication approximation, and retain a number of salient physical properties: Gravitational forces are incorporated in the diffusion term ST 1 (H)H x and are responsible for making the flow susceptible to the classical Couette-Poiseuille instability that appears when a heavy fluid lies above a lighter fluid beside the difference between basic velocities. This scenario emerges only when > 1 , in which case the diffusion term destabilizes the waves motion.

Growth Rate
In this section we proceed to investigate the stability properties of the flow as a result to small disturbances of small amplitudes superimposed on the basic state solution as is a small deflection from the flat interface solution H = 1, and Γ p (x, +t) is a small deviation from the initial case of uniform surfactant. Substituting H, Γ into evolution Eqs. (14), (15), and then linearizing with respect to H p , Γ p to have where the coefficients C i are given in the Appendix.
We seek traveling wave solutions for Eqs. (17) and (18) in the form where over bar denotes the complex conjugate, k is the positive real wave number, = r + i i is the complex angular frequency while H ,Γ stand for the initial constant disturbances. The symbols i = √ −1. Here r is the growth wave while i ∕k represents the phase speed. Insertion these solutions into Eqs. (17) and (18) yields the dispersion relation where In view of the dispersion relation (20), what we're interested in here with the real part r of the angular frequency that governing the stability criteria of the linearized surface waves. It is clear that when the growth rate r is negative, the disturbance of wave elevation will oscillates periodically or decay with time, while the positive values of r admit s growing or a nonuniform disturbance for such waves. Each curve grows monotonically as the wave number increases until it reaches its maximum value (peak value). After reaching that value, the curve gradually begins to decrease, which is attributed to the stabilization effect that played by the surface tension. In Fig. 2, we have plotted the growth rate as a function of the wave number at various the velocity of the upper plate V. In our model calculations, we have considered the fixed values of the parameters Ca = 1000, Z = 1.3, = .1, = 0.8, = 0.5, q = −.005, S = 1, Pe = 600, Ma = 100. We observe that the growth rate as well as the cut-off wave number separating the stable region about the unstable one is monotonically increases with the increase in the values of V as pointed by Frenkel and Halpern (2017) and Frenkel et al. (2019). This enhances the destabilizing impact of the parameter V. Numerical simulations are plotted in Fig. 3 to study the variation of r with respect to the wave number for = 0.1, 1, 3, while the values of the remaining physical parameters are held constant, as given in caption of Fig. 3. It is to be noted that both the growth rate and the instability window are extending due to the increase in , which means that this parameter also tends to grow the disturbance of the interfacial waves profile. As for Fig. 4, we have displayed the evolution of the growth rate with the wave number as a result of the increase in the values of q. It is evident that the parameter q has acquired a tenuous destabilizing effect on stability behavior of the liquid films where the vales of the maximum growth rate and the cut-off wave number posses a slight increase owing to the increase in values of q. Unlike to the effects of V, q, the numerical computations displayed in Fig. 5 indicate that the Marangoni number Ma has a considerable stabilizing effect, where the cut-off wave number, growth rate and the stability window have a significant contract with respect to larger values of Ma. This clarifies that the stability property is very sensitive to the variation in the values of Ma. Such a result is in a well agreement with several of previous studies such as (Frenkel et al. 2019;Sirwah and Alkharashi 2021). In Fig. 6, we proceed to investigate the influence of the location of the upper plate, specified by the parameter Z, on the temporal growth rate versus the waves number. We distinctly note that the cut-off wave number kc, the maximum growth rate m and the corresponding wave number km will diminish with increasing Z. This confirms that the more the upper plate is repositioned, the more stable the system becomes where the growth rate as well as well as the instability window is constantly decreasing. Figure 7 illustrates the effect of the density ratio on stability picture of the linearized system where the temporal growth rate of the interfacial waves is plotted as a function of the Marangoni number at several values of for the sake of comparison. Demonstration of the graphs displayed in Fig. 7 show that as = 0.1 there are two modes of the flow pattern, one is unstable region specified by Ma < 70 where the growth rate has negative values and the other is the stable for Ma > 70. It is also noticed that when is increasing to the value to.5, there is a noticeable contraction in the stability region to become limited in range Ma > 150. With a further increase in the value of the density ratio = 0.7, the shrinkage continues further in the wave surface motion stability region separating the two liquid films to be characterized by Ma > 267. Accordingly, that is, the more dense the upper fluid is than the lower liquid, the less stable the movement of the interfacial surface separating them, and it also weakens the stabilizing role of the Marangoni number and these results, also, is similar to resluts of Frenkel and Halpern (2017), taking into account the difference in the hypotheses of the two problems.

Neutral Stability Curves
We now proceed to obtain the equations of the transition curves that separate the regions of linearizes stable flow from those of destabilizing flow. To do so, we substitute for by r + i i into Eq. (20) and then split the real and imaginary parts to have As we eliminate i from these two equations, a fourthorder polynomial equation of r will be yielded  The required marginal curves that separate the stable from unstable regions will be specified as we set r = 0, which in turn gives which can be rewritten in the form 2 1 3 + 1 2 4 − 2 4 = 0. (24) Equation (24) has in general two complex conjugate roots and a third real root which utilized to control stability criterion of the linearized flow. This real root is expressed as The wave number k max corresponding to the maximum growth rate rmax will be determined by simultaneously solving (22) and the equation designated by differentiating (22) with respect to k, taking into consideration that r k = 0 at the maximum growth rate. The last equation is expressed as The numerical results displayed in Figs. 8, 9, 10 and 11 represent the transition curves that separate regions of stability, specified by S, from regions of instability, specified by U, in the plane k − Ma In Fig. 8, the critical curves are plotted at three different values of the travel velocity of the upper plate. It is to be observed that there is a noticeable increase in the regions of instability, and in turn a large contraction in the regions of stability as a result of increasing the velocity values V, particularly at small values of the wavelength. Hence, these results reinforce the unstable role of the upper plate velocity, especially for long waves. As for Fig. 9, the neutral curves are computed for various values of parameter . . We also note that the variable tends to generate growing waves at the interface separating the two liquids, but this effect appears to be moderate and very weak for short waves when compared to the effect of the traveling velocity of the upper plate, as was observed in Fig. 8 and as it appeared in the results of Bhat and Samanta (2018) and Samanta (2020). The curves presented in Fig. 10 clearly show that there is a monotonic increasing in the long wave stabilization regions with increasing the values of Z , while this increase grows at a smaller rate for short waves. This means that the stability of the wave motion can be controlled by increasing the thickness of the fluid layers. Figure 9 clarifies that the parameter q has a similar destabilizing effect as well as that of the parameters Z, and the unstable role of q is much more than of that of and less than the instability influence of V.
To gain a deeper understanding of flow behavior, In order to better understand the flow behavior, Fig. 12 clarifies the structures of the dimensionless growth rate in the plane k − by changing the values of q corresponding to the sections of this graph. In Fig. 12a, where the value −0.2 is assigned for q, the zero-value contour corresponds to the neutral case that separate the stable regions specified by counters of negative values from the unstable regions represented by counters of positive values. When increasing the value of q to −0.005, represented in Fig. 12b, we notice the zero-value contour will be shifted up results in an increase in the k − region corresponding to the contours with positive values and a contraction of the domain associated with the counters with positive values. Further increase in values of q to 0.01 as performed in Fig. 12c causes more enlargement in the k − domain of unstable flow pattern and more shrinking of negative-values contours admitting stable flow. In order to examine the influence of variation of the parameter Z on the stability criteria of the linearized model, we have plotted the critical wave number against the viscosity ratio for several values of Z in Fig. 13. In Fig. 13a, where Z = 1.1, we have some contours with positive values and others of negative values and then according to the specified parameters two modes of an unstable and one stable system are observed. In Fig. 13b and c, we display graphically the numerical results of the same system investigated in Fig. 13a with increasing Z values to 1.2 and 1.4, respectively. An investigation of contours of the growth rate curves of Fig. 13 shows that the greater the thickness of the fluid layer represented by the increase in Z, the more significant shrinkage occurs in the regions of instability. It is expected that at large enough values for Z we will only have a stable flow behavior pattern. The effect of increasing the velocity of the upper plate in the k − Ma plane on the stability picture can easily be observed through the comparison Fig. 14a-c. It is worth noting by investigating these graphs that as V increases there will be an expansion in the unstable region along the Ma-axis and it will also increase the critical Marangoni number Ma c .

Limit of Long Waves
The study of special case of long waves can yield some analytical asymptotic formulas to calculate the wave growth rate as well as its maximum value. The solutions of Eq. (20) are then given by . Therefore, keeping the four leading terms in the series to have Therefore Then the approximated wave growth rate is given by where the expressions of H 1 and H 2 are defined in the Appendix.
To determine the wave number corresponding to the maximum growth rate k max , we must solve the equation r k = 0 . Then we have and thereby the maximum growth rate is given by or with and where the formula of H 3 and H 4 are given in the Appendix.
In Figs. 15, 16 and 17, we investigate the variation of maximum growth rmax and its corresponding wave number k max with respect to the Marangoni number Ma, Z and . The solid curves display such variations as calculated from the general Fig. 13 Investigates the dimensionless growth rate features in the plane k − for the same system of parameter' values used in Fig. 12 and q = −.005 when: a Z = 1.1, b Z = 1.2 Z = 1.4 Page 13 of 32 30  (20), while the dashed curves are corresponding to these evaluated from the long-wave approximation model represented by Eq. (32) and (33). We observe that there is a good agreement between the two models. The other parameter values are held constant, such as given in the caption of these figures. It is evident from Fig. 15 that there is an increase in both rmax and k rmax with an increase in Ma up to the value of Ma = 24.7. The opposite occurs as Ma > 24.7 where increasing Ma leads to a monotonic decrease in values of rmax , k max . Figure 16 clarifies the evolution in both rmax and k rmax with . We notice that the values of rmax and k m grow gradually with increasing which confirms the initializing role of . Figure 17 indicates that increasing the velocity of the upper plate in either direction will increase the values of rmax , k max and thus the growing waves are dominated.

Weakly Nonlinear Analysis of Surface Resonated Waves
In this section, we will focus our efforts to study the weakly nonlinear dynamical behavior of interacting surface waves. To accomplish this, we will resort to using Taylor's expansion about the unperturbed flat interface h = 1 and initial uniform concentration Γ = 1 to obtain the evolution equations while keeping the terms in h and Γ to the second degree, and then solve the resulting system of equations using the multiple scales method, as will be explained later. Such weakly nonlinear equations read  No, we proceed to utilize the method of multiple scales to obtain second-order convergent solutions of Eqs. (37) and (38) describing the wave motion of small but finite amplitude. To do so we introduce the new slow variables t 1 and x 1 as well the two variables t 0 and x 0 that correspond to rapid variations in such away that Nayfeh and Balachandran (2004) where is a small dimensionless parameters expressing the wave steepness ratio. In this analysis the unknown functions of wave interfacial elevation and surfactant concentration are expressed in powers of in the form Fig. 16 The variation of rmax and k max with Z according to the general growth rate (20), and the long-wave approximation model given by Eqs. (30)   The coefficients i , i are very lengthy and not included here and are ready on request from the authors.

The Problem of First-Order ( O())
In what follows, we proceed to discuss the dynamical behavior of the weak;y nonlinear surface waves owing to the interacting between two-to-one internal waves. Such case of resonance has been demonstrated in several works, for example, Henderson and Miles (1991), Zakaria (2008) and Sirwah (2014). Henderson & Miles studied both theoretically and experimentally the behavior of surface waves arising from resonance between internal waves that are preserved against the effect of the weakly damped viscosity inside a closed basin subjected to periodic vertical vibrational motion. Zakaria (2008) studied Faraday waves with three modes, on the interface between superposed liquids in an infinite boxed basin subjected to vertical excitation. Sirwah (2014) considered the chaotic behavior of the three-dimensional resonated standing surface waves propagated at the interface of a weakly viscous, incompressible magnetic liquid within a rectangular basin as the system is subjected to a normal time-periodic magnetic field as well as external vertical excitation. The solutions of the first order problem are sought in the form (43) Here, A 1 x 1 , t 1 , A 2 x 1 , t 1 are two complex unknown functions to be governed by the solvability conditions of the second order problem and 1 , 2 are the roots of the dispersion relation (20). The next step is to insert the first order solutions in the second order equation. Then, the resulting equations will posses uniform convergent solutions are obtained when the probable secular terms are omitted. Such process will result the so called the solvability conditions. We must distinguish between two basic cases: the first one is nonresonance case in which there is no nearness between the two frequencies of the internal waves, while the second one occurs when the frequency 2 approaches the frequency 1 .
In the nonresonance case, the solvability condition read We are interested in the resonance case specified by 2 ∼ 2 1 . To express this nearness, we introduce the detuning parameter according to Such case of resonance will admit the following solvability conditions where the coefficients J ij , K ij are given in the Appendix.
It is convenient to use the transformation to remove the step mismatch of the non-linear terms. Consequently, the solvability conditions (50) and (51) take the form (46) Γ 1 = fA 1 x 1 , t 1 e i(k 1 x 0 −t 0 1) + gA 2 x 1 , t 1 e i(k 2 x 0 −t 0 2) , (49) 2 = 2 1 + . (50) The solution to the modal dynamical system of these equations can be assumed to be in the form The two amplitudes 1 , 2 are, therefore, related via the equations In order to determine the possible steady-state solutions of the dynamic system represented by the Eqs. (56) and (57), it is very appropriate to make use of the polar form where = t 1 ; and i , i stand for the amplitudes and phases of the two modes, which are presumed to be real. Substituting (58) into Eqs. (56) and (57) and then separating the real and imaginary parts, we get where (54) The steady-state solution of the dynamical system (59)- (61) are the values of a 1 , a 2 , 0 that correspond to the conditions i = = 0. Therefore, such solution can be calculated by solving the following equations Thus, we have the possible solutions or with B 31 = 1 K 2 11 + K 2 12 J 2 11 + J 2 12 K 2 11 (J 12 J 12 − 2 J 14 − 2 J 13 J 11 + J 2 11 )+ K 12 K 12 −2 J 13 J 11 − 2 J 12 J 14 + J 2 11 + J 2 12 + 2 K 14 J 2 11 + J 2 12 + 2 K 13 K 11 J 2 11 + J 2 12 , = 2 1 − 2 .
(68) (72) (J 12 −2J 15 p 2 + 2J 16 q 2 +2 J 14 + J 12 + 2J 11 (J 16 p 2 + J 15 q 2 + L 13 ) + L 2 11 )p 1 − 2 L 12 J 16 p 2 + J 15 q 2 − J 13 +J 11 J 15 p 2 − J 16 q 2 + J 14 q 1 , The first part of this section's computations aims to portray graphical representation for the evolutions of the assigned fixed points or the equilibria of the autonomous dynamical system controlling the modified amplitudes of the interacting resonance waves of the weakly nonlinear system in the presence of insoluble interfacial surfactant. Secondly, we proceed to discuss several two-dimensional projections of the phase portraits of the modulation dynamical system onto some planes as well as the time histories of the involved parameters. The linear stability of the assigned equilibria will be demonstrated via the eigenvalues of the Jacobian of the dynamical system governing the perturbed quantities superimposed on the fixed points. Figure 18 shows the effect of the change in the variable Z (78) K 12 K 15 − K 11 K 16 p 2 1 − 2 K 12 K 14 + K 12 + K 13 K 11 + K 2 11 p 2 −2 K 11 K 15 + K 12 K 16 p 1 q 1 + K 11 K 16 − K 12 K 15 q 2 1 +2 K 11 K 14 − K 12 K 13 q 2 .

Fig. 19
Bifurcation diagram in the plane − Z for the same system investigated in Fig. 18 on the behavior of the specified fixed points. These points are drawn as functions of the variable sigma corresponding to the values Z = 1.2, 1.4, 1.5 in the parts Fig. 18a-c, respectively. In our model calculations, the values of other parameters are given in the caption of figure. Having noted the curves displayed in Fig. 18a, it is to be noted that the dynamical system has always three pairs of fixed points whatever the value of sigma, the zero fixed point which is stable in the region specified by the relation > −0.39 and is unstable as < −0.39. The other two nontrivial fixed points (a 1 , a 2 ) are found to be stable as < −0.39 and > 0.61, while these equilibria become unstable as −0.39 < < 0.61. Note that the stable equilibrium is presented by solid line whereas unstable one presented by dashed line. We observe that at the point = −0.39, the three equilibria exchange stability which in turn suggests that = −0.39 represents a transcritical bifurcation point. Nevertheless, there are several scientific positions in which this type of bifurcation may occur, for example, in the case of logistic equation and models for the growth of a single species where there exists a fixed point at zero population despite the value of the growth rate. At the point = −0.61, the characteristic equation admits two pairs of imaginary conjugate eigenvalues. Also, it is found that as decreases below the value 0.61, the equilibria lose their stabilities with pairs of complex conjugate eigen values of positive real parts. Hence, the system is subject to Hopf bifurcation at = 0.61 and thereby results in creating a limit-cycle solution near the point = 0.61 of amplitude proportional to = −0.61 as it will be enhanced by the numerical results illustrated later in the rest of figures. Limit cycle solutions appear in such systems that fluctuate even if there are no periodic forces exerted on the system or equivalently the self-sustained oscillating systems. These solutions are of increasing scientific importance because of their many applications in various fields. Among them, we mention the systems of heart beat, spontaneous oscillating chemical reactions, dangerous self-vibrations in bridges and wings of aircraft. The conditions for Hopf bifurcation, are given in Lee and Mei (1996) and Usha and Uma (2000) by 1 > 0, 2 > 0 and 1 2 = 3 .

(a)
An investigation of Fig. 18b clearly shows that the qualitative behavior of equilibria to variation of the detuning parameter are similar to those of Fig. 18a. However, the increase in Z leads to the displacement of the Hopf bifurcation point to the left to occur at = 0.05, while the transcritical bifurcation point shifts to the right to become at = −0.3, which indicates that there is a contraction of the sigma range corresponding to the unstable solutions of the system. Consequently, such behavior enhances the stabilizing role of Z, which is consistent with the results of the linear case studied above. As expected, another increase in Z will lead to a further contraction of the unstable fixed solution regions without a qualitative change in the system behavior as shown in Fig. 18c. The hopf bifurcation point has occurred at = −0.004, and the transcritical bifurcation point came at = −0.17. Fig. 19 shows the points of both transcritical (TC) and Hopf bifurcations (HF) in the (Z − ) plane.
The influence of variation of the parameter on the evolution of the steady-state solutions of system (61)-(63) with the detuning parameter is shown in Fig. 20. Figure 20a is drawn as = 0.1, while Fig. 20b and c are plotted at = .3 and 0.7, respectively. The frequency response curves of Fig. 20a show that there are three possible fixed point solution. The trivial steady-state solution is found to be stable (sink) as > −0.07 and unstable (saddle) when < −0.07. As for the nontrivial stationary solutions, they are unstable in the region specified by −0.07 < < 0.06 while otherwise are stable solutions. Moving to Fig. 20b-c, the numerical calculations show that there are no essential differences regarding the qualitative behavior of these solutions, except that there is an expansion in the region of stability of trivial solutions, and on the contrary, there is a contraction in the domain of stability of non-trivial steady state solutions. These results reinforce the dual role of stability that the variable plays within the nonlinear system, in contrast to the unstable role of the Beta, which was previously observed in the framework of linear theory. It is shown that the system exhibits transcritical bifurcation at = −0.07 and Hopf bifurcation at = 0.06 corresponding to the case of = 0.3 ; while these types of bifurcations are occurred at the points = −0.2 and = 0.22 when = 0.7. Then an increase in the values will cause the points of the two types of bifurcations to move away from each other, as shown in Fig. 21.
The numerical results presented through Fig. 22 were accomplished by solving the modified equations for the dynamical system using the fourth-order Runge Kutta fourth -order method. In these figures, the phase portraits have been drawn in the planes p 1 − q 1 , p 2 − q 2 , as well as time history of the modified amplitudes p 1 , p 2 , q 1 , q 2 are demonstrated. The graphics of Two-dimensional projections of the phase portraits onto the p 2 − q 2 , plane and the longtime histories of p 2 , q 2 , for the same system studied in Fig. 24 (a)

(b)
are presented by stable spiral and the period of oscillation will decrease as the time increases. Such a spiral or helical trajectory may occur as the harmonic oscillator were slightly damped. In this case, the path fails to close, because the oscillator loses very little energy per cycle. The path of the helical tip is an important characteristic distinctive of a spiral wave. It can be round in shape, among other things. The spiral is rigid as the path is circular. Otherwise it slaloms. In exciting systems, the helical meander caused by a bifurcation that resembles a Hopf bifurcation. The spiral paths can be modified by external influence fashions (Zykov and Mtiller 1996;Hu et al. 2010).
In Fig. 23, we have selected the parametric values, which are not satisfying the stability conditions (72). The graphs of As the detuning parameter is increased to − 0.44, the system undergoes a quadrable bifurcation period demonstrated. Additional increases in are likely to lead to a series of a more doubling period solutions with some noise. Such results are in a well agreements of those of Arafat et al. (1998), Sirwah (2014. When has the value −0.4, the system behavior is no longer periodic and it exhibits some noise. The phase portraits evaluated at = −0.33 show that the system behaves in a chaotic manner where chaotic attractor will merge. Moreover, the time histories expresses an attractor that is expected to grow and deform as is increased further. Consequently, the system behavior is no longer periodic or quasi-periodic as a result of increasing the values of the variable to −0.33 and then an irregular movement of the surface waves at the surface separating the two fluid layers.

Conclusion
In this study, the behavior of surface waves propagating at the separating interface between two layers of viscous Newtonian fluids has been studied taking into account the influence (c) Fig. 25 (continued) of insoluble surface surfactant. The mathematical equations describing the problem under study are communicated with the Navier-Stokes equations, continuity equation for incompressible fluids as well as the evolution equation for surface surfactant. We have reformulated the equations of motion and the associated boundary conditions in the dimensional form. We have applied long-wave technique for obtaining analytical solutions of the mathematical system describing for the physical model.
The two fluid layers are confined between a non-permeable top plate and a permeable one located at the bottom. The analysis is performed under the slippery state experiments that designate the condition necessary to determine the tangential component of the velocity at the interface between the lower fluid layer and the permeable plate. Consequently, we have derived two nonlinear coupled evolution equations governing the wave elevation at the free surface between the two films and the concentration of the surface surfactant. The multiple scales method is exploited to specify the solvability conditions necessary for obtaining asymptotic analytical solutions. In the frame work of the linear analysis, the technique of the normal mode is exploited to obtain the dispersion relation that relates the wave growth rate to the wave number and other parameters. Thereby the stability criteria of the linear system via numerical simulations are developed to demonstrate the effect of some parameters on the stability picture. Moreover, we have given a focusing to study the special case of small wave number limit which allows to obtain an analytical formula for the wave growth rate. The results of numerical simulation of the linear problem show that the dimensionless growth rate and the peak of wave number increase small with increasing the parameter q, indicating that q has a weak destabilizing influence on the linear stability of the flow. It is also found that the growth rate as well as the peak wave number increases significantly, revealing that theses parameters further enhance the surface instability. In contrast, it has been shown that the instability region, maximum growth rate and the corresponding wave number decreases pointedly, manifesting that such parameters suppress the flow instability significantly. In the nonlinear stage, we have devoted our attention to discuss the instability analysis of weakly nonlinear interfacial resonant waves. Such waves are structured by the interaction between the first and second modes or the case of two-to-one internal resonance. We have applied the modulation concept to the non-secular conditions of the resonate nonlinear waves and then derived a third-order autonomous dynamic system that governs the modified amplitudes and phases of the interacting waves. The steady-state solutions of this system are thereby obtained and their stability analysis is discussed. Numerical calculations are achieved by means of Runge-Kutta fourth-order method and their results are graphically illustrated to reveal the response of the system behavior to variation in some selected parameters. These influences are investigated through some diagrams include the frequency-response curves, phase portraits on some planes and time histories of the modulated amplitudes. Numerical results have also shown that the system undergoes two types of bifurcations with the detuning parameter as a control parameter: the transcritical and Hopf bifurcations. It worth noting that increasing the parameter Z leads to contracting the instability region of nontrivial stationary solutions while increase the instability region of the trivial fixed point solutions. In other words, the stabilizing or destabilizing role of Z on the fluid flow depends on the type of the assigned steady-state solution. Whereas, we have concluded that the increase in values causes an enlargement in instability region of the nontrivial fixed points and shrinks this of instability of the zero fixed point solutions. Moreover, it is remarked that system response does not qualitatively depend on the variation in Z or . The projections of the modulation equations as well as the time histories of modulated amplitudes in the vicinity of the transcritical bifurcation points are illustrated which explained that the behavior of the system varies between the centered stable and unstable spiral depending on the type of the stationary solution, whether it is stable or unstable. In the vicinity of the Hopf bifurcation points, the numerical computations admit periodic solutions and thereby a limit cycle is born. the limit cycle develops, distorts, and experiences reduplicated period-doubling bifurcations and in turn transmitted to a chaotic behavior.