Air-Aided Shear on a Thin Film Subjected to a Transverse Magnetic Field of Constant Strength: Stability and Dynamics

Abstract

The effect of air shear on the hydromagnetic instability is studied through (i) linear stability, (ii) weakly nonlinear theory, (iii) sideband stability of the filtered wave, and (iv) numerical integration of the nonlinear equation. Additionally, a discussion on the equilibria of a truncated bimodal dynamical system is performed. While the linear and weakly nonlinear analyses demonstrate the stabilizing (destabilizing) tendency of the uphill (downhill) shear, the numerics confirm the stability predictions. They show that (a) the downhill shear destabilizes the flow, (b) the time taken for the amplitudes corresponding to the uphill shear to be dominated by the one corresponding to the zero shear increases with magnetic fields strength, and (c) among the uphill shear-induced flows, it takes a long time for the wave amplitude corresponding to small shear values to become smaller than the one corresponding to large shear values when the magnetic field intensity increases. Simulations show that the streamwise and transverse velocities increase when the downhill shear acts in favor of inertial force to destabilize the flow mechanism. However, the uphill shear acts oppositely. It supports the hydrostatic pressure and magnetic field in enhancing films stability. Consequently, reduced constant flow rates and uniform velocities are observed.

1. Introduction

A nonlinear fourth-order degenerate parabolic differential equation of the form: ℎ_𝑡+𝒜(ℎ)ℎ_𝑥+𝜕ℬ(ℎ)ℎ𝜕𝑥^𝑥+𝒞(ℎ)ℎ_𝑥𝑥𝑥=0, (1) where 𝒜, ℬ, and 𝒞 are arbitrary continuous functions of the interfacial thickness ℎ(𝑥,𝑡), represents a scalar conservation law associated to the flow of a thin viscous layer on an incline under different conditions. The stability and dynamics of equations of the type of (1) is a subject of major interest [1–18] because of their robustness in regimes where viscosity dominates inertia [19]. Such studies have focused attention primarily on the isothermal and nonisothermal instability analysis, mainly for nonconducting fluids.

Since the investigation of Chandrasekar [20] on the stability of a flow between coaxial rotating cylinders in the presence of a magnetic field held in the axial direction, the laminar flow of an electrically conducting fluid under the presence of a magnetic field has been studied extensively. For instance, Stuart [21] has reported on the stability of a pressure flow between parallel plates under the application of a parallel magnetic field. Among other earlier investigations, Lock [22] examined the stability when the magnetic field is applied perpendicular to the flow direction and to the boundary planes. Hsieh [23] found that the magnetic field stabilizes the flow through Hartmann number when the electrically conducting fluid is exposed to a transverse magnetic field, provided that the surface-tension effects are negligible in a horizontal film. Ladikov [24] studied a problem in the presence of longitudinal and transverse magnetic fields. The author observed that the longitudinal magnetic field plays a stabilizing role and that the effect of instability at small wave numbers could be removed if the longitudinal magnetic field satisfies certain conditions. Lu and Sarma [25] investigated the transverse effects of the magnetic field in magnetohydrodynamic gravity-capillary waves. The flow of an electrically conducting fluid over a horizontal plane in the presence of tangential electric and magnetic fields was reported by Gordeev and Murzenko [26]. They found that the flow suffers from instability not due to the Reynolds number but due to the strength of the external electric field.

In applications such as magnetic-field-controlled material-processing systems, aeronautics, plasma engineering, MEMS technology, and magnetorheological lubrication technologies, the hydromagnetic effects are important. Also, liquid metal film flows are used to protect the solid structures from thermonuclear plasma in magnetic confinement fusion reactors, and this application requires a better understanding of the instability mechanism arising in a thin magnetohydrodynamic flow over planar substrates [27]. Furthermore, the presence of an external magnetic field regulates the thickness of a coating film. It prevents any form of direct electrical or mechanical contact with the fluid thereby reducing the risk of contamination [28]. Renardy and Sun [29] pointed out that the magnetic fluid is effective in controlling the flow of ordinary fluids and in reducing hydraulic resistance. The reason is that the magnetic fluids can be easily controlled with external magnetic fields and that coating streamlined bodies with a layer of less-viscous magnetic fluid significantly reduces the shear stress in flow boundaries. Magnetohydrodynamic flow can be a viable option for transporting weakly conducting fluids in microscale systems such as flows inside a micro-channel network of a lab-on-a-chip device [30, 31]. In all of the above applications, considering the associated stability problem is important because it gives guidance in choosing the flow parameters for practical purposes. In this regard, in the past two decades, the emerging studies on the hydromagnetic effects have focused their attention on stability problems and on analyzing the flow characteristics [28, 29, 32–43].

The shearing effect of the surrounding air on the fluid in realistic situations induces stress tangentially on the interfacial surface. The hydrodynamic instability in thin films in the presence of an external air stream attracted attention in the mid 1960s, which led Craik [44] to conduct laboratory research and study theory. He found that the instability occurs regardless of the magnitude of the air stream when the film is sufficiently thin. Tuck and Vanden-Broeck [45] reported on the effect of air stream in industrial applications of thin films of infinite extent in coating technology. Sheintuch and Dukler [46] did phase plane and bifurcation analyses of thin wavy films subjected to shear from countercurrent gas flows and found satisfactory agreement of their results connected to the wave velocity along the flooding curve with the experimental results of Zabaras [47]. Although the experimental results related to the substrate thickness did not match exactly with the theoretical predictions, still the results gave qualitative information about the model. Thin liquid layer supported by steady air-flow surface-traction was reported by King and Tuck [48]. Their study models the surface-traction-supported fluid drops observed on the windscreen of a moving car on a rainy day. Incorporating the superficial shear offered by the air, Pascal [49] studied a problem which models the mechanism of wind-aided spreading of oil on the sea and found quantitative information regarding the maximum upwind spread of the gravity current. Wilson and Duffy [50] studied the steady unidirectional flow of a thin rivulet on a vertical substrate subjected to a prescribed uniform longitudinal shear stress on the free surface. They categorized the possible flow patterns and found that the direction of the prescribed shear stress affects the velocity in the entire rivulet. The generation of roll waves on the free boundary of a non-Newtonian liquid was numerically assessed using a finite volume method by Pascal and D’Alessio [51], revealing the significant effect of air shear on the evolution of the flow. Their study showed that the instability criteria was conditional and depends on the directionally induced shear. Kalpathy et al. [52] investigated an idealized model suitable for lithographic printing by examining the shear-induced suppression of two stratified thin liquid films confined between parallel plates taking into account the van der Waals force. A film thickness equation for the liquid-liquid interface was derived in their study using lubrication approximation. They found that the effect of shear affects the imaginary part of the growth rate, indicating the existence of traveling waves. Furthermore, they also observed a critical shear rate value beyond which the rupture mechanism could be suppressed. This study motivated Davis et al. [53] to consider the effect of unidirectional air shear on a single fluid layer. For a two-dimensional ultrathin liquid, Davis et al. [53] showed that the rupture mechanism induced by the London van der Waals force could be suppressed when the magnitude of the wind shear exceeds a critical value, as observed by Kalpathy et al. [52]. Recently, Uma [18] measured numerically the profound effect of unidirectional wind stress on the stability of a condensate/evaporating power-law liquid flowing down an incline.

In the present investigation, the effects of downhill and uphill air shear on a thin falling film in the presence of a transverse magnetic field are studied. Such an investigation will illustrate the realistic influence of the natural environment acting upon the flows. Or, it may illustrate the need to control the flow mechanism through artificial techniques, which blow air when the hydromagnetic effects are considered. The outline of this paper is organized as follows. Section 1 presents the introduction. In Section 2, mathematical equations governing the physical problem are presented. Section 3 discusses the long-wave Benney-type equation. In Section 4, linear stability analysis, weakly nonlinear stability analysis, and the instability arising due to sideband disturbances are analyzed. The equilibria of a truncated bimodal dynamical system is mathematically presented in Section 5. While the results in Section 6 discuss nonlinear simulations of the film thickness evolution, Section 7 highlights the main conclusions of the study and includes future perspectives.

2. Problem Description

A thin Newtonian liquid layer of an infinite extent falling freely over a plane under the influence of gravitational acceleration, 𝑔, is considered. The flow is oriented towards the 𝑥-axis, and the plane makes an angle 𝛽 with the horizon. Properties of the fluid like density (𝜌), viscosity (𝜇), and surface-tension (𝜎) are constants. The magnetic flux density is defined by the vector 𝐁=𝐵₀𝑗, where 𝐵₀ is the magnitude of the magnetic field imposed along the 𝑦-direction (Figure 1). It is assumed that there is no exchange of heat between the liquid and the surrounding air, but an air flow (either in the uphill or in the downhill direction) induces a constant stress of magnitude 𝜏_𝑎 on the interface and moves tangentially along the surface. The 𝑦-axis is perpendicular to the planar substrate such that, at any instant of time 𝑡, ℎ(𝑥,𝑡) measures the film thickness.

Figure 1 

Sketch of the flow configuration. When 𝜏𝑎>0, the air shears the surface along the downhill direction (in the direction of the arrow). However, if 𝜏𝑎<0, the air shears the surface along the uphill direction (in this case, the arrow would point to the opposite direction). The air shear effect is zero when 𝜏𝑎=0.

The magnetohydrodynamic phenomena can be modeled by the following equations, which express the momentum and mass balance: 𝜌𝜕𝐔𝜕𝑡+𝐔⋅∇𝐔=𝜌𝐆−∇𝑝+𝜇∇²𝐔+𝐅,∇⋅𝐔=0. (2)(3) The last term in (2) arises due to the contribution of Lorenz body force, based on Maxwell’s generalized electromagnetic field equations [28, 54, 55].

