Evolution of a thin film down an incline: A new perspective

Evolution of a thin film down an incline: A new perspective Cite as: Phys. Fluids 32, 013603 (2020); https://doi.org/10.1063/1.5127815 Submitted: 16 September 2019 . Accepted: 26 December 2019 . Published Online: 10 January 2020 Usha Ranganathan, Geetanjali Chattopadhyay , and Naveen Tiwari Phys. Fluids 32, 013603 (2020); https://doi.org/10.1063/1.5127815 © 2020 Author(s). 32, 013603 Physics of Fluids ARTICLE scitation.org/journal/phf Evolution of a thin film down an incline: A new perspective Cite as: Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Submitted: 16 September 2019 • Accepted: 26 December 2019 • Published Online: 10 January 2020 Usha Ranganathan,1,a) Geetanjali Chattopadhyay,2 and Naveen Tiwari2 AFFILIATIONS 1 Department of Mathematics, Indian Institute of Technology Madras, Chennai 600036, Tamilnadu, India 2 Department of Chemical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, Uttar Pradesh, India a) Electronic mail: [email protected] ABSTRACT A new model which accounts for energy balance while describing the evolution of a thin viscous, Newtonian film down an incline at high Reynolds numbers and moderate Weber numbers has been derived. With a goal to improve the predictions by the model in inertia dominated regimes, the study employs the Energy Integral Method with ellipse profile EIM(E) as a weight function and is motivated by the success of EIM in effectively and accurately predicting the squeeze film force in squeeze flow problems and in predicting the inertial effects on the performance of squeeze film dampers [Y. Han and R. J. Rogers, “Squeeze film force modeling for large amplitude motion using an elliptical velocity profile,” J. Tribol. 118(3), 687–697 (1996)]. The focus in the present study is to assess the performance of the model in predicting the instability threshold, the model successfully predicts the linear instability threshold accurately, and it agrees well with the classical result [T. Benjamin, “Wave formation in laminar flow down an inclined plane,” J. Fluid Mech. 2, 554–573 (1957)] and the experiments by Liu et al. [“Measurements of the primary instabilities of film flows,” J. Fluid Mech. 250, 69–101 (1993)]. The choice of the ellipse profile allows us to have a free parameter that is related to the eccentricity of the ellipse, which helps in refining the velocity profile, and the results indicate that as this parameter is increased, there is a significant improvement in the inertia dominated regimes. Furthermore, the full numerical solutions to the coupled nonlinear evolution equations are computed through approximations using the finite element method. Based on a measure {used by Tiwari and Davis [“Nonmodal and nonlinear dynamics of a volatile liquid film flowing over a locally heated surface,” Phys. Fluids 21, 102101 (2009)]} of the temporal growth rate of perturbations, a comparison of the slope of the nonlinear growth rate with the linear growth rate is obtained and the results show an excellent agreement. This confirms that the present model’s performance is as good as the other existing models, weighted residual integral boundary layer (WRIBL) by Ruyer-Quil and Manneville [“Improved modeling of flows down inclined planes,” Eur. Phys. J: B 15, 357–369 (2000)] and energy integral method with parabolic profile [EIM(P)] by Usha and Uma [“Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method,” Phys. Fluids 16, 2679–2696 (2004)]. Furthermore, for any fixed inclination θ of the substrate, 0 < θ < π/2, there is no significant difference between the EIM(E) and EIM(P) results for weaker inertial effects, but when the inertial effects become stronger, the EIM(E) results for energy contribution from inertial terms to the perturbation at any streamwise location is enhanced. More detailed investigation on the model’s performance due to this enhancement in energy contribution and the assessment of the model as compared to the other existing theoretical models, experimental observations, and numerical simulations, in the inertia dominated regimes, will be reported in a future study. Published under license by AIP Publishing. https://doi.org/10.1063/1.5127815., s I. INTRODUCTION A wide range of dynamic phenomena is displayed by flow of liquid films down an inclined plane in the presence of a free surface. Many simplified model equations have been developed starting from the full Navier-Stokes equations by using a long-wave Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing approximation and by assuming various order-of-magnitude of the dimensionless parameters (Reynolds numbers and Weber numbers). These models have been developed with a goal to capture main physical features of the flow. The fascinating dynamics and features of gravity driven thin film flow have been captured and understood using developed models for a wide range of flow systems.1–5 32, 013603-1 Physics of Fluids Following Kapitza and Kapitza,6 who performed experimental as well as theoretical studies, Benjamin7 has obtained the instability threshold using a power series expansion of the Orr-Sommerfeld equations. Benjamin7 and later Yih8 have shown that the Nusselt flat film stationary solution with a parabolic velocity profile is unstable to long waves, if the Reynolds number (Rec ) exceeds a critical value, that is, Rec ≙ 56 cot θ, where θ is the inclination of the inclined substrate with the horizontal. Few numerical solutions of the full Navier-Stokes equations have been considered, and they are restricted to small Reynolds numbers and two-dimensions. Although this is difficult and involves high computational cost due to the presence of a free surface, there are studies9–14 which have performed numerical scrutiny of the governing equations and boundary conditions. Chang et al.15 and Demekhin et al.16 have performed the analytical scrutiny of the reduced Navier-Stokes equations having the same dimensionality obtained by assuming that the streamwise variations are small as compared to crosswise variations (boundary layer equations). This has initiated the development of two-equation models for thickness of the film and flow rate. Shkadov17 was the first to propose such a two equation model starting from the boundary layer equations15 [an Integral Boundary Layer (IBL) model]. He has employed the averaging method of Kárman-Pohlhausen with Nusselt solution parabolic velocity as the downstream velocity and integrated the streamwise momentum equation, and this has yielded a two-equation model for film thickness and local flow rate. Although this model has reproduced solitary waves observed in experiments, it has failed to predict the correct instability threshold; the model equation was not consistent with the Benney equation, and the parabolic velocity profile did not satisfy the dynamic boundary condition at the second order. However, this model is capable of describing the dynamics of waves in thin film flows but is not consistent with the long-wave analysis of the Orr-Sommerfeld equations. Note that the correct critical Reynolds number is recovered by this model for films down vertical substrates, in the long-wave limit, but for moderate inclinations of the inclined plane, the predicted instability threshold is not correct. Even if the Nusselt flow is unstable, the self-similar parabolic profile of the Nusselt flow is employed to derive the depth averaged Navier-Stokes equations. After Benney,18 many consistent models in the limit of longwave asymptotic approximation have been derived for the film thickness evolution. Such a first order scalar hyperbolic conservation equation displays finite time discontinuity near the critical Reynolds number. Also, there is a blow up of solutions of the Benney type evolution equation, when linear stability conditions are violated.19–21 However, the Benney equation resembles the Kuramoto-Sivashinsky equation22,23 in the limit of small-amplitude modulations. Takeshi24 resolved this deficiency by employing a Padé-like regularization technique; but the model failed to yield accurate wave solutions away from the critical Reynolds number. The above investigations with one equation model for film thickness clearly indicate that such a model is inadequate to provide accurate and complete details of wave development in a thin film down an incline. In the subsequent investigations, the focus has been on deriving (i) models of low dimensionality so as to get a reasonable computational cost and (ii) models with a desired mathematical Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing ARTICLE scitation.org/journal/phf structure so as to employ robust numerical schemes which will yield reliable results. The models are based on long-wave asymptotic expansion as the flow remains close to Nusselt flow in most cases. They are of lower dimensions since the depth average of the Navier-Stokes equations is performed. The computations become faster and easier due to inclusion of the boundary conditions in the model. Motivated by the experimental observations that solitary waves in thin film flows display a deviation from the parabolic profile,25,26 Nguyen and Balakotaiah27 have employed a general parabolic profile and have obtained a three-equation model for film thickness, flow rate, and wall shear stress, which has predicted physically meaningful results and quantitative predictions of large amplitude waves. Ruyer-Quil and Manneville28 have attributed the failure of Shkadov’s model to the lack of flexibility of the assumed velocity profile. Their model, derived using a Galerkin method, is based on a long-wave expansion of the Navier-Stokes equations and consists of four evolution equations for the local flow depth, local flow rate, and two other quantities measuring the departure from the parabolic velocity profile. The Nusselt solution and three more polynomials appearing in the Benney equation derivation have been used. A simplified two-equation model [Weighted Residual Integral Boundary Layer (WRIBL) model] for local flow depth and local flow rate has been obtained from the above four-equation model, and this is consistent with the Benney equation of the second order and predicts the instability threshold accurately. The deficiency in the Benney equation which displayed unphysical solutions that grow indefinitely in time19,29–32 and the failure of Shkadov’s17 equation to predict the critical condition for instability correctly could be overcome by the Weighted Residual Integral Boundary Layer (WRIBL) model.28,33 Their first and second order models, for large Kapitza numbers, have predicted the linear instability threshold accurately. In a wide window of Reynolds numbers, their solutions remain bounded. When the Reynolds number slightly exceeds the critical Reynolds number, then the above model fails to capture the solitary wave solutions accurately. Following the success of this model, several investigations involving both confined geometries and free surface flows34,35 have employed WRIBL and have analyzed and predicted reliable and accurate stability features. It is important to note that the results obtained from the Weighted Residual Integral Boundary Layer (WRIBL) method extended to three-dimensional flows show a very good agreement with the experimental results.36 Mudunuri and Balakotaiah37 have assumed a self-similar parabolic velocity profile and have obtained a two-mode model similar to the two-equation model by Ruyer-Quil and Manneville,28 using a Galerkin projection, and this describes the dynamics of film down an incline/vertical substrate for Weber numbers >1. The above studies on lowdimensional models for thin films are based on a closure assumption for the self-similar parabolic velocity profile in the streamwise direction.38 Teshukov39 has proposed a method not involving the notion of a profile for turbulent flows. Motivated by its success in describing turbulent roll waves and in the study of hydraulic jump,40,41 Richard et al.42 have generalized it to laminar viscous film flows down an incline with a long-wavelength asymptotic expansion. 32, 013603-2 Physics of Fluids They have derived a three equation model for the film thickness, average velocity, and entropy related to the variance of the velocity, which describes the flow of a thin viscous film down an incline. The structure of the resulting equation is the same as that of Euler equations governing compressible fluids. The equation includes the effects of diffusivity and capillarity in addition to relaxation of source terms. Classical, reliable, and robust numerical schemes have been employed to solve these equations. The model results and the experimental observations43,44 agree well. Amaouche et al.45 have included third order terms in order to capture the effects of the small Weber number and high Reynolds numbers and have developed a model for describing the evolution of waves on a film flow down an incline using WRIBL. They have attempted to modify or refine the velocity profile through a free parameter and have performed a linear stability analysis, and the results are shown to be in good agreement with OrrSommerfeld equations for all Weber and Kapitza numbers. A substantial improvement in the predictions is observed in the inertia dominated regimes. Lee and Mei46 have developed a two-equation Momentum Integral Model (MIM) for local flow depth and local flow rate, valid for moderate and high Reynolds numbers, and is based on a parabolic velocity profile. The model equations captured the Hopf bifurcation features but failed to predict the correct critical Reynolds number. The dynamics and stability of a thin film flowing down an inclined substrate has been examined by Usha and Uma47 using the Energy Integral Method (EIM) with parabolic profile as the weight function. The two equation model represents a second order approximation to the Navier-Stokes equations and the boundary conditions and captures the effects for the large Reynolds number and moderate surface tension. Their study has been motivated by evolution equations based on EIM, employed by Elkouh48 and Crandall and Ei-Shafei49 for squeeze film flows and the investigations on the computation of normal force on a squeeze film for large inertial effects by Han and Rogers50 using EIM. The above studies on squeeze film flow have revealed that the inclusion of energy approximations has shown to be more effective in understanding the effects of inertia on the squeeze-film damper performance. The exact solutions of squeeze-film damper flows49 and the predictions by EIM agree well for Reynolds numbers of order one. At the low Reynolds number, there is a 20% error at the first order approximation, thereby indicating the capability of EIM as compared to the Momentum Integral Method (MIM) in the prediction of inertial effects on the squeezefilm damper performance. Elkouh’s48 accurate prediction on pressure distribution by EIM and the available experimental observations51 agree well. There are also other evidences which suggest the use of EIM as the predicted results are close to either experimental observations or numerical simulations.50,52–54 It is important to observe that the squeeze film force has been accurately predicted by the EIM for small to moderate inertial effects and for large amplitude motion. Note that the simplified two equation model, deduced from a four equation model,28,33 is identical with the two equation model for local film depth and flow rate developed by EIM with parabolic profile by Usha and Uma.47 The EIM model with parabolic velocity profile by Usha and Uma47 has predicted the instability threshold correctly in the linear regime (as given in the work of Benjamin7 and Yih8 ). They have also presented qualitative and quantitative Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing ARTICLE scitation.org/journal/phf features of traveling wave solutions by deriving a dynamical system by considering the model equations in a stationary frame of reference. Abderrahmane and Vatistas55 have pointed out that the second order models by Lee and Mei46 and Usha and Uma47 have not accounted for the tangential component of the dynamic interface condition up to the second order as their choice of the parabolic velocity profile cancels the shear stress at the free surface. This drawback has been removed in the model developed by them by incorporating the shear stress at the interface as variable in the velocity profile. Their new model thus is a refinement of Usha and Uma’s model and is consistent to the third order of accuracy. This approach has been effectively employed in subsequent investigations by Ruyer-Quil et al.56 in their two-equation model that describes film flow on a vertical cylinder to the first and second orders of approximation to the governing Navier-Stokes equations. They have considered a linear stability analysis of the model equations and have investigated the traveling wave solutions. Novbari and Oron,57,58 in their study on nonlinear dynamics of an axisymmetric thin film flowing down a vertical cylinder, have derived their model equation based on EIM proposed by Usha and Uma,47 approximating the Navier-Stokes equations to the first order of a small aspect ratio parameter. Their study has demonstrated the success of EIM in predicting results which are “sufficiently good to fit several experimental observations better than other theories of higher order of approximation existing in the literature in terms of wave shapes and amplitudes.” Furthermore, they have shown that EIM admits both traveling wave and nonstationary wave flows. Luchini and Charru59 have noted that the kinetic-energy equation derived by Usha and Uma47 gives the right stability threshold of long waves in a falling film. Furthermore, they have pointed out its generality as a consistent nonlinear formulation and have presented the method of consistent section-averaged equations using the kinetic-energy balance appropriate to the boundary layer approximation in their investigation on shallow-water flow over a perturbed bottom. Using EIM, Sadiq60 has investigated thin film flow over a topography having sinusoidally varying longitudinal furrows. Sadiq60 has examined the model equations numerically by posing it as an Initial Value Problem (IVP) on a periodic domain for a configuration resembling a vertically shifted topography with a shift value equal to the amplitude value. The film flow is stabilized as compared to a topography which does not undergo any vertical shift. Han and Rogers50 have demonstrated the advantage of employing EIM with ellipse velocity profile over EIM either with parabolic profile/MIM or with parabolic/ellipse profile, while exploring squeeze film force modeling for large amplitude motion. This study clearly shows the efficiency of the EIM with ellipse profile to predict the effect of inertia on the squeeze film damper’s performance. Motivated by the above study and by the curiosity to develop more and more accurate models capable of providing quantitative information about the fully developed regime far away from the inlet where nonlinearities are expected to dominate the wave profiles, in this article, a new two-equation model for evolutions of 32, 013603-3 Physics of Fluids ARTICLE local flow depth and local flow rate is developed. The focus in this study is only on the development of a new model and examining how correctly it predicts the instability threshold. Furthermore, in order to highlight the model’s capability to capture the features in the nonlinear regime, the full numerical solutions to the coupled nonlinear evolution equations are obtained through computations based on the finite element method (FEM). A comparison of the slope of the nonlinear growth rate with the linear growth rate based on a measure61 of temporal growth rate of perturbations reveals an excellent agreement. This confirms how accurately the present EIM(E) model involving a free parameter, related to the eccentricity of the ellipse profile chosen as the weight function, also predicts the solutions as good as the other existing models such as those by Ruyer-Quil and Manneville,33 Amaouche et al.,45 and Usha and Uma.47 The present study presents a perspective of the developed model and a detailed rigorous study on the model’s capability to capture the other features such as traveling-wave and solitary-wave solutions in the strongly nonlinear regime, and the assessment of the model’s performance as compared to the other theoretical models through experimental observations and numerical simulations will be considered in a future study. density, p is the pressure, and μ is the dynamic viscosity of the fluid. The conditions at the boundaries are ⃗ ≙⃗ v 0 at ht + uhx ≙ v ⃗ ⋅ n̂ ≙ 0 t̂ ⋅ T ⃗ ≙ 0, ∇⋅v (1) y ≙ η(x, t), (4) y ≙ η(x, t), (5) at y ≙ η(x, t), (6) g sin θh20 2y y2 [ − 2 ], v ≙ 0, 2ν h0 h0 where h0 is the depth of the primary flow. The corresponding depth averaged velocity is ū0 ≙ h0 g sin θh20 1 u0 dy ≙ , ∫ h0 0 3ν (7) and the flow rate is Q̄0 ≙ u¯0 h0 . The dimensionless variables are defined by x∗ ≙ (2) in a domain Ω(t) ≙ {(x, y)∣x ∈ R, 0 < y < η(x, t)} occupied at ⃗ ≙ (u, v) is the velocity vector, ρ denotes time t by the fluid, where v (3) which correspond to the no-slip condition at the bottom, the kinematic free surface condition, and the balance of normal and ⃗ is the stress tensor and shear stresses at the free surface. Here, T (−ηx , 1) n̂ ≙ √ is the unit outward normal vector to the free bound1 + η2x ary. The dimensionless form of the above equations and the boundary conditions (1)–(6) are obtained by choosing scales based on the time dependent Nusselt solution for the system (1)–(6) which corresponds to solution of the system (1)–(6) for uniform flow given by II. MATHEMATICAL FORMULATION ⃗ + g, ρ∥⃗ v t + (⃗ v ⋅ ∇)⃗ v ∥ ≙ −∇p + μ∇2 v at y ≙ 0, ⃗ ⋅ n̂ ≙ −σ0 (∇ ⃗ ⋅ n̂) at n̂ ⋅ T u0 ≙ A viscous, incompressible, Newtonian thin film fluid flowing down an incline making an angle θ (0 < θ < π2 ) with the horizontal under the action of gravity and surface tension σ 0 is considered. A Cartesian coordinate system with the x-axis along the inclined plane and the y-axis normal to it is the reference frame for the above two-dimensional motion [Fig. 1(a)]. The free-surface is located at y = η(x, t), and the governing equations are the continuity and Navier-Stokes equation, given by scitation.org/journal/phf η Lv x ∗ y u , , y ≙ , H ≙ , u∗ ≙ , v ∗ ≙ ū0 L h0 h0 h0 ū0 t∗ ≙ h0 tū0 ∗ p − pa ,p ≙ ,ϵ ≙ , L L ρū20 (8) where L is the length scale in the direction of flow and it is of the order of the wave length. This study is pursued under the assumption that the ratio of layer depth to the wave length is very small. The corresponding dimensionless equations (omitting ∗) are ux + vy ≙ 0, ut + uux + vuy ≙ −px + ϵ2 ∥vt + uvx + vvy ∥ ≙ −py − ϵ 1 3 + uxx + uyy , ϵRe Re ϵRe (10) 3 cot θ ϵ3 ϵ + vxx + vyy , Re Re Re (11) y ≙ 0, (12) u ≙ 0, v ≙ 0 FIG. 1. (a) Schematic diagram of a thin viscous film flow down an incline and (b) the ellipse profiles [Eq. (19), 0 < y < 1] for different values of the free parameter A with corresponding values of eccentricity (e listed in Table I); the insets show the zoom-in view of the velocity profiles at different y locations. Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing p+ (9) at −1 2ϵ [uy Hx + ϵ2 vx Hx − ϵ2 ux Hx2 − vy ][1 + ϵ2 Hx2 ] Re + ϵ2 We Hxx [1 + ϵ2 Hx2 ] − 23 ≙0 at y ≙ H, (13) 32, 013603-4 Physics of Fluids ARTICLE (uy + ϵ2 vx )[1 − ϵ2 Hx2 ] + 2ϵ2 ∥vy − ux ∥Hx ≙ 0 v ≙ Ht + uHx at y ≙ H, y ≙ H, at This in turn gives u as (14) (15) g sin θh30 ū0 h0 ≙ is the Reynolds number and where Re ≙ ν 3ν2 σ0 We ≙ is the Weber number. ρū20 h0 Our goal is to derive a model that describes the temporal and spatial evolution of the local flow rate Q and local flow depth H(x, t) from equations (9)–(15). This is achieved by using the energy integral method47 with an appropriate ellipse type profile which is derived below. u (x, y, t) + B1 (x, t)y + C1 (x, t)u + D1 (x, t)y + E1 (x, t) ≙ 0 u ≙ 0 at ∂u ≙ 0 ∂y at y ≙ 0, (17) y ≙ H. (18) Equation (17) implies that E1 (x, t) ≡ 0 and (18) gives D1 (x, t) = −2HB1 (t). This reduces (16) to (u + C21 )2 (y − H)2 + ≙ 1, a21 b21 where C1 ≙ (19) AH 2a1 √ 2 , b1 − H 2 , b1 ≙ b1 2 where a1 and b1 are semimajor and semiminor axes of the ellipse profile and are related to B1 and C1 by a21 ≙ C12 a2 + B1 H 2 , B1 ≙ 21 , 4 b1 and A is the free parameter that allows us to refine the ellipse profile. Simplification of (19) yields the ellipse type velocity profile as √ a1 √ 2 2 (20) [ A H − 4(y − H)2 − A2 H 2 − 4H 2 ]. u(x, y, t) ≙ AH The associated flow rate Q(x, t) ≙ ∫H 0 u(x, y, t) dy is obtained as Q(x, t) ≙ a1 [ 2 H√ 2 AH sin−1 ( ) − A − 4]. 4 A 2A Now, a1 is expressed in terms of Q, H, and A as a1 ≙ Q AH 4 −1 sin H ( A2 ) − 2A . √ A2 − 4 Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing (22) Kp ≙ −1 √ A2 − 4. (24) 1 A2 4 sin i.e ( A2 ) − 21 ,B ≙ √ A2 − 4 (25) (u + Kp B)2 (y − 1)2 + ≙ 1, A2 /4 Kp2 A2 (26) A 2 ( ) ≙ Kp2 A2 (1 − e2 ), 2 (27) where the major axis is K p A and the minor axis is A/2. Therefore, and eccentricity (e) of the ellipse is e2 ≙ 1 − √ 2 [A2 sin−1 ( A2 ) − 2 A2 − 4] 64 . (28) As A increases, eccentricity (e) increases and approaches the value 1 (Table I). Thus, the parameter A modifies the ellipse profile [Fig. 1(b)]. Figure 1(b) shows the ellipse profiles [used as weight functions; Eq. (19), 0 < y < 1] for different values of the free parameter A, which is related to the eccentricity of the ellipse profile. The velocity exhibits a nonmonotonic trend as the parameter A increases. The velocity increases with an increase in A up to A = 40, but beyond that there is a decrease in the velocity. The velocity for A = 40 is larger than that for A = 20 or A = 60. It is for this reason that in the computations that follow, the A value has been fixed as 40. A more detailed analysis in future will focus on what is the optimal value of this parameter A as the surface tension effects become larger. TABLE I. The eccentricity (e) of the ellipse as parameter A increases. A (21) (23) √ u ≙ Kp [ A2 − 4(y − 1)2 − Kp B], (16) and satisfying the no-slip condition Kp Q √ 2 2 [ A H − 4(y − H)2 − BH], H2 The free parameter A serves to refine the velocity profile and thus acts on the flexibility of the velocity profile which is derived by satisfying the boundary conditions for flow of a flat smooth uniform Nusselt film down an inclined plane. The dimensionless velocity profile is An ellipse type velocity profile (as the weight function) is obtained by taking its equation as 2 u(x, y, t) ≙ where A. Ellipse profile 2 scitation.org/journal/phf 20 30 40 60 80 100 Kp e 14.9549 22.4699 29.9775 44.9850 59.9887 74.9910 0.999 441 0.999 752 0.999 861 0.999 938 0.999 965 0.999 978 32, 013603-5 Physics of Fluids ARTICLE The continuity Eq. (9) then gives the transverse velocity v(x, y, t) as ⎡ ⎫ √ ⎪ ⎪⎤ ⎢ QHx ⎧ ⎥ ⎪ 2 H 2 − 4(y − H)2 ⎪ A v ≙ −Kp ⎢ ⎨ BHy − y ⎬⎥ ⎢ H3 ⎪ ⎥ ⎪ ⎢ ⎪ ⎪ ⎩ ⎭⎥ ⎣ ⎦ − Kp [ T4 ≙ − ∫ (29) (30) 3 cot θ 2ϵ Q Hx Qx (y − H) + (A − B)Kp [ 2 − ] − ϵ2 We Hxx , Re Re H H (31) from which px is computed and substituted in (32). The momentum Eq. (10) along the x-direction is multiplied by the weight function u and is integrated across the film thickness. This gives [u(ut + uux + vuy )dy + upx − 3u ϵ 1 − uuxx − uuyy] dy ≙ 0. ϵRe Re ϵRe (32) Using the ellipse velocity profile in (32) and performing integration of each term in (32) gives T 1 + T 2 + T 3 + T 4 = 0, where 0 T1 ≙ ∫ H 0 ≙ u( ∂u ∂u ∂u +u + v )dy ∂t ∂x ∂y Q2 Ht Kp2 QQt Kp2 2 4 A2 B −1 2 [A − − sin ]+ H 3 2 A H2 + Published under license by AIP Publishing C1E (35) 1 ∂2u 2 2 Q2 A−2 u 2 dy ≙ − Kp [A ln∣ ∣ + 2]. ϵRe ∂y ϵRe H 3 B (36) Qt QHt QQx Q2 Hx 3Hx cot θ + C2E 2 + C3E 2 + C4E + H H H H3 Re − ϵ2 We Hxxx − C1E ≙ Kp2 [A2 − C2E ≙ Kp2 [− 3 C5E Q + ≙ 0, ϵRe ϵRe H 3 (33) (37) 2 4 1 2 − A B sin−1 ( )], 3 2 A 1 2 3 A2 8 + + AB + A2 B sin−1 ( )], 2 3 4 A 1 37 A4 A3 B 2 C3E ≙ Kp3 [( − 8 A2 + ) sin−1 ( ) 2 16 2 A + (B − 37 A2 B − A3 + 4 A)], 8 C4E ≙ Kp3 [(3 A2 − C5E ≙ −2Kp2 [A ln∣ 2 15 A2 B 15 A4 ) sin−1 ( ) + − B], 16 A 8 (A − 2) ∣ + 2]. B The first four terms in Eq. (37) contribute to the inertial effects in the developed model; the next two terms describe the pressure and surface tension effects, and the last but one term and the last term arise due to gravitational acceleration and shear stress at the wall. Equations (30) and (37) give the new two-equation model for the spatial and temporal evolutions of Q(x, t) and H(x, t), and it is based on energy balance with the ellipse type profile. The energy integral method (EIM) with parabolic velocity profile (as weight function) formulated by Usha and Uma47 [EIM(P)] reads as follows, if one takes terms accurate up to the first order, + 37 4 AB 2 37 A − 8A2 + } sin−1 + B − A2 B − A3 + 4A] 16 2 A 8 Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 3 3 udy ≙ − Q, ϵRe ϵRe H (p) 3 Q3 Hx Kp3 2 15 15 [{− A4 + 3A2 } sin−1 + A2 B − B], H3 16 A 8 (34) (p) 3Ht Q(p) 81Q(p) Qx 6 Qt − + 5 H 5H 2 35H 2 Q2 Qx Kp3 3 8 A2 B −1 2 × [− A2 + + AB + sin ]+ 2 3 4 A 2H 2 × [{ ∂p 3 cot θ dy ≙ Hx Q − ϵ2 We QHxxx , ∂x Re The evolution model for Q and H [EIM(E)] is given by p≙− H H 0 The first-order accurate y-momentum Eq. (11) yields ∫ u 0 We now employ a Galerkin based formulation with the weight function as the ellipse type velocity profile u(x, y, t) satisfying Nusselt film boundary conditions and derive the two-equation model for local flow rate Q(x, t) ≙ ∫H 0 u(x, y, t) dy and local flow depth H(x, t) by requiring that u, v, p also satisfy the weighted depth averaged system (10) and (11) along with (9) and (12)–(15). As mentioned in the Introduction, a first-order accurate EIM is developed based on the confidence of the results obtained using EIM in other configurations,47,55,57,58 when both ϵRe and We are of order one. Now, the depth average of the continuity Eq. (9) along with the use of the kinematic free surface condition (15) gives Ht + Qx ≙ 0. H 0 T3 ≙ − ∫ 2(y − H) (y − H) Qx A2 H 2 { sin−1 + H2 4 AH 2 √ A2 H 2 2 sin−1 × A2 H 2 − 4(y − H)2 + 4 A √ H2 √ 2 A − 4 − A2 − 4 Hy}]. + 2 T2 ≙ ∫ scitation.org/journal/phf 2 54Hx Q(p) 3 − − 35H 3 ϵRe 1 3Q(p) 3 cot θ + Hx − ϵ2 We Hxxx ≙ 0, ϵRe H 3 Re (38) (p) where Q(p) ≙ ∫H dy, with 0U U (p) ≙ 3Q(p) 2y y2 ( − ). 2H H H 2 (39) 32, 013603-6 Physics of Fluids ARTICLE We observe that the new model (37) contains exactly the same terms but with different coefficients depending on the parameter A. Note that the constants K p and B, in Eq. (24), also depend on A. Following the above computations, the following two equation model also has been developed for comparison purposes. In the present study. This is based on the momentum integral method and employs an ellipse profile (23) as a weight function. The model equations [MIM(E)] valid up to the fist order are Ht + Qx ≙ 0, − C5M Q 3 + ≙ 0, ϵRe ϵRe H 3 C1M C3M ≙ C5M ≙ P1E + iQ1E ≙ 4 A2 B −1 2 + sin ( )], 3 2 A iC5E 1 [C2E − C3E + ], C1E ϵRe 3iC5E 1 3 cot θ − C4E − ϵ2 We − ]. [− C1E Re ϵ Re −Q1E ± N1 −P1E ± M1 and ci ≙ , 2 2 √ 2 − Q2 − 4P ) + i(2P Q − 4Q ). Writwhere M1 + iN1 ≙ (P1E 2E 1E 1E 2E 1E ing P1E , Q1E , P2E , Q2E in terms of C1E , C2E , C3E , C4E , and C5E , we have cr ≙ Also, Lee and Mei46 have used a parabolic profile (39) as a weight function and have developed a model based on the Momentum Integral Method [MIM(P)], and it reads as ≙ 0, P1E ≙ (42) 2 (p) Qt 3 Q(p) Ht 9 Q(p) Qx 6 Q(p) Hx 3Hx cot θ − + − + H 2 H2 10 H2 5 H3 Re P2E ≙ C5E C2E − C3E , , Q1E ≙ C1E C1E ϵ Re 1 3 cot θ C5E 1 [− − C4E − ϵ2 We], Q2E ≙ −3 . C1E Re C1E ϵ Re 3 Q(p) 3 + ≙ 0. (43) ϵRe ϵRe H 3 In what follows, we perform a linear stability analysis of the base Nusselt flow and obtain the critical condition for instability threshold predicted by the new model (37) and compare the results with the other models [Eqs. (38), (41), and (43)]. For ci = 0, the flow is neutrally stable and this yields III. LINEAR STABILITY OF THE BASE FLOW This can be rewritten in terms of the free parameter A as − ϵ2 We Hxxx − A. Instability threshold predicted by EIM(E) Corresponding to infinitesimal two-dimensional disturbances (η̂) of the basic uniform flow for the free surface deflection and local flow rate given by H ≙ 1 + η̂ and Q ≙ 1 + Q̂, the evolution equations for Q and H [(30) and (37)] are linearized; Q̂ is then eliminated. This gives C1E η̂tt + (C3E − C2E )η̂xt − C4E η̂xx − 3 cot θ η̂xx Re C5E (η̂t + 3η̂x ) ≙ 0, ϵRe with an accuracy up to the first order. + ϵ2 We η̂xxxx + Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing (45) The sign of ci determines whether the primary Nusselt film flow is stable or unstable. For ci > 0, it is unstable, and for ci < 0, it is stable. The real and imaginary parts of two distinct solutions of Eq. (45) are given by 4Kp . B (p) c2 + c(P1E + iQ1E ) + (P2E + iQ2E ) ≙ 0, P2E + iQ2E ≙ 2 8 A − − A2 B sin−1 ( )], 3 A (p) Ht + Q x where c = cr + ici is the complex phase velocity, cr is the phase speed, and ci is the linear growth rate, and substituting in (44), a quadratic equation in c is obtained and is given by (41) 2 C4M ≙ Kp2 [−A2 + η̂ ≙ exp{i(x − ct)}, where Kp 2 −1 2 ≙ [A sin ( ) − 2B], 4 A Kp2 [2 Taking a wavelike disturbance for η̂ as (40) QQx Q2 Hx 3Hx cot θ Qt + C3M 2 + C4M − ϵ2 We Hxxx + C1M H H H3 Re scitation.org/journal/phf (44) Q22E − P1E Q1E Q2E + P2E Q1E 2 ≙ 0, which simplifies to cot θ C4E ϵ2 We ≙− + 3C1E + (C2E − C3E ) − . Re 3 3 ⎤ ⎡ ⎫ ⎪ ⎪ ⎥ ϵ2 ⎢ −1 2 ⎧ cot θ ⎪ 5 ⎪ 4 3 sin ( )⎨ − A2 B⎬ − + A2 + AB⎥ ≙ Kp2 ⎢ ⎥ − 3 We ⎢ ⎪ Re A ⎪ 4 3 2 ⎥ ⎢ ⎪ ⎪ ⎩ ⎭ ⎦ ⎣ ⎧ ⎫ 3⎡ 3 ⎪ ⎪ Kp ⎢ −1 2 ⎪ ⎢sin ( )⎨ 81 A4 − 18A2 + 3A B ⎪ − ⎬ ⎪ 6 ⎢ A ⎪ 16 2 ⎢ ⎪ ⎪ ⎩ ⎭ ⎣ ⎤ ⎧ ⎫⎥ ⎪ ⎪ ⎪ 81 ⎪ + ⎨ − A2 B + B − 3A3 + 12A⎬⎥ ⎥, ⎪ ⎪ 8 ⎪ ⎪ ⎩ ⎭⎥ ⎦ which gives the instability threshold predicted by EIM with ellipse type profile as the weight function. 32, 013603-7 Physics of Fluids ARTICLE scitation.org/journal/phf D. Instability threshold based on MIM(P) cot θ ≙ 1. Re (48) IV. RESULTS AND DISCUSSIONS cot θ as a funcRe cot θ approaches the tion of the parameter A. We observe that Re value 1.2 predicted by the other models [EIM(P),43 WRIBL33 ], thereby demonstrating its capability of predicting the correct threshold for instability for a thin film flow down an incline. The inertial coefficients Q/H 3 , Q2 H x /H 3 , QQx /H 2 in the models MIM(E) [Eq. (41), present study], MIM(P) [Eq. (43)],46 EIM(P) [Eq. (38)]47 are compared with the coefficients of the same terms in EIM(E) [Eq. (37), present study] and are presented in Table II for A = 20, 30, and 40. We observe that as A increases, the improvement is significant in the inertia dominated regime as is evident by the result in Table II. In addition, the corresponding coefficients in the second order model developed by Ruyer-Quil and Manneville33,67 and Scheid et al.68,69 are computed, and it is observed that they match well with the corresponding coefficients in EIM(P) [Eq. (38)].47 It is also observed that all the inertial coefficients given by EIM(E) lie between those predicted by either MIM(E) or MIM(P) and EIM(P) or WRIBL. This suggests that the present model’s [Eq. (37) EIM(E)] predictions in the nonlinear regime (that is, inertia dominated regions) would be as close to reality as the other three models. In order to facilitate the comparison of our results with the existing models, computations are performed with a wavelike disturbance for η̂ as η̂ ≙ exp∥i(kx − ct)∥, and the phase speed, linear growth rate, and the critical condition on threshold for instability are obtained as functions of wave number k. Figure 3 displays neutral stability curves of the EIM(E) model for different values of the parameter A and EIM(P),47 in the wave numberWeber number regime, for a vertically falling film. The neutral curve for A = 40 is bounded between that for A = 30 and the EIM(P) model. Figure 2 shows the instability threshold cot θ predicted by the two-equation model based Re on EIM(E) [Eq. (37)] and MIM(E) [Eq. (41)] as a function of the free parameter A. FIG. 2. The instability threshold B. Instability threshold predicted by MIM(E) Proceeding as above, the instability threshold predicted by the model based on MIM(E) with ellipse profile [MIM(E)] as the weight function [Eq. (41)] is obtained and is given by cot(θ) −5 A2 20 5 2 2 ≙ Kp2 [ + + A B sin−1 ( )] Re 3 9 6 A 3 2 3 ϵ2 We + Kp [ A2 sin−1 ( ) − B] − . 4 A 2 3 (46) In Subsections III C and III D, we present the instability threshold based on EIM(P) [Eq. (38)]47 and MIM(P) [Eq. (43)]46 for the sake of completeness, for the corresponding first order model. C. Instability threshold based on EIM(P) cot θ 6 ≙ . Re 5 (47) TABLE II. Comparison of coefficients in the evolution equations [Eqs. (37), (38), (41), and (43)] from different models. Note that the superscript (p) is suppressed in the coefficients appearing in MIM(P) [Eq. (43)] and EIM(P) [Eq. (38)]. The coefficients of the terms in the table are evaluated at A = 20, 30, and 40 for MIM(E) [Eq. (41)] and EIM(E) [Eq. (37)]. Coefficient of terms in the evolution equation MIM(E) A = 20 Eq. (41) MIM(E) A = 30 Eq. (41) MIM(E) A = 40 Eq. (41) MIM(P) Eq. (43) EIM(E) A = 20 Eq. (37) EIM(E) A = 30 Eq. (37) EIM(E) A = 40 Eq. (37) EIM(P) or WRIBL Eq. (38) Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing QQx H2 Q2 Hx H3 Q H3 2.399 311 2.399 695 2.399 828 2.4 1.188 865 1.189 099 2.428 388 2.428 571 −1.199 655 −1.199 847 1.199 914 1.2 −1.541 824 −1.542 399 1.285 592 1.285 714 3.006 047 3.002 676 3.001 503 3 3.000 005 3.000 001 2.5 2.5 32, 013603-8 Physics of Fluids FIG. 3. Neutral stability curves for a vertically falling film for different values of the parameter A in the EIM(E) model and the EIM(P) model. The inset shows the zoom-in view of results near 1/(ϵ2 We) = 9. Comparison between wave speeds from the experimental work of Jones and Whitaker,62 Stainthrop and Allen,63 Strobel and Whitaker64 and from EIM(P) and present EIM(E) models is shown in Fig. 4(a) (A = 40) and Fig. 4(b) (A = 20, 30, 40) for vertically falling FIG. 4. The nondimensional phase speed of the perturbation down a vertically falling water film, as a function of the Reynolds number (Re). (a) and (b) comparison of the present results [EIM(E); solid lines] with the experimental results by Jones and Whitaker62 (◽), Stainthrop and Allen63 (◇), Strobel and Whitaker64 (△), and EIM(P) (dashed curve), the zoom-in view in (b) shows that there is no significant difference between the phase speed values for A = 40 and A = 60. (c) and (d) comparison of the model EIM(E) with A = 40 and EIM(P) results with approximate solutions of the Orr-Sommerfeld equations of Krantz and Goren65 and direct numerical solutions of Orr-Sommerfeld equations by Pierson and Whitaker.66 Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing ARTICLE scitation.org/journal/phf FIG. 5. Dimensionless spatial growth rate of linear waves as a function of the wavenumber, for glycerine-water films with θ = 4.6○ , Re = 15.33, and We = 12.13. The solid line corresponds to the present model [EIM(E)] predictions for A = 40, and the dashed and dashed-dotted lines correspond to the EIM(P) model (Usha and Uma47 ) and the WRIBL model (Ruyer-Quil and Manneville67 ). The ◽ symbol represents experimental observations by Liu et al.70 . water films for Re ≤ 100. The results reveal that both EIM(P) and the present EIM(E) models follow the trend predicted by experimental work62–64 and that the values for EIM(E) are closer to the experimental predictions for moderate Reynolds numbers. The above conclusion is also confirmed by the approximate solutions from the momentum integral form of the Orr-Sommerfeld equation by Krantz and Goren65 [Fig. 4(c)] and the direct numerical solutions of the Orr-Sommerfeld equation by Pierson and Whitaker66 [Fig. 4(d)]. The close agreement of the EIM(E) model results with experimental work clearly supports the validity of the linearized theory in the parameter regime Re ≈ O(1/ϵ). The dimensionless growth rate predicted by the EIM(E) model along with those of EIM(P) (Usha and Uma47 ) and the WRIBL model (Ruyer-Quil and Manneville67 ) is presented in Fig. 5 when θ = 4.6○ , Re = 15.33, and We = 12.13. Note that growth rate values are rescaled to account for the scales employed in the work of Liu et al.70 The present EIM(E) model with A = 40 performs as good as the other models such as EIM(P) and WRIBL and follows the trend observed in the experimental study by Liu et al.70 for flow of glycerine-water film down an inclined plane for small wave numbers. For each wave number k in a window of wavenumbers considered, EIM(E) predictions are bounded between EIM(P) (Usha and Uma47 or Sadiq60 for the planar surface) and the WRIBL model (Ruyer-Quil and Manneville67 ). The comparison demonstrates that the model equation [EIM(E)] is valid for very long wavelengths. Furthermore, the results confirm the validity of the approximations incorporated in the linear stability theory. V. NONLINEAR EVOLUTION OF THE FREE SURFACE From the results in Sec. IV, we observed the present EIM(E) model with ellipse velocity profile as the weight function gives satisfactory results for the threshold for instability and compares well with the characteristics of the waves such as the growth rate 32, 013603-9 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 6. (a) A measure of the temporal growth rate of perturbations is plotted as time (t) varies using ξ = ln[∥H(x, t) − 1∥/ξ 0 ], where ∥⋅∥ is the L2 -norm defined over the computational domain and ξ 0 = ∥H(x, t = 0) − 1∥ is the magnitude of the perturbation given to the primary flow. (b) Comparison of the slope of the nonlinear growth rate with the linear growth (linear fit; black squares) for different values of cot(θ) with A = 40. The linear stability (Table III) seems to Re predict the growth accurately at the initial time. TABLE III. The comparison of the growth rate from the two models using ellipse and parabolic profiles. Three different parameter values for A are chosen. The entries within bracket shows the linear growth rate recovered from the nonlinear evolution of the perturbation. cot(θ) Re EIM(P) (ci ) EIM(E) A = 20 EIM(E) A = 40 EIM(E) A = 60 0.1 0.5 0.9 1.1 1.2 1.5 0.442 23 0.264 88 0.105 21 0.033 616 −0.000 000 −0.092 19 0.442 41 (0.426) 0.265 04 (0.247) 0.105 35 (0.091) 0.033 737 (0.023) 0.000 113 (−0.011) −0.092 10 (−0.115) 0.442 27 (0.432) 0.264 92 (0.268) 0.105 24 (0.105) 0.033 646 (0.023) 0.000 028 (−0.000) −0.092 17 (−0.115) 0.442 249 (0.426) 0.264 900 (0.246) 0.105 228 (0.091) 0.033 630 (0.023) 0.000 113 (−0.010) −0.092 179 (−0.115) and phase speed predicted by laboratory experiments and numerical solutions of the Orr-Sommerfeld system. In order to assess the validity of our model in the nonlinear regime, in this section, we try to understand the existence and properties of the nonlinear twodimensional waves that the model generates beyond threshold. This is achieved by obtaining the full numerical solutions to the coupled nonlinear evolution equations in a computation domain of length 2π with uniform base state of unit thickness. The computations are performed with 1000 mesh-points and with finer mesh at the periodic boundaries. The growth of an initial sinusoidal perturbation to the base state (H = 1) is studied by using the finite element method (FEM) to numerically approximate the actual solutions of the coupled nonlinear PDEs through COMSOL 5.1. The growth of the perturbations is obtained using ξ = ln[∥H(x, t) − 1∥]/ξ 0 , here ∥⋅∥ is the L2 -norm defined over the computational domain and ξ 0 = ∥H(x, t = 0) − 1∥ is the magnitude of the perturbation given to the primary flow. An absolute tolerance of 10−7 is defined to ensure the accuracy of the solutions to the time dependent evolution equations. A periodic type boundary condition is specified at the edges of the domain to maintain the continuation in the phases of the waves propagating outside and entering inside the computational domain. Figure 6(a) presents the nonlinear growth rate for different valθ . The comparison of the slope of the nonlinear growth ues of cot Re rate with the linear growth rate (linear fit; black squares), using ξ = ln[∥H(x, t) − 1∥/ξ 0 ] we infer from Table III and Fig. 6(b) that linear stability predicts the growth accurately at the initial time. The results reveal that the growth rate as predicted by nonlinear simulation and linear stability theory agrees very well. Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing There is a natural curiosity to understand the performance of the EIM(E) model as compared to the EIM(P) model for moderate to high Reynolds numbers. Note that both the models give satisfactory results for instability threshold. The basic difference in the models is that while a parabolic profile is employed as a weight function in EIM(P), EIM(E) considers an ellipse profile as a weight function and there is a free parameter A related to eccentricity in this model. With this in view, the energy contributions from the inertial terms [T 1 in Eq. (33), present study] in the EIM(E) model and [t 1 in Eq. (26) in the article by Usha and Uma47 ] in EIM(P) are plotted as a function of the streamwise coordinate (x) in Figs. 7(a) θ . It is interesting to note that and 7(b) for two different values of cot Re FIG. 7. The energy contribution from the inertial terms in EIM with parabolic (blue = 0.1 and (b) circles) and ellipse velocity (black solid line) profiles. (a) cot(θ) Re cot(θ) Re = 0.5, when We = 1 and A = 40. 32, 013603-10 Physics of Fluids ARTICLE scitation.org/journal/phf FIG. 8. Film thickness as a function of the streamwise coordinate x, for We = 1, Re = 1/ϵ, and ϵ ≪ 1. (a) and (c) cot(θ) Re ≙ 1.1. The ≙ 0.1 and (b) and (d) cot(θ) Re wave profiles are displayed at nondimensional times with step-size Δt = 0.1 in (a) and (b) and Δt = 1 in (c) and (d). contribution to the absolute value of energy from inertial terms in the EIM(E) model is more than that from the EIM(P) model. For any inclination θ of the substrate with 0 < θ < π/2, there is no significant difference between the EIM(E) and EIM(P) results for a lower Re [Fig. 7(b)], but at a higher Re [Fig. 7(a)], when inertial effects become stronger, EIM(E) results for energy contribution to the perturbation at any streamwise location is enhanced. It is of interest to understand if this enhancement in energy contribution helps in improving the performance of the EIM(E) model as compared to other models in the inertia dominated regimes, and this will be pursued in a future study. Table III presents the comparison of the growth rate from the two models with ellipse [EIM(E)] and parabolic [EIM(P)] profiles. The results are shown for three different values of A (A = 20, 40, and 60). The terms within brackets in each entry show the linear growth rate values recovered from the nonlinear evolution of the perturbation. The results also suggest that the choice of A = 40 at this We is justified. Figure 8 shows the film thickness as the disturbance waves evolve as a function of the stream-wise coordinate, at different θ ≙ 0.1 and Figs. 8(b) nondimensional times [Figs. 8(a) and 8(c): cot Re cot θ and 8(d): Re ≙ 1.1], when We = 1 and Re = 1/ϵ, ϵ << 1. As time elapses, the amplitude of the perturbations decays and the thickness of the film approaches the base state. As remarked by Liu et al.,70 this instability is of the convective type; as a result, the perturbation leaves the computational domain soon after it is introduced. The transition to absolute instability is not possible to be captured in the laboratory frame of reference using long wave70 analysis and hence remains beyond the scope of this study. Figure 9 presents the maximum and minimum amplitude of the film thickness as a function of nondimensional time for differθ . There is a decrease in the absolute deviation from ent values of cot Re θ increases indicating flow stabilization due to the base state as cot Re decay in the imposed perturbation. This behavior was also reflected in Fig. 8 as time progresses. VI. CONCLUSION FIG. 9. Maximum [Hmax , (a)] and minimum [Hmin , (b)] amplitude of the film thickness as a function of dimensionless time, when We = 1, Re = 1/ϵ, ϵ ≪ 1, and A = 40. Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing A new model for describing the dynamics of a thin film down an incline at high Reynolds numbers and moderate Weber numbers has been derived using the energy integral method with ellipse profile as a weight function. A free parameter appearing in the ellipse profile is related to the eccentricity of ellipse, and it allows one to refine the velocity profile which in turn modifies the performance of the model in the inertia dominated regimes. The goal is to assess the model’s capacity in predicting the linear instability threshold, and the results reveal that the model successfully predicts the instability threshold accurately and it agrees with earlier theoretical7,8 32, 013603-11 Physics of Fluids and experimental predictions.70 Furthermore, the remark by Luchini and Charru59 that the kinetic-energy balance equation is a consistent nonlinear formulation gives us confidence in the performance of the current model EIM(E) in predicting the nonlinear dynamics. Furthermore, the success of EIM with ellipse profile in predicting accurately the normal force on a squeeze film for large inertial effects and for large amplitude motion enhances our expectation about the performance of the present EIM(E) model, in the inertia dominated regimes. As a first step toward assessing the validity of the present EIM(E) model in the nonlinear regime, the numerical solution of the nonlinear evolution Eq. (37) describing the propagation of nonlinear waves is performed, and the results reveal that the amplitude of the perturbations decay as time progresses and this is in line with similar results in the work of Sadiq and Usha.71 The growth rate predicted by the linear stability theory matches well with that extracted from the nonlinear evolution, thereby indicating that the eigenvalue analysis is able to provide physical insights about the instability in such flows. In summary, the analysis in the present study shows that the current EIM(E) model’s performance is satisfactory both in predicting the linear instability threshold and characteristics of the nonlinear waves at the moderate distance from the threshold. The value of the free parameter A that is intrinsic in the ellipse profile taken as the weight function considered in the present study (namely, A = 40) is based on the closeness of the linear growth rate and phase speed for this value of A to the available experimental observations. This choice has also been justified by the values of linear growth rate derived from nonlinear evolution (as is evident from Table III). The model’s performance is as good as the other available models such as EIM(P)47 and WRIBL.33 A detailed and thorough investigation on the performance of the model predicting wave shapes and amplitudes as compared to other existing models is postponed to a future study. ACKNOWLEDGMENTS The authors sincerely thank the referees for their very valuable comments and suggestions. These have helped to improve the quality and content of the manuscript. The authors also thank Professor A. Jeffrey Giacomin, Editor-in-chief, for his very encouraging and useful remarks. REFERENCES 1 S. Veremieiev and D. H. Wacks, “Modelling gravity-driven film flow on inclined corrugated substrate using a high fidelity weighted residual integral boundarylayer method,” Phys. Fluids 31, 022101 (2019). 2 E. Ellaban, J. Pascal, and S. D’Alessio, “Instability of a binary liquid film flowing down a slippery heated plate,” Phys. Fluids 29, 092105 (2017). 3 S. Ghosh and R. Usha, “Stability of viscosity stratified flows down an incline: Role of miscibility and wall slip,” Phys. Fluids 28, 104101 (2016). 4 M. Dauth and N. Aksel, “Breaking of waves on thin films over topographies,” Phys. Fluids 30, 082113 (2018). 5 S. D’Alessio, C. Seth, and J. Pascal, “The effects of variable fluid properties on thin film stability,” Phys. Fluids 26, 122105 (2014). 6 P. Kapitza and S. Kapitza, “Wave flow of thin layers of viscous liquids. Part III. Experimental research of a wave flow regime,” Zh. Eksp. Teor. Fiz. 19, 105–120 (1949). 7 T. B. Benjamin, “Wave formation in laminar flow down an inclined plane,” J. Fluid Mech. 2, 554–573 (1957). Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing ARTICLE scitation.org/journal/phf 8 C.-S. Yih, “Stability of liquid flow down an inclined plane,” Phys. Fluids 6, 321–334 (1963). 9 P. Bach and J. Villadsen, “Simulation of the vertical flow of a thin, wavy film using a finite-element method,” Int. J. Heat Mass Transfer 27, 815–827 (1984). 10 T. R. Salamon, R. C. Armstrong, and R. A. Brown, “Traveling waves on vertical films: Numerical analysis using the finite element method,” Phys. Fluids 6, 2202–2220 (1994). 11 B. Ramaswamy, S. Chippada, and S. Joo, “A full-scale numerical study of interfacial instabilities in thin-film flows,” J. Fluid Mech. 325, 163–194 (1996). 12 Y. Y. Trifonov, “Stability and bifurcations of the wavy film flow down a vertical plate: The results of integral approaches and full-scale computations,” Fluid Dyn. Res. 44, 031418 (2012). 13 N. Malamataris, M. Vlachogiannis, and V. Bontozoglou, “Solitary waves on inclined films: Flow structure and binary interactions,” Phys. Fluids 14, 1082–1094 (2002). 14 T. Nosoko and A. Miyara, “The evolution and subsequent dynamics of waves on a vertically falling liquid film,” Phys. Fluids 16, 1118–1126 (2004). 15 H.-C. Chang, E. Demekhin, and D. Kopelevich, “Nonlinear evolution of waves on a vertically falling film,” J. Fluid Mech. 250, 433–480 (1993). 16 E. Demekhin, I. Demekhin, and V. Y. Shkadov, “Solitons in flowing layer of a viscous fluid,” Izv. Akad. Nauk. SSSR, Mekh. Zhid. I Gaza 4, 9–16 (1983). 17 V. Y. Shkadov, “Wave conditions in the flow of a thin layer of a viscous liquid under the action of gravity,” Izv. Akad. Nauk. SSSR, Mekh. Zhidk. I Gaza 1, 43–51 (1967). 18 D. Benney, “Long waves on liquid films,” J. Math. Phys. 45, 150–155 (1966). 19 A. Pumir, P. Manneville, and Y. Pomeau, “On solitary waves running down an inclined plane,” J. Fluid Mech. 135, 27–50 (1983). 20 S. Joo, S. Davis, and S. Bankoff, “Long-wave instabilities of heated falling films: Two-dimensional theory of uniform layers,” J. Fluid Mech. 230, 117–146 (1991). 21 A. Oron, S. H. Davis, and S. G. Bankoff, “Long-scale evolution of thin liquid films,” Rev. Mod. Phys. 69, 931–980 (1997). 22 Y. Kuramoto and T. Tsuzuki, “Persistent propagation of concentration waves in dissipative media far from thermal equilibrium,” Prog. Theor. Phys. 55, 356–369 (1976). 23 G. Sivashinsky, “Diffusional-thermal theory of cellular flames,” Combust. Sci. Technol. 15, 137–145 (1977). 24 O. Takeshi, “Surface equation of falling film flows with moderate Reynolds number and large but finite Weber number,” Phys. Fluids 11, 3247–3269 (1999). 25 S. Alekseenko, V. Y. Nakoryakov, and B. Pokusaev, “Wave formation on a vertical falling liquid film,” AIChE J. 31, 1446–1460 (1985). 26 V. Nakoryakov, B. Pokusaev, S. Alekseenko, and V. Orlov, “Instantaneous velocity profile in a wavy fluid film,” J. Eng. Phys. Thermophys. 33, 1012–1016 (1977). 27 L. T. Nguyen and V. Balakotaiah, “Modeling and experimental studies of wave evolution on free falling viscous films,” Phys. Fluids 12, 2236–2256 (2000). 28 C. Ruyer-Quil and P. Manneville, “Further accuracy and convergence results on the modeling of flows down inclined planes by weighted-residual approximations,” Phys. Fluids 14, 170–183 (2002). 29 P. Rosenau, A. Oron, and J. Hyman, “Bounded and unbounded patterns of the Benney equation,” Phys. Fluids A 4, 1102–1104 (1992). 30 A. Oron and O. Gottlieb, “Nonlinear dynamics of temporally excited falling liquid films,” Phys. Fluids 14, 2622–2636 (2002). 31 B. Scheid, C. Ruyer-Quil, U. Thiele, O. A. Kabov, J. C. Legros, and P. Colinet, “Validity domain of the Benney equation including the Marangoni effect for closed and open flows,” J. Fluid Mech. 527, 303–335 (2005). 32 O. Gottlieb and A. Oron, “Stability and bifurcations of parametrically excited thin liquid films,” Int. J. Bifurcation Chaos 14, 4117–4141 (2004). 33 C. Ruyer-Quil and P. Manneville, “Improved modeling of flows down inclined planes,” Eur. Phys. J. B 15, 357–369 (2000). 34 A. Samanta, “Role of slip on the linear stability of a liquid flow through a porous channel,” Phys. Fluids 29, 094103 (2017). 35 A. Oron and C. Heining, “Weighted-residual integral boundary-layer model for the nonlinear dynamics of thin liquid films falling on an undulating vertical wall,” Phys. Fluids 20, 082102 (2008). 32, 013603-12 Physics of Fluids 36 C. Park and T. Nosoko, “Three-dimensional wave dynamics on a falling film and associated mass transfer,” AIChE J. 49, 2715–2727 (2003). 37 R. R. Mudunuri and V. Balakotaiah, “Solitary waves on thin falling films in the very low forcing frequency limit,” AIChE J. 52, 3995–4003 (2006). 38 S. Kalliadasis, C. Ruyer-Quil, B. Scheid, and M. G. Velarde, Falling Liquid Films (Springer Science & Business Media, New York, 2011). 39 V. Teshukov, “Gas-dynamic analogy for vortex free-boundary flows,” J. Appl. Mech. Tech. Phys. 48, 303–309 (2007). 40 G. L. Richard and S. L. Gavrilyuk, “A new model of roll waves: Comparison with Brock’s experiments,” J. Fluid Mech. 698, 374–405 (2012). 41 G. L. Richard and S. L. Gavrilyuk, “The classical hydraulic jump in a model of shear shallow-water flows,” J. Fluid Mech. 725, 492–521 (2013). 42 G. Richard, C. Ruyer-Quil, and J. Vila, “A three-equation model for thin films down an inclined plane,” J. Fluid Mech. 804, 162–200 (2016). 43 J. Liu and J. P. Gollub, “Solitary wave dynamics of film flows,” Phys. Fluids 6, 1702–1712 (1994). 44 J. Tihon, K. Serifi, K. Argyriadi, and V. Bontozoglou, “Solitary waves on inclined films: Their characteristics and the effects on wall shear stress,” Exp. Fluids 41, 79–89 (2006). 45 M. Amaouche, N. Mehidi, and N. Amatousse, “An accurate modeling of thin film flows down an incline for inertia dominated regimes,” Eur. J. Mech. B: Fluids 24, 49–70 (2005). 46 J.-J. Lee and C. C. Mei, “Stationary waves on an inclined sheet of viscous fluid at high Reynolds and moderate Weber numbers,” J. Fluid Mech. 307, 191–229 (1996). 47 R. Usha and B. Uma, “Modeling of stationary waves on a thin viscous film down an inclined plane at high Reynolds numbers and moderate Weber numbers using energy integral method,” Phys. Fluids 16, 2679–2696 (2004). 48 A. Elkouh, “Inertia effect in laminar radial flow between parallel plates,” Int. J. Mech. Sci. 9, 253–255 (1967). 49 S. Crandall and A. Ei-Shafei, “Momentum and energy approximations for elementary squeeze-film damper flows,” J. Appl. Mech. 60, 728–736 (1993). 50 Y. Han and R. Rogers, “Squeeze film force modeling for large amplitude motion using an elliptical velocity profile,” J. Tribol. 118, 687–692 (1996). 51 D. C. Kuzma, “Fluid inertia effects in squeeze films,” Appl. Sci. Res. 18, 15–20 (1968). 52 V. Kapur and K. Verma, “Energy integral approach for MHD hydrostatic thrust bearing,” J. Lubr. Technol. 97, 647–649 (1975). 53 S. Turns, “Annular squeeze films with inertial effects,” J. Lubr. Technol. 105, 361–363 (1983). 54 R. Usha and P. Vimala, “Squeeze film force using an elliptical velocity profile,” J. Appl. Mech. 70, 137–142 (2003). Phys. Fluids 32, 013603 (2020); doi: 10.1063/1.5127815 Published under license by AIP Publishing ARTICLE scitation.org/journal/phf 55 H. A. Abderrahmane and G. Vatistas, “Improved two-equation model for thin layer fluid flowing down an inclined plane problem,” Phys. Fluids 19, 098106 (2007). 56 C. Ruyer-Quil, P. Treveleyan, F. Giorgiutti-Dauphiné, C. Duprat, and S. Kalliadasis, “Modelling film flows down a fibre,” J. Fluid Mech. 603, 431–462 (2008). 57 E. Novbari and A. Oron, “Energy integral method model for the nonlinear dynamics of an axisymmetric thin liquid film falling on a vertical cylinder,” Phys. Fluids 21, 062107 (2009). 58 E. Novbari and A. Oron, “Analysis of time-dependent nonlinear dynamics of the axisymmetric liquid film on a vertical circular cylinder: Energy integral model,” Phys. Fluids 23, 012105 (2011). 59 P. Luchini and F. Charru, “The phase lead of shear stress in shallow-water flow over a perturbed bottom,” J. Fluid Mech. 665, 516–539 (2010). 60 I. M. R. Sadiq, “First-order energy-integral model for thin Newtonian liquids falling along sinusoidal furrows,” Phys. Rev. E 85, 036309 (2012). 61 N. Tiwari and J. M. Davis, “Nonmodal and nonlinear dynamics of a volatile liquid film flowing over a locally heated surface,” Phys. Fluids 21(10), 102101 (2009). 62 L. O. Jones and S. Whitaker, “An experimental study of falling liquid films,” AIChE J. 12, 525–529 (1966). 63 F. Stainthorp and J. Allen, “The development of ripples on the surface of liquid film flowing inside a vertical tube,” Trans. Am. Inst. Chem. Eng. 43, 85 (1965). 64 W. Strobel and S. Whitaker, “The effect of surfactants on the flow characteristics of falling liquid films,” AIChE J. 15, 527–532 (1969). 65 W. B. Krantz and S. Goren, “Stability of thin liquid films flowing down a plane,” Ind. Eng. Chem. Fundam. 10, 91–101 (1971). 66 F. W. Pierson and S. Whitaker, “Some theoretical and experimental observations of the wave structure of falling liquid films,” Ind. Eng. Chem. Fundam. 16, 401–408 (1977). 67 C. Ruyer-Quil and P. Manneville, “Modeling film flows down inclined planes,” Eur. Phys. J. B 6, 277–292 (1998). 68 B. Scheid, S. Kalliadasis, C. Ruyer-Quil, and P. Colinet, “Interaction of threedimensional hydrodynamic and thermocapillary instabilities in film flows,” Phys. Rev. E 78, 066311 (2008). 69 B. Scheid, C. Ruyer-Quil, and P. Manneville, “Wave patterns in film flows: Modelling and three-dimensional waves,” J. Fluid Mech. 562, 183–222 (2006). 70 J. Liu, J. D. Paul, and J. P. Gollub, “Measurements of the primary instabilities of film flows,” J. Fluid Mech. 250, 69–101 (1993). 71 I. M. R. Sadiq and R. Usha, “Thin Newtonian film flow down a porous inclined plane: Stability analysis,” Phys. Fluids 20, 022105 (2008). 32, 013603-13