Formation and evolution of roll waves in a shallow free surface flow of a power-law fluid down an inclined plane
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b Formation and evolution of roll waves in a shallow free surface flowof a power-law fluid down an inclined plane
Alexander Chesnokov a,b a Lavrentyev Institute of Hydrodynamics SB RAS, 15 Lavrentyev Ave., Novosibirsk 630090, Russia b Novosibirsk State University, 1 Pirogova Str., Novosibirsk 630090, Russia
Abstract
A laminar flow of a thin layer of mud down an inclined plane under the action of gravityis considered. The instability of a film flow and the formation of finite amplitude wavesare studied in the framework of both two-dimensional governing equations of a power-lawfluid and its depth-averaged hyperbolic simplification. Conditions of roll waves existence forthese models are formulated in terms of Whitham criterion. Numerical calculations of thefree surface evolution and the roll waves development are performed. It is shown that theroll waves amplitude obtained by the 2D equations is slightly larger than for the 1D model.Moreover, for certain flow parameters, small perturbations of the basic solution grow forthe 2D equations and decay for the depth-averaged model. A two-parameter class of exactpiecewise-smooth solutions of the 1D model is obtained and a comparison with a numericalsolution is made. In the region of these parameters, diagrams of the roll waves existence areconstructed.
Keywords: power-law fluid, thin films, roll waves, hyperbolic equations
1. Introduction
The gravity-driven flow of a power-law liquid film down an inclined plane is an importantproblem with many industrial applications (chemical engineering, coating processes) andnatural hazards (mud flows, debris flow impact). Steady flows of thin layers of Newtonian ornon-Newtonian fluid can become unstable for sufficiently high slopes against long-wavelengthinfinitesimal perturbations. Such long-wave instability frequently develops into a seriesof progressing bores, known as roll waves. The control of this instability concerns manytechnologies. For instance, a uniform flow thickness is required in coating processes and,consequently, the instability should be avoided. On the other hand, in chemical engineeringsuch waves can be useful since they enhance heat and mass exchange. In environmentalhydraulics, muddy intermittent flows are frequently encountered in mountainous regions,especially after torrential rains. Moreover, these type of waves can also be observed as anevent following volcano eruptions.Since the remarkable work of Kapitza [17], wave evolution of liquid films flowing downan inclined plane has attracted a considerable number of studies [1, 8]. Direct numerical
Email address: [email protected] (Alexander Chesnokov)
Preprint submitted to Elsevier March 1, 2021 imulation of such flows is difficult and requires large computational costs. Therefore, depth-averaged models are widely used to describe and study the film flow evolution. Existingtheoretical models of thin layer motion of Newtonian and non-Newtonian fluid are based onsimplifications of the governing equations obtained from vertical averaging. The procedure ismuch like the von K´arm´an–Pohlhausen power integral method used in the classical theory ofboundary layers, and proceeds by first prescribing a fixed vertical structure to the flow. Thenthe governing equations can be vertically integrated to remove entirely one of the spatialvariables. In this case the shear stress is replaced by a drag term as in the Saint-Venantmodel used in hydraulics. For the same reasons, one uncovers roll-wave instability, themathematical theory of which was developed by Dressler [15] for shallow flows in inclinedchannels. Shkadov [28] was the first who applied this approach to obtain a two-equationmodel for laminar film flows of a viscous liquid. Many other researchers adopted his idea tomore complex geometries or non-Newtonian fluids. Ng and Mei [22] investigated the laminarflow of a shallow layer of mud by using a shear thinning power-law rheological model andapplied a theory of roll waves for such flows. The behaviour and stability of a thin layer ofnon-Newtonian fluids accounting for the presence of yield stress was studied by Liu and Mei[19] for the Bingham rheology. Roll waves in mud modelled as a viscoplastic fluid with theHerschel–Bulkley constitutive law was considered by Balmforth and Liu [2].Numerous subsequent works are aimed at applying the theory of roll waves for varioushydrodynamic problems, as well as constructing depth-averaged models, their theoreticalanalysis and use for modelling the evolution of finite amplitude waves. For instance, in [11] afluid flow in a long tube with elastic walls is considered and shown that for certain parametersof the flow small perturbations at the inlet section of the tube give rise to roll waves. Thestability of two-layer flows in open and closed channels, as well as in a Hele–Shaw cell, isconsidered in [6, 3] and [12], where the possibility of the roll waves formation in such systemsis shown. In [4, 5] the modulation equations for free surface shallow flows are derived andthe stability criterion of roll waves is formulated in terms of hyperbolicity of these equations.Starting with an expansion of the velocity field in terms of polynomial test functions andusing a weighted-residual technique based on a classical long-wavelength expansion morecomplex depth-averaged models for Newtonian and power-law fluid are derived and verifiedin [27, 26]. Consistent shallow water equations for the flow of thin films of power-law fluid aredeveloped in [20] through a second-order asymptotic expansion of the fluid velocity field andfluid strain. Roll waves in flows of thin films of a viscous liquid and shear shallow water flowsin inclined channels are studied in works [25] and [16] in the framework of advanced three-equation models. It is shown that these hyperbolic governing equations more accuratelyreproduce the experimental data in comparison with the widely used two-equation models.Last years, considerable attention has been paid to the study of the stability and rollwaves development in shallow flows of non-Newtonian fluids. Spatial evolution of a smalldisturbance in an open-channel flow of a power-law fluid at non-uniform initial conditions, upto the occurrence of roll waves in mild and steep slope channels is investigated in [7]. Withreference to the propagation of mud flows described by means of a power-law rheology, thepotentialities of different approximations of the linearised version of the two-equation model2nd the simplified kinematic-wave model are studied in [13, 14]. The linear stability of shear-thinning fluid down an inclined plane in the framework of the generalized Orr-Sommerfeldequation is performed in [23]. The influence of a prescribed superficial shear stress on thegeneration and structure of roll waves developing from infinitesimal disturbances on thesurface of a power-law fluid layer is investigated in [24]. The effect of oblique perturbationsthat propagate at an arbitrary angle to the velocity of the undisturbed fluid flow moving downan incline plane under gravity is studied in [32]. It is shown that under certain conditions,oblique perturbations can lead to the flow instability.In the present paper we consider the occurrence of roll waves in a free surface shallowflow of a power-law fluid using both the two-dimensional boundary layer equations and thedepth-averaged hyperbolic model proposed by Ng and Mei [22]. A generalized theory ofcharacteristics and the notion of hyperbolicity for integro-differential equations introducedby Teshukov [29] (see also [10, 9]) allow one to apply Whitham’s stability criterion toboth considered models. Numerical simulations are carried out in order to confirm theachievements of the theoretical analysis as well as to compare the results of calculationsusing fully nonlinear models of different levels of complexity. We also construct and study atwo-parameter class of exact solutions for the 1D model describing the roll waves regime. Itis shown that with an appropriate choice of the critical depth and amplitude of the wave, thissolution is consistent with the results of numerical calculations. Diagrams of the roll wavesexistence that bound all possible values of the critical depth and fluid depth are obtained.The paper is structured as follows. In Section 2 the governing equations are reported.Then the Whitham stability criterion is formulated in Section 3 for both 1D and 2D models.For numerical treatment of the 2D equations in Section 4, a multilayer approximation isderived in the form of 1D balance laws. Numerical simulation of the development of roll wavesinstability using the considered models is performed in Section 5. A two-parameter class ofpiecewise-smooth travelling wave solutions is studied in Section 6. Finally, conclusions aredrawn in Section 7.
2. Mathematical model
We consider a power-law fluid flowing down an inclined plane under the action of gravity.The flow is assumed to be incompressible and the fluid density ρ , the angle of inclination θ and the gravity acceleration g are constant. An ( x, z ) coordinate system is defined with the x -axis along and the z -axis normal to the plane bed. The velocity field and the total pressureare denoted by v = ( u, w ) T and p , respectively. The dimensional governing equations read ∂u∂t + u ∂u∂x + w ∂u∂z + 1 ρ ∂p∂x = 1 ρ (cid:16) ∂τ xx ∂x + ∂τ xz ∂z (cid:17) + g sin θ,∂w∂t + u ∂w∂x + w ∂w∂z + 1 ρ ∂p∂z = 1 ρ (cid:16) ∂τ zx ∂x + ∂τ zz ∂z (cid:17) − g cos θ,∂u∂x + ∂w∂z = 0 , τ ij = 2 µ n (2 D kl D kl ) ( n − / D ij , (1)where D ij are the components of strain-rate tensor D = ( ∇ v + ( ∇ v ) ∗ ) /
2, positive constants n and µ n are the power-law index and the fluid consistency respectively. The sum convention3s employed in the last formula in Eq. (1) and, consequently, the terms τ ij have the form τ xx = 2 τ n u x , τ xz = τ zx = τ n ( u z + w x ) , τ zz = 2 τ n w z ,τ n = µ n (cid:0) u x + 2 w z + ( u z + w x ) (cid:1) ( n − / . If the fluid behaviour is shear-thinning then n <
1, whereas n > n = 1, the Newtonian case is recovered, with the consistency µ n coinciding with the dynamic viscosity. We note that power-law fluids serve as a simple butan efficient model for a broad variety of non-Newtonian fluids. The system of equations iscompleted by the boundary conditions at the wall z = 0 and at the free surface z = h : u (cid:12)(cid:12) z =0 = 0 , w (cid:12)(cid:12) z =0 = 0 , ∂h∂t + u ∂h∂x − w (cid:12)(cid:12)(cid:12) z = h = 0 , ( τ xx − p ) ∂h∂x − τ xz (cid:12)(cid:12)(cid:12) z = h = 0 , τ xz ∂h∂x + τ zz − p (cid:12)(cid:12)(cid:12) z = h = 0 . (2)The surface tension is neglected here.System (1), (2) admits the Nusselt-type solution u ( z ) = f ( z ) U , f ( z ) = 1 + 2 n n (cid:18) − (cid:16) − zh (cid:17) /n (cid:19) , (3)where U = n n (cid:18) ρg sin θµ n (cid:19) /n H /n , (4)and w = 0, p = ( h − z ) g cos θ , h = H = const. Indeed, direct calculations show that theseformulae give an exact solution to equations (1) with boundary conditions (2). Velocityprofile (3) and relation (4) will be used further.Let us now write this system in a non-dimensional form and in the shallow-water scaling.We assume that the characteristic scales of the film depth, its velocity and longitudinalwavelength are H , U and L = ( g sin θ ) − U , respectively. The shallowness of the flowis determined by the aspect ratio ε = H /L ≪
1. In addition, it is assumed that thecharacteristic velocity U and the flow depth H are related by formula (4). Then we scalethe variables as follows x = L x ′ , ( h, z ) = H ( h ′ , z ′ ) , t = L U t ′ , u = U u ′ , w = εU w ′ ,p = ρU p ′ , ( τ xx , τ zz ) = εµ n (cid:18) U H (cid:19) n ( τ ′ xx , τ ′ zz ) , τ xz = µ n (cid:18) U H (cid:19) n τ ′ xz . Following [22], we defined the non-dimensional Froude and Reynolds numbers as
F r = U gH = sin θε , Re = U ρµ n (cid:18) H U (cid:19) n = 1 εK , where K = ( n/ (1 + 2 n )) n . 4hen, by dropping the primes, one may rewrite the momentum equations and thedynamic boundary conditions in the form ∂u∂t + u ∂u∂x + w ∂u∂z + ∂p∂x = K (cid:16) ε ∂τ xx ∂x + ∂τ xz ∂z (cid:17) + 1 ,ε (cid:16) ∂w∂t + u ∂w∂x + w ∂w∂z (cid:17) + ∂p∂z = ε K (cid:16) ∂τ xz ∂x + ∂τ zz ∂z (cid:17) − ε cot θ ; ε K (cid:0) τ xx − p (cid:1) ∂h∂x − Kτ xz (cid:12)(cid:12)(cid:12) z = h = 0 , ε K (cid:16) τ xz ∂h∂x + τ zz (cid:17) − p (cid:12)(cid:12)(cid:12) z = h = 0 . (5)The rest equations in formulae (1) and (2) are unchanged. Neglecting terms of order ε inEq. (5), we obtain a simplified version of the governing equations ∂u∂t + u ∂u∂x + w ∂u∂z + ∂p∂x = K ∂τ xz ∂z + 1 , ∂p∂z = − ε cot θ, ∂u∂x + ∂w∂z = 0; u (cid:12)(cid:12) z =0 = 0 , w (cid:12)(cid:12) z =0 = 0 , h t + uh x − w (cid:12)(cid:12)(cid:12) z = h = 0 , τ xz (cid:12)(cid:12) z = h = 0 , p (cid:12)(cid:12) z = h = 0 . (6)In this approximation the fluid pressure obeys the hydrostatic law p = ( h − z ) ε cot θ and thedimensionless stress tensor component τ xz has the form τ xz = | u z | n − u z . Below we assumethat u z ≥
0. Taking this into account, Eq. (6) can be written as ∂u∂t + u ∂u∂x + w ∂u∂z + α ∂h∂x = K ∂∂z (cid:18) ∂u∂z (cid:19) n + 1 ,∂h∂t + ∂∂x (cid:18) Z h u dz (cid:19) = 0 , u (cid:12)(cid:12) z =0 = 0 , ∂u∂z (cid:12)(cid:12)(cid:12) z = h = 0 . (7)where α = ε cot θ, w = − Z z u ( t, x, z ′ ) dz ′ . Let us integrate the first equation of system (7) with respect to z under assumption thatthe velocity profile has the form u = f ( z )¯ u ( t, x ), where f ( z ) is given by formula (3) with h = h ( t, x ). Profiles of f ( z ) for different n are plotted in Fig. 1. Using other equations ofsystem (7), we obtain the following depth-averaged model [22]:(¯ uh ) t + (cid:0) β ¯ u h + αh / (cid:1) x = h − (cid:18) ¯ uh (cid:19) n , h t + (¯ uh ) x = 0 . (8)Here β = 1 h Z h f ( y ) dy = 2(1 + 2 n )2 + 3 n > , unknowns ¯ u ( t, x ) and h ( t, x ) > α ≥ n > .5 1.0 1.50.20.40.60.81.00 n =1.00.20.6 z / h u u / - l l l -+ < h l h * Figure 1: Nusselt-type velocity profiles u/ ¯ u = f ( z )given by formula (3) for power-law index n = 1, 0.6and 0.2. Figure 2: The ‘frozen’ λ ± ( h ) and ‘equilibrium’ˆ λ ( h ) characteristic velocities for α = 1, n = 0 . c ± ( h ) of Eq. (7) are shown bydashed curves.
3. Whitham criterion
It is easy to verify that system (8) is hyperbolic. There are two families of characteristicswith local slopes λ ± = β ¯ u ± p β ( β − u + αh. (9)Together with Eq. (8) we consider more simple kinematic model where inertial effects are nottaken into account. Then the first equation in (8) is replaced by the condition ¯ u = h /n .Substitution of this expression in the second equation of (8) yields h t + ( ϕ ( h )) x = 0 , ϕ ( h ) = h /n . (10)The velocity of characteristic of the kinematic-wave model (10) isˆ λ = ϕ ′ ( h ) = 2 n + 1 n h /n . Eq. (10) is used for formulation a criterion of the roll waves existing for system (8). Thecriterion was first described by Whitham in [31] for the shallow water flows in inclined openchannels. This criterion is based on comparison of ‘frozen’ characteristic velocities λ ± ( h )of system (8) with ‘equilibrium’ characteristic speed ˆ λ ( h ) of Eq. (10). The variables λ ± ( h )are obtained by substitution of ¯ u = ϕ ( h ) /h into the right-hand side of (9). According toWhitham criterion [31] the roll waves exist in the intervals where one of inequalities is validˆ λ ( h ) > λ + ( h ) or ˆ λ ( h ) < λ − ( h ) . (11)A typical behaviour of the ‘frozen’ λ ± ( h ) and ‘equilibrium’ ˆ λ ( h ) characteristic velocities areshown in Fig. 2 by solid curves. The figure is obtained for α = 1 and n = 0 . λ = λ + ( h ) and λ = ˆ λ ( h ) intersect at the point h = h ∗ = (cid:18) αn n (cid:19) n/ (2+ n ) . (12)6herefore, in accordance with condition (11) and formula (12), the constant solution h = h ,¯ u = h /n is stable if h < h ∗ and unstable otherwise. Without loss of generality, one canchoose h = 1. Then the stability condition takes the form α > α n = 1 + 2 nn . (13)Inequality (13) coincides with the condition obtained in [22] based on the linear stabilityanalysis of the constant solution h = 1, ¯ u = 1. As far as we know, linear analysis andWhitham’s criterion always give the same result for hyperbolic two-equation systems (see[11, 12] for instance).A generalized theory of characteristics for integro-differential equations was introduced byTeshukov [29] (see also [30, 10, 9]), who developed new mathematical tools for the qualitativestudy of such systems. As it follows from the cited works, characteristic velocities c = c ( t, x )of the integro-differential system (7) on a solution u ( t, x, z ), h ( t, x ) are determined by theequation χ ( c ) = 1 − α Z h dz ( u − c ) = 0 . (14)It is known that for flows with a monotonic velocity profile (e.g. 0 ≤ u z < ∞ ) this equationhas exactly two real roots c = c − < u ( t, x,
0) and c = c + > u ( t, x, h ). Let us calculate thecharacteristic velocities c ± on the Nusselt-type solution u = f ( z ) h /n , h = const of Eq. (7).Substitution of this solution into Eq. (14) gives an implicit equation of the form1 − αh (1 + n )( c − u s ) c (cid:0) n + F (1 , κ ; u s /c ) (cid:1) = 0 , (cid:18) u s = κh /n , κ = 1 + 2 n n (cid:19) (15)which determines the characteristic velocities c = c ± ( h ). Here F ( a, b ; c ; z ) is the Gaussianhypergeometric function. For α = 1 and n = 0 . c = c ± ( h ) are shown inFig. 2 by dashed lines. As we can see in Fig. 2, the curves λ = λ + ( h ) and λ = c + ( h ) arelocated quite close. The same is true for other values of the parameters α and n . We alsonote that the ‘equilibrium’ simplification of system (7) coincides with the kinematic-waveequation (10). This allows one to apply Whitham’s criterion to study the stability of theNusselt-type solution of integro-differential system (7). To do this, we can use formulae (11),where λ ± ( h ) should be replaced by c ± ( h ). Since the function f ( z ) is convex (its secondderivative does not change sign on the interval z ∈ (0 , h ) for n = 1), according to [9] thecharacteristic equation (14) has no complex roots (such that Im ( c ) = 0).Let us note, that substitution of h = 1 and c = ˆ λ (1) = 2 + 1 /n in Eq. (15) yields ananalogue of the stability condition (13) in the form α > α ∗ n = (1 + 2 n ) α n n + F (1 , κ ; n/ (1 + n )) , (16)where α n is defined in (13). Thus, the Nusselt-type solution u = f ( z ), h = 1 of Eq. (7)is stable if condition (16) is satisfied. We note that α ∗ n > α n for all n >
0. For instance, α n ≈ . α ∗ n ≈ .
47 at n = 0 .
7, and α n = 3, α ∗ n ≈ .
50 at n = 1. This circumstance7llows us to construct examples in which small perturbations of the Nusselt-type solution ofsystem (7) lead to the formation of finite amplitude waves, while the corresponding constantsolution of the depth-averaged model (8) turns out to be stable. Further, we carry outnumerical calculations based on both considered models and demonstrate the Whithamcriterion applicability.
4. Multilayer approximation
The construction of a numerical solution to the hyperbolic system of balance laws (8)does not cause difficulties and can be performed according to Godunov-type schemes. Thisapproach is not directly applicable to integro-differential model (7). However, a multilayerapproximation of shear flow equations proposed in [30, 10] allows one to treat model (7) asa system of one-dimensional balance laws.As it was shown in [18, 10], the evolution of the smooth solution of system (7) caninvolve a gradient catastrophe. The further description of the solution is possible onlyin the class of discontinuous functions. It leads to the necessity to formulate this modelin the form of balance laws. To do this, we present Eq. (7) in conservative form usingsemi-Lagrangian coordinates. In these variables Eq. (7) are reduced to a quasilinear integro-differential system [30, 10] u t + uu x + α Z H x dξ = 1 + KH ∂∂ξ (cid:18) u ξ H (cid:19) n , H t + ( uH ) x = 0 . (17)with boundary conditions u | ξ =0 = 0 and u ξ | ξ =1 = 0. Here ξ ∈ [0 ,
1] is Lagrangian coordinate, H = Φ ξ > z = Φ( t, x, ξ ), function Φ is a solutionof the Cauchy problemΦ t + u ( t, x, Φ)Φ x = w ( t, x, Φ) , Φ | t =0 = ξh (0 , x ) . System (17) can be rewritten in an equivalent conservative form (see [30, 10] for details) H t + ( uH ) x = 0 , ( ωH ) t + ( uωH ) x = K ∂∂ξ (cid:18) H ∂ω n ∂ξ (cid:19) ,∂∂t Z uH dξ + ∂∂x (cid:18) Z u H dξ + αh (cid:19) = h − Kω n (cid:12)(cid:12) ξ =0 , ω (cid:12)(cid:12) ξ =1 = 0 , (18)where h = Z H dξ, u = Z ξ ωH dξ ′ . The first equation in (18) is the local mass conservation law, the second and last equationscorrespond to the local and total horizontal momentum balance laws.Let us divide the interval [0 ,
1] into M subintervals 0 = ξ < ξ < ... < ξ M = 1 andintroduce new variables z i = Φ( t, x, ξ i ) , u i = u ( t, x, ξ i ) , h i = z i − z i − ,ω i = u i − u i − h i , ¯ u i = u i + u i − , Q = M X i =1 ¯ u i h i . (19)8hen we integrate Eq. (18) with respect to ξ over the intervals ( ξ i − , ξ i ), i = 1 , ..., M . Doingso, we take into account the equality Hdξ = dz , boundary conditions u = 0, ω M = 0 andassume a piecewise linear approximation for horizontal velocity u ( t, x, z ) ≈ ω i ( t, x )( z − z i − ) + u i − , z ∈ [ z i − , z i ] . As a result, for 2 M − h , ..., h M , ω , ..., ω M − , Q ), which are layers’depths h i , vorticities in the layers ω i , and a total fluid rate Q , we obtain the followingsystem of balance laws: ∂h i ∂t + ∂∂x (cid:0) ¯ u i h i (cid:1) = 0 , ( i = 1 , ..., M ) ∂∂t (cid:0) ω i h i (cid:1) + ∂∂x (cid:0) ¯ u i ω i h i (cid:1) = K (cid:18) ω ni +1 − ω ni h i − ω ni − ω ni − h i − (cid:19) , ( i = 2 , ..., M − ∂Q∂t + ∂∂x (cid:18) M X i =1 (cid:16) ¯ u i h i + ω i h i (cid:17) + αh (cid:19) = h − Kω n , (20)where h = M X i =1 h i , ¯ u i = − ω i h i i X j =1 ω j h j ,ω = 2(2 h − h ) h (cid:18) Q + M X i =2 h i (cid:16) ω i h i − i X j =2 ω j h j (cid:17)(cid:19) . To solve the differential conservation laws (20) numerically, one can apply standard methodsbased on various modifications of Godunov’s scheme. In our case, due to a large numberof equations in system (20), it is convenient to use central schemes, which do not requireexact or approximate solution of the Riemann problem. In this work we implement theNessyahu–Tadmor second-order central scheme [21].
5. Numerical simulation of roll waves
In this Section we perform numerical modelling of the finite amplitude waves developmentfrom infinitesimal perturbation using both multilayer equations (20) and depth-averagedmodel (8). Let at the initial time t = 0 the fluid depth h = 1 and depth-averaged velocity¯ u = 1. Note that this is a solution to Eq. (8). For multilayer model (20) at t = 0 we choose Q = 1, z i = i/M , u i = f ( z i ) ( i = 0 , ..., M ) and according to formulae (19) one can determinevalues of h i and ω i . We recall that f ( z ) is the Nusselt-type profile given by (3). We take thenumber of layers M = 20. On the left boundary x = 0 the constant flow rate is disturbedas follows Q | x =0 = ¯ uh | x =0 = 1 + a p sin(Ω t ) (21)with a p = 0 .
005 and Ω = 2 π . For unknown variables U we set the following conditions U N = U N − on the right boundary of the computational domain x = L . Here U j is thevalue of vector-function U in the nodal point x j . An uniform grid with respect to the9
10 20 30 40 50 600.80.91.01.11.21.3 x h x h a b Figure 3: Free surface z = h at t = 35, obtained by Eq. (8) (solid curves) and Eq. (20) (dashed curves) for α = 1. Power-law index n = 1 (a) and n = 0 . n = 1 n = 0 . n = 0 . n = 0 . n = 0 . n = 0 . δa ×
10 2 .
644 2 .
476 2 .
295 2 .
103 1 .
997 1 . δl × .
07 3 .
03 2 .
85 2 .
52 2 .
24 2 . x is used for calculations, the number of nodes N = 8000. The time step isdetermined by the Courant condition.Let us take α = 1 and carry out calculation for different values of the power-law index n ∈ [0 . , α < α n < α ∗ n hold, since α n ≥ α ∗ n ≥ . n , perturbations develop more intensively. Fig. 4corresponds to part of Fig. 3 (b) and demonstrates the form of the developed roll waves thatconsist of continuous flow regions connected by strong discontinuities. In this case the wavesmove with some positive velocity. In the continuity domain of the solution, a transitionfrom a subcritical to a supercritical flow occurs in the frame moving with wave velocity.The dash-dotted curve in Fig. 4 corresponds to the exact solution of Eq. (8) in the classof travelling waves that will be constructed in the next section. Note that the exact andnumerical solutions almost completely coincide for the developed roll waves.The wave amplitude obtained by Eq. (20) turns out to be higher than that calculatedby model (8), while the wavelengths practically coincide. The relative differences in theamplitudes δa ( n ) = | a ∗ − a | /a ∗ and wavelengths δl ( n ) for different values of n are shown inTable 1. Here ( a ∗ , l ∗ ) and ( a, l ) are the amplitude and wavelength of the roll waves accordingto Eq. (20) and (8), respectively. In both cases, small perturbations are given in the form (21)10 x h x h Figure 4: Part of Fig. 3 (b) on the interval x ∈ (50 , z = h at t = 90 obtained byEq. (8) (solid curve) and Eq. (20) (dashed curve) for n = 1 and α = 3 . with a p = 0 .
005 and Ω = 2 π . As follows from Table 1, the differences between these wavesbecome smaller with decreasing of the power-law index n .The probable reason for this is that for small values of n , the velocity profile weaklydepends on the variable z everywhere, with the exception of a small neighbourhood of thebottom (see Fig. 1). That is why the calculation results using depth-averaged equation (8)and multilayer approximation (20) become closer for small n . We also note that the indicatedfeatures of the velocity profile complicate the calculations by Eq. (20) for highly non-Newtonian fluid. This is due to the fact that the values of ω i are close to zero in asignificant number of layers and to perform calculations it is necessary to apply a moreaccurate resolution in the variable x .As noted above, there is a difference between conditions (16) and (13), which ensure thestability of a Nusselt-type solution within the framework of models (8) and (7), respectively.In particular, α n = 3 and α ∗ n ≈ . n = 1. Let us take α n < α < α ∗ n and carryout numerical simulation of the evolution of small unsteady perturbations given on the leftboundary based on Eq. (8) and multilayer system (20). We choose α = 3 . n = 1.The calculation results at t = 90 for perturbations (21) with parameters a p = 0 . . π are shown in Fig. 5. As we can see in the figure, small perturbations of the basicsolution are reduced when we use Eq. (8). On the contrary, the same perturbations growif we apply Eq. (20). Thus, if the parameter α belongs to the interval ( α n , α ∗ n ), then theresults of calculations using two-dimensional equation (7) and depth-averaged model (8) arequalitatively different.
6. Travelling waves
The construction of a periodic piecewise-smooth solution of Eq. (8), called roll waves,in the class of travelling waves was carried out in [22]. Here we focus on the fact thatthis family of solutions is two-parameter and construct diagrams of roll waves existence,11s well as compare the exact and numerical solutions. It should be noted that obtainingsimilar piecewise-smooth solutions for more general models (7) or (20) is a rather complicatedproblem which requires a separate consideration.The solutions to system (8) in the class of travelling waves are determined from theequations (cid:0) (¯ u − D ) h (cid:1) ′ = 0 , (cid:18) ( β ¯ u − D )¯ uh + αh (cid:19) ′ = h − (cid:18) ¯ uh (cid:19) n , (22)where the prime denotes differentiation with respect to the variable ζ , and D is the constantvelocity of the travelling wave. Consider the case of 0 < ¯ u < D . Then (¯ u − D ) h = − m ,where m is a positive constant which is the discharge rate as seen by an observer moving atthe wave speed D . Let us introduce the notation G ( h ) = (cid:0) ( β − Dh − βm (cid:1)(cid:18) D − mh (cid:19) + αh , ∆( h ) = dGdh = ( β − D − βm h + αh, F ( h ) = h − h n (cid:18) D − mh (cid:19) n . (23)Then system (22) reduces to the ordinary differential equation dhdζ = F ( h )∆( h ) . (24)A characteristic feature of roll waves is the presence of a critical depth y such that∆( y ) = 0. For h = y , a smooth transition occurs from a supercritical (∆( h ) <
0) to asubcritical (∆( h ) >
0) flow in the system of coordinates moving with the wave velocity D .To guarantee a smooth transition, one needs the condition F ( y ) = 0 to be satisfied. In viewof (23) we obtain the following representations for the parameters m and D in terms of thecritical depth y : m = Dy − y /n , D = βy /n + q β ( β − y /n ) + αy . (25)Note that the wave speed D coincides with the ‘frozen’ characteristic velocity D = λ + (¯ u ( h ) , h )for h = y .Mutual arrangement of the curves F ( h ) and ∆( h ) from the right hand side of Eq. (24)is shown in Fig. 6 (a). Such a qualitative arrangement of the curves takes place when theinequality α < nn y /n (26)is fulfilled. Indeed, ∆( h ) is a monotonically increasing function (∆ ′ ( h ) > h >
0) and F ( h ) has two real roots h = y and h = h a for h > m/D . In the interval between these roots,the function F ( h ) is negative. Obviously, h a < y (as shown in the figure) if F ′ ( y ) > h a > y if F ′ ( y ) <
0. The derivative of the function F ( h ) at the point h = y is F ′ ( y ) = 1 + n (cid:0) − Dy − − /n (cid:1) . In view of the second formula in (25), the inequality F ′ ( y ) > .7 0.8 0.9 1.0 1.1-0.2-0.10.1 z z= h D ( )/15 z=F h ( ) a h a C hG b A BC
D<0 D>0 hh yh a Figure 6: Mutual arrangement of the curves z = F ( h ) and z = ∆( h ) /
15 (a); graph of the function G ( h ) (b).The graphs are obtained for n = 0 . α = 1 and y = 1 . For given parameters α and n of the model, one can construct a two-parameter family ofperiodic piecewise-smooth solutions as follows. First, we arbitrarily choose the critical depth y > h fromthe interval ( h a , y ). These are points C and A in Fig. 6 (b). Instead of h one can choosethe maximum fluid depth h or the wave amplitude a = h − h . At h = y the ratio F/ ∆has the finite value because both the functions F and ∆ have the first order zeros at thispoint. For h > h the right-hand side of Eq. (24) is positive, so the function h ( ζ ) increasesuntil the maximum value of h is reached (point B in the graph). The Hugoniot conditionsat the shock front ( D − ¯ u ) h = ( D − ¯ u ) h = m, G ( h ) = G ( h ) . (27)follow from the balance laws (8). We denote the values of the functions on the right and leftsides of the discontinuity by adding subscripts 1 and 2, respectively. The second relationin (27) allows one to determine the maximum depth h > h . Thus, the periodic solutionconsists of a continuous interval ACB and a shock BA as shown in Fig. 6 (b). In this casethe shock stability condition λ + (¯ u , h ) < D < λ + (¯ u , h ) (28)is satisfied.A periodic piecewise-smooth solution of Eq. (8) constructed according to the abovealgorithm for n = 0 . α = 1, y = 1 .
01 and h = 0 .
84 is shown Fig. 4 (dash-dotted curve).Using formulae (25) and (27), we find the maximum depth h ≈ .
21, the wave velocity D ≈ .
30 and m ≈ .
29. It can be seen from Fig. 4 that the constructed exact solution almostcompletely coincides with the non-stationary calculation performed according to Eq. (8) forsinusoidal perturbations of a constant solution with frequency of Ω = 2 π .The construction of roll waves is impossible if F ′ ( y ) <
0. As mentioned above, in thiscase F ( h ) >
0, ∆( h ) < h ∈ ( m/D, y ) and F ( h ) <
0, ∆( h ) > h ∈ ( y, h a ).Therefore, the right-hand side in (24) is negative, i.e. the function h ( ζ ) is decreasing. Let ustake h from the interval ( y, h a ). Then we find h < h from the Hugoniot condition (27).Eq. (24) allows one to construct a continuous solution h = h ( ζ ) connecting the values of h .0 0.5 1.0 1.5 2.00.00.51.01.52.0 0.0 0.5 1.0 1.5 2.00.00.51.01.52.0 WW h * yh yh h * h * h * a b +- W - W + y * Figure 7: Roll wave diagrams for the parameters α = 1, n = 0 . n = 0 . F ( h, y ) = 0 ( h = y ), curve 2 is conjugate. Dashed line is the diagonal h = y . and h . However, it is impossible to pass from h = h to h = h using a shock wave, since thecondition (28) is violated. Thus, in the case F ′ ( y ) <
0, there is no periodic piecewise-smoothsolution (roll waves).The solution constructed above contains two parameters: the critical depth and thewave amplitude. Here we present a diagram of roll waves, which defines all possible pointson the ( y, h )-plane for which there exist roll waves. Let us take α = 1 (for other values α the qualitative behaviour of the diagrams does not change) and consider two values of n ,namely n = 0 . > n ∗ and n = 0 . < n ∗ , where n ∗ = 1 / √
2. This choice is explained bythe fact that [22] noted some differences in the construction of roll waves for a power-lawfluid with n > n ∗ and n < n ∗ . To construct the diagram, we use the functions G , F , and∆ defined by formulae (23) and (25). These functions depend on variable h and parameter y in such a way that ∆( y, y ) = 0 and F ( y, y ) = 0. We find numerically the dependence h = s ( y ) ( s = y ) at which the equation F ( s ( y ) , y ) = 0 is valid. This dependence is shownby curve 1 in Fig. 7. Curve 2 in the same figure is conjugate to curve 1 and is obtainedfrom the condition G ( s ( y ) , y ) = G ( h, y ) for h = s ( y ). The equation ∆( h, y ) = 0 has onlyone branch of solutions h = y and does not impose additional restrictions on the range ofadmissible parameters. In the region bounded by curves 1 and 2, the ratio F/ ∆ does notchange sign. However, for y < h ∗ in the region Ω − this ratio is negative and, as shown above,the construction of a piecewise-smooth periodic solution is impossible. The value h ∗ , definesthe point of intersection of the ‘equilibrium’ and ‘frozen’ characteristic velocities, is given byformula (12). In the domain Ω + , the ratio F/ ∆ is positive. Thus, Ω + defines all admissiblevalues of the parameters of roll waves on the ( y, h )-plane.For a fluid of flow index above n ∗ , the function s ( y ) (such that F ( s ( y ) , y ) = 0) increasesmonotonically and s ( y ) > h ∗ in the domain Ω + . For n < n ∗ , i.e. highly non-Newtonianfluid, the function s ( y ) is non-monotonic and s ( y ) < h ∗ on some interval y ∈ ( h ∗ , y ∗ ). Thismeans that it is possible to construct roll waves for which the minimum flow depth h = h
14s less than the stability threshold h = h ∗ . However, the critical depth h = y is higher thanthe specified threshold.
7. Conclusion
In the paper we perform a nonlinear stability analysis on the film flow of a power-lawfluid down an inclined plane by applying Whitham’s criterion to the governing equations. Tostudy the flow evolution, we use both 2D equations of the boundary layer theory (7) and 1Ddepth-averaged model (8). In dimensionless variables, these models contain two constants n and α , which are the power-law index and the flow parameter. It turn out that the velocitiesof the characteristics on a steady solution with a velocity profile of the Nusselt-type formodel (7) and the corresponding constant solution of equations (8) are somewhat different(see Fig. 2). Therefore, the stability criteria (13) and (16) for a steady flow of constantdepth for the models under consideration do not completely coincide. This circumstancemakes it possible to construct examples of flows for which small perturbations grow withinthe framework of Eq. (7) and decay when we use depth-averaged Eq. (8). However, thisrange of parameters in the ( n, α )-plane is rather narrow and, as a rule, both models givesimilar results of calculating the fluid depth and flow rate.The given numerical calculations confirm and illustrate the indicated above theoreticalresults. For the numerical simulation of the development of finite amplitude waves frominfinitesimal disturbances on the surface, we use standard methods developed for solvingsystems of hyperbolic balance laws. The application of this approach to Eq. (8) presentsno difficulties. To solve long-wave Eq. (7) by the same method, we derive 1D multilayerequations (20) that approximate the original 2D system. The calculation results of thedevelopment of roll waves are shown in Fig. 3 and 4. As we can see, the amplitude of rollwaves obtained by Eq. (20) is slightly higher than calculated by Eq. (8), while the wavelengthpractically coincides. The qualitatively different evolution of small perturbations for the 2Dand 1D models is presented in Fig. 5. We also consider the construction of exact piecewise-smooth solutions of equations (8) in the class of travelling waves. We focus on the fact thatthis class of solutions is two-parameter. In the plane of parameters, the critical depth – fluiddepth, we construct diagrams of the region of roll waves existence for the given values of n and α (see Fig. 7).The study shows that qualitative differences in the evolution of small perturbations onthe solutions of the 2D and 1D models are manifested only in a narrow range of n and α .Outside this range of parameters, the use of the depth-averaged equations is justified andexpedient. Acknowledgements
This work is supported by Russian Foundation for Basic Research (Project No. 19-01-00498). The author thanks S.L. Gavrilyuk, V.Yu. Liapidevskii and I.V. Stepanova forfruitful discussions. 15 eferenceseferences