In a realistic situation corresponding to a three-dimensional flow, the total current flow can be defined using Ohm’s law as follows: 𝐉=Σ(𝐄+𝐔×𝐁), (4) where 𝐉, 𝐄, 𝐔, and Σ represent the current density, electric field, velocity vector, and electrical conductivity, respectively. The Lorenz force acting on the liquid is defined as 𝐅=𝐉×𝐁. In the rest of the analysis, a short circuited system corresponding to a two-dimensional problem is considered by assuming that 𝐄=0. This assumption simplifies the last term corresponding to the electromagnetic contribution in (2). In this case, the pondermotive force acting on the flow (the last term in (2)) has only one nonvanishing term in the 𝑥-direction; therefore, 𝐅=𝐹_𝑥,𝐹_𝑦=−Σ𝐵²₀𝑢,0, (5) where 𝐆=(𝑔cos𝛽,−𝑔sin𝛽) and 𝑝 is the pressure.

Boundary conditions on the planar surface and the interface are added to complete the problem definition of (2)-(3). On the solid substrate, the no-slip and the no-penetration conditions are imposed, which read as 𝐔=(𝑢,𝑣)=0. (6) The jump in the normal component of the surface-traction across the interface is balanced by the capillary pressure (product of the mean surface-tension coefficient and the local curvature of the interface), which is expressed as 𝑝_𝑎−𝑝+𝜇∇𝐔+∇𝐔^𝑇(ℎ)(),⋅𝐧⋅𝐧=2𝜎Πon𝑦=ℎ𝑥,𝑡 (7) where Π(ℎ)=−(1/2)∇⋅𝐧 is the mean film curvature and 𝑝_𝑎 is the pressure afforded by the surrounding air. The unit outward normal vector at any point on the free surface is 𝐧=∇(𝑦−ℎ)/|∇(𝑦−ℎ)|, and 𝐭 represents the unit vector along the tangential direction at that point such that 𝐧⋅𝐭=0. The tangential component of the surface-traction is influenced by the air stress 𝜏_𝑎 and reads as 𝜇∇𝐔+∇𝐔^𝑇⋅𝐭⋅𝐧=𝜏_𝑎on𝑦=ℎ(𝑥,𝑡). (8) The location of the interface can be tracked through the following kinematic condition: 𝜕ℎ𝜕𝑡=𝑣−𝑢𝜕ℎ().𝜕𝑥on𝑦=ℎ𝑥,𝑡 (9)

In order to remove the units associated with the model (2)–(9) involving the physical variables, reference scales must be prescribed. In principle, one can nondimensionalize the system based on the nature of the problem by choosing one of the following scales [9, 16]: (a) kinematic viscosity based scales, (b) gravitational acceleration as the flow agent based scales, (c) mean surface-tension based scales, and (d) Marangoni effect as the flow agent based scales. However, for very thin falling films, the main characteristic time is the viscous one [7–9, 16]. An advantage of choosing the viscous scale is that either both the large and the small inclination angles could be considered by maintaining the sine of the angle and the Galileo number as separate entities [7, 16] or the Reynolds number as a product of Galileo number and sine of the inclination angle could be defined as a single entity [9]. Also, such a scale plays a neutral role while comparing the action of gravity and Marangoni effect in nonisothermal problems [16]. Choosing the viscous scale, the horizontal distance is scaled by 𝑙, vertical distance by ℎ₀, streamwise velocity by 𝜈/ℎ₀, transverse velocity by 𝜈/𝑙, pressure by 𝜌𝜈²/ℎ²₀, time by 𝑙ℎ₀/𝜈, and, finally, the shear stress offered by the wind by 𝜌𝜈²/ℎ²₀. In addition, the slenderness parameter (𝜀=ℎ₀/𝑙) is considered small, and a gradient expansion of the dependent variables is done [7–9]. The horizontal length scale, 𝑙, is chosen such that 𝑙=𝜆, where 𝜆 is a typical wavelength larger than the film thickness.

The dimensionless system is presented in Appendix A (the same symbols have been used to avoid new notations). The set of nondimensional parameters arising during nondimensionalization procedure are Re=Gasin𝛽 (the Reynolds number) [9], Ga=𝑔ℎ³₀/𝜈² (the Galileo number), Ha=Σ𝐵²₀ℎ²₀/𝜇 (the Hartmann number, which measures the relative importance of the drag force resulting from magnetic induction to the viscous force arising in the flow), 𝒮=𝜎ℎ₀/𝜌𝜈² (the surface-tension parameter), and 𝜏=𝜏_𝑎ℎ²₀/𝜌𝜈² (the shear stress parameter). The surface-tension parameter is usually large; therefore, it is rescaled as 𝜀²𝒮=𝑆 and set as 𝒪(1) in accordance to the waves observed in laboratory experiments. All of the other quantities are considered 𝒪(1). The long-wave equation is derived in the next step.

3. Long-Wave Equation

The dependent variables are asymptotically expanded in terms of the slenderness parameter 𝜀 to derive the long-wave equation. Using the symbolic math toolbox available in MATLAB, the zeroth and the first-order systems are solved. These solutions are then substituted in the kinematic condition (A.5) to derive a Benney-type model accurate up to 𝒪(𝜀) of the form (1) as ℎ_𝑡+𝐴(ℎ)ℎ_𝑥+𝜀𝐵(ℎ)ℎ_𝑥+𝐶(ℎ)ℎ_{𝑥𝑥𝑥𝑥}+𝒪𝜀²=0. (10) The standard procedure for the derivation [7, 9, 14, 15, 28] is skipped here. The expressions for 𝐴(ℎ), 𝐵(ℎ), and 𝐶(ℎ) are, the following: 𝐴(ℎ)=ReHatanh²√+𝜏Haℎ√√HatanhHaℎsech√Haℎ, (11)𝐵(ℎ)=Recot𝛽Ha^3/2tanh√Haℎ−Recot𝛽ℎ+HaRe²ℎ2Ha^5/2tanh√Haℎsech⁴√+HaℎRe²ℎHa^5/2tanh√Haℎsech²√−Haℎ3Re²2Ha³tanh²√Haℎsech²√Haℎ+𝜏3Re2Ha^5/2tanh√Haℎsech√−Haℎ3ReHa^5/2tanh√Haℎsech³√+HaℎReℎHa²sech⁵√+3Haℎ2sech³√√Haℎ−sechHaℎ+𝜏²ℎtanh√HaℎHa^3/2−sech²√−1Haℎ2sech⁴√+3Haℎ2Ha²tanh²√Hasech²√,𝐶(ℎ)=Haℎ𝑆ℎ−𝑆HaHa^3/2tanh√Haℎ. (12)(13) It should be remarked that, when 𝜏=0, the above terms agree with the long-wave equation derived by Tsai et al. [28] when the phase-change effects on the interface are neglected. The dimensionless parameters differ from the nondimensional set presented in Tsai et al. [28] because viscous scales are employed here. Although one can consider different scales, in principle, the structure of the evolution equation remains unchanged regardless of the dimensionless parameters appearing in the problem. Furthermore, the error in the second term associated with 𝐵(ℎ) in Tsai et al. [28] is corrected here, which as per the convention followed in their paper should read as (4𝛼Reℎ/𝑚⁵)sech²𝑚ℎtanh𝑚ℎ. Also, the case corresponding to Ha=0 can be recovered from (11)–(13) when Ha→0. In this case, the functions in (11)–(13) read as 𝐴(ℎ)=Reℎ²𝐵(ℎ)=2+𝜏ℎ,15Reℎ⁵()−1Reℎ+𝜏3Recot𝛽ℎ³,𝐶(ℎ)ℎ=𝑆³3, (14) which match with the evolution equation derived by Miladinova et al. [9] in the absence of Marangoni and air shearing effect. The effect of magnetic field and air shear affects the leading order solution through 𝐴(ℎ) term in (11). This term contributes towards wave propagation and steepening mechanism. The effect of hydrostatic pressure is measured by the terms within the first flower bracket in 𝐵(ℎ) in (12). The rest of the terms in 𝐵(ℎ) affect the mean flow due to inertial, air shear and magnetic field contributions. The function 𝐶(ℎ) corresponds to the mean surface-tension effect. Although the effect of 𝜏 cannot be properly judged based on its appearance in (11)–(13), it is obvious that in the absence of the magnetic field both 𝐴(ℎ) and 𝐵(ℎ) increase (decrease) when 𝜏>0(<0). The stability of the long-wave model (10) subject to (11)–(13) is investigated next.

4. Stability Analysis

The Nusselt solution corresponding to the problem isℎ₀𝑢=1,₀=sinh√Ha𝑦√HaRe√√HatanhHa+𝜏sech√+HaRe√Ha1−cosh,𝑣Ha𝑦₀𝑝=0,₀=Recot𝛽1−𝑦. (15a)(15b)(15c)(15d)

