Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1323 https://internationalpubls.com Nonlinear Hydrodynamic Instability Control on a Shear-Imposed Two-Layer Film Flow by Means of Velocity Slip and Insoluble Surfactant 1Muhammad Sani, 2S. Karthick*, 1Iliyasu Bello, 1James Andrawus, 3G. Jayaraman . 1 Department of Mathematics, Faculty of Physical Sciences, Federal University Dutse-7156, Jigawa State, Nigeria. msdjunior9@gmail.com, iliyasubello13@gmail.com, jamesandy7772@gmail.com 2* Department of Mathematics, College of Engineering and Technology, SRM Institute of Science and Technology, Kattankulathur-603203, Tamilnadu, India. karthickmaths007@gmail.com 3 Department of Mathematics, Vels Institute of Science, Technology and Advanced Studies, Chennai, Tamilnadu, India. jayaram07maths@gmail.com Corresponding Author: karthickmaths007@gmail.com Article History: Received: 12-01-2025 Revised: 15-02-2025 Accepted: 01-03-2025 Abstract: The impact of an insoluble surfactant on the free surface of a two-layer fluid flowing down an inclined slippery plane is analyzed under external shear stress. The flow dynamics in both layers are governed by the Navier-Stokes equations, with shear stress applied to the upper layer. Using perturbation and normal-mode analysis, the governing equations and boundary conditions are transformed into the Orr-Sommerfeld system and solved via the spectral collocation method. The study reveals the presence of an unstable Yih mode due to the free surface. This instability is suppressed by the Marangoni effect, which arises from the surfactant-induced surface tension gradient. The Marangoni mode stabilizes the interface by counteracting perturbations. However, as the external shear stress increases, the stabilizing effect of the Marangoni mode weakens. Beyond a certain threshold, the Marangoni mode destabilizes, leading to interfacial disturbances. This study highlights the interplay between surfactants, shear stress, and hydrodynamic stability. The findings have applications in coating processes, microfluidics, and geophysical fluid dynamics. Keywords: Hydrodynamic stability, Two-layer fluid, Insoluble surfactant, Orr- Sommerfeld, Marangoni number, Shear stress. AMS Subject Classification: 76E05, 76E06, 76E09, 76E25, 76E30 1. Introduction The hydrodynamic instability of falling flow models is a topic of much interest because of its wide range of applications in the chemical and coating industries [6, 26], as well as in countless other industrial and natural settings [27, 31]. In order to increase the heat and mass transfer rates in process equipment, the surface wave is essential. Different dynamical features, referred to as free surface instabilities, are observed on the surface of a liquid film [19]. In recent decades, numerous studies have been carried out to reveal the underlying physical mechanisms of these instabilities and to offer various explicit and implicit management strategies to either enhance or mitigate single- or double-layer free surface flow instabilities, depending on the applications [24, 25]. Controlling the surface wave instability of fluid flows is required in a number of real-world situations. Controlling the surface wave Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1324 https://internationalpubls.com instability of fluid flows is required in a number of real-world situations. These situations can be addressed either by surface-active or surface-inactive agents depending on the use [1, 4, 17]. For instance, consistent film flows with superposed layers are essential in coating technology to guarantee the quality of photographic films [26]. The literature contains a vast number of studies offering different control tactics for surface wave instability of single- or multi-layer laminar, Newtonian, or non-Newtonian fluid flows over an inclined bed [20, 28]. Adding an insoluble surfactant to the fluid surface as a surface-active agent is one method of managing the instability [1]. By resolving the associated Orr-Sommerfeld boundary value problem, Blyth and Pozrikidis [4] discovered that the local variation in surface tension caused by the surfactant agent is the source of the surfactant mode, also known as the Marangoni mode, connected with the surfactant. Faruq Ahmad Bhat et al. [2] demonstrated the influence of surfactants on the stability of the Yih mode using an insoluble surfactant at the single-layer fluid surface. For a number of applications, this is highly significant in the scientific community [8]. In addition to the instability management method mentioned above, another control technique to stabilize the most unstable modes is to insert a floating flexible membrane or plate at the top surface. In accordance with the aforementioned concepts, Selvan et al. [28] examined the use of a floating, flexible plate to control hydrodynamic instability in both linear and nonlinear regimes at the surface of a gravity-viscous liquid film over a rigid substrate. The authors observed that the instability mechanism of the floating elastic membrane is the same as the Marangoni mechanism of a single-layer falling liquid film polluted by surfactants. Their research verified that, in situations where the compressive force is either absent or present, the flexible plate on a falling fluid surface might be employed as a substitute for passive or active control. The Orr-Sommerfeld eigenvalue problem obtained by the authors was solved numerically using the Chebyshev collocation method, and it was established that the floating elastic plate suppresses the unstable Yih mode [32], which is connected to the disturbance wave superimposed on the liquid surface. Physically, by decreasing the flow velocity in the boundary layer beneath the plate, the elastic membrane prevents the laminar flow from transitioning to turbulence. Moreover, imposing an external shear at the free surface is another method for managing the energy of surface waves [20]. Physically, this procedure can be used in the airway occlusion process [16]. Another active mechanism for controlling the Yih mode different from the ones mentioned above is by imposing shear at the free surface, the influence of the shear helps in the dissipation of energy at the free surface. Because of this, Samanta [21] investigated how external shear affected the falling film caused by gravity down a sharp slope and flow-directed external force which lowers the essential Reynolds number in order to encourage long wave instability. Subsequently, an external force was applied to the top surface of a surfactant-laden single-layer fluid film flowing down a rigid bed, as studied by Samanta [20]. The authors examined how wave energy created by external shear was attenuated using an insoluble surfactant at the fluid surface. In addition to the surface and surfactant modes, they discovered the existence of a third mode, called the shear mode. The onset of the shear mode is primarily caused by the inclination angle and large inertia force. It was determined that strong stabilization of the surface mode is possible if the applied force acts in the reverse flow direction with a stronger Marangoni force [18]. Recently, researchers have shown increasing interest in studying interfacial instability in double- or multi-layer liquid films influenced by kinetic factors. In a single-layer falling film over an inclined wall, instability mainly arises due to inertial effects, leading to surface wave formation. In multi-layer films, additional factors such as viscosity contrast, density differences, and interfacial tension play a significant role. These complexities can either amplify or suppress disturbances, depending on the flow conditions. Understanding these instabilities is crucial for applications in coating technologies, microfluidics, and industrial fluid processes [11, 12]. Concurrently, due to its widespread applications in various technical and chemical processes, such as liquid-liquid extraction, food manufacturing, and liquid coating, several researchers have shown interest in studying instability control mechanisms in double-layer fluid flow systems [10, 29, 30]. In Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1325 https://internationalpubls.com coating processes and related applications, interfacial and free surface instabilities can arise under specific flow conditions, potentially compromising the quality of the final product. Effective control of these instabilities is therefore crucial for ensuring consistency, precision, and overall product reliability. Gao and Lu [14] examined the effect of an insoluble surfactant on the free surface and interface in the inertialess limit while studying a two-layer system. They identified four modes but primarily focused on the surface and interface modes. Their study revealed that longwave region inertialess instability could be either suppressed or amplified by surface or interfacial surfactants. Furthermore, two-layer flows are widely relevant as they occur in various natural settings and play a crucial role in biological and technical systems, including the plastics and petroleum industries [7, 13, 15]. By analyzing the effects of depth ratio, viscosity ratio, and density ratio on the fundamental instability produced by surface and interface modes, Kao [13], first introduced the concept of double- layered fluid flow. Subsequently, Chen [5] investigated the linear stability of a double-layered falling film at low Reynolds numbers. The author observed that viscosity stratification plays a crucial role in the interface mode and confirmed that instability occurs when the viscosity of the lower layer is greater than that of the overlying layer. Samanta [19] added inertia force up to moderate Reynolds number values, extending the work of Gao and Lu [14] even further, investigated how the free surface and interface mode-induced primary instability was affected by insoluble surfactant. It was demonstrated that when the top layer is less viscous than the bottom layer, the interfacial instability weakens due to the presence of the interfacial surfactant. Bhat and Samanta[3] further examined the effects of surfactants on two-layer falling films over a slippery incline, highlighting the dual role of the viscosity ratio in surface mode instability. Additionally, studies with various flow configurations have explored the use of external shear as a surface-active agent to control instabilities caused by free surface waves [9, 22, 23]. It is worth noting that several applications of heterogeneous fluid mostly deal with liquid- liquid extraction. However, there are few applications related to large airways and bolus-dispersal surfactant replacement therapy. In viewing the respiratory tract at the microscopic level, there is a thin layer of membrane that lies between the loosely bound connective tissue and mucus layer, which behaves like a two-layer fluid in the presence of insoluble surfactant. In the present study, the instability control mechanism by insoluble surfactant on the two-layer fluid due to the imposed shear on the airway track is analyzed. The Orr-Sommerfeld system obtained for the governing equations and boundary conditions is solved numerically using the spectral collocation method. The numerical results such as dispersion and neutral stability curves are analyzed extensively for different surfactant concentration, shear, and slip parameters for smaller Reynolds numbers. The paper is organized as follows: the mathematical formulation along with the solution of the physical problem is discussed in Section 2, the numerical results are analyzed in Section 3, and the concluding remark is highlighted in Section 4 2. Mathematical Formulation A two-dimensional shear-imposed flow of immiscible, incompressible and Newtonian fluids having different density and viscosity in an inclined slippery plane is considered as in Fig. 1. The upper layer is denoted as fluid-1 with density 𝜌1 and viscosity πœ‡1 , whereas the lower layer is denoted as fluid-2 with density 𝜌2 and viscosity πœ‡2. The distance between the slippery bottom and mean free surface is 𝐻 and width of the upper layer fluid is 𝑑. The origin is assumed at the mean free surface as 𝑂 and the positive 𝑦-axis is pointing downward perpendicular to the flow direction. 𝜏 represent the external share pointing in the opposite direction of the fluid flow. Moreover, β„Ž1(π‘₯, 𝑑) and β„Ž2(π‘₯, 𝑑) are denoted as the elevation of upper and lower layers, respectively. Further, the angle between the flat and tilted bottom is denoted as πœƒ. The governing equations of fluid flow in both layers are given by the two-dimensional Navier-Stokes equations, Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1326 https://internationalpubls.com πœ•π‘₯𝑒𝑗 + πœ•π‘¦π‘£π‘— = 0, (2.1) πœ•π‘‘π‘’π‘— + π‘’π‘—πœ•π‘₯𝑒𝑗 + π‘£π‘—πœ•π‘¦π‘’π‘— = βˆ’ 1 πœŒπ‘— πœ•π‘₯𝑝𝑗 + πœˆπ‘—[πœ•π‘₯π‘₯ 2 𝑒𝑗 + πœ•π‘¦π‘¦ 2 𝑒𝑗] + 𝑔sinπœƒ, (2.2) πœ•π‘‘π‘£π‘— + π‘’π‘—πœ•π‘₯𝑣𝑗 + π‘£π‘—πœ•π‘¦π‘£π‘— = βˆ’ 1 πœŒπ‘— πœ•π‘¦π‘π‘— + πœˆπ‘—[πœ•π‘₯π‘₯ 2 𝑣𝑗 + πœ•π‘¦π‘¦ 2 𝑣𝑗] βˆ’ 𝑔cosπœƒ, (2.3) Figure 1: Schematic diagram for a two-layer thin film flow down in an inclined slippery plane in the presence of an insoluble surfactant at the free surface. where 𝑗 = 1 and 𝑗 = 2 denote the upper and lower layer fluids respectively, and 𝑔 be the acceleration due to gravity. In the upper layer, the surface is covered by the monolayer insoluble surfactant with concentration Ξ“(π‘₯, 𝑑), which diffuses on the free surface. But, there is no much penetration in the bulk fluid. Now, the local surface tension on the free surface (𝜎1) varies with respect to the local surfactant concentration (πœŽπ‘). The relation between the surface concentration and tension (Gibbs relation) is given by, πœŽπ‘ βˆ’ 𝜎1 = Ξ“ 𝑅𝑇, (2.4) where and are the ideal gas constant and absolute temperature respectively. Thus, the linear equation is given by 𝜎1 = πœŽπ‘ (1 βˆ’ 𝛾 Ξ“ Ξ“0 ), (2.5) where 𝛾 = Ξ“0𝑅𝑇/πœŽπ‘, which is related to the surface elasticity 𝐸 = π›ΎπœŽπ‘/Ξ“0. The reference surfactant concentration and surface tension are denoted as Ξ“0 and 𝜎0, respectively. Using the Taylor series expansion about Ξ“0, the local surface tension is given as 𝜎1 = 𝜎0 βˆ’ 𝐸(Ξ“ βˆ’ Ξ“0). (2.6) The boundary conditions associated with the upper layer fluid such as kinematic condition, tangential and normal stress balances, are given by 𝑣1 = πœ•π‘‘β„Ž1 + 𝑒1πœ•π‘₯β„Ž1 at 𝑦 = β„Ž 1(π‘₯, 𝑑), (2.7) πœ‡1[βˆ’4πœ•π‘₯𝑒1πœ•π‘₯β„Ž1 + [πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1][1 βˆ’ (πœ•π‘₯β„Ž1)2]] = [πœ•π‘₯𝜎1 + 𝜏] √1 + (πœ•π‘₯β„Ž1)2 at 𝑦 = β„Ž1(π‘₯, 𝑑), (2.8) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1327 https://internationalpubls.com 𝑝1 = 2πœ‡1 [1+(πœ•π‘₯β„Ž1)2] [πœ•π‘₯ 𝑒1 (πœ•π‘₯ β„Ž1)2 βˆ’ (πœ•π‘¦ 𝑒1 + πœ•π‘₯𝑣1) πœ•π‘₯ β„Ž1 + πœ•π‘¦ 𝑣1] βˆ’ 𝜎1 πœ•π‘₯π‘₯ 2 β„Ž1 [1+πœ•π‘₯ β„Ž1]3/2 at 𝑦 = β„Ž1(π‘₯, 𝑑), (2.9) where 𝜏 is the external shear imposed on the fluid flow. In the case of lower layer fluid, the kinematic condition along with the balance of velocities, tangential and normal stresses, are given by, 𝑒1 = 𝑒2 and 𝑣1 = 𝑣2 at 𝑦 = β„Ž2(π‘₯, 𝑑), (2.10) 𝑣2 = πœ•π‘‘ β„Ž2 + 𝑒2 πœ•π‘₯ β„Ž2 at 𝑦 = β„Ž 2(π‘₯, 𝑑), (2.11) πœ‡2 [(πœ•π‘¦ 𝑒2 + πœ•π‘₯ 𝑣2)(1 βˆ’ (πœ•π‘₯ β„Ž2)2) βˆ’ 4 πœ•π‘₯𝑒2 πœ•π‘₯β„Ž2] = πœ‡1 [(πœ•π‘¦ 𝑒1 + πœ•π‘₯ 𝑣1)(1 βˆ’ (πœ•π‘₯ β„Ž1)2) βˆ’ 4πœ•π‘₯ 𝑒1πœ•π‘₯ β„Ž1] at 𝑦 = β„Ž2(π‘₯, 𝑑), (2.12) 𝑝1 + 2 πœ‡1 [1+(πœ•π‘₯β„Ž1)2] [πœ•π‘₯𝑒1(1 βˆ’ (πœ•π‘₯β„Ž1)2) + πœ•π‘₯β„Ž1 (πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1)] + 𝜎2 πœ•π‘₯π‘₯ 2 β„Ž1 [1+(πœ•π‘₯β„Ž1)2]3/2 = 𝑝2 + 2 πœ‡2 [1+(πœ•π‘₯β„Ž2)2] [πœ•π‘₯𝑒2(1 βˆ’ (πœ•π‘₯β„Ž2)2) + (πœ•π‘¦π‘’2 + πœ•π‘₯𝑣2)] at 𝑦 = β„Ž2(π‘₯, 𝑑). (2.13) Where 𝜎2 represents the interfacial tension between the two layers. The boundary conditions at slippery bottom are given by, 𝑒2 = 𝛽 πœ•π‘¦ 𝑒2 and 𝑣2 = 0 at 𝑦 = 𝐻, (2.14) where 𝛽 be the slip parameter having range as 0 ≀ 𝛽 ≀ 1. The reference length and velocity scale are chosen as 𝐻 and 𝑉, respectively. Further, the chosen scale of velocity is the maximum value of velocity for an uniform Nusselt film over the rigid substrate, which is defined as 𝑉 = 𝑔𝐻2sinπœƒ/2𝜈1, with 𝜈1 being the kinematic viscosity of upper layer fluid. Now, the governing equations along with the boundary conditions are made dimensionless using the following non-dimensional variables, are given as οΏ½Μ…οΏ½ = π‘₯ 𝐻 , οΏ½Μ…οΏ½ = 𝑦 𝐻 , �̅�𝑗 = 𝑒𝑗 𝑉 , �̅�𝑗 = 𝑣𝑗 𝑉 , 𝑑̅ = 𝑑𝑉 𝐻 , �̅�𝑗 = 𝑝𝑗 πœŒπ‘—π‘‰2 , β„ŽΜ…π‘— = β„Žπ‘— 𝐻 , οΏ½Μ…οΏ½ = 𝑑 𝐻 , Ξ“Μ… = Ξ“ Ξ“0 , �̅�𝑗 = πœŽπ‘— 𝜎0 . The set of non-dimensional form of governing equations and boundary conditions for the two-layer fluid system are given by πœ•π‘₯𝑒𝑗 + πœ•π‘¦π‘£π‘— = 0, (2.15) πœ•π‘‘π‘’π‘— + 𝑒𝑗 πœ•π‘₯𝑒𝑗 + 𝑣𝑗 πœ•π‘¦π‘’π‘— = βˆ’πœ•π‘₯𝑝𝑗 + 1 𝑅𝑒𝑗 (πœ•π‘₯π‘₯ 2 𝑒𝑗 + πœ•π‘¦π‘¦ 2 𝑒𝑗) + 𝐺, (2.16) πœ•π‘‘π‘£π‘— + π‘’π‘—πœ•π‘₯𝑣𝑗 + π‘£π‘—πœ•π‘¦π‘£π‘— = βˆ’πœ•π‘₯𝑝𝑗 + 1 𝑅𝑒𝑗 (πœ•π‘₯π‘₯ 2 𝑣𝑗 + πœ•π‘¦π‘¦ 2 𝑣𝑗) βˆ’ 𝐺 cotπœƒ, (2.17) 𝑣𝑗 = πœ•π‘‘β„Žπ‘— + π‘’π‘—πœ•π‘₯β„Žπ‘— at 𝑦 = β„Žπ‘— (π‘₯, 𝑑), (2.18) βˆ’4πœ•π‘₯𝑒1πœ•π‘₯β„Ž1 + (πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1)(1 βˆ’ (πœ•π‘₯β„Ž1)2) = [βˆ’ π‘€π‘Ž πΆπ‘Ž πœ•π‘₯Ξ“ 𝜏] √1 + (πœ•π‘₯β„Ž1)2 at 𝑦 = β„Ž1(π‘₯, 𝑑), (2.19) 𝑅𝑒 𝑝1 = 2 (1+(πœ•π‘₯β„Ž1)2) [πœ•π‘₯𝑒1 (πœ•π‘₯β„Ž1)2 βˆ’ (πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1) πœ•π‘₯β„Ž1 + πœ•π‘₯𝑣1] βˆ’ [1βˆ’π‘€π‘Ž(Ξ“βˆ’1)] πœ•π‘₯π‘₯ 2 β„Ž1 πΆπ‘Ž [1+(πœ•π‘₯β„Ž)2]3/2 at 𝑦 = β„Ž1(π‘₯, 𝑑), (2.20) πœ•π‘‘Ξ“ + 𝑒1 πœ•π‘‘Ξ“ + Ξ“πœ•π‘₯𝑒1 = 1 𝛼 πΆπ‘Ž πœ•π‘₯π‘₯ 2 Ξ“ at 𝑦 = β„Ž1(π‘₯, 𝑑), (2.21) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1328 https://internationalpubls.com 𝑒1 = 𝑒2 and 𝑣1 = 𝑣2 at 𝑦 = β„Ž2(π‘₯, 𝑑), (2.22) βˆ’4 πœ•π‘₯𝑒1 πœ•π‘₯β„Ž1 + (πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1)(1 βˆ’ (πœ•π‘₯β„Ž1)2) = π‘š [βˆ’4 πœ•π‘₯𝑒2 πœ•π‘₯β„Ž2 + (πœ•π‘¦π‘’2 + πœ•π‘₯𝑣2)(1 βˆ’ (πœ•π‘₯β„Ž2)2)] at 𝑦 = β„Ž2(π‘₯, 𝑑), (2.23) 𝑝1 + 2 𝑅𝑒1 (1+(πœ•π‘₯β„Ž1)2) [(πœ•π‘¦π‘’1 + πœ•π‘₯𝑣1) πœ•π‘₯β„Ž1 βˆ’ πœ•π‘₯𝑒1 (πœ•π‘₯β„Ž1)2 βˆ’ πœ•π‘¦π‘£1] = π‘Ÿ 𝑝2 + 2 𝑅𝑒1 (1+(πœ•π‘₯β„Ž2)2) [(πœ•π‘¦π‘’2 + πœ•π‘₯𝑣2) πœ•π‘₯β„Ž2 βˆ’ πœ•π‘₯𝑒2 (πœ•π‘₯β„Ž2)2 βˆ’ πœ•π‘¦π‘£2] βˆ’ π‘Šπ‘’ πœ•π‘₯π‘₯ 2 β„Ž2 [1+(πœ•π‘₯β„Ž2)2] at 𝑦 = β„Ž2(π‘₯, 𝑑) (2.24) 𝑒2 = 𝛽 πœ•π‘¦π‘’2 and 𝑣2 = 0 at 𝑦 = 𝐻, (2.25) where π‘š = πœ‡2/πœ‡1 and π‘Ÿ = 𝜌2/𝜌1 be the viscosity and density ratios. The other parameters are given as follows: Marangoninumber π‘€π‘Ž = 𝐸Γ0 𝜎0 , Reynoldsnumber 𝑅𝑒1 = 𝑉𝐻 𝜈1 , Capillarynumber πΆπ‘Ž = π‘‰πœ‡1 𝜎0 , Surfactantsurfaceconstant 𝛼 = 𝜎0𝐻 πœ‡1 , Webernumber π‘Šπ‘’ = 𝜎2 𝜌 𝑔𝐻2. 2.1 Base solution and stability equation The base state solution corresponding to the two-layer fluid system in the case of velocity and pressure are obtained by assuming 𝑒1 = (π‘ˆ1𝐡(𝑦),0), 𝑒2 = (π‘ˆ2𝐡(𝑦),0), 𝑝1 = (𝑃1𝐡(𝑦),0) and 𝑝2 = (𝑃2𝐡(𝑦),0), where π‘ˆ1𝐡(𝑦), 𝑃1𝐡(𝑦), π‘ˆ2𝐡(𝑦), and 𝑃2𝐡(𝑦) denote the base state velocity and pressure in upper and lower layers. By applying these assumptions in the above system of equations (2.15)–(2.25) yield the following solutions π‘ˆ1𝐡(𝑦) = βˆ’π‘¦2 + 𝜏 𝑦 + [(2𝑑2 βˆ’ 𝜏 𝑑 βˆ’ 2𝑑) (1 βˆ’ 1 π‘š ) + (1 βˆ’ 𝜏 π‘š )] + 𝛽 [2𝑑 (1 βˆ’ 1 π‘š ) + 𝜏 π‘š βˆ’ 2], (2.26) π‘ˆ2𝐡(𝑦) = βˆ’π‘¦2 + [2𝑑 (1 βˆ’ 1 π‘š ) + 𝜏 π‘š ] 𝑦 + [1 βˆ’ 2𝑑 (1 βˆ’ 1 π‘š ) βˆ’ 𝜏 π‘š ] +𝛽 [2𝑑(1 βˆ’ 1 π‘š ) + 𝜏 π‘š βˆ’ 2], (2.27) 𝑃1𝐡(𝑦) = βˆ’ cosπœƒ πΉπ‘Ÿ2 𝑦, (2.28) 𝑃2𝐡(𝑦) = βˆ’ cosπœƒ πΉπ‘Ÿ2 𝑦 + ( π‘Ÿβˆ’1 π‘Ÿ ) cosπœƒ πΉπ‘Ÿ2 𝑑. (2.29) It is worth to be noted that the behavior of base flow velocity and pressure do not alter due to the presence of insoluble surfactant. However, the base velocity depends on the slip parameter and shear stress. By considering 𝛽 = 0 and 𝜏 = 0, the base profiles are similar to that in the case of two-layer film flow down the inclined plane in the absence of shear stress for rigid bottom. On the other hand, for 𝑑 = 0 and π‘š = 1, the base velocity profile in Eq. (2.26) is same as the velocity profile for homogeneous fluid flow in the presence of external shear stress. Further, an infinitesimal disturbance is applied on the system and the physical quantities can be written as the sum of initial and perturbed state. Substituting 𝑒𝑗(π‘₯, 𝑦, 𝑑) = π‘ˆπ‘—π΅(𝑦) + �̃�𝑗(π‘₯, 𝑦, 𝑑), 𝑣𝑗(π‘₯, 𝑦, 𝑑) = �̃�𝑗(π‘₯, 𝑦, 𝑑), 𝑝𝑗(π‘₯, 𝑦, 𝑑) = 𝑃𝑗𝐡(𝑦) + 𝑝𝑗(π‘₯, 𝑦, 𝑑), β„Ž1(π‘₯, 𝑑) = β„ŽΜƒ1(π‘₯, 𝑑), β„Ž2(π‘₯, 𝑑) = 𝑑 Β± β„ŽΜƒ2(π‘₯, 𝑑) and Ξ“(π‘₯, 𝑑) = 1 + Ξ“Μƒ in the governing equations along with the boundary conditions. Let the normal solution for the two-dimensional two-layer fluid flows are given as Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1329 https://internationalpubls.com �̃�𝑗(π‘₯, 𝑦, 𝑑) = πœ™π‘—(𝑦) eiπ‘˜(π‘₯βˆ’π‘π‘‘), (2.30) Ξ“Μƒ(π‘₯, 𝑦, 𝑑) = Ξ“(𝑦) eiπ‘˜(π‘₯βˆ’π‘π‘‘), (2.31) β„ŽΜƒπ‘—(π‘₯, 𝑦, 𝑑) = πœ‚π‘—(𝑦) eiπ‘˜(π‘₯βˆ’π‘π‘‘), (2.32) where π‘˜, 𝑐 and �̃�𝑗 are the wavenumber, complex wave speed and stream function, respectively and i ≑ βˆšβˆ’1. Using Eqs. (2.30)–(2.32) in the governing equations and boundary conditions (2.15)–(2.25) yield the following Orr-Sommerfeld system for both upper and lower layers, (𝐷2 βˆ’ π‘˜2)2 πœ™π‘— βˆ’ 𝑅𝑒𝑗[i π‘˜ (π‘ˆπ‘—π΅ βˆ’ 𝑐) (𝐷2 βˆ’ π‘˜2) πœ™π‘— βˆ’ iπ‘˜ π‘ˆβ€²β€²π‘—π΅ πœ™π‘—] = 0, (2.33) πœ™π‘— + πœ‚π‘— (π‘ˆπ‘—π΅ βˆ’ 𝑐) = 0 at 𝑦 = 0 (for 𝑗 = 1) and 𝑦 = 𝑑 (for 𝑗 = 2), (2.34) π·πœ™1 + Ξ“ [(π‘ˆ1𝐡 βˆ’ 𝑐) βˆ’ iπ‘˜ π›ΌπΆπ‘Ž ] = 0 at 𝑦 = 0, (2.35) π‘ˆβ€²β€²1𝐡 πœ‚1 + (𝐷2 + π‘˜2) πœ™1 + iπ‘˜π‘€π‘Ž πΆπ‘Ž Ξ“ = 0 at 𝑦 = 0, (2.36) 𝐷3πœ™1 + [𝑅𝑒1iπ‘˜(𝑐 βˆ’ π‘ˆ1𝐡) βˆ’ 3π‘˜2]π·πœ™1 βˆ’ iπ‘˜π‘…π‘’1π‘ˆβ€² 1π΅πœ™1 βˆ’[2iπ‘˜cotπœƒ + 2π‘ˆβ€²1𝐡 π‘˜2 + iπ‘˜3 πΆπ‘Ž ] πœ‚1 = 0 at 𝑦 = 0, (2.37) (𝐷2 + π‘˜2) πœ™1 = (π‘ˆβ€²β€²1𝐡 + π‘š π‘ˆβ€²β€²2𝐡) πœ‚2 + π‘š (𝐷2 + π‘˜2) πœ™2 at 𝑦 = 𝑑, (2.38) π·πœ™1 = [π‘ˆβ€²1𝐡 + π‘ˆβ€²2𝐡] πœ‚2 + π·πœ™2 at 𝑦 = 𝑑, (2.39) πœ™1 = πœ™2 at 𝑦 = 𝑑, (2.40) 2π‘˜2 π‘ˆβ€² 1𝐡 πœ‚1 + [3π‘˜2𝐷 βˆ’ 𝐷3 + iπ‘˜π‘…π‘’1((π‘ˆ1𝐡 βˆ’ 𝑐)𝐷 βˆ’ π‘ˆβ€² 1𝐡)] = (2icotπœƒ (π‘Ÿ + 1) + 2π‘šπ‘˜2π‘ˆβ€² 2𝐡 βˆ’ Iπ‘˜3 π‘Šπ‘’ 𝑅𝑒1) πœ‚2 + [π‘š (3π‘˜2𝐷 βˆ’ 𝐷3) + π‘Ÿπ‘…π‘’1 (iπ‘˜(π‘ˆ2𝐡 βˆ’ 𝑐) 𝐷 βˆ’ iπ‘˜π‘ˆβ€² 2𝐡)] πœ™2 at 𝑦 = 𝑑, (2.41) πœ™2 = 0 and π·πœ™2 = 𝛽 𝐷2πœ™2 at 𝑦 = 1, (2.42) where 𝐷 denote the derivative with respect to 𝑦. 3 Numerical Results Here, the spectral collocation method [33] is employed to solve the Orr-Sommerfeld system (Eqs.(2.33) – (2.42)) by MATLAB (2014a). The entire domain splits into the two sub-domain, which are discretized with 𝑛-points using the Chebyshev polynomials as given by 𝑦𝑗 = 0.5 (cos[πœ‹ (π‘—βˆ’1) (π‘›βˆ’1) ] + 1) where 𝑗 = 1,2,3, … 𝑛. In order to implement this method, the sub-domains such as [0, 𝑑] and [𝑑, 1] are transformed to the computational interval [βˆ’1,1] ([10]). Further, the complex phase velocity 𝑐 = π‘π‘Ÿ + i𝑐𝑖 is obtained from the Orr-Sommerfeld system, from which the non-dimensional scaled parameter known as growth rate (πœ”π‘– = π‘˜π‘π‘–) is derived. Using the special process, the extra eigenvalues are eliminated, which, in turns, increases the accuracy of obtained eigenvalues. In Fig. 2, the present result is validated with the result of Chin et al. [34] in the case of homogeneous fluid (𝑑 = 1, π‘š = 1, πΆπ‘Ž = 0.0169, 𝜏 = 0, 𝛽 = 0, πœƒ = 4∘ and 𝑅𝑒1 = 2500). The result depicts that the dispersion curve for single layer fluid (π‘˜π‘π‘– vs π‘˜) obtained from the present theory, matches precisely with the existing result of Chin et al. [34]. Moreover, the nature of growth rate is two-fold and becomes more stable for the wavenumber π‘˜ β‰₯ 1.5. In the case of the two-layer fluid, the dispersion curves for different interface position (𝑑) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1330 https://internationalpubls.com and viscosity ratio (π‘š) are plotted in Figs. 3(a) and (b), respectively. The patterns are two-fold, which are similar to the single-layer fluid flow as in Fig. 2. However, the long wave stability is more as Figure 2: Validation of free surface single layer dispersion curve with the existing theory of [34] by considering 𝒅 = 𝟏, π’Ž = 𝟏, π‘ͺ𝒂 = 𝟎. πŸŽπŸπŸ”πŸ—, 𝝉 = 𝟎, 𝜷 = 𝟎 and π‘Ήπ’†πŸ = πŸπŸ“πŸŽπŸŽ. compared to that for single-layer fluid i.e., the growth rate of dominant disturbance increases and attains maximum, then befall for the smaller values of π‘˜. While increasing the interface position (𝑑) as in Fig. 3(a), the growth rate increases. This result in more perturbed wave when the interface is closer to the slippery bottom. On the other hand, the growth rate increases for increase in the values of π‘š in Fig. 3(b) due to the presence of less viscous fluid in upper layer fluid, which destabilizes the flow. (a) (b) Figure 3: Growth rate against the wavenumber for different values of (a) interface position 𝒅 with π’Ž = 𝟎. πŸ– and (b) viscosity ratio π’Ž with 𝒅 = 𝟎. πŸ“, and the other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, 𝑴𝒂 = 𝟎. 𝟏, π‘Ήπ’†πŸ = πŸ“πŸŽ, 𝝉 = 𝟎. πŸ” and 𝜷 = 𝟎. πŸ“. The effect of insoluble surfactant concentration (π‘€π‘Ž) on the dispersion curve is discussed in Fig. 4(a). The patterns resemble as in Fig. 3, and it is observed that the presence of insoluble surfactant Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1331 https://internationalpubls.com decreases the instability. Further, the stability enhances for the increase in the concentration of insoluble surfactant π‘€π‘Ž. In the case of Fig. 4(b), the dispersion curves are plotted for different values of density ratio π‘Ÿ. As π‘Ÿ increases, the instability increases due to increase in the value of growth rate till π‘˜ ≀ 0.09. After that the growth rate follows the reverse pattern for increase in the values of π‘Ÿ. This is due to the rapid decline of the instability for higher wavenumbers π‘˜ > 0.09, but the rate of befall has no major deviations for larger values of π‘Ÿ. (a) (b) Figure 4: Growth rate against the wavenumber for different values of (a) Marangoni number 𝑴𝒂 with 𝒓 = 𝟎. πŸ” and (b) density ratio 𝒓 with 𝑴𝒂 = 𝟎. 𝟏, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, π‘Ήπ’†πŸ = πŸ“πŸŽ, 𝒅 = 𝟎. πŸ“, 𝝉 = 𝟎. πŸ” and 𝜷 = 𝟎. πŸ“. ( (a) (b) Figure 5: Growth rate against the wavenumber for different values of (a) shear parameter 𝝉 with 𝜷 = 𝟎. πŸ“ and (b) slip parameter 𝜷 with 𝝉 = 𝟎. πŸ”, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, π‘Ήπ’†πŸ = πŸ“πŸŽ, 𝒅 = 𝟎. πŸ“, 𝑴𝒂 = 𝟎. 𝟏 and 𝒓 = 𝟎. πŸ”. The effect of shear parameter (𝜏) and slip parameter (𝛽) on the dispersion curve for two-layer film flow are illustrated in Figs. 5(a) and (b), respectively. As the shear imposed on the upper layer Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1332 https://internationalpubls.com increases as in Fig. 5(a), the growth rate increases for wavenumber 0 ≀ π‘˜ ≀ 0.038 and it decreases for the value of π‘˜ > 0.038. This bifurcation is due to the fact that the perturbed wave amplitude exceeds the applied force for 0 ≀ π‘˜ ≀ 0.038 and it follows the reverse trends for π‘˜ > 0.038. On the other hand, in the case of Fig. 5(b), it is observed that the fluid flow destabilizes for larger values of 𝛽. Moreover, there is a negligible effect for wavenumber π‘˜ β‰₯ 0.01. (a) (b) Figure 6: Growth rate against the Reynolds number for different values of (a) Marangoni number 𝑴𝒂 with 𝒅 = 𝟎. πŸ“ and (b) interface position 𝒅 with 𝑴𝒂 = 𝟎. 𝟏, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, 𝒅 = 𝟎. πŸ“, 𝑴𝒂 = 𝟎. 𝟏 and 𝒓 = 𝟎. πŸ”. (a) (b) Figure 7: Growth rate against the Reynolds number for different values of (a) shear parameter 𝝉 with 𝜷 = 𝟎. πŸ“ and (b) slip parameter 𝜷 with 𝝉 = 𝟎. πŸ”, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, 𝒅 = 𝟎. πŸ“, 𝑴𝒂 = 𝟎. 𝟏 and 𝒓 = 𝟎. πŸ”. The inertial effect on the maximum growth rate is plotted for the different values of Marangoni number (π‘€π‘Ž) and interface position (𝑑) in Figs. 6(a) and (b), respectively. From both the figures, it is noticed that the growth rate increases and attains maximum, then decreases for increase in the values of 𝑅𝑒1. The fluid flow becomes stable for smaller values of 𝑅𝑒1, whereas it becomes unstable gradually for larger values of 𝑅𝑒1. This is due to the fact that the viscous force dominates the inertial force for Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1333 https://internationalpubls.com smaller values of 𝑅𝑒1 and vice versa. The maximum value of growth rate is due to the balance between the inertial and viscous force. Moreover, the inertial force eventually overcomes the viscous force for larger values of 𝑅𝑒1. From Fig. 6(a), it is observed that the growth rate decreases for increase in the concentration of surfactant. This implies that the higher concentration of surfactant enhances the stability of flow system. In the case of Fig. 6(b), the growth rate follows reverse pattern for increase in the values of 𝑑 against the Reynolds number 𝑅𝑒1. This is due to the fact that there is more instability as the interface approach closer to the slippery bottom. In Figs.7(a) and (b), the inertial effects on the maximum growth rate are plotted for different values of shear parameter (𝜏) and slip parameter (𝛽), respectively. The patterns are similar to Fig. 6. As the value of 𝜏 increases, the instability of two-layer fluid increases for 𝑅𝑒1 < 6 and it follows a reverse pattern for 𝑅𝑒1 > 6. This may be due to the stabilizing effect of shear stress for increase in the inertial force between the layers of fluid. In the case of Fig. 7(b), the instability of fluid flow increases due to increase in the velocity gradient at the bottom for larger values of 𝛽. (a) (b) Figure 8: Contours of growth rate for different values of (a) shear parameter 𝝉 with 𝜷 = 𝟎. πŸ“ and (b) slip parameter 𝜷 with 𝝉 = 𝟎. πŸ”, against the Reynolds number π‘Ήπ’†πŸ, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, 𝒅 = 𝟎. πŸ“, 𝑴𝒂 = 𝟎. 𝟏 and 𝒓 = 𝟎. πŸ”. Figs. 8(a) and (b) depict the contours of growth rate for different values of shear parameter (𝜏) and slip parameter (𝛽) as a function of 𝑅𝑒1, respectively. In both the figures, the stokes flow is validated as 𝑅𝑒1 β†’ 0, where the flow is completely stabilized. As the 𝑅𝑒1 increases, the growth rate increases and destabilizes the flow due to the increase in the inertial effect. In the case of Fig. 8(a), the shear imposed Marangoni instability increases as the value of 𝜏 increases. However, the variation of slip parameters (Fig. 8(b)) destabilizes the fluid and it decreases gradually for the larger values of 𝑅𝑒1. In Figs. 9(a), (b) and (c) depict the marginal stability curve for different values of shear parameter, slip parameter and Marangoni number, respectively. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1334 https://internationalpubls.com (a) (b) (c) Figure 9: Neutral stability curve for different values of (a) shear parameter 𝝉 with 𝜷 = 𝟎. πŸ“ and 𝑴𝒂 = 𝟎. 𝟏, (b) slip parameter 𝜷 with 𝝉 = 𝟎. πŸ” and 𝑴𝒂 = 𝟎. 𝟏 and (c) Marangoni number 𝑴𝒂 with 𝜷 = 𝟎. πŸ“ and 𝝉 = 𝟎. πŸ”, and other parameters are π‘ͺ𝒂 = 𝟎. 𝟎𝟐, π’Ž = 𝟎. πŸ–, 𝒅 = 𝟎. πŸ“ and 𝒓 = 𝟎. πŸ”. In all figures, it is noticed that the right hand side of each curve is known as unstable mode bandwidth, while the left hand side is known as stable mode bandwidth. The unstable mode bandwidth known as shear mode, which increases for increase in the values of shear parameters and slip parameters as in Figs. 9(a), (b) and (c), respectively. Further, the critical Reynolds number shift towards left for larger values of 𝜏 and 𝛽. This increases the unstable shear mode, which destabilizes the fluid flow. On the other hand, the unstable mode bandwidth decreases for increase in the values of π‘€π‘Ž. As a consequence, the presence of insoluble surfactant stabilizes the fluid flow with increase in its concentration, which is evident from Fig. 9(c). Hence, the insoluble surfactant present on the heterogeneous fluid, stabilizes the uniform external shear-imposed flow. 4 Conclusion The present study deals with the effect of insoluble surfactant on the shear-imposed fluid flow having the slippery bottom using normal mode analysis. The resultant Orr-Sommerfeld system is derived and solved by using the spectral collocation method, where the domain is discritized by Chebyshev collocation points. The effect of surfactant concentration, shear and slip parameters on the spatio- temporal growth rate are analyzed for smaller Reynolds number. The dispersion curves shows that the Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1335 https://internationalpubls.com instability enhances for the larger values of shear parameters, whereas it decreases for the larger surfactant concentration. Also, the instability increases when the interface is closer to the slippery bottom and the presence of high viscous fluid on the upper layer. Further, the shear-imposed instability decreases as a consequence of insoluble surfactant. It is concluded that the presence of surfactant stabilizes the both the gravity-driven and shear-driven heterogeneous flows. This study can be helpful in understanding the instabilities occurred due to the shear imposed on the respiratory airways. References [1] A. Anjalaiah, R. Usha, S. Millet, Thin film flow down a porous substrate in the presence of an insoluble surfactant: Stability analysis, Phys. Fluids., 25.2(2013). [2] F.A. Bhat, A. Samanta, Linear stability analysis of a surfactant-laden shear-imposed falling film, Phys. Fluids., 31.5(2019). [3] F.A. Bhat, A. Samanta, Linear stability for surfactant-laden two-layer film flows down a slippery inclined plane, Chem. Eng. Sci., 220(2020), 115611. [4] M.G. Blyth, C. Pozrikidis, Effect of surfactant on the stability of film flow down an inclined plane, J. Fluid. Mech., 521(2004), 241-250. [5] K.P. Chen, Wave formation in the gravity driven lower Reynolds number flow of two liquid films down an inclined plane, Phys. Fluids A: Fluid Dyn., 5.12(1993), 3038-3048. [6] P.G. Drazin, Introduction to hydrodynamic stability, Cambridge university press, 32 (2002). [7] F.F. Espinosa, R.D. Kamm, Bolus dispersal through the lungs in surfactant replacement therapy, J. Appl. Phys., 86.1(1999), 391-410. [8] A.E. Hosoi, L. Mahadevan, Peeling, healing, and bursting in a lubricated elastic sheet, Physical review letters, 93.13(2004), 137802. [9] M.M. Hossain, S. Ghosh, H. Behera, Linear instability of a surfactant-laden shear imposed falling film over an inclined porous bed, Physics of Fluids., 34.8 (2022). [10] M. Sani, SA Selvan, S. Ghosh, H. Behera Effect of imposed shear on the dynamics of a contaminated two-layer film flow down a slippery inclined. Physics of Fluids (2020), 32(10). [11] T.W. Kao, Role of the interface in the stability of stratified flow down an inclined plane, Phys. Fluid., 8.12(1965), 2190-2194. [12] T.W. Kao, Role of viscosity stratification in the stability of two-layer flow down an incline, J. Fluid. Mech., 33.3(1968), 561-572. [13] T.W. Kao, Stability of two-layer viscous stratified flow down an inclined plane, Technical report, NASA-CR-57878, (1964). [14] P. Gao, X.Y. Lu, Effect of surfactants on the inertialess instability of a two-layer film flow, J. Fluid, Mech., 591(2007), 495-507. [15] N.A. Malamataris, M. Vlachogiannis, V. Bontozoglou, Solitary waves on inclined films: Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1336 https://internationalpubls.com Flow structure and binary interactions, Phys. Fluids., 14.3(2002), 1082-1094. [16] D.R. Otis, M. Johnson, T.J. Pedley, R.D. Kamm, Role of pulmonary surfactant in airway closure: a computational study, J. Appl. Phys., 75.3(1993), 1323-1333. [17] X. Li, C. Pozrikidis, The effect of surfactants on drop deformation and on the rheology of dilute emulsions in Stokes flow, J. Fluid. Mech., 341(1997), 165-194. [18] I. Sadiq, R. Usha, Thin Newtonian film flow down a porous inclined plane: Stability analysis, Phys. Fluids., 20.2(2008). [19] A. Samanta, Effect of surfactants on the instability of a two-layer film flow down an inclined plane, Phys. Fluids., 26.9 (2014). [20] A. Samanta, Shear-imposed falling film, J. Fluid. Mach., 753(2014), 131-149. [21] A. Samanta, Linear stability of a plane CouetteΓ’β‚¬β€œPoiseuille flow overlying a porous layer, Int. J. Multi. Flow., 123(2020), 103160. [22] A. Samanta, Effect of surfactant on the linear stability of a shear-imposed fluid flowing down a compliant substrate, J. Fluid. Mech., 920(2021), A23. [23] S. Pal, A. Samanta, Linear stability analysis of a shear-imposed oldroyd-b liquid flowing down a slippery inclined plane, Bulletin of the American Phys. Soc., 67(2022). [24] M. Sani, H. Behera, S. Ghosh, Stability analysis of a film flow down an incline in the presence of a floating flexible membrane, Math. Model. Sci. Comput. Appl., 19 (2020), 253-263. [25] M. Sani, S. Ghosh, H. Behera, Effect of a floating elastic membrane for stabilizing the film flow down a porous inclined plane, J. Appl. Math. Mach., 101.2, (2021), e201900246. [26] I.R. Siqueira, R.L. Thompson, M.S. Carvalho, P.R. de Souza Mendes, Slot coating of viscoplastic materials: A computational study of the effects of viscoplasticity on the flow dynamics and low- flow limit, J. Non-Newtonian Fluid Mech., 327 (2024), 105222. [27] M. Sivapuratharasu, S. Hibberd, M.E. Hubbard, H. Power, Inertial effects on thin-film wave structures with imposed surface shear on an inclined plane, Phys. D: Nonlinear Phenom.,325 (2016), 86-97. [28] S.A. Selvan, S. Ghosh, H. Behera, M.H. Meylan, Hydroelastic response of a floating plate on the falling film: A stability analysis, Wave Motion, 104 (2021), 102749. [29] Y.Y. Suy,The stability of superposed polymeric fluids flowing down an inclined plane and its application to coating processes, Polymer-Plastics Tech. Eng., 34.2 (1995), 243-269. [30] E.O. Tuck, Continuous coating with gravity and jet stripping, Phys. Fluid., 26.9 (1983), 2352- 2358. [31] C.Y. Wang,Thin film flowing down a curved surface, J. Appl. Math. Phys., 35 (1984), 532-544. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 1337 https://internationalpubls.com [32] C.S. Yih, Stability of liquid flow down an inclined plane, World Scientific, (1991), 357-370. [33] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T.A. Zang, Spectral methods: evolution to complex geometries and applications to fluid dynamics. Springer Science and Business Media, 2007. [34] R.W. Chin, F.H. Abernath, and J.R. Bertschy, Gravity and shear wave stability of free surface flows. Part 1. Numerical calculations. J. Fluid Mech. 168 (1986), 501–513.