For a parallel shear flow, (10) subject to (11)–(13) admits normal mode solutions of the form ℎ=1+𝐻(𝑥,𝑡), where 𝐻(𝑥,𝑡) is the unsteady part of the film thickness representing the disturbance component such that 𝐻≪1 [56, 57]. Inserting ℎ=1+𝐻(𝑥,𝑡) in (10) and invoking a Taylor series expansion about ℎ=1, the unsteady nonlinear equation representing a slight perturbation to the free surface is obtained as 𝐻_𝑡+𝐴₁𝐻_𝑥+𝜀𝐵₁𝐻_𝑥𝑥+𝜀𝐶₁𝐻_{𝑥𝑥𝑥𝑥}=−𝐴^′₁𝐴𝐻+^1′′2𝐻²+𝐴₁^′′′6𝐻³𝐻_𝑥−𝜀𝐵^′₁𝐵𝐻+^1′′2𝐻²𝐻_𝑥𝑥−𝜀𝐶^′₁𝐶𝐻+^1′′2𝐻²𝐻_{𝑥𝑥𝑥𝑥}−𝜀𝐵^′₁+𝐵₁^′′𝐻𝐻_𝑥²−𝜀𝐶^′₁+𝐶₁^′′𝐻𝐻_𝑥𝐻_𝑥𝑥𝑥+𝒪𝜀𝐻⁴,𝐻⁵,𝜀²𝐻, (16) where (a prime denotes the order of the derivative with respect to ℎ) 𝐴₁=𝐴(ℎ=1),𝐴^′₁=𝐴^′(),𝐴ℎ=1_1′′=𝐴^′′(),𝐵ℎ=1₁=𝐵(ℎ=1),𝐵^′₁=𝐵^′(),𝐵ℎ=1_1′′=𝐵^′′(),𝐶ℎ=1₁=𝐶(ℎ=1),𝐶^′₁=𝐶^′(),𝐶ℎ=1_1′′=𝐶^′′(ℎ=1). (17) It should be remarked that, while expanding the Taylor series, there are two small parameters, namely, 𝜀≪1 and 𝐻≪1 whose orders of magnitude should be considered such that 𝜀≪𝐻≪1. When 𝒪(𝜀𝐻³) terms are retained and since 𝜀𝐻³≪𝐻⁴, (𝐴₁^′′′/6)𝐻³𝐻_𝑥 appears as a unique contribution of order 𝐻⁴ in (16). This term, although present in the unsteady equation (16), does not contribute when a multiple-scale analysis is done (refer to Section 4.2.1 and Appendix B), where equations only up to 𝒪(𝛼³) are considered while deriving a complex Ginzburg-Landau-type equation [12, 14, 17, 56, 57]. In addition, such a term did not appear in earlier studies [12, 14, 17, 56, 57] because 𝐴(ℎ) was a mathematical function of second degree in ℎ, whose higher-order derivatives are zero.

Equation (16) forms the starting point for the linear stability analysis and describes the behavior of finite-amplitude disturbances of the film. Such an equation predicts the evolution of timewise behavior of an initially sinusoidal disturbance given to the film. It is important to note that the constant film thickness approximation with long-wave perturbations is a reasonable approximation only for certain segments of flow and implies that (16) is only locally valid.

4.1. Linear Stability

To assess the linear stability, the linear terms in (16) are considered. The unsteady part of the film thickness is decomposed as (a tilde denotes the complex conjugate) 𝐻(𝑥,𝑡)=𝜂𝑒^{𝑖(𝑘𝑥−𝐶𝐿𝑡)+𝐶𝑅𝑡}+𝜂𝑒^{−𝑖(𝑘𝑥−𝐶𝐿𝑡)−𝐶𝑅𝑡}, (18) where 𝜂(≪1) is a complex disturbance amplitude independent of 𝑥 and 𝑡. The complex eigenvalue is given by 𝜆=𝐶_𝐿+𝑖𝐶_𝑅 such that 𝑘∈[0,1] represents the streamwise wavenumber. The linear wave velocity and the linear growth rate (amplification rate) of the disturbance are 𝐶_𝐿/𝑘, 𝐶_𝑅∈(−∞,∞), respectively. Explicitly, they are found as 𝐶_𝐿𝑉=𝐶_𝐿𝑘=𝐴₁,𝐶_𝑅=𝜀𝑘²𝐵₁−𝑘²𝐶₁. (19)

The disturbances grow (decay) when 𝐶_𝑅>0 (𝐶_𝑅<0). However, when 𝐶_𝑅=0, the curve 𝑘=0 and the positive branch of 𝑘²=𝐵₁/𝐶₁ represent the neutral stability curves. Identifying the positive branch of 𝐵₁−𝑘²𝐶₁=0 as 𝑘_𝑐=(𝐵₁/𝐶₁) (𝑘_𝑐 is the critical wavenumber), the wavenumber 𝑘_𝑚 corresponding to the maximal growth rate is obtained from (𝜕/𝜕𝑘)𝐶_𝑅=0. This gives 𝐵₁−2𝑘²_𝑚𝐶₁=0 such that 𝑘_𝑚=𝑘_𝑐/√2. The linear amplification of the most unstable mode is calculated from 𝐶_𝑅|_{𝑘=𝑘𝑚}.

A parametric study considering the elements of the set 𝕊={Re,Ha,𝑆,𝛽,𝜏} is done in order to trace the neutral stability and linear amplification curves by assuming the slenderness parameter to be 0.1. Only those curves which are relevant in drawing an opinion are presented.

The influence of the magnetic field on the critical Reynolds number (Re_𝑐) varying as a function of the shear parameter is presented in Figure 2. For each Ha, there is a Re_𝑐 below which the flow is stable. The critical Reynolds number decreases when the angle of inclination increases, and, therefore, the flow destabilizes. The stabilizing effect of the magnetic field is also seen when Ha increases. There exists a certain 𝜏>0 such that the flow remains unstable beyond it.

(a)

(b)

Figure 2 

Critical Reynolds number as a function of 𝜏 when 𝑆=0: (a) 𝛽=45∘; (b) 𝛽=80∘.

Figures 3 and 4 display the neutral stability curves, which divide the 𝑘−Re and 𝑘−𝜏 planes into regions of stable and unstable domains. On the other hand, Figure 5 shows the linear amplification curves. The shear stress offered by the air destabilizes the flow when it flows along the downhill direction (𝜏>0) and increases the instability threshold compared to the case corresponding to 𝜏=0 (Figures 3 and 4). However, the flow mechanism is better stabilized when the applied shear stress offered by the air is in the uphill direction (𝜏<0) than when 𝜏=0 (Figures 3 and 4). As seen from Figure 3, the portion of the axis corresponding to the unstable Reynolds numbers increases and extends towards the left when the angle of inclination increases, thereby reducing the stabilizing effect offered by the hydrostatic pressure at small inclination angles. From curves 2 and 3 corresponding to Figure 3, it is observed that the force of surface-tension stabilizes the flow mechanism. Figure 4 supports the information available from Figure 3. When the magnitude of the Hartmann number increases, the instability region decreases because the value of the critical wavenumber decreases. Comparing Figure 4(a) with Figure 4(c), it is also observed that the inertial force destabilizes the flow mechanism. The growth rate curves (Figure 5) agree with the results offered by the neutral stability curves (Figures 3 and 4).

(a)

(b)

Figure 3 

Neutral stability curves. The symbols ST and US indicate the stable and unstable regions. Here, (1) Ha=0.001 and 𝑆=1; (2) Ha=0.2 and 𝑆=1; (3) Ha=0.2 and 𝑆=10. (a) 𝛽=30∘ and (b) 𝛽=60∘.

(a)

(b)

(c)

Figure 4 

Neutral stability curves in the 𝑘−𝜏 plane for 𝑆=1: (1) Ha=0; (2) Ha=0.05; (3) Ha=0.1; and (4) Ha=0.15. (a) 𝛽=60∘ and Re=1; (b) 𝛽=75∘ and Re=1; and (c) 𝛽=60∘ and Re=2.5.

(a)

(b)

Figure 5 

Growth rate curves when We=1: (a) 𝛽=75∘ and (b) 𝛽=90∘.

The linearly increasing graphs of 𝐶_𝐿 and the decreasing plots of 𝐶_𝐿𝑉 with respect to 𝑘 and Ha, respectively, are presented in Figure 6. The effect of inertia doubles the linear wave speed, 𝐶_𝐿𝑉. But when the magnitude of Ha, increases, the linear wave speed decreases. The downhill effect of the shear stress on the interface makes the linear wave speed larger than the cases corresponding to 𝜏=0 and 𝜏<0.

Figure 6 

Variation of the linear wave velocity, 𝐶𝐿𝑉, with respect to Ha, and 𝐶𝐿 versus 𝑘.

The linear stability results give only a firsthand information about the stability mechanism. The influence of air-induced shear on the stability of the flow under the application of a transverse magnetic field will be better understood only when the nonlinear effects are additionally considered. To analyze and illustrate the nonlinear effects on the stability threshold, a weakly nonlinear study is performed in the next step.

4.2. Multiple-Scale Analysis

4.2.1. Weakly Nonlinear Theory

In order to do multiple-scale analysis, the following slow scales are introduced following Sadiq and Usha [14] (the justification for stretching the scales is provided in Lin [56] and in Krishna and Lin [57]): 𝑡₁=𝛼𝑡,𝑡₂=𝛼²𝑡,𝑥₁=𝛼𝑥, (20) such that 𝜕⟶𝜕𝜕𝑡𝜕𝜕𝑡+𝛼𝜕𝑡₁+𝛼²𝜕𝜕𝑡₂,𝜕⟶𝜕𝜕𝑥𝜕𝜕𝑥+𝛼𝜕𝑥₁. (21) Here, 𝛼 is a small parameter independent of 𝜀 and measures the distance from criticality such that 𝐶_𝑅∼𝒪(𝛼²)(or,𝐵₁−𝑘²𝐶₁∼𝒪(𝛼²)). The motivation behind such a study lies in deriving the complex Ginzburg-Landau equation (CGLE), which describes the evolution of amplitudes of unstable modes for any process exhibiting a Hopf bifurcation. Using such an analysis, it is possible to examine whether the nonlinear waves in the vicinity of criticality attain a finite height and remain stable or continue to grow in time and eventually become unstable. Refer to Appendix B for the detailed derivations of the threshold amplitude, 𝛼Γ₀, and the nonlinear wave speed, 𝑁_𝑐.

The threshold amplitude subdivides the flow domain according to the signs of 𝐶_𝑅 and 𝐽₂ [14]. If 𝐶_𝑅<0 and 𝐽₂<0, the flow exhibits subcritical instability. However, when 𝐶_𝑅>0 and 𝐽₂>0, the Landau state is supercritically stable. If, on the other hand, 𝐶_𝑅<0 and 𝐽₂>0, the flow is subcritically stable. A blow-up supercritical explosive state is observed when 𝐶_𝑅>0 and 𝐽₂<0.

Different regions of the instability threshold obtained through multiple-scale analysis are illustrated in Figure 7. The subcritical unstable region is affected due to the variation of 𝜏. For 𝜏<0, such a region is larger than the ones corresponding to 𝜏≥0. In addition, the subcritical unstable and the stable regions increase when Ha increases. The explosive state region (also called the nonsaturation zone) decreases when 𝜏<0 than when 𝜏≥0. The strip enclosing asterisks is the supercritical stable region where the flow, although linearly unstable, exhibits a finite-amplitude behavior and saturates as time progresses [8, 9, 14]. The bottom line of the strip is the curve 𝑘_𝑠=𝑘_𝑐/2 which separates the supercritical stable region from the explosive region [56, 58, 59]. Within the strip 𝑘_𝑠<𝑘<𝑘_𝑐, different possible shapes of the waves exist [7, 9, 14].

Figure 7 

Variation of 𝑘 with respect to Re when 𝛽=50∘ and 𝑆=5: (1) 𝐶𝑅<0 and 𝐽2<0; (2) 𝐶𝑅<0 and 𝐽2>0; (3) 𝐶𝑅>0 and 𝐽2<0. The strip enclosing asterisks indicates that 𝐶𝑅>0 and 𝐽2>0. The broken line within the strip represents 𝑘𝑚=𝑘𝑐/√2.

Tables 1 and 2 show the explosive and equilibration state values for different flow parameters. The wavenumber value corresponding to the explosive state increases when the inertial effects increase. This increases the unstable region 0<𝑘<𝑘_𝑠 (Table 1). Such a wavenumber decreases either when the hydromagnetic effect is increased or when the air shear is in the uphill direction. It is evident from Table 1 that the explosive state occurs at small values of Re and 𝛽, which is also true for large values of Ha when 𝜏<0. Also, it is observed that either when the Hartmann number or the value of the uphill shear is increased, the critical wavenumber becomes zero (Table 2). Therefore, 𝑘_𝑠=0 (Table 1).

Table 1 

Parameter values for which the nonlinear wave explodes.

Table 2 

Parameter values for which the nonlinear wave attains a finite amplitude.

The threshold amplitude (𝛼Γ₀) profiles display an asymmetric structure (Figure 8), increasing up to a certain wavenumber (>𝑘_𝑚) and then decreasing beyond it in the supercritical stable region. The 𝜏>0 induced amplitudes show larger peak amplitude than the cases corresponding to 𝜏=0 and 𝜏<0. The peak amplitude value decreases when Ha increases. The nonlinear wave speed (𝑁_𝑐) curves represent a 90^∘ counterclockwise rotation of the mirror image of the alphabet 𝐿 when Ha is small. The nonlinear speed decreases to a particular value in the vertical direction and thereafter traces a constant value beyond it as the wavenumber increases in the supercritical stable region. However, when the effect of the magnetic field is increased, the nonlinear wave speed decreases in magnitude and sketches an almost linear constant profile. The magnitude of 𝑁_𝑐 in the supercritical stable region corresponding to 𝜏=0 remains inbetween the values corresponding to 𝜏<0 and 𝜏>0.

Figure 8 

Threshold amplitude and nonlinear wave speed in the supercritical stable region for 𝛽=50∘ and 𝑆=5: (1) 𝜏=−1, (2) 𝜏=0, and (3)𝜏=1.

4.2.2. Sideband Instability

Let us consider the quasi-monochromatic wave of (B.4) exhibiting no spatial modulation of the form 𝜂=Γ_𝑠𝑡₂=Γ₀𝑒^{−𝑖𝑏𝑡2}, (22) where Γ₀ and 𝑏 are defined by (B.10) and (B.11) such that 𝜔_𝑠+𝑏 is the wave frequency and 𝑘(>0) is the modulation wavenumber.

If one considers a band of frequencies centered around 𝜔_𝑠, the interaction of one side-mode with the second harmonic would be resonant with the other side-mode causing the frequency to amplify. This leads to an instability known as sideband instability [60, 61].

To investigate such an instability, Γ_𝑠(𝑡₂) is subjugated to sideband disturbances of bandwidth 𝑘𝛼 [12, 56, 57]. The explicit expression for the eigenvalues is found as (refer to Appendix C) 1𝜗=2tr(ℳ)±tr²(ℳ)−4det(ℳ)=−𝑐_𝑖+𝑘²𝐽_1𝑟±124𝑐²_𝑖+𝑘²𝑣²−8𝑖𝐽₄𝑐_𝑖𝑣𝑘𝐽₂=𝜛₁±𝜛₂, (23) where tr(ℳ)=2𝑐_𝑖+2𝑘²𝐽_1𝑟−4𝐽₂Γ_𝑠²=−2𝑐_𝑖+2𝑘²𝐽_1𝑟,det(ℳ)=𝐽²₁𝑘⁴−𝑣²+2𝑐_𝑖𝐽₁𝑘²−2𝑖𝑐_𝑖𝑣𝐽₄𝐽₂𝑘. (24)

It should be remarked that the above expression for 𝜗 is true only if 𝑣∼𝒪(𝛼). However, when 𝑣=0, it is easily seen from (23) that 𝜗₁=𝜛₁+𝜛₂=𝑘²𝐽_1𝑟<0 and 𝜗₂=𝜛₁−𝜛₂=−2𝑐_𝑖+𝑘²𝐽_1𝑟<0 (when 𝑐_𝑖>0), implying that the system is stable to the sideband disturbances as 𝑡₂→∞. For nonzero 𝑣, the eigenvalues depend on the dimensionless flow parameters. If 𝜛₂>0 and is less than the absolute value of 𝜛₁, the sideband modes stabilize the system as 𝑡₂→∞. However, if 𝜛₂>0 and is greater than the absolute value of 𝜛₁, only one of the modes is sideband stable. On the other hand, if (𝑐²_𝑖+𝑘²𝑣²)−(8𝑖𝐽₄𝑐_𝑖𝑣𝑘⧵𝐽₂)<0, again, one of the modes is sideband stable.

5. Equilibria of a Bimodal Dynamical System

Considering the initial thickness of the amplitude to be one, Gjevik [58] analyzed the amplitude equations by representing the amplitude using a truncated Fourier series and by imposing restrictions on its coefficients. The velocity along the mean flow direction and the corresponding surface deflection moving with this velocity were deduced by posing the problem into a dynamical system. Gottlieb and Oron [62] and Dandapat and Samanta [12] also expanded the evolution equation using a truncated Fourier series to derive a modal dynamical system. The results in Gottlieb and Oron [62] showed that a two-mode model was found to coincide with the numerical solution along the Hopf bifurcation curve. Based on this confidence, the stability of the bimodal dynamical system was assessed in Dandapat and Samanta [12].

In this section, the stability of a truncated bimodal dynamical system is analyzed using the approach followed by the above authors, but using the assumptions considered in Gjevik [58]. The coupled dynamical system and its entries are listed in Appendix D.

It should be remarked that in the 𝑘−Re plane the equations 𝑏₁₁₁=0 and 𝜕[𝑘b₁₁₁]/𝜕𝑘=0 give the neutral stability curve and the curve corresponding to the maximum rate of amplification for linear disturbances. The expression 𝑣_𝑐=𝑑𝜃𝑑𝑡¹(𝑡)=𝑐₁₁₁+𝑐₁₂₁cos𝜙(𝑡)−𝑏₁₂₁sin𝜙(𝑡)×𝑎₂(𝑡)+𝑐₁₃₁𝑎²₁(𝑡)+𝑐₁₃₂𝑎²₂(𝑡) (25) measures the velocity along the direction of the mean flow of a steady finite-amplitude wave. The surface deflection, ℎ(𝑥,𝑡), moving with velocity 𝑣_𝑐 along the mean flow direction in a coordinate system 𝑥^′=𝑥+𝜃₁ is calculated from (D.1) and reads as ℎ(𝑥,𝑡)=1+2𝐵₁cos𝑥^′+𝐵₂cos2𝑥^′−𝜙. (26)

The steady solutions of the system (D.2a) and (D.2b) correspond to the fixed points of (D.4a)–(D.4c). In addition to the trivial solution 𝑎₁(𝑡)=0, 𝑎₂(𝑡)=0, and 𝜙(𝑡)=0 (𝑧₁(𝑡)=0;𝑧₂(𝑡)=0), system (D.4a)–(D.4c) offers nontrivial fixed points [62–64]. These fixed points can be classified as pure-mode fixed points (where one of the fixed points is zero and the others are nonzero), mixed-mode fixed points (nonzero fixed points with a zero or nonzero phase difference), and traveling waves with nonzero constant phase difference.

5.1. Pure-Mode

The fixed points in this case are obtained by setting 𝑓_𝑖(𝑎₁,𝑎₂,𝜙)=0, 𝑖=1,2,3, in (D.4a)–(D.4c). This gives 𝑎₁(𝑡)=0 and 𝑎²₂(𝑡)=−𝑏₂₁₁/𝑏₂₃₂. The solution exists if sgn(𝑏₂₁₁𝑏₂₃₂)=−1. For convenience, 𝑎²₂(𝑡)=𝑎² is set. From (D.4c), 𝑐₁₂₁cos𝜙(𝑡)−𝑏₁₂₁sin𝜙(𝑡)=(𝑘/2)𝐴₁^′′𝑎 is obtained. From this condition, the fixed point for 𝜙(𝑡) is derived as 𝜙(𝑡)=𝜙(𝑡)=tan⁻¹4𝑏₁₂₁𝑐₁₂₁±𝑘𝐴₁^′′𝑎×4𝑏²₁₂₁+4𝑐²₁₂₁−𝑘²𝐴₁^′′2𝑎^21/2×4𝑏²₁₂₁−𝑘²𝐴₁^′′2𝑎²⁻¹, (27) provided that the quantity within the square root is real and positive. The stability of the nonlinear dynamical system (D.4a)–(D.4c) can be locally evaluated using the eigenvalues of the matrix obtained after linearizing the system around the fixed points. The linear approximation of the dynamical system (D.4a)–(D.4c) can be represented in matrix notation as 𝑑𝑎₁(𝑡)𝑑𝑡𝑑𝑎₂(𝑡)(𝑡)𝑑𝑡𝑑𝜙=𝑝𝑑𝑡¹000𝑝₂00𝑝₃𝑝₄𝑎₁(𝑡)𝑎₂(𝑡)𝜙(𝑡)+𝑞₁𝑞₂𝑞₃, (28) such that 𝑝₁=𝑏₁₁₁+𝑏₁₂₁cos𝜙+𝑐₁₂₁sin𝜙𝑎+𝑏₁₃₂𝑎²;𝑝₂=𝑏₂₁₁+3𝑏₂₃₂𝑎²,𝑝₃=−2𝑏₁₂₁sin𝜙+2𝑐₁₂₁cos𝜙−2𝑘𝐴₁^′′𝑝𝑎;₄=−2𝑏₁₂₁cos𝜙+2𝑐₁₂₁𝜙𝑞sin𝑎,₁=0;𝑞₂=−2𝑏₂₃₂𝑎³;𝑞₃=𝑘𝐴₁^′′𝑎²+2𝑏₁₂₁cos𝜙+2𝑐₁₂₁sin𝜙𝜙𝑎. (29)

The stability of the linear system (28) depends on the eigenvalues of the coefficient matrix with entries 𝑝_𝑖. The eigenvalues are found as 𝜆₁=𝑏₂₁₁+3𝑎²𝑏₂₃₂=−2𝑏₂₁₁,𝜆₂=−2𝑎𝑏₁₂₁cos𝜙+𝑐₁₂₁𝜙,𝜆sin₃=𝑏₁₁₁+𝑎𝑏₁₂₁cos𝜙+𝑐₁₂₁sin𝜙+𝑏₁₃₂𝑎². (30)

The Hopf bifurcation at the critical threshold is defined by 𝐵₁=𝑘²𝐶₁, and this yields 𝜆₁=24𝜀𝑘⁴𝐶₁>0. This eigenvalue being independent of 𝜏 increases when the effects of surface-tension and Hartmann number increase. Therefore, at the critical threshold, the fixed points corresponding to the pure-mode are unstable. Beyond the neutral stability limit where the flow is linearly stable (𝐶_𝑅<0) and where the stable wavenumber region decreases when the direction of the shear offered by the wind changes from uphill to downhill direction, the eigenvalue may remain positive and still destabilize the system as shown in Figure 9. The stable wavenumbers corresponding to the linear stability threshold (refer to Figure 5 with Ha=0.3 and Re=2) are considered to plot the eigenvalue 𝜆₁. Although the magnitude of 𝜆₁ remains larger for 𝜏<0 than for 𝜏≥0, being positive, it plays a destabilizing role.

Figure 9 

Eigenvalue 𝜆1 as a function of 𝑘 when We=1 and 𝛽=75∘.

5.2. Mixed-Mode

When 𝑎₁(𝑡), 𝑎₂(𝑡)≠0, the fixed points in this case correspond to mixed-mode. Considering 𝜙(𝑡)≠0 and 𝑎^∗₁, 𝑎^∗₂, and 𝜙^∗ to be the fixed points of (D.4a)–(D.4c), the nonlinear system can be linearized around the fixed points. Then, the stability (instability) of the fixed points demands all of the eigenvalues of the linearized Jacobian matrix to be negative (at least one of them to be positive). The nine entries 𝑡_𝑖𝑗(𝑖,𝑗=1,2,3) of the 3 × 3 Jacobian matrix arising due to linearization are the following: 𝑡₁₁=𝑏₁₁₁+𝑏₁₂₁cos𝜙^∗+𝑐₁₂₁sin𝜙^∗𝑎^∗₂+3𝑏₁₃₁𝑎₁^∗2+𝑏₁₃₂𝑎₂^∗2,𝑡₁₂=𝑏₁₂₁cos𝜙^∗+𝑐₁₂₁sin𝜙^∗𝑎^∗₁+2𝑏₁₃₂𝑎^∗₁𝑎^∗₂,𝑡₁₃=−𝑏₁₂₁sin𝜙^∗+𝑐₁₂₁cos𝜙^∗𝑎^∗₁𝑎^∗₂,𝑡₂₁=2𝑏₂₂₁cos𝜙^∗−𝑐₂₂₁sin𝜙^∗𝑎^∗₁+2𝑏₂₃₁𝑎^∗₁𝑎^∗₂,𝑡₂₂=𝑏₂₁₁+𝑏₂₃₁𝑎₁^∗2+3𝑏₂₃₂𝑎₂^∗2,𝑡₂₃=−𝑏₂₂₁sin𝜙^∗+𝑐₂₂₁cos𝜙^∗𝑎₁^∗2,𝑡₃₁=2𝑘𝐴₁^′′𝑎^∗₁−2𝑎^∗₁𝑎^∗₂𝑏₂₂₁sin𝜙^∗+𝑐₂₂₁cos𝜙^∗,𝑡₃₂=−2𝑘𝐴₁^′′𝑎^∗₂+2𝑐₁₂₁cos𝜙^∗−𝑏₁₂₁sin𝜙^∗+𝑎₁^∗2𝑎₂^∗2𝑏₂₂₁sin𝜙^∗+𝑐₂₂₁cos𝜙^∗,𝑡₃₃=−2𝑐₁₂₁sin𝜙^∗+𝑏₁₂₁cos𝜙^∗𝑎^∗₂−𝑏₂₂₁cos𝜙^∗−𝑐₂₂₁sin𝜙^∗𝑎₁^∗2𝑎^∗₂. (31)

It should be remarked that Samanta [65] also discussed the mixed-mode for the flow of a thin film on a nonuniformly heated vertical wall. However, it should be noted that the Jacobian was computed not by considering a zero phase difference with 𝑎₁(𝑡),𝑎₂(𝑡)≠0, but by evaluating the Jacobian first by imposing 𝑎₁(𝑡),𝑎₂(𝑡),𝜙(𝑡)≠0 and then by substituting 𝜙=0 in the computed Jacobian. If one considers 𝑎₁(𝑡),𝑎₂(𝑡),𝜙(𝑡)=0, an overdetermined system is obtained.

5.3. Traveling Waves

Fixed points of (D.4a)–(D.4c) with 𝜙 being a nonzero constant correspond to traveling waves. The modal amplitudes are considered small in the neighborhood of the neutral stability limit. After rescaling 𝑎₁(𝑡)→𝛼𝑚₁(𝑡) and 𝑎₂(𝑡)→𝛼²𝑚₂(𝑡), the following system is obtained:𝑑𝑚₁(𝑡)𝑑𝑡=𝑏₁₁₁𝑚₁(𝑡)+𝛼²×𝑏₁₂₁cos𝜙(𝑡)+𝑐₁₂₁sin𝜙(𝑡)𝑚₁(𝑡)𝑚₂(𝑡)+𝛼²𝑏₁₃₁𝑚³₁(𝑡)+𝒪𝛼⁴,𝑑𝑚₂(𝑡)𝑑𝑡=𝑏₂₁₁𝑚₂(𝑡)+𝑏₂₂₁cos𝜙(𝑡)−𝑐₂₂₁sin𝜙(𝑡)𝑚²₁(𝑡)+𝛼²𝑏₂₃₁𝑚²₁(𝑡)𝑚₂(𝑡)+𝒪𝛼⁴,𝑑𝜙(𝑡)𝑏𝑑𝑡=−₂₂₁sin𝜙(𝑡)+𝑐₂₂₁(𝑡𝑚cos𝜙)²¹(𝑡)𝑚₂(𝑡)+𝛼²𝑘𝐴₁^′′𝑚²₁(𝑡)+2𝛼²×𝑐₁₂₁cos𝜙(𝑡)−𝑏₁₂₁sin𝜙(𝑡)𝑚₂(𝑡)+𝒪𝛼⁴. (32a)(32b)(32c)

To the 𝒪(𝛼²), there are no nontrivial fixed points in the above system. The phase evolution is governed by equating the right-hand side of (32c) to zero by considering terms up to 𝒪(𝛼²): tan𝜙(𝑡)=−𝑙₁𝑙₂±𝑙₃𝑙²₁+𝑙²₂−𝑙²₃^1/2𝑙²₁−𝑙²₃,, (33) where 𝑙₁=(1+2𝛼²𝑏₁₂₁𝑚²₂(𝑡)/𝑏₂₂₁𝑚²₁(𝑡)), 𝑙₂=(𝑐₂₂₁/𝑏₂₂₁)−(2𝛼²𝑐₁₂₁𝑚²₂(𝑡)/𝑏₂₂₁𝑚²₁(𝑡)), and 𝑙₃=𝛼²𝑘𝐴₁^′′𝑚₂(𝑡)/𝑏₂₂₁. In the above equation, if only the leading order effect is considered, the fixed point is found as 𝜙^∗(𝑡)=tan⁻¹−𝑐₂₂₁𝑏₂₂₁+𝒪𝛼². (34)

Equating (32a) to zero and using (34), it is found that 𝛼²𝑏₁₃₁𝑚²₁(𝑡)=−𝑏₁₁₁−𝛼²𝜑₁𝑚₂(𝑡). (35) Such a solution exists if [𝑏₁₁₁+𝛼²𝜑₁𝑚₂(𝑡)]𝑏₁₃₁<0, where 𝜑₁=(𝑏₁₂₁cos𝜙^∗+𝑐₁₂₁sin𝜙^∗). A quadratic equation in 𝑚₂(𝑡) is obtained by considering (32b) as 𝛼⁴𝑏₂₃₁𝜑₁𝑚²₂(𝑡)+𝛼²𝑏₁₁₁𝑏₂₃₁−𝑏₂₁₁𝑏₁₃₁+𝜑₁𝜑₂×𝑚₂(𝑡)+𝜑₂𝑏₁₁₁=0, (36) where 𝜑₂=𝑏₂₂₁cos𝜙^∗−𝑐₂₂₁sin𝜙^∗. For this equation, there are two solutions. The existence of such solutions to be a real number demands (𝑏₁₁₁𝑏₂₃₁−𝑏₂₁₁𝑏₁₃₁+𝜑₁𝜑₂)²−4𝛼⁴𝜑₁𝜑₂𝑏₁₁₁𝑏₂₃₁>0. Since 𝑏₁₁₁=0 in the neutral stability limit, (36) gives the amplitude of the nonzero traveling wave as 𝑚₂(𝑡)=𝑏₂₁₁𝑏₁₃₁−𝜑₁𝜑₂𝛼²𝑏₂₃₁𝜑₁=3𝑘²𝐶₁𝐽₂𝛼²𝐵₁^′′−4𝑘²𝐶₁^′′𝜑₁. (37) When 𝐽₂=0,  𝑚₂(𝑡) is zero, and this corresponds to finding the critical Reynolds number in the neutral stability limit [12]. If 𝑚^∗₁ and 𝑚^∗₂ are the fixed points obtained from (35) and (36), the stability or instability of the fixed points corresponding to traveling waves depends on the sign of the eigenvalues of the following matrix:𝐽=2𝛼²𝑏₁₃₁𝑚₁^∗2𝛼²𝜑₁𝑚^∗₁−𝛼²sin𝜙^∗𝑏₁₂₁+𝑐₁₂₁𝑏₂₂₁𝑐₂₂₁𝑚^∗₁𝑚^∗₂2𝜑₂+𝛼²𝑏₂₃₁𝑚^∗₂𝑚^∗₁𝛼²𝑏₂₁₁𝑚₁^∗202𝛼²𝑘𝐴₁^′′𝑚^∗₁−2𝛼²𝑏₁₂₁+𝑐₁₂₁𝑏₂₂₁𝑐₂₂₁sin𝜙^∗−2𝛼²𝜑₁𝑚^∗₂+𝜑₂𝑚₁^∗2𝑚^∗₂. (38)

6. Nonlinear Development of the Interfacial Surface

The nonlinear interactions are studied numerically beyond the linear stability threshold in a periodic domain, 𝔻={(𝑥,𝑡)∶𝑥∈(−𝜋/𝑘_𝑚,𝜋/𝑘_𝑚),𝑡∈[0,∞)}, by solving (10) subject to the initial condition ℎ(𝑥,0)=1−0.1cos(𝑘_𝑚𝑥). A central difference scheme, which is second-order accurate in space, and a backward Euler method, which is implicit in forward time, are adopted to solve the problem in MATLAB. The local truncation error for this numerical scheme is 𝒯=𝒪(Δ𝑡,Δ𝑥²). The computations are performed with a small time increment, Δ𝑡=10⁻³. The numerical scheme is always stable (since a backward Euler method is used) and does not build up errors. The number of nodes along the spatial direction is 𝑁=800 such that the spatial step length is Δ𝑥=2𝜋/𝑁𝑘_𝑚. An error tolerance of 10⁻¹⁰ is set, and the simulations are stopped once the absolute value of the error becomes smaller than this value. Also, the numerical simulations show no particular deviation from the results obtained neither when the spatial grid points are doubled nor when the time step is further reduced. Figures 5(b), 7(a), and 9(a) available in Joo et al. [7] are reproduced in Figure 10 to show the correctness of the numerical scheme. This gives confidence in applying it to the evolution equation considered here.

(a)

(b)

(c)

Figure 10 

(a), (b), and (c) represent the surface-wave instability curves reproduced from Joo et al. [7] corresponding to Figures 5(b), 7(a), and 9(a), respectively, in their study.

The effect of the transverse magnetic field is illustrated in Figure 11. The wave structure when the magnetic field strength is zero displays a steep curvier wave than when Ha>0. The surface-wave instability decreases when the strength of the applied magnetic field increases. This reveals the stabilizing mechanism of the transversely applied magnetic field.

Figure 11 

Free surface profiles at Re=5, 𝛽=45∘, and 𝑆=5 when 𝜏=0 at time 𝑡=250.

Figure 12 displays maximum amplitude profiles of the sheared flow for two different values of Ha on a semi-inclined plane. The maximum amplitude initially increases and then decreases. Due to saturation of the nonlinear interactions, the initial perturbation of the free surface is damped after a long time. For a short time, the amplitude profile corresponding to 𝜏>0 remains smaller than the ones corresponding to 𝜏≤0. Also, for different Ha, the amplitudes corresponding to 𝜏<0 remain larger than those corresponding to 𝜏=0 up to a certain time (see also Figure 13). However, these trends change over time. When the magnitude of the Hartmann number increases, the time required for the amplitudes corresponding to 𝜏<0 to decrease below the amplitude corresponding to 𝜏=0 increases (Figures 12(a), 12(b), and 13). In addition, it is also observed that the time needed for the amplitude profiles corresponding to 𝜏=−1 to dominate the amplitude profiles corresponding to 𝜏=−1.5 increases when the value of Ha increases (Figures 12(a) and 12(b)).

(a)

(b)

Figure 12 

Maximum film thickness profiles at Re=5, 𝛽=45∘, and 𝑆=5: (a) Ha=0.05 and (b) Ha=0.2.

Figure 13 

Maximum and minimum amplitude profiles corresponding to Re=5, 𝛽=45∘, 𝑆=5, and Ha=0.3.

The surface waves emerging at the free surface are captured and presented in Figure 14. For 𝜏=0, a wave which has a stretched front is observed after a long time. When 𝜏=−1, a one-humped solitary-like wave is formed when a large amplitude wave and a small capillary ripple coalesce together (Figure 14(b)) at the time of saturation (a certain time after which all of the waves have the same structure and shape). Figure 14 demonstrates that the instability measured by the wave height decreases when the shear is induced along the uphill direction.

(a)

(b)

Figure 14 

Evolution of the surface-wave instability at Re=5, 𝛽=45∘, 𝑆=5, and Ha=0.3: (a) 𝜏=0 and (b) 𝜏=−1. The waves are plotted at a gap of 0.1 time units in a periodic domain.

Tsai et al. [28] pointed out that, when the intensity of the magnetic field increases, the flow retards and stabilizes the system. The interfacial surface subjected to air shear affects the velocity and flow rate, 𝑞(𝑥,𝑡)=∫^ℎ₀𝑢(𝑥,𝑡)𝑑𝑦, for different values of 𝜏 and Ha (Figure 15). The numerical interactions reveal that the velocity and the flow rate can either increase or decrease depending on the direction the air shears the deformable free interface. In conjunction to the nonlinear wave speed traced in Figure 8, the velocity profiles traced in Figure 15 show a similar response. For Ha=0 and 𝜏>0, the streamwise velocity and the flow rate increase. Also, the transverse velocity is larger for Ha=0 than when Ha>0. Furthermore, the velocity and the flow rate decrease when the strength of the applied magnetic field increases. When the magnitude of the shear induced in the uphill direction is increased (by considering a small negative 𝜏 value) for large values of Ha, constant velocity and flow rate profiles are observed.

Figure 15 

(1) Streamwise velocity, (2) flow rate, and (3) transverse velocity profiles at Re=5, 𝛽=45∘, and 𝑆=5 in a periodic domain.

It is inferred from the nonlinear simulations that the effect of magnetic field and wind shear greatly affects the thickness of the interfacial free surface. When the strength of the magnetic field increases, the flow has a uniform velocity and constant discharge. For Ha>0, the transverse velocity shows a small growth in the positive direction. The growth of the unwanted development of the free surface is retarded by the magnetic field, when it acts opposite to the flow direction. Such a growth-reducing mechanism of the amplitude, in particular, can be further enhanced by applying uphill shear on the free surface. This inclusion causes further reduction in the values of streamwise and transverse velocities, thereby leading to constant flux situation and improved stability.

7. Conclusions and Perspectives

Linear and nonlinear analyses on the stability of a thin film subjected to air shear on a free surface in the presence of a transversely applied external magnetic field of constant strength have been studied. Instead of obtaining the solutions at different orders of approximation using a power series approach [37–39, 43], the solutions of equations arising at various orders have been straightforwardly solved [28].

The linear stability results gave firsthand information about the stability mechanism. The modal interaction phenomenon was studied by deriving a complex Ginzburg-Landau-type equation using the method of multiple scales. The stability regimes identified by the linear theory were further categorized using the weakly nonlinear theory by considering a filtered wave to be the solution of the complex Ginzburg-Landau equation. When the nonlinear amplification rate (𝐽₂) is positive, an infinitesimal disturbance in the linearly unstable region attains a finite-amplitude equilibrium state. The threshold amplitude and the nonlinear wave speed exist in the supercritical stable region. Stability of a filtered wave subject to sideband disturbances was considered, and the conditions under which the eigenvalues decay in time were mathematically identified. Considering the initial thickness to be one and imposing restrictions upon the amplitude coefficients, a truncated Fourier series was used to derive a bimodal dynamical system. The stability of the pure-mode, mixed-mode, and the traveling wave solution has been mathematically discussed.

Although there are several investigations based on weakly nonlinear theory [28, 37–40, 43], the complete nonlinear evolution equation was not studied numerically beyond the linear stability threshold in the presence of an external magnetic field, especially in a short-circuited system. In this regard, the nonlinear interactions have been numerically assessed using finite-difference technique with an implicit time-stepping procedure. Using such a scheme, the evolution of disturbances to a small monochromatic perturbation was analyzed. The destabilizing mechanism of the downhill shear was identified with the help of computer simulation. Also, it takes a long time for the amplitude profiles corresponding to the zero-shear-induced effect to dominate those corresponding to the uphill shear, when the strength of the magnetic field increases. Furthermore, among the uphill shear-induced flows, the smaller the uphill shear-induced effect is, the longer it takes for the respective amplitude to be taken over by the amplitude corresponding to large values of uphill shear, provided that the magnetic fields intensity increases.

The velocity and the flow rate profiles gave a clear understanding of the physical mechanism involved in the process. When the magnetic field effect and the air shear are not considered, the gravitational acceleration, the inertial force, and the hydrostatic pressure determine the mechanism of long-wave instability, by competing against each other. These forces trigger the flow, amplify the film thickness, and prevent the wave formation, respectively [66]. When the magnetic field effect is included, it competes with other forces to decide the stability threshold. The applied magnetic field retards the flow considerably by reducing the transverse velocity. This phenomenon reduces the wave thickness, but it favors hydrostatic pressure and surface-tension to promote stability on a semi-inclined plane. However, when the air shear is considered, it can either stabilize or destabilize the system depending on the direction it shears the free surface. For downhill shear, the transverse velocity increases resulting in increased accumulation of the fluid under a typical cusp, thereby favoring inertial force. Overall, this mechanism causes the film thickness to increase. For uphill shear, the results are opposite.

The present investigation suggests the inclusion of superficial shear stress on the interfacial surface to enhance the stability of the film, or on the other hand, shows the natural effect of air shear on the films stability under natural conditions, in particular to those films subjected to a transversely applied uniform external magnetic field in a short-circuited system. Future research activities may include (i) assessing nonisothermal effects, (ii) studying non-Newtonian films, and (iii) considering the contribution of electric field. The last point would demand considering Poisson equation to describe electric potential together with appropriate boundary conditions.

Appendices

A. Dimensionless Equations

The dimensionless equations are given by 𝑢_𝑥+𝑣_𝑦𝜀𝑢=0,_𝑡+𝑢𝑢_𝑥+𝑣𝑢_𝑦=−𝜀𝑝_𝑥+𝜀²𝑢_𝑥𝑥+𝑢_𝑦𝑦𝜀+Re−Ha𝑢,₂𝑣_𝑡+𝑢𝑣_𝑥+𝑣𝑣_𝑦=−𝑝_𝑦+𝜀³𝑣_𝑥𝑥+𝜀𝑣_𝑦𝑦−Recot𝛽, (A.1) and the respective boundary conditions are 𝑢=0,𝑣=0on𝑦=0,𝑝+2𝜀1−𝜀²ℎ²_𝑥1+𝜀²ℎ²_𝑥𝑢_𝑥+2𝜀ℎ_𝑥1+𝜀²ℎ²_𝑥𝑢_𝑦+𝜀²𝑣_𝑥𝜀=−²𝒮ℎ_𝑥𝑥1+𝜀²ℎ²_𝑥^3/2on𝑦=ℎ,1−𝜀²ℎ²_𝑥1+𝜀²ℎ²_𝑥𝑢_𝑦+𝜀²𝑣_𝑥−4𝜀²ℎ_𝑥1+𝜀²ℎ²_𝑥𝑢_𝑥=𝜏on𝑦=ℎ,𝑣=ℎ_𝑡+𝑢ℎ_𝑥on𝑦=ℎ. (A.2)(A.3)(A.4)(A.5)

B. Multiple Scale Analysis

The film thickness, 𝐻, is expanded as 𝐻𝜀,𝑥,𝑥₁,𝑡,𝑡₁,𝑡₂=^∞_𝑖=1𝛼^𝑖𝐻_𝑖𝜀,𝑥,𝑥₁,𝑡,𝑡₁,𝑡₂. (B.1) Exploiting (16) using (20) and (B.1), the following equation is obtained: 𝐿₀+𝛼𝐿₁+𝛼²𝐿₂𝛼𝐻₁+𝛼²𝐻₂+𝛼³𝐻₃=−𝛼²𝑁₂−𝛼³𝑁₃, (B.2) where the operators 𝐿₀, 𝐿₁, and 𝐿₂ and the quantities 𝑁₂ and 𝑁₃ are given as follows𝐿₀=𝜕𝜕𝑡+𝐴₁𝜕𝜕𝑥+𝜀𝐵₁𝜕²𝜕𝑥²+𝜀𝐶₁𝜕⁴𝜕𝑥⁴,𝐿₁=𝜕𝜕𝑡₁+𝐴₁𝜕𝜕𝑥₁+2𝜀𝐵₁𝜕²𝜕𝑥𝜕𝑥₁+4𝜀𝐶₁𝜕⁴𝜕𝑥³𝜕𝑥₁,𝐿₂=𝜕𝜕𝑡₂+𝜀𝐵₁𝜕²𝜕𝑥²₁+6𝜀𝐶₁𝜕⁴𝜕𝑥²𝜕𝑥²₁,𝑁₂=𝐴^′₁𝐻₁𝐻_1𝑥+𝜀𝐵^′₁𝐻₁𝐻_1𝑥𝑥+𝜀𝐶^′₁𝐻₁𝐻_{1𝑥𝑥𝑥𝑥}+𝜀𝐵^′₁𝐻²_1𝑥+𝜀𝐶^′₁𝐻_1𝑥𝐻_{1𝑥𝑥𝑥},𝑁₃=𝐴^′₁𝐻₁𝐻_2𝑥+𝐻₁𝐻_1𝑥1+𝐻₂𝐻_1𝑥+𝜀𝐵^′₁𝐻₁𝐻_2𝑥𝑥+2𝐻₁𝐻_1𝑥𝑥1+𝐻₂𝐻_1𝑥𝑥+𝜀𝐶^′₁𝐻₁𝐻_{2𝑥𝑥𝑥𝑥}+4𝐻₁𝐻_{1𝑥𝑥𝑥𝑥1}+𝐻₂𝐻_{1𝑥𝑥𝑥𝑥}+𝜀𝐵^′₁2𝐻_1𝑥𝐻_2𝑥+2𝐻_1𝑥𝐻_1𝑥1+𝜀𝐶^′₁𝐻_{2𝑥𝑥𝑥}𝐻_1𝑥+3𝐻_{1𝑥𝑥𝑥1}𝐻_1𝑥+𝐻_{1𝑥𝑥𝑥}𝐻_2𝑥+𝐻_{1𝑥𝑥𝑥}𝐻_1𝑥1+𝐴₁^′′2𝐻²₁𝐻_1𝑥𝐵+𝜀^1′′2𝐻²₁𝐻_1𝑥𝑥𝐶+𝜀^1′′2𝐻²₁𝐻_{1𝑥𝑥𝑥𝑥}+𝜀𝐵₁^′′𝐻₁𝐻²_1𝑥+𝜀𝐶₁^′′𝐻₁𝐻_1𝑥𝐻_{1𝑥𝑥𝑥}. (B.3a)(B.3b)(B.3c)(B.3d)(B.3e)The equations are solved order-by-order up to 𝑂(𝛼³) [12, 14, 56, 57] to arrive at an equation related to the complex Ginzburg-Landau type for the perturbation amplitude 𝜂(𝑥₁,𝑡₁,𝑡₂): 𝜕𝜂𝜕𝑡₂−𝛼⁻²𝐶_𝑅𝜂+𝜀𝐵₁−6𝑘²𝐶₁𝜕²𝜂𝜕𝑥²₁+𝐽₂+𝑖𝐽₄𝜂²𝜂=0, (B.4) where 𝑒_𝑟=2𝐵^′₁−𝑘²𝐶^′₁−4𝐵₁+16𝑘²𝐶₁,𝑒_𝑖=−𝐴^′₁𝜀−4𝑘𝐵₁+16𝑘³𝐶₁,𝐽₂=𝜀7𝑘⁴𝑒_𝑟𝐶^′₁−𝑘²𝑒_𝑟𝐵^′₁−𝑘²2𝐵₁^′′+𝑘⁴2𝐶₁^′′−𝑘𝑒_𝑖𝐴^′₁,𝐽₄=𝜀7𝑘⁴𝑒_𝑖𝐶^′₁−𝑘²𝑒_𝑖𝐵^′₁+𝑘𝑒_𝑟𝐴^′₁𝐴+𝑘^1′′2. (B.5) It should be remarked that in (B.4) an imaginary, apparently convective term of the form 𝑖𝑣𝜕𝜂𝜕𝑥₁=2𝑖𝑘𝜀𝐵₁−2𝑘²𝐶₁𝛼⁻¹𝜕𝜂𝜕𝑥₁ (B.6) could be additionally considered on the left-hand side of (B.4) as done in Dandapat and Samanta [12] and Samanta [65], provided that such a term is of 𝒪(𝛼) and 𝑣<0. Such a term arises in the secular condition at 𝑂(𝛼³) as a contribution arising from one of the terms from 𝐿₁𝐻₁ at 𝑂(𝛼²). The solution of (B.4) for a filtered wave in which the spatial modulation does not exist and the diffusion term in (B.4) becomes zero is obtained by considering 𝜂=Γ₀(𝑡₂)𝑒^{[−𝑖𝑏(𝑡2)𝑡2]}. This leads to a nonlinear ordinary differential equation for Γ₀, namely, 𝑑Γ₀𝑑𝑡₂−𝑖Γ₀𝑑𝑑𝑡₂𝑏𝑡₂𝑡₂−𝛼⁻²𝐶_𝑅Γ₀+𝐽₂+𝑖𝐽₄Γ³₀=0, (B.7) where the real and imaginary parts, when separated from (37), give 𝑑Γ₀𝑡₂𝑑𝑡₂=𝛼⁻²𝐶_𝑅−𝐽₂Γ²₀Γ₀,𝑑𝑏𝑡₂𝑡₂𝑑𝑡₂=𝐽₄Γ²₀. (B.8)(B.9)

If 𝐽₂=0 in (B.8), the Ginzburg-Landau equation is reduced to a linear ordinary differential equation for the amplitude of a filtered wave. The second term on the right-hand side in (B.8) induced by the effect of nonlinearity can either accelerate or decelerate the exponential growth of the linear disturbance depending on the signs of 𝐶_𝑅 and 𝐽₂. The perturbed wave speed caused by the infinitesimal disturbances appearing in the nonlinear system can be modified using (B.8). The threshold amplitude 𝛼Γ₀ from (B.8) is obtained as (when Γ₀ is nonzero and independent of 𝑡₂) 𝛼Γ₀=𝐶_𝑅𝐽₂. (B.10) Equation (B.9) with the use of (B.10) gives 𝐶𝑏=^𝑅𝐽₄𝛼²𝐽₂. (B.11) The speed 𝑁_𝑐 of the nonlinear wave is now obtained from (26) by substituting 𝜂=Γ₀𝑒^{[−𝑖𝑏𝑡2]} and by using 𝑡₂=𝛼²𝑡. This gives the nonlinear wave speed as 𝑁_𝑐=𝐶_𝐿𝑉+𝛼²𝑏=𝐶_𝐿𝑉+𝐶_𝑅𝐽₄𝐽₂. (B.12)

C. Sideband Stability

The sideband instability is analyzed by subjecting Γ_𝑠(𝑡₂) to sideband disturbances in the form of (Δ≪1) [12, 56, 57] 𝜂=Γ_𝑠𝑡₂+ΔΓ₊𝑡₂𝑒^𝑖𝑘𝑥+ΔΓ₋𝑡₂𝑒^{−𝑖𝑘𝑥}𝑒^{−𝑖𝑏𝑡2}, (C.1) in (B.4). Neglecting the nonlinear terms and separating the coefficients of Δ𝑒^{𝑖(𝑘𝑥−𝑏𝑡2)} and Δ𝑒^{−𝑖(𝑘𝑥−𝑏𝑡2)}, the following set of equations are derived (𝑣=(2𝜀𝑘/𝛼)(𝐵₁−2𝑘²𝐶₁)<0, 𝑐_𝑖=𝛼⁻²𝐶_𝑅, and 𝐽_1𝑟=𝜀(𝐵₁−6𝑘²𝐶₁)<0 since 𝐵₁−𝑘²𝐶₁=𝒪(𝛼²)):𝑑Γ₊𝑡₂𝑑𝑡₂=𝑖𝒫+𝑘𝑣+𝑐_𝑖+𝑘²𝐽_1𝑟−2𝐽₂+𝑖𝐽₄Γ_𝑠²×Γ₊𝑡₂−𝐽₂+𝑖𝐽₄Γ_𝑠²Γ₋𝑡₂,𝑑Γ₋𝑡₂𝑑𝑡₂=𝑖𝒫−𝑘𝑣+𝑐_𝑖+𝑘²𝐽_1𝑟−2𝐽₂+𝑖𝐽₄Γ_𝑠²×Γ₋𝑡₂−𝐽₂+𝑖𝐽₄Γ_𝑠²Γ₊𝑡₂, (C.2a)(C.2b)where the barred quantities are the complex conjugates corresponding to their counterparts. For convenience, the above system is represented as 𝑑Γ₊𝑡₂𝑑𝑡₂𝑑Γ₋𝑡₂𝑑𝑡₂Γ=ℳ⁺𝑡₂Γ₋𝑡₂, (C.3) where ℳ is a 2 × 2 matrix whose entries are 𝑀₁₁=𝑖𝒫+𝑘𝑣+𝑐_𝑖+𝑘²𝐽_1𝑟−2𝐽₂+𝑖𝐽₄Γ_𝑠²,𝑀₂₂=−𝑖𝒫−𝑘𝑣+𝑐_𝑖+𝑘²𝐽_1𝑟−2𝐽₂−𝑖𝐽₄Γ_𝑠²,𝑀₁₂=−𝐽₂+𝑖𝐽₄Γ_𝑠²,𝑀₂₁=𝑀₁₂. (C.4)

Setting Γ₊(𝑡₂)∝𝑒^𝜗𝑡2 and Γ₋(𝑡₂)∝𝑒^𝜗𝑡2, the eigenvalues are obtained from |ℳ−𝜗𝐼|=0, and the stability of (22) subject to sideband disturbances as 𝑡₂→∞ demands that 𝜗<0. The eigenvalues are found as follows 1𝜗=2tr(ℳ)±tr²(ℳ)−4det(ℳ). (C.5)

D. Bimodal Dynamical System

In order to analyze the stability of the bimodal dynamical system, the film thickness is expanded as a truncated Fourier series (𝑧_𝑛(𝑡) is complex, and a bar above it designates its complex conjugate) ℎ(𝑥,𝑡)=1+²_𝑛=1𝑧_𝑛(𝑡)𝑒^{𝑛𝑖𝑘𝑥}+𝑧_𝑛(𝑡)𝑒^{−𝑛𝑖𝑘𝑥} (D.1) and is substituted in the evolution equation (10). A Taylor series expansion is then sought to be about “1”, and the following coupled dynamical system is obtained from the coefficients of 𝑒^𝑖𝑘𝑥 and 𝑒^{2𝑖𝑘𝑥} terms:𝑑𝑧𝑑𝑡¹(𝑡)=𝑎₁₁₁𝑧₁(𝑡)+𝑎₁₂₁𝑧₁(𝑡)𝑧₂(𝑡)+𝑧₁(𝑡)𝑎₁₃₁𝑧₁²+𝑎₁₃₂𝑧₂²,𝑑𝑧𝑑𝑡²(𝑡)=𝑎₂₁₁𝑧₂(𝑡)+𝑎₂₂₁𝑧²₁(𝑡)+𝑧₂(𝑡)𝑎₂₃₁𝑧₁²+𝑎₂₃₂𝑧₂². (D.2a)(D.2b)It should be remarked that the higher-order terms arising while deducing (D.2a) and (D.2b) are neglected by assuming that |𝑧₁|≪1 and |𝑧₂|≪|𝑧₁|. The coefficients of the above system are given by 𝑎₁₁₁=𝜀𝑘²𝐵₁−𝑘⁴𝐶₁−𝑖𝑘𝐴₁,𝑎₂₁₁=𝜀4𝑘²𝐵₁−16𝑘⁴𝐶₁−2𝑖𝑘𝐴₁,𝑎₁₂₁=𝜀𝑘²𝐵^′₁−7𝑘⁴𝐶^′₁−𝑖𝑘𝐴^′₁,𝑎₂₂₁=𝜀2𝑘²𝐵^′₁−2𝑘⁴𝐶^′₁−𝑖𝑘𝐴^′₁,𝑎₁₃₁=𝜀2𝑘²𝐵₁^′′−𝑘⁴𝐶₁^′′−𝑖2𝑘𝐴₁^′′,𝑎₂₃₁=𝜀4𝑘²𝐵₁^′′−16𝑘⁴𝐶₁^′′−2𝑖𝑘𝐴₁^′′,𝑎₁₃₂=𝜀𝑘²𝐵₁^′′−𝑘⁴𝐶₁^′′−𝑖𝑘𝐴₁^′′,𝑎₂₃₂=𝜀2𝑘²𝐵₁^′′−8𝑘⁴𝐶₁^′′−𝑖𝑘𝐴₁^′′. (D.3)

The complex conjugates occurring in system (D.3) can be avoided by representing 𝑧_𝑛(𝑡)=𝑎_𝑛(𝑡)𝑒^{𝑖𝜃𝑛(𝑡)} in polar form, where 𝑎_𝑛(𝑡) and 𝜃_𝑛(𝑡) are real functions such that 𝑛=1,2. Using polar notation, the following three equations are obtained:𝑑𝑎₁(𝑡)𝑑𝑡=𝑓₁𝑎₁(𝑡),𝑎₂(𝑡),𝜙(𝑡)=𝑏₁₁₁𝑎₁(𝑡)+𝑏₁₂₁cos𝜙(𝑡)+𝑐₁₂₁sin𝜙(𝑡)×𝑎₁(𝑡)𝑎₂(𝑡)+𝑎₁(𝑡)𝑏₁₃₁𝑎²₁(𝑡)+𝑏₁₃₂𝑎²₂(𝑡),𝑑𝑎₂(𝑡)𝑑𝑡=𝑓₂𝑎₁(𝑡),𝑎₂(𝑡),𝜙(𝑡)=𝑏₂₁₁𝑎₂(𝑡)+𝑏₂₂₁cos𝜙(𝑡)−𝑐₂₂₁sin𝜙(𝑡)×𝑎²₁(𝑡)+𝑎₂(𝑡)𝑏₂₃₁𝑎²₁(𝑡)+𝑏₂₃₂𝑎²₂(𝑡,(𝑡))𝑑𝜙𝑑𝑡=𝑓₃𝑎₁(𝑡),𝑎₂(𝑡),𝜙(𝑡)=𝑘𝐴₁^′′𝑎²₁−𝑎²₂+2𝑐₁₂₁cos𝜙(𝑡)−𝑏₁₂₁sin𝜙(𝑡)×𝑎₂(𝑡)−𝑏₂₂₁sin𝜙(𝑡)+𝑐₂₂₁(𝑡𝑎cos𝜙)²¹(𝑡)𝑎₂(𝑡). (D.4a)(D.4b)(D.4c)

While the coefficient of 𝑒^{𝑖𝜃𝑛(𝑡)} term yields (D.4a) and (D.4b), (D.4c) is obtained from the two equations which result as the coefficient of 𝑖𝑒^{𝑖𝜃𝑛(𝑡)} term with the imposition that the phase relationship is defined as 𝜙(𝑡)=2𝜃₁(𝑡)−𝜃₂(𝑡). The terms 𝑏_𝑖𝑗𝑘 and 𝑐_𝑖𝑗𝑘 are real numbers and represent the real and imaginary parts of a typical 𝑎_𝑖𝑗𝑘 term in (D.3) such that 𝑖,𝑘=1,2 and 𝑗=1,2,3.