The paper theoretically investigates the stability properties of the dam-break wave of a fluid with power-law rheology. Assuming the long-wave approximation, a depth-averaged flow model is considered. The linear stability analysis of the wave is carried out to individuate the marginal stability conditions. To this aim, the multiple-scale technique is applied with reference to the kinematic wave solution, which formally limits the validity of the theoretical achievements to relatively long time scales. Both shear-thinning and shear-thickening fluids are considered. Similarly to the case with uniform conditions, the analysis indicates that stable conditions can be associated with a marginal value of the Froude number. However, differently from the uniform conditions, the marginal Froude number is shown to be a function not only of the power-law index but also of the streamwise gradient of the base flow velocity and of the disturbance wavelength. The critical Froude number is found to be larger than the corresponding one in uniform conditions. Numerical solutions of the full model confirmed the outcomes of the linear stability analysis for both shear-thinning and shear-thickening fluids.

A fluid layer flowing on a sloping bed may exhibit free surface instabilities, which are often referred to as roll-waves. Their presence may cause overflows and intermittent forces when impacting against obstacles, increasing the risk of infrastructure damages.1 The presence of roll-waves has been observed in water but also in non-Newtonian fluids, such as mud, debris, and lava flows.2 The majority of literature works studied surface instabilities under parallel-flow conditions, even if they may develop also over non-parallel and transient base flows, such as in dam-break phenomenon, represented by a rapid release of a mass of fluid, which can be water or exhibit a non-Newtonian behavior.3 The analysis of the kinematics of dam-break still has a great interest.4 

This study is focused on the definition of the unstable conditions of a non-parallel time-varying base flow of a non-Newtonian fluid, represented by a dam-break wave generated by the rapid release of a mass of power-law fluid. The power-law model is suitable in describing shear-thickening (with power-law index n > 1) and shear-thinning (n < 1) behavior of different polymeric liquids, such as colloids and suspensions,5 but also of geophysical flows involving mud or lava.6,7

The formation of roll-waves in non-Newtonian fluids has been widely studied within the framework of the long-wave approximation and using depth-integrated models. This approach offers the advantage that the dynamics of the flowing layer can be analyzed without requiring the detailed knowledge of the internal structure of the flow.8,9 In this regard and considering the power-law rheology, the instability development has been early studied performing the linear stability analysis of a parallel (or uniform) flow over an inclined plane.10 Berezin et al.11 investigated the effect of the inertial terms on the instability of a thin film of a fluid obeying to a power-law constitutive relation of the Ostwald–de Waele type through a linear stability analysis and a non-linear numerical solution. Lin and Hwang12 approximated the governing equations in the long-wave regime and applied the method of perturbation with multiple scales to derive the non-linear stability conditions of a power-law film. The results showed that the stability conditions are strongly affected by the power-law exponent n with the system becoming more unstable as n decreases. Dandapat and Mukhopadhyay,13,14 using the method of integral relations, derived an evolution equation revealing the presence of both kinematic and dynamic wave processes, which may either act together or singularly dominate, depending on the value of the problem parameters. Miladinova et al.,15 using the Benney long-wave model, studied the non-linear evolution of falling uniform films of a power-law fluid flowing over an inclined plane. The results showed that a saturation of non-linear interactions occurs, resulting in a finite-amplitude permanent wave. Amaouch et al.16 developed a method that combines the lubrication theory and the weighted residual approach to describe the non-linear behavior a thin film of a power-law fluid both far from and close to the instability threshold in the presence of a small-to-moderate Reynolds number. Ng and Mei10 and Longo,17 using the simplified von Kármán's momentum integral approach, studied finite-amplitude permanent roll-waves in a shallow layer of shear-thinning and shear-thickening fluids, respectively. The conditions for the existence of a periodic discontinuous solution, constituted by smooth profiles with a monotonically increasing flow depth between periodic shocks, have been individuated for both fluids. Zayko and Eglit18 investigated the effect of oblique perturbations that propagate with an arbitrary angle to the velocity of the undisturbed fluid flow moving down an incline plane, showing that, for a power-law fluid and under certain conditions, oblique perturbations can lead to flow instability. Chesnokov19 performed a non-linear stability analysis on a film of a power-law fluid down an inclined plane using both two-dimensional governing equations and a depth-averaged hyperbolic simplification. The comparison of the results of the two models indicates that the amplitude of roll-waves obtained by the 2D equations is slightly larger than that for the 1D model and for certain flow parameters the small perturbations of the basic solution grow for the 2D equations and decay for the depth-averaged model. Other studies have addressed the stability conditions of power-law films including the effects of wind stress,20 bed porosity,21–25 electric and magnetic fields,26,27 thermocapillarity,28,29 and oscillation of the inclined plane.30 

All of the above cited works address instabilities developing over an unperturbed state represented by the uniform flow of a power-law fluid. The stability problem of the steady solutions of the flow equations different from the uniform one has been investigated for several gradually varying profiles31 and in the presence of a wavy impermeable32 and porous33 bottom wall.

Focusing the attention on transient flow conditions, roll-waves have been observed in clear-water flows such as in the dam-break experiments carried out by Logan and Iverson.34 It has been found that instabilities develop after some time, i.e., at some distance from the dam where the bulk of the flow reaches a quasi-uniform and quasi-steady state. Motivated also by the above-mentioned experimental evidence, Bohorquez35 carried out a linear stability analysis to predict the occurrence of surface instabilities superposed to a dam-break water wave. A multiple-scale technique36,37 was applied to the 1D shallow-water model considering a space and time varying base flow. The study assumes the kinematic approximation to be valid. The validity of such an assumption, although restricted to relatively long time scales, has been shown by the author himself and confirmed by several further studies.38–41 Bohorquez35 concludes that the non-parallel and time-dependent characteristics of the wave play a fundamental role in causing the growth/decay of surface disturbances, as the result of the competition between kinematic and dynamic waves. Numerical simulations carried out using the fully dynamic 1D model confirmed the theoretical achievements.

Similarly to the clear-water case, surface instabilities have been experimentally observed even in association with dam-break waves of non-Newtonian fluids, e.g., Cochard and Ancey42 and Iverson et al.43 However, to the authors' knowledge, their description in terms of an instability mechanism has never been attempted.

The present paper aims to investigate the occurrence of unstable conditions in dam-break waves of power-law fluid according to the wave competition dynamics. To this aim, the depth-averaged 1D model of Ng and Mei10 is considered. Similarly to Bohorquez,35 the multiple-scale technique has been applied assuming as base flow the wave predicted by the kinematic model, widely used also with non-Newtonian fluids.3,44,45 The applicability of the kinematic model is first investigated. Then, the linear stability analysis is performed, assuming as base flow the space/time varying solution of the kinematic model. To verify the theoretical achievements, the predictions of the linear stability analysis are finally compared with the outcomes of numerical simulations of the full model. The study furnishes different original outcomes. First, the multi-scale technique is applied for the first time to a non-Newtonian fluid to investigate the development of hydrodynamic instabilities in dam-break waves, which constitutes a non-modal stability problem of great interest since the non-parallel time-varying characteristics of the base flow may affect the stability conditions. Moreover, the results permit to pinpoint the intrinsic differences between water and non-Newtonian power-law rheology behavior of free surface small disturbances. Finally, it furnishes indications on the development of roll-waves in geophysical flows such as those involving mud or lava, which is an aspect of practical and technical interest for the assessment of the correlated risk and for the design of appropriate countermeasures.46–48 

The paper is organized as it follows. Section II presents the formulation of the problem, while in Sec. III the applicability of the kinematic approximation is discussed. In Sec. IV, the linear stability analysis is performed and the results are illustrated. In Sec. V, theoretical predictions of the linear stability analysis are compared with the numerical results of the fully non-linear model. Finally, conclusions are given in Sec. VI. The full expression of the key parameters required for the evaluation of the stability conditions of the flow is given in the  Appendix.

We consider a two-dimensional homogeneous layer of an incompressible power-law fluid flowing over a fixed bed, with a constant slope [θ, with tan θ O ( 1 )] respect to the horizontal plane, without lateral inflows or outflows. Neglecting the surface tension and considering a laminar gradually varied flow where spatial variations occur over length scales larger than the flow depth, the dimensional depth-averaged mass and momentum conservation equations are as follows:10,
h ̃ t ̃ + u ̃ h ̃ x ̃ = 0 ,
(1)
h ̃ u ̃ t ̃ + x ̃ ( β u ̃ 2 h ̃ ) + g cos θ x ̃ ( h ̃ 2 2 ) = g h ̃ sin θ τ ̃ b ρ ,
(2)
where t ̃ is the time, x ̃ is the coordinate along the bed, h ̃ is the flow depth measured along the y ̃ axis orthogonal to the bed, u ̃ is the depth-averaged velocity component along x ̃, and g and ρ are the gravity and the fluid density, respectively. The expressions of the momentum correction coefficient β and of the bottom shear stress τ ̃ b in laminar conditions read10 
β = 2 2 n + 1 3 n + 2 ,
(3)
τ ̃ b = μ n ( 2 n + 1 n u ̃ h ̃ ) n ,
(4)
where μn and n represent the consistency and the rheological index of the power-law fluid, respectively. For shear-thinning fluids the rheological index is smaller than one, while values larger than one represent shear-thickening fluids.
The wave arising from an idealized dam-break in which the fluid at rest occupies a triangular wedge is studied (see Fig. 1). Therefore, the spatial distributions of the flow depth and the velocity at t ̃ = 0 are as follows:
h ̃ ( x ̃ , 0 ) = x ̃ tan θ 0 x ̃ h ̃ 0 tan θ , h ̃ ( x ̃ , 0 ) = 0 x ̃ > h ̃ 0 tan θ , u ̃ ( x ̃ , 0 ) = 0 x ̃ .
(5)
with h ̃ 0 the flow depth at the dam. Indeed, the initial dimensional volume (per unit width) of fluid is as follows:
A ̃ = h ̃ 0 2 2 tan θ .
(6)
FIG. 1.

Sketch of the dam-break problem.

FIG. 1.

Sketch of the dam-break problem.

Close modal
Denoting with η ̃ 0 a reference length scale related to the initial flow depth ( 0 < η ̃ 0 < h ̃ 0), let us introduce the following dimensionless variables:
x = x ̃ η ̃ 0 , t = t ̃ g η ̃ 0 , h = h ̃ η ̃ 0 , u = u ̃ g η ̃ 0 .
(7)
The choice of the actual value of η ̃ 0 will be detailed in Sec. III. Taking into account for (4) and (7), the dimensionless form of (1) and (2) reads
h t + u h x + h u x = 0 ,
(8)
u t + ( 2 β 1 ) u u x + ( β 1 ) u 2 h h x + cos θ h x = sin θ B d u n h n + 1 ,
(9)
Bd being the basal drag coefficient:
B d = μ n ρ ( 2 n + 1 n ) n g n 2 η ̃ 0 n + 2 .
(10)
The initial conditions in dimensionless variables read
h ( x , 0 ) = x tan θ 0 x h 0 tan θ , h ( x , 0 ) = 0 x > h 0 tan θ , u ( x , 0 ) = 0 x .
(11)
Next, the first-order outer solution of the dam-break problem (11) governed by Eqs. (8) and (9) is sought for through the method of matched asymptotic expansion. This outer solution corresponds to the kinematic approximation of the governing equations (8) and (9).8 Let us denote with l ̃ b the length scale of the flood wave (base flow) in the x direction and with ε b = η ̃ 0 / l ̃ b the corresponding shallowness parameter. To describe the flow at sufficiently far distance from the dam, i.e., for ε b 1, it is convenient to introduce the following slower space and time variables:49 
x ̂ = ε b x , t ̂ = ε b t .
(12)
Taking into account for Eq. (12), the governing equations (8) and (9) in the slower variables read
h t ̂ + u h x ̂ + h u x ̂ = 0 ,
(13)
ε b [ u t ̂ + ( 2 β 1 ) u u x ̂ + ( β 1 ) u 2 h h x ̂ + cos θ h x ̂ ] = sin θ B d u n h n + 1 .
(14)
The kinematic solution h kin ( x ̂ , t ̂ ) and u kin ( x ̂ , t ̂ ) of (8) and (9)3,44,50 can be viewed as the solution at the zeroth order of ε b of the problem (13) and (14)
u kin t ̂ + 2 n + 1 n u kin u kin x ̂ = 0 ,
(15)
h kin = ( B d sin θ ) 1 n + 1 u kin n n + 1 .
(16)
Moreover, the initial condition (11) is schematized as a point-wise source at x ̂ = 0 with a volume per unit of width A * = ε b A.35 To be consistent with volume conservation, the shallowness parameter ε b is assumed to be51,
ε b = 1 A = η ̃ 0 2 A ̃ .
(17)
It therefore results A * = 1. In the present dimensionless variables, the solution of Eqs. (15) and (16) with a Dirac delta function at x ̂ = 0 as the initial condition reads52,
h kin ( x ̂ , t ̂ ) = ( B d sin θ ) 1 n + 1 ( n 2 n + 1 ) n n + 1 ( x ̂ t ̂ ) n n + 1 x ̂ x ̂ s ( t ̂ ) , h kin ( x ̂ , t ̂ ) = 0 x ̂ > x ̂ s ( t ̂ ) ,
(18)
u kin ( x ̂ , t ̂ ) = n 2 n + 1 x ̂ t ̂ x ̂ x ̂ s ( t ̂ ) , u kin ( x ̂ , t ̂ ) = 0 x ̂ > x ̂ s ( t ̂ ) ,
(19)
in which the distance x ̂ s ( t ̂ ) represents the position reached by the downstream front of the dam-break wave at time t ̂, given by52,
x ̂ s ( t ̂ ) = ( 2 n + 1 ) ( t ̂ n ) n 2 n + 1 ( B d sin θ ) 1 2 n + 1 ( 1 n + 1 ) n + 1 2 n + 1 .
(20)
The influence of the rheological parameters on the dam-break wave shape is deeply discussed by Ganguly et al.52 

The conditions for the applicability of the kinematic approximation are addressed in Sec. III B.

The system constituted by Eqs. (15) and (16) would represent an appropriate approximation of Eqs. (13) and (14) provided that
ε b | u kin t ̂ + ( 2 β 1 ) u kin u kin x ̂ + ( β 1 ) u kin 2 h h kin x ̂ + cos θ h kin x ̂ | | sin θ | .
(21)
Next, the applicability conditions of the kinematic model are discussed imposing the fulfillment of (21). To this aim, the following measure of the spatial velocity variation, i.e., of the non-parallelism of the flow, is introduced:
ϕ = ε b u kin sin θ u kin x ̂ .
(22)
Taking into account for definition (22), inequality (21) may be rewritten as follows:
ϕ u kin | [ 2 n + 1 n + ( 2 β 1 ) ] u kin + [ ( β 1 ) u kin 2 h kin + cos θ ] d h kin d u kin | 1 ,
(23)
where
d h kin d u kin = n n + 1 ( B d sin θ ) 1 n + 1 u kin 1 n + 1 .
(24)
Taking into account for Eqs. (3) and (24), inequality (23) reads
ϕ | n n + 1 cos θ ( sin θ ) 1 / n ( B d u kin n + 2 ) 1 n + 1 2 n + 1 n ( n + 1 ) | 1.
(25)
Denoting with
F = u cos θ h ,
(26)
the local Froude number, the validity condition (25) for the kinematic approximation becomes
ϕ ϕ * with ϕ * = n n + 1 F * 2 | ( F * F ) 2 1 |
(27)
in which
F * = n 2 n + 1
(28)
represents the value of the Froude number in the marginal (critical) stability condition of the uniform flow.10,31

Independently of the n value, Eq. (27) reveals that the kinematic approximation loses its validity for small values of the Froude number, i.e., ϕ * 0 when F 0. Conversely, the limiting value ϕ * diverges when the Froude number equals the F * value and reduces when F increases. Figure 2 shows, for three different values of n, namely, n = 0.5 , 1.0 , 1.5, the limiting value ϕ * as a function of the Froude number. Considering that the kinematic wave model may be applied whenever ϕ ϕ *, i.e., in the region below the curves, Fig. 2 indicates that the increase in the power-law index favors the applicability of the kinematic approximation, consistently with the theoretical findings of Hogg and Pritchard.53 Indeed, Hogg and Pritchard53 have shown that, for sufficiently high values of n, i.e., n > 0.26, a reduction of the power-law exponent implies a bottom drag increase, which induces a reduction of the dam-break wave front celerity and therefore of the length of the flood wave l ̃ b.

FIG. 2.

Dependence of ϕ * on F for three different values of n.

FIG. 2.

Dependence of ϕ * on F for three different values of n.

Close modal
Assuming the kinematic solution validity, Eqs. (19) and (22) show that, at each instant t ̂, the maximum of ϕ takes place at the downstream front, i.e., at x ̂ = x ̂ s. Let ϕ s ( t ̂ ) be the maximum value of ϕ at t = t ̂. The expression of ϕ s ( t ̂ ) may be deduced by substituting Eq. (19), along with its derivative with respect to x ̂, in Eq. (22), setting x ̂ = x ̂ s and accounting for Eq. (20). The following expression is indeed deduced for ϕ s ( t ̂ ):
ϕ s ( t ̂ ) = ε b 2 n + 1 [ 1 B d ( n t ̂ ) 3 n + 2 ( 1 n + 1 ) n + 1 ( 1 sin θ ) 2 n ] 1 2 n + 1 .
(29)
Equation (29) shows that ϕ s ( t ̂ ) is a decreasing function of the slower time. Starting from this observation the applicability condition (27) may be reformulated in terms of the slower time. In particular, the kinematic model is applicable for instants larger than t ̂ c, with t ̂ c such that
ϕ s ( t ̂ c ) ϕ * = Ψ with Ψ 1.
(30)
To find the t ̂ c value, the reference length scale η ̃ 0 is selected imposing that
h kin ( x ̂ s ( t ̂ c ) , t ̂ c ) = 1.
(31)
Substituting (18) into (31), and accounting for (20), the following relation is found:
t ̂ c = n n + 1 ( B d sin θ ) 1 n .
(32)
The reference length scale η ̃ 0 and the corresponding actual value of t ̂ c may then be evaluated by substituting (32) into (29) and finally enforcing the condition (30). By further accounting for the definition of Bd, Eq. (10), the following non-linear equation in η ̃ 0 is obtained:
( n + 1 ) 2 n ( 2 n + 1 ) 2 η ̃ 0 2 A ̃ sin θ a ( η ̃ 0 ) | a ( η ̃ 0 ) ( 2 n + 1 ) cos θ 1 | = Ψ
(33)
in which
a ( η ̃ 0 ) = [ ( μ n ρ ) 2 g n 2 η 0 ̃ n + 2 sin 2 θ ] 1 n .
(34)
Therefore, fixed the tolerance Ψ and assigned the values of the ( μ n , ρ , n) triplet, the bottom slope angle θ and the released fluid volume A ̃, the value of the reference length scale η ̃ 0 is evaluated by solving Eq. (33), and therefore, the values of Bd, ε b, and t ̂ c [see Eqs. (10), (17), and (32), respectively] are known. Finally, the value of the minimum time for the kinematic model applicability in fast coordinate follows from t c = t ̂ c / ε b.

In Sec. IV, the linear stability analysis of (8) and (9) will be carried out assuming as the base flow the solution of the kinematic model, i.e., Eqs. (18) and (19), along with (20).

This section presents the linear stability analysis of the dam-break wave of a power-law fluid, by applying the multi-scale technique This technique has been fruitfully employed to analyze the stability conditions of two-dimensional laminar boundary layer54–56 and recently of water dam-break wave in turbulent conditions.35 

We want to analyze the stability properties of a base flow slowly varying in the space and in the time, i.e., h b ( x ̂ , t ̂ ) and u b ( x ̂ , t ̂ ). In addition to the characteristic length scales η ̃ 0 and l ̃ b (with η ̃ 0 l ̃ b), the perturbation length l ̃ p has to be introduced as a third length scale. We further assume that not only the characteristic length scale of the base flow, l ̃ b, but also the perturbation length scale, l ̃ p, is much larger than η ̃ 0.57 

Following the multiple-scale technique,35 we expand the flow depth h ( x , t ) and the velocity u ( x , t ) with respect to the perturbation parameter ε p = l ̃ p / l ̃ b 1 in the form:
h ( x , t ; ε p ) = h b ( x ̂ , t ̂ ) + ε p h ( x ̂ , t ̂ , x , t ; ε p ) ,
(35)
u ( x , t ; ε p ) = u b ( x ̂ , t ̂ ) + ε p u ( x ̂ , t ̂ , x , t ; ε p ) .
(36)
As follows, we assume ε p = ε b 1, similarly to Saric and Nayfeh54 and Citro and Luchini57,58 in the derivation of the stability bounds for the steady and unsteady boundary layer, respectively, and to Bohorquez35 with reference to the water dam-break problem.
Moreover, the stability analysis is carried out for times t larger than tc. Consequently, we assume
h b ( x ̂ , t ̂ ) = h kin ( x ̂ , t ̂ ) , u b ( x ̂ , t ̂ ) = u kin ( x ̂ , t ̂ ) with t ̂ > t ̂ c ,
(37)
with h kin ( x ̂ , t ̂ ) and u kin ( x ̂ , t ̂ ) given by (18) and (19), respectively. Substituting (35) and (36) into (8) and (9), accounting for Eqs. (37), (15), and (16), at the first order of ε p, the following equation in terms of the perturbation u T = [ h , u ] is found:
u t + A u x + ( B + ϕ C ) u = f .
(38)
In deducing Eq. (38), only the two terms O ( ε p 2 ) involving the spatial derivative of the base flow have been retained, since, accounting for Eq. (22) and under the assumption ε p = ε b, their order becomes O ( ε p ). Keeping these two terms, the measure of the velocity spatial variation of the base flow, i.e., ϕ, explicitly appears in the perturbation Eq. (38).
The expressions of the A, B, and C matrices and of the source term f, all of which depend on x ̂ and t ̂ only, in (38) are as follows:
A = ( u kin h kin ( β 1 ) u kin 2 h kin + cos θ ( 2 β 1 ) u kin ) ,
(39)
B = ( 0 0 ( n + 1 ) B d u kin n h kin n + 2 B d n u kin n 1 h kin n + 1 ) ,
(40)
C = sin θ u kin ( 1 d h kin d u kin 1 h kin [ ( 2 β 1 2 n + 1 n ) u kin + d h kin d u kin cos θ ] ( 2 β 1 ) + 2 ( β 1 ) u kin h kin d h kin d u kin ) ,
(41)
f = u kin x ̂ ( 0 ( 2 β 1 2 n + 1 n ) u kin + [ u kin 2 h kin ( β 1 ) + cos θ ] d h kin d u kin ) .
(42)
Owing to the independence of A, B, C, and f on the (x, t) pair, the Laplace transform is applied to Eq. (38) for x > 0, so that the following o.d.e. is deduced:
d L { u } ( s , t ) d t + D L { u } ( s , t ) = f s ,
(43)
where
L { u } ( s , t ) = 0 u ( x , t ) e s x d x
(44)
is the Laplace transform of u in x [ 0 , + ) and
D = A s + B + ϕ C .
(45)
Let us suppose that the −D matrix can be diagonalized and let us denote with Λ and V the corresponding eigenvalues and eigenvectors matrices, respectively, i.e., D = V Λ V 1. Under these assumptions, it is easy to verify that Eq. (43) admits the following solution:
L { u } ( s , t ) = V e Λ t V 1 L { u } 0 + 1 s V Λ 1 ( e Λ t I ) V 1 f ,
(46)
with L { u } 0 the Laplace transform of u at t = 0 and I the identity matrix.
The stability of the base flow, i.e., the condition in which a time asymptotic decay of the perturbation occurs, is governed by the real part of the eigenvalues of the diagonal matrix Λ.35,59 Next, with the aim of performing the temporal stability analysis of the flow model (38), the complex number s = s r + i s i is assumed to have null real part, i.e., sr = 0, and an imaginary part s i = 2 π / Λ 0, with Λ 0 the disturbance wavelength. Moreover, it is convenient to introduce the local wavenumber α defined as follows:
α = Λ 0 2 π sin θ u kin 2
(47)
and, under the validity of the kinematic approximation, to express the Froude number as
F = u kin n + 2 2 ( n + 1 ) 1 cos θ ( sin θ B d ) 1 n + 1 .
(48)
Taking into account for (48), it is easy to verify that the matrix −D admits two distinct eigenvalues λ ̂ ± = 2 n ( 3 n + 2 ) ( n + 1 ) u kin sin θ λ ±, with
λ ± = [ n 2 ( n + 1 ) ( 3 n + 2 ) + 2 n ( 5 n 2 + 6 n + 2 ) ϕ + 4 i α n ( n + 1 ) ( 2 n + 1 ) ] ± { n [ 4 n ( n 2 ( 3 n + 2 ) 2 F 2 + n 4 16 n 3 32 n 2 20 n 4 ) ϕ 2 ± 4 n 2 ( n + 1 ) ( 3 n + 2 ) ( n 2 + 4 n + 2 ) ϕ + 4 α 2 ( 3 n + 2 ) 2 F 2 + 8 n ( 2 n + 1 ) α 2 + 4 ( 3 n + 2 ) ( 2 n + 1 ) ( n + 2 ) α n 2 ( 3 n + 2 ) 2 + i ( n + 1 ) [ 8 α n 2 ( 3 n + 2 ) 2 F 2 + 4 α ( 2 n + 1 ) ( 2 n 3 9 n 2 12 n 4 ) ] ϕ 4 i α n ( 3 n + 2 ) ( 2 n + 1 ) ( n + 2 ) ( n + 1 ) 2 ] } 1 2 .
(49)
Unstable conditions, i.e., the temporal growth of the perturbation, occur whenever ( λ ) > 0. The real part of λ ± reads
( λ ± ) = 2 n ( 5 n 2 + 6 n + 2 ) ϕ n 2 ( n + 1 ) ( 3 n + 2 ) ± a 2 + b 2 + a 2 ,
(50)
with
a = 4 n 2 F 2 ( 3 n + 2 ) 2 [ n 2 ϕ 2 α 2 ( 1 + n ) 2 ] + 4 n 2 ( n 4 16 n 3 32 n 2 20 n 4 ) ϕ 2 4 n 3 ( n + 1 ) ( 3 n + 2 ) ( n 2 + 4 n + 2 ) ϕ n 3 ( n + 1 ) 2 [ 8 ( 2 n + 1 ) α 2 n ( 3 n + 2 ) 2 ] ,
(51)
b = 4 α n ( n + 1 ) [ 2 n 2 ϕ F 2 ( 3 n + 2 ) 2 + ( 2 n + 1 ) ϕ n 2 ( 3 n + 2 ) ( 2 n + 1 ) ( n + 2 ) ( n + 1 ) ] .
(52)

Owing to the negative definiteness of ( λ ) , λ ̂ does not play any role in defining the stability conditions; i.e., only a temporal decay of the perturbation is associated with the λ ̂ eigenvalue.

Conversely, ( λ + ) may be positive or negative. Therefore, either a growth or a decay of the perturbation may take place. Setting ( λ + ) = 0 in (50) and solving for F provides the expression of the threshold stability Froude number, Fm (marginal Froude number), as a function of the rheological exponent n, the parameter ϕ, and the disturbance wavenumber α. The lengthy expression of Fm is given in the  Appendix, along with a more concise approximation, accurate for moderate values of the ϕ parameter.

Figure 3 reports in the ( F , α) plane the contour lines of ( λ + ) for n = 0.5 (a), n = 1.0 (b), n = 1.5 and (c), assuming ϕ = 0. Independently of the n values, Fig. 3 shows that the neutral stability curves, i.e., the ( F , α) pairs such that ( λ + ) = 0, correspond to a vertical line in the ( F , α) plane. Therefore, the marginal stability condition of a parallel flow is independent of the disturbance wavenumber and the critical Froude number Fc, defined as the minimum value of Fm, is a function of the power-law index only. Moreover, for F close to Fc the growth rate is so small that a remarkably long channel would be needed for permitting the development of unstable waves up to a perceivable threshold. Finally, Fig. 3 suggests that Fc coincides with F * (red lines), in agreement with the results of normal mode10 and near-front expansion31 analyses.

FIG. 3.

Contour lines of the constant growth/decay rate ( λ + ) in the ( F , α) plane for ϕ = 0 and three different n values. Solid line: ( λ + ) = 0 , 1 , 2 , 3 , 4 , 5 , 6. Dashed line: ( λ + ) = 1. Red lines: marginal stability of normal mode analysis F c = F *.

FIG. 3.

Contour lines of the constant growth/decay rate ( λ + ) in the ( F , α) plane for ϕ = 0 and three different n values. Solid line: ( λ + ) = 0 , 1 , 2 , 3 , 4 , 5 , 6. Dashed line: ( λ + ) = 1. Red lines: marginal stability of normal mode analysis F c = F *.

Close modal

Figure 4, which is the counterpart of Fig. 3, refers to the ϕ = 10 3 case. Figure 4 puts in evidence that the presence of a velocity gradient strongly modifies the linear stability conditions. Independently of the n value, the neutral stability condition in the ( F , α) plane does not show up as a straight vertical line. Indeed, Fm is found to be a function, not only of the n value as in the case of parallel flow, but also of the disturbance wave number α. At very high values of α, Fm tends to a critical value Fc, which represents the lower bound in terms of the Froude number for the instability. A reduction of α leads to an increase in Fm. Moreover, there exists a minimum value of α below which the base flow is always stable, for all Froude number values (horizontal asymptote, α c , ). A similar result has been found by Bohorquez35 in analyzing the stability conditions of a dam-break water wave in turbulent regime. Finally, Fig. 4 shows that the increase in n induces an increase in the α c , value and then an enlargement of the stability region. Similarly to the uniform flow conditions, the stabilizing effect due to the increase in the power-law exponent may be mainly ascribed to the increase in the effective viscosity of the film flows with n.12 

FIG. 4.

Contour lines of the constant growth/decay rate ( λ + ) in the ( F , α) plane for ϕ = 10 3 and three different n values. Solid line: ( λ + ) = 0 , 1 , 2 , 3 , 4 , 5 , 6. Dashed line: ( λ + ) = 1. Red lines: marginal stability of normal mode analysis F c = F * ( ϕ = 0).

FIG. 4.

Contour lines of the constant growth/decay rate ( λ + ) in the ( F , α) plane for ϕ = 10 3 and three different n values. Solid line: ( λ + ) = 0 , 1 , 2 , 3 , 4 , 5 , 6. Dashed line: ( λ + ) = 1. Red lines: marginal stability of normal mode analysis F c = F * ( ϕ = 0).

Close modal

With the aim of deeply analyzing the marginal stability conditions in the presence of a spatial velocity variation, Fig. 5 reports the marginal Froude number as function of α, for different values of ϕ and n.

FIG. 5.

Marginal stability curves for different values of ϕ and three different n values.

FIG. 5.

Marginal stability curves for different values of ϕ and three different n values.

Close modal

Independently of n, Fig. 5 shows that two branches exist. The former refers to very low values of the marginal Froude number and it is particularly evident at high ϕ values, i.e., ϕ > 10 2. However, this branch does not belong to the region in which the kinematic model is applicable [see Eq. (27) and Fig. 2], and therefore, its interest is marginal. Conversely, the second branch develops always for F > F *, and therefore, it refers to conditions in which the kinematic model can confidently approximate the full model.

Independently of the n values, Fig. 5 confirms the dependence of the marginal condition on ϕ and the existence of the vertical (Fc) and horizontal ( α c , ) asymptotes, both depending on the ( ϕ , n) pair. The analytical expressions of both bounds Fc and α c , , reported in the  Appendix, have been found starting from Eq. (50), accounting for (51) and (52), setting ( λ + ) = 0 and taking the limit for α + and F + , respectively.

Figure 6 highlights, for three different n values, the dependence of both Fc [Fig. 6(a)] and α c , [Fig. 6(b)] on ϕ. Figure 6(a) shows that Fc is always larger than F * and it tends to F * as ϕ 0. Therefore, the unsteady non-parallel nature of the base flow leads to increase the stability region expressed in terms of the Froude number.

FIG. 6.

F c / F * (a) and α c , (b) as function of ϕ for three different n values.

FIG. 6.

F c / F * (a) and α c , (b) as function of ϕ for three different n values.

Close modal

Moreover, also the minimum value of the wave number below which the flow is stable for all Froude number values increases with ϕ starting from 0 [Fig. 6(b)]. It is also evident that the F c / F * ratio decreases as the rheological index n increases, while the opposite is observed for α c , .

Present results indicate that similarly to the turbulent clear-water flow,35 also for the power-law model, the temporal evolution of small disturbances superimposed on a kinematic wave strongly differs from the corresponding one in the presence of a plane-parallel flow. Indeed, there exists a minimum value of the disturbance wavenumber below which the disturbance does not grow; i.e., a cutoff wavenumber for the spectrum of the unstable perturbations occurs.

On the other hand, due to the considered rheological model, the kinematic approximation of the dam-break wave exhibits at a given time, a variable Froude number along the wave. This marks a difference with the turbulent clear-water flows, characterized by kinematic wave with a constant value of the Froude number.35 In fact, by substituting Eq. (19) in Eq. (48), it is easily verified that the Froude number monotonically increases along the flow direction, attaining its maximum value, Fs, at the wave front, i.e.,: for x ̂ = x ̂ s, which reads
F s = ( n n + 1 1 t ̂ ) n + 2 2 n + 1 ( sin θ B d ) 3 2 n + 1 1 cos θ .
(53)

Moreover, as stated before also ϕ increases in the flow direction until the maximum value ϕ s expressed by Eq. (29) is reached at the downstream front. Since the increase in F and ϕ has opposite effect in terms of stability, in the present case, the detection of unstable conditions requires, in general, the analysis of the local fulfillment of the stability criterion along the dam-break wave.

Under unstable conditions, it is expected that the instability predicted by the above analysis is observed even in the solution of the fully non-linear governing equations, provided that a suitable disturbance is applied to the unperturbed flow. Section V is devoted to verify this issue. The spatial/temporal evolution of small perturbations superimposed on a flood wave is analyzed by numerically solving Eqs. (8) and (9).

As follows, the results of the stability analysis developed in the previous section are compared with the outcomes of numerical solution of the fully non-linear model. Indeed, the triggering of instability as the consequence of the suitable choice of the disturbance under linearly stable and unstable conditions is verified. To this aim, the propagation on an inclined plane with a constant slope θ of the wave arising from the sudden release of a mass of power-law fluid initially confined in a triangular wedge (Fig. 1) is simulated. The initial conditions in dimensionless variables are given by Eq. (11). The non-linear hyperbolic system (8) and (9) has been numerically solved with a finite volume scheme, which is second order in both time and space.31 The solution of the associated Riemann problem, i.e., the evaluation of the convective fluxes at the interface between adjacent finite volumes, is performed with the HLL approximation, whereas the temporal integration is performed with a Runge–Kutta scheme.60,61 Local absorbent boundary conditions have been imposed at the upstream and downstream limits of the channel.62 Additional details about the numerical method can be found in Campomaggiore et al.31 

Two power-law fluid rheologies are considered. The former describes a lava (herein denoted as LF) with a shear-thinning behavior characterized by the following parameters: n = 0.763, ρ = 2780 kg / m 3 , μ n = 51.88 Pa s n.63 The latter refers to a shear-thickening mud ( denoted as M F ) with the following rheology: n = 1.53, ρ = 1200 kg / m 3 , μ n = 0.13 Pa s n.64 In both fluids, the initial dimensional volume (per unit width) A ̃ is assumed equal to 10 m 2. The bottom inclination angles are as follows: θ = 10 ° (LF) and θ = 3 ° (MF). Therefore, the initial depths are h ̃ 0 = 1.88 m and h ̃ 0 = 1.02 m, for the LF and MF, respectively.

The stability analysis developed in Sec. IV assumes that the base profile may be approximated (with a tolerance Ψ) by the solution of the kinematic model when the critical time tc is overwhelmed. Indeed, preliminary tests have been carried out in order to numerically verify this conclusion, along with the theoretical estimate of tc given in Sec. III A.

As first step, for both fluids, the tolerance has been preliminary fixed using the following Ψ values: 0.1 and 0.01, for the LF and MF tests, respectively. The solution of Eq. (33) therefore gives the value of the reference length scale η ̃ 0 as η ̃ 0 = 0.15 m (LF) and η ̃ 0 = 0.06 m (MF). For both tests, Table I reports the computed values of the Basal drag coefficient Bd, of the shallowness parameter ε b and of the minimum applicability time tc, deduced through Eqs. (10), (17), and (32), respectively.

TABLE I.

Values of the governing parameters for the lava and mud fluids.

Fluid Bd ε b tc
LF  1.56 × 10 1  2.24 × 10 3  160 
MF  4.23 × 10 2  3.42 × 10 4  1600 
Fluid Bd ε b tc
LF  1.56 × 10 1  2.24 × 10 3  160 
MF  4.23 × 10 2  3.42 × 10 4  1600 

All the simulations have been carried out keeping fixed the following values of both mesh spacing and time step: Δ x = 0.25 , Δ t = 5 × 10 4. The results have been verified to be mesh-independent upon repetitions of the simulations with smaller values of both mesh spacing and time step.

Figure 7 compares the velocity profiles predicted by the fully dynamic (8)–(11) and the kinematic (18)–(20) models at different instants for both fluids. As theoretically expected, the agreement between the two solutions strongly improves as the time proceeds. Indeed, at t = tc, the deviation between the slopes of the two velocity profiles is lower than 5%, in agreement with the results of Bohorquez.35 Moreover, it has been verified that the difference on the peak value of the velocity is lower than 3% and 4% and the difference on the front position is lower than 2% and 3% for LF and MF, respectively. Therefore, present results confirm the validity the theoretical estimate of Sec. III A. Indeed, it may be concluded that in both tests the kinematic model can be confidently applied for t t c.

FIG. 7.

Velocity profiles at several instants. Fully dynamic model: solid lines; kinematic model: dashed lines for lava (a) and mud (b) tests.

FIG. 7.

Velocity profiles at several instants. Fully dynamic model: solid lines; kinematic model: dashed lines for lava (a) and mud (b) tests.

Close modal
With the aim of verifying the theoretical achievements of Sec. IV and accounting for the modest differences at t = tc between the two models, simulations have been carried out assuming as initial condition the solution of the kinematic model at t = tc,
h ( x , 0 ) = h kin ( x , t c ) , u ( x , 0 ) = u kin ( x , t c ) x ,
(54)
where h kin ( x , t c ) and u kin ( x , t c ) are given by (18) and (19) (with t = tc), respectively. The values of x s ( t c ), given by Eq. (20), are as follows: x s ( t c ) = 630 (LF) and x s ( t c ) = 4700 (MF).

For a given fluid, the linear stability analysis suggests that, at given Froude number F and the non-parallelism parameter ϕ values, a small disturbance superimposed to the base flow may trigger or not the instability depending on the local wavenumber α and therefore on the disturbance wavelength Λ 0.

In order to confirm this result, a small perturbation in the region of the wave between xmin and xmax, with xmin > 0 and x max < x s ( t c ) has been superposed to the longitudinal velocity profile at t = tc. The imposed velocity disturbance is as follows:
u p ( x ) = U 0 sin ( 2 π x x min Λ 0 ) , for x min x x max ,
(55)
with U0 the disturbance amplitude, assumed as U 0 = 5 × 10 3 for both fluids.

At t = 0 and for all abscissas in the range x min x x max, the values of F and ϕ are known from Eqs. (48) and (22), respectively. Moreover, for a given value of Λ 0, also the α values are known [Eq. (47)]. Indeed, a proper choice of the wavelength Λ 0 allows to define perturbed initial conditions, i.e., the ( F , ϕ , α ) spatial distribution at t = 0, for which the linear theory predicts stable or unstable conditions.

For each fluid, two values of the perturbation wavelength have been considered, namely, Λ 0 , s and Λ 0 , u. The former has been chosen to simulate linearly stable flow conditions, while the latter refers to unstable ones. The extents of the perturbed region, i.e., the values of xmin and xmax, along with the wavelengths of disturbance Λ 0 , s and Λ 0 , u assumed in the numerical simulations are reported in Table II.

TABLE II.

Perturbation parameters.

Fluid xmin xmax Λ 0 , s Λ 0 , u
LF  100  500  400  10 
MF  2000  4000  1000  100 
Fluid xmin xmax Λ 0 , s Λ 0 , u
LF  100  500  400  10 
MF  2000  4000  1000  100 

Figure 8 represents the spatial distribution of ( λ + ) for the perturbed initial conditions obtained assuming the Λ 0 values given in Table II. As shown in Fig. 8, ( λ + ) corresponding to the perturbation Λ 0 , s is everywhere negative for both fluids, and therefore, the linear theory predicts stable conditions. Conversely, the spatial distribution of ( λ + ), evaluated with reference to the Λ 0 , u values, shows the presence of ( λ + ) positive values and therefore of unstable flow conditions. Thus, a perturbation growth (respectively decay) is expected in response to the disturbance characterized by wavelength Λ 0 , u (respectively Λ 0 , s).

FIG. 8.

Real part of the most unstable eigenvalue for the lava (a) and mud (b) tests.

FIG. 8.

Real part of the most unstable eigenvalue for the lava (a) and mud (b) tests.

Close modal

Figure 9 reports the flow depth [Fig. 9(a)] and the velocity profiles [Fig. 9(b)] for the test LF at three different times for Λ 0 , s, confirming that the imposed perturbation does not grow, nor in time neither in space. Based on the numerical results for Λ 0 , u, Fig. 10 depicts the flow depth [Fig. 10(a)] and the velocity profiles [Fig. 10(b)] at three different times. The results indicate that the instabilities become apparent at time t = 50, when small amplitude waves are observed on the free surface. Present results, therefore, confirm the predictions of the stability analysis for the lava fluid. Consistently with the linear analysis, the amplitude of the perturbation increases from t = 50 to t = 200. In the same time interval, the disturbances are observed to descend the dam-break wave up to reach the downstream front x = xs. Moreover, a non-linear wave dynamics is clearly observed between t = 50 and t = 200, with a progressive coarsening of the perturbations and their evolution toward the form of a shock-wave train. This dynamics is typical of roll-waves superposed to uniform conditions.65 

FIG. 9.

Flow depth (a) and velocity (b) profiles at different times for the test LF with Λ 0 , s.

FIG. 9.

Flow depth (a) and velocity (b) profiles at different times for the test LF with Λ 0 , s.

Close modal
FIG. 10.

Flow depth (a) and velocity (b) profiles at different times for the test LF with Λ 0 , u.

FIG. 10.

Flow depth (a) and velocity (b) profiles at different times for the test LF with Λ 0 , u.

Close modal

With reference to the numerical simulations of the mud fluid with Λ 0 , s, Fig. 11, which reports the flow depth [Fig. 11(a)] and the velocity profiles [Fig. 11(b)] at three different instants, indicates the absence of disturbance growth. Similarly, Fig. 12 shows the flow depth [Fig. 12(a)] and the velocity profiles [Fig. 12(b)] at three different instants with reference to Λ 0 , u. The growth of the disturbance on the free surface is confirmed, as also indicated from the velocity evolution [Fig. 12(b)]. The amplitude of the perturbation grows in time while descending the dam-break wave, as it is evident from the profiles at t = 1000, where the perturbation has reached the downstream front. Also, in this case, the fully non-linear analysis shows the progressive coarsening and steepening of the initially sinusoidal disturbance, peculiar of the roll-waves.65 

FIG. 11.

Flow depth (a) and velocity (b) profiles at different times for the test MF with Λ 0 , s.

FIG. 11.

Flow depth (a) and velocity (b) profiles at different times for the test MF with Λ 0 , s.

Close modal
FIG. 12.

Flow depth (a) and velocity (b) profiles at different times for the test MF with Λ 0 , u.

FIG. 12.

Flow depth (a) and velocity (b) profiles at different times for the test MF with Λ 0 , u.

Close modal

In conclusion for both shear-thinning and shear-thickening fluids, the predictions of the linear stability analysis are confirmed by the numerical simulations performed with the fully non-linear model.

The present paper has investigated the stability of dam-break waves of a non-Newtonian fluid propagating on an inclined plane. A depth-averaged shallow-water model with a power-law rheology has been adopted and both shear-thinning and shear-thickening fluids have been analyzed. The linear stability of non-uniform time-dependent conditions has been studied by applying the multiple-scale technique, widely used to investigate the stability conditions of two-dimensional laminar steady and unsteady boundary layer. The analysis has been performed assuming as initial condition the solution of the kinematic model, i.e., the outer solution of matched asymptotic expansion of the governing equations. Therefore, present results are formally restricted to relatively long time-scales, i.e., t > t c with tc the critical time, for which the kinematic approximation is valid within a prescribed tolerance. It has been found that the critical time tc strongly depends on the power-law index, increasing as n decreases. For both shear-thinning and shear-thickening fluids, the linear stability analysis has shown that the non-uniform time-dependent character stabilizes the kinematic wave and increases the marginally stable Froude number Fm required for instability. Moreover, the marginal Froude number Fm is not only a function of the power index n, as in uniform conditions. In fact, in the kinematic wave case, Fm depends also on the local wavenumber α, which involves the disturbance wavelength, and on the measure of the spatial velocity variation ϕ. The analysis concludes that there exists a lower bound of the local wavenumber α c , , which depends on both n and ϕ, below which the kinematic wave is stable for all Froude numbers. The α c , value increases as the spatial velocity variation and the rheological index increase. Furthermore, the critical Froude number Fc is larger than the value in uniform condition. The closed expressions of both Fc and α c , have been provided. Numerical simulations have been carried out solving the full model using a finite volume scheme, with second-order accuracy both in time and space. Both shear-thinning and shear-thickening fluids have been considered for analyzing the spatial/temporal evolution of small perturbations superposed to a kinematic flood wave in its applicability conditions. Numerical results have confirmed the theoretical achievements of the linear stability analysis.

Future research will be devoted to extend the range of validity of the assumed base flow, including higher-order terms in the solution.35 Furthermore, the definition of the convective/absolute character of the stability, along with the upset of finite-amplitude disturbances, will be worthy of exploration. As far as the former issue is concerned, the Briggs criterion could be applied,36 while for the latter the weakly non-linear stability analysis of the problem66,67 could be carried out. Moreover, a possible transient growth caused by the non-modal character of the spatial differential operators would deserve investigation, as well.68,69 Finally, additional efforts could be devoted to extend the presented study to a two-dimensional steep-slope flow.70 

The authors are indebted to Professor Patricio Bohorquez of Universidad de Jaén for fruitful discussions and suggestions. The role of the anonymous reviewers, which allowed to improve the quality of the paper, is also gratefully acknowledged.

The paper was supported by the SEND intra-university project, financed through the V:ALERE 2019 program of the University of Campania “Luigi Vanvitelli” (Grant No. B68D19001880005).

The authors have no conflicts to disclose.

Cristana Di Cristo: Data curation (equal); Validation (equal); Writing – original draft (equal). Michele Iervolino: Data curation (equal); Funding acquisition (equal); Validation (equal); Writing – original draft (equal). Andrea Vacca: Conceptualization (equal); Methodology (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

As follows, the expression of the marginal Froude number predicted by the stability analysis is provided. As already pointed out, this value may be considered as a function of the rheological exponent of the power-law fluid, n, of the non-parallelism parameter, ϕ, and of the spatial wavenumber of the perturbation α. A relatively concise expression may be worked out as follows:
F m ( n , ϕ , α ) = n 2 2 f 1 ( n , ϕ , α ) f 2 ( n , ϕ , α )
(A1)
in which
f 1 ( n , ϕ , α ) = ( 3 n + 2 ) 2 [ 10 ϕ 2 + 4 ( 3 n + 1 ) ϕ + n ( 3 n + 2 ) ( n + 1 ) ] { 4 n 4 ( 5 n 2 + 6 n + 2 ) ϕ 6 + 4 n 5 ( 3 n + 2 ) ( n + 1 ) ( 5 n 2 + 6 n + 2 ) ϕ 5 + [ 8 ( 27 n 4 + 76 n 3 + 87 n 2 + 44 n + 8 ) α 2 + n 4 ( 3 n + 2 ) 2 ] ϕ 4 + 8 n 3 ( 5 n + 4 ) ( 3 n + 2 ) ( n + 1 ) 4 α 2 ϕ 3 + [ 2 ( 33 n 4 + 100 n 3 + 122 n 2 + 64 n + 12 ) α 2 n 4 ( 3 n + 2 ) 2 ] ϕ 2 + 4 n ( 3 n + 2 ) ( 9 n 2 + 16 n + 6 ) ( n + 1 ) 5 α 4 ϕ + n 2 ( 3 n + 2 ) 2 ( n + 1 ) 6 α 4 } 1 2 8 n 2 ( 3 n + 2 ) ( 5 n 2 + 6 n + 2 ) 2 ϕ 4 8 n 3 ( 3 n + 2 ) 2 ( 5 n 2 + 6 n + 2 ) ( n + 1 ) ϕ 3 + [ 8 ( 3 n + 2 ) ( n + 1 ) ( 21 n 4 + 76 n 3 + 89 n 2 + 44 n + 8 ) α 2 2 ( 3 n + 2 ) 3 ( n + 1 ) 2 n 4 ] ϕ 2 + 8 n ( 7 n + 4 ) ( 3 n + 2 ) 2 ( n + 1 ) 4 α 2 ϕ + 2 n 2 ( 3 n + 2 ) 3 ( n + 1 ) 4 α 2
(A2)
and
f 2 ( n , ϕ , α ) = ( n + 1 ) 2 { 16 n 2 ( 2 n + 1 ) ( 5 n 2 + 6 n + 2 ) 2 ϕ 4 + 8 n 3 ( 3 n + 2 ) 2 ( 5 n 2 + 6 n + 2 ) ϕ 3 + [ ( 2 n + 1 ) ( 64 n 6 + 140 n 5 + 46 n 4 111 n 3 126 n 2 52 n 8 ) α 2 + 4 n 4 ( n + 1 ) ( 12 n 2 + 15 n + 5 ) ( 3 n + 2 ) 2 ] ϕ 2 + [ ( 2 n + 1 ) ( 8 n 4 19 n 2 16 n 4 ) α 2 + n 4 ( n + 1 ) ( 3 n + 2 ) 2 ] ϕ n 2 ( 2 n + 1 ) ( n + 1 ) 2 ( 3 n + 2 ) 3 α 2 } .
(A3)
As it is apparent from Eqs. (A2) and (A3), the expression of Fm involves up to the eighth power of the ϕ parameter. However, it is easy to verify that in both the presented test-cases, the values of ϕ are far below the unity [e.g., ϕ s ( t c ) = 1.2 × 10 2 and ϕ s ( t c ) = 4.1 × 10 3 for LF and MF, respectively]. Under such conditions, assuming the non-parallelism parameter ϕ 1 and developing in series Eqs. (A2) and (A3) around ϕ = 0, the following more handy form is obtained:
F m = F * + ( 4 n + 3 ) ( n + 2 ) 2 ( 2 n + 1 ) 2 α 2 + n 4 ( n + 1 ) ( 3 n + 2 ) ( 9 n 3 + 21 n 2 + 16 n + 4 ) ( 2 n + 1 ) 3 2 α 2 ϕ + O ( ϕ 2 ) .
(A4)
For the sake of completeness, we also provide the expression of the critical Froude number for the onset of instability, Fc,
F c = F * g 1 ( n , ϕ ) g 2 ( n , ϕ )
(A5)
with
g 1 ( n , ϕ ) = 2 ( 3 n + 2 ) { [ 16 ( 33 n 4 + 100 n 3 + 122 n 2 + 64 n + 12 ) ( 5 n 2 + 6 n + 2 ) 2 ϕ 4 + 32 n ( 39 n 4 + 117 n 3 + 133 n 2 + 66 n + 12 ) ( 3 n + 2 ) ( n + 1 ) ϕ 3 + 8 n 2 ( 7 n + 4 ) ( 17 n 2 + 23 n + 8 ) ( 3 n + 2 ) 2 ( n + 1 ) 3 ϕ 2 + 8 n 3 ( 7 n + 4 ) ( 3 n + 2 ) 3 ( n + 1 ) 4 ϕ + n 4 ( 3 n + 2 ) 4 ( n + 1 ) 4 ] 1 2 + 8 ( 21 n 4 + 76 n 3 + 89 n 2 + 44 n + 8 ) ϕ 2 + 8 n ( 7 n + 4 ) ( 3 n + 2 ) ( n + 1 ) 2 ϕ + 2 n 2 ( 3 n + 2 ) 2 ( n + 1 ) 2 } ,
(A6)
g 2 ( n , ϕ ) = ( 64 n 6 + 140 n 5 + 46 n 4 111 n 3 126 n 2 52 n + 8 ) ϕ 2 + 2 n ( 3 n + 2 ) ( n + 1 ) ( 8 n 4 19 n 2 16 n 4 ) ϕ n 2 ( 3 n + 2 ) 3 ( n + 1 ) 2 .
(A7)
The corresponding approximation for small ϕ values reads
F c = F * [ 1 + ( 4 n + 3 ) ( n + 2 ) 2 ( 2 n + 1 ) n ( 9 n 3 + 21 n 2 + 16 n + 4 ) ϕ ] + O ( ϕ 2 ) ,
(A8)
which clearly puts in evidence that the spatio-temporal variation of the base flow leads to stabilize the base flow.
Finally, the analytical expression of the minimum value of the wavenumber below which the flow is stable for any Froude number, α c , , reads as follows:
α c , = ϕ g 3 ( n , ϕ ) g 4 ( n , ϕ )
(A9)
in which
g 3 ( n , ϕ ) = 16 ( 2 n + 1 ) ( 5 n 2 + 6 n + 2 ) 2 ϕ 3 + 8 n ( 5 n 2 + 6 n + 2 ) ( 3 n + 2 ) 3 ϕ 2 + 4 n 2 ( n + 1 ) ( 12 n 2 + 15 n + 5 ) ( 3 n + 2 ) 2 ϕ + 2 n 3 ( n + 1 ) 2 ( 3 n + 2 ) 3 ,
(A10)
g 4 ( n , ϕ ) = ( 2 n + 1 ) ( 64 n 6 + 140 n 5 + 46 n 4 111 n 3 126 n 2 52 n 8 ) ϕ 2 + 2 n ( 3 n + 2 ) ( 1 + 2 n ) ( n + 1 ) ( 8 n 4 19 n 2 16 n 4 ) ϕ n 2 ( 2 n + 1 ) ( n + 1 ) 2 ( 3 n + 2 ) 3 .
(A11)
Expanding in series around ϕ = 0, the following approximation is found:
α c , = 2 ϕ 2 n + 1 n 3 2 + O ( ϕ 3 2 ) .
(A12)

Both Eqs. (A8) and (A12) may be employed to estimate the asymptotic boundaries for sufficiently small values of the ϕ parameter.

Figure 13 compares the exact (lines) and the approximated (lines plus circles) values of F c / F * [Fig. 13(a)] and α c , [Fig. 13(b)] for the same n values reported in Fig. 6. In the investigated range of n values, Fig. 13 shows that the provided approximations perform with reasonable accuracy (5% of the exact value) up to about ϕ = 0.05 [Eq. (A8)] and up to approximately ϕ = 0.0125 [Eq. (A12)].

FIG. 13.

Approximation of the critical Fc and the asymptotic α c , values for small values of ϕ. Exact expressions: lines; approximated expressions: lines plus circles.

FIG. 13.

Approximation of the critical Fc and the asymptotic α c , values for small values of ϕ. Exact expressions: lines; approximated expressions: lines plus circles.

Close modal
1.
B.
Yu
and
V.
Chu
, “
Impact force of roll waves against obstacles
,”
J. Fluid Mech.
969
,
A31
(
2023
).
2.
R. V.
Craster
and
O. K.
Matar
, “
Dynamics and stability of thin liquid films
,”
Rev. Mod. Phys.
81
,
1131
1198
(
2009
).
3.
C.
Ancey
,
N.
Andreini
, and
G.
Epely-Chauvin
, “
Viscoplastic dambreak waves: Review of simple computational approaches and comparison with experiments
,”
Adv. Water Resources
48
,
79
91
(
2012
).
4.
T.
Tan
,
Y.
Ma
,
J.
Zhang
,
X.
Niu
, and
K.-A.
Chang
, “
Experimental study on flow kinematics of dam-break induced surge impacting onto a vertical wall
,”
Phys. Fluids
35
,
025127
(
2023
).
5.
C.
Routa
,
P.
Balaji
, and
A. K.
Sahu
, “
Numerical investigation of effect of surface pattern and rotation on power-law fluid flow and heat transfer around a cylinder in laminar flow regime
,”
Phys. Fluids
35
,
073101
(
2023
).
6.
X.
Zhang
,
Y.
Bai
, and
C.
Ng
, “
Rheological properties of some marine muds dredged from China coasts
,” in
Proceedings of the 28th International Offshore and Polar Engineering Conference, Beijing, China
(
2010
).
7.
A.
Tallarico
,
M.
Dragoni
,
M.
Filippucci
,
A.
Piombo
,
S.
Santini
, and
A.
Valerio
, “
Cooling a channeled lava flow with non-Newtonian rheology: Crust formation and surface radiance
,”
Ann. Geophys.
54
,
1
(
2011
).
8.
X.
Huang
and
M. H.
Garcia
, “
A Herschel-Bulkley model for mud flow down a slope
,”
J. Fluid Mech.
374
,
305
333
(
1998
).
9.
E. D.
Fernandez-Nieto
,
P.
Noble
, and
J. P.
Vila
, “
Shallow water equations for non-Newtonian fluids
,”
J. Non-Newtonian Fluid Mech.
165
,
712
732
(
2010
).
10.
C. O.
Ng
and
C. C.
Mei
, “
Roll waves on a shallow layer of mud modelled as a power–law fluid
,”
J. Fluid Mech.
263
,
151
184
(
1994
).
11.
Y.
Berezin
,
K.
Hutter
, and
L.
Spodareva
, “
Stability analysis of gravity driven shear flows with free surface for power-law fluids
,”
Arch. Appl. Mech.
68
,
169
178
(
1998
).
12.
J. S.
Lin
and
C. C.
Hwang
, “
Finite amplitude long-wave instability of power-law liquid films
,”
Int. J. Non-Linear Mech.
35
,
769
777
(
2000
).
13.
B. S.
Dandapat
and
A.
Mukhopadhyay
, “
Waves on a film of power–law fluid flowing down an inclined plane at moderate Reynolds number
,”
Fluid Dyn. Res.
29
,
199
220
(
2001
).
14.
B.
Dandapat
and
A.
Mukhopadhyay
, “
Waves on the surface of a falling power-law fluid film
,”
Int. J. Non-Linear Mech.
38
,
21
38
(
2003
).
15.
S.
Miladinova
,
G.
Lebon
, and
E.
Toshev
, “
Thin–film flow of a power-law liquid falling down an inclined plat
,”
J. Non-Newtonian Fluid Mech.
122
,
69
78
(
2004
).
16.
M.
Amaouche
,
A.
Djema
, and
L.
Bourdache
, “
A modified Shkadov's model for thin film flow of a power law fluid over an inclined surface
,”
C. R. Mec.
337
,
48
52
(
2009
).
17.
S.
Longo
, “
Roll waves on a shallow layer of a dilatant fluid
,”
Eur. J. Mech. B
30
,
57
67
(
2011
).
18.
J.
Zayko
and
M.
Eglit
, “
Stability of downslope flows to two-dimensional perturbations
,”
Phys. Fluids
31
,
086601
(
2019
).
19.
A.
Chesnokov
, “
Formation and evolution of roll waves in a shallow free surface flow of a power-law fluid down an inclined plane
,”
Wave Motion
106
,
102799
(
2021
).
20.
J. P.
Pascal
and
J. D.
D'Alessio
, “
Instability of power-law fluid flows down an incline subjected to wind stress
,”
Appl. Math. Modell.
31
,
1229
1248
(
2007
).
21.
J. P.
Pascal
, “
Linear stability of fluid flow down a porous inclined plane
,”
J. Phys. D: Appl. Phys.
32
,
417
422
(
1999
).
22.
R.
Usha
,
S.
Millet
,
H. B.
Hadid
, and
F.
Rousset
, “
Shear–thinning film on a porous substrate: Stability analysis of a one–sided model
,”
Chem. Eng. Sci.
66
,
5614
5627
(
2011
).
23.
C.
Di Cristo
,
M.
Iervolino
, and
A.
Vacca
, “
Gravity-driven flow of a shear–thinning power–law fluid over a permeable plane
,”
Appl. Math. Sci.
7
,
1623
1641
(
2013
).
24.
M.
Iervolino
,
J.
Pascal
, and
A.
Vacca
, “
Instabilities of a power–law film over an inclined permeable plane: A two-sided model
,”
J. Non-Newtonian Fluid Mech.
259
,
111
124
(
2018
).
25.
S.
Chakraborty
,
T.
Sheu
, and
S.
Ghosh
, “
Dynamics and stability of a power-law film flowing down a slippery slope
,”
Phys. Fluids
31
,
013102
(
2019
).
26.
K.
Zakaria
,
M. A.
Sirwah
, and
S. A.
Alkharashi
, “
Non-linear analysis of creeping flow on the inclined permeable substrate plane subjected to an electric field
,”
Int. J. Non-Linear Mech.
47
,
577
598
(
2012
).
27.
Y.
Gamiel
,
W.
Zahra
, and
M.
El-Behairy
, “
Stability criteria of streaming conducting fluids through porous media under the influence of a uniform normal magnetic field
,”
Int. J. Math. Trends Technol.
41
,
147
154
(
2017
).
28.
J. P.
Pascal
,
J. D.
D'alessio
, and
M.
Hasan
, “
Instability of gravity-driven flow of a heated power-law fluid with temperature dependent consistency
,”
AIP Adv.
8
,
105215
(
2018
).
29.
M.
Iervolino
,
J. P.
Pascal
, and
A.
Vacca
, “
Thermocapillary instabilities of a shear–thinning fluid falling over a porous layer
,”
J. Non-Newtonian Fluid Mech.
270
,
36
50
(
2019
).
30.
E.
Mogilevskiy
and
R.
Vakhitova
, “
Falling film of power-law fluid on a high-frequency oscillating inclined plane
,”
J. Non-Newtonian Fluid Mech.
269
,
28
36
(
2019
).
31.
F.
Campomaggiore
,
C.
Di Cristo
,
M.
Iervolino
, and
A.
Vacca
, “
Development of roll-waves in power-law fluids with non-uniform initial conditions
,”
J. Hydraul. Res.
54
,
289
306
(
2016
).
32.
C.
Heining
and
N.
Aksel
, “
Effects of inertia and surface tension on a power-law fluid flowing down a wavy incline
,”
Int. J. Multiphase Flow
36
,
847
857
(
2010
).
33.
J.
Pascal
and
A.
Vacca
, “
Instabilities of a shear-thinning fluid falling over an undulating porous layer
,”
J. Non-Newtonian Fluid Mech.
298
,
104693
(
2021
).
34.
M.
Logan
and
R. M.
Iverson
, “
Video documentation of experiments at the USGS debris flow flume 1992–2006
,”
Technical report Open-File Report 2007-1315
(
US Geological Survey
,
2007
), http://pubs.usgs.gov/of/2007/1315
35.
P.
Bohorquez
, “
Competition between kinematic and dynamic waves in floods on steep slopes
,”
J. Fluid Mech.
645
,
375
409
(
2010
).
36.
P.
Huerre
and
P. A.
Monkewitz
, “
Local and global instabilities in spatially developing flows
,”
Annu. Rev. Fluid Mech.
22
,
473
537
(
1990
).
37.
P. J.
Schmid
and
D. S.
Henningson
,
Stability and Transition in Shear Flows
(
Springer
,
2001
).
38.
H.
Capart
, “
Analytical solutions for gradual dam breaching and downstream river flooding
,”
Water Resources Res.
49
,
1968
1987
(
2013
).
39.
F.
Stilmant
,
M.
Pirotton
,
P.
Archambeau
,
S.
Erpicum
, and
B.
Dewals
, “
Determination of dam releases to generate warning waves in a mountain stream: Performance of an analytical kinematic wave model
,”
J. Hydraul. Eng.
144
,
05017006
(
2018
).
40.
T.-Y. K.
Chen
and
H.
Capart
, “
Kinematic wave solutions for dam-break floods in non-uniform valley
,”
J. Hydrol.
582
,
124381
(
2020
).
41.
P.
Bohorquez
, “
Discussion of ‘Computing nonhydrostatic shallow-water flow over steep terrain’ by Roger P. Denlinger and Daniel R. H. O'Connel
,”
J. Hydraul. Eng.
137
,
140
141
(
2011
).
42.
S.
Cochard
and
C.
Ancey
, “
Experimental investigation of the spreading of viscoplastic fluids on inclined planes
,”
J. Non–Newtonian Fluid Mech.
158
,
73
84
(
2009
).
43.
R. M.
Iverson
,
M.
Logan
,
R. G.
LaHusen
, and
M.
Berti
, “
The perfect debris flow? aggregated results from 28 large–scale experiments
,”
J. Geophys. Res.
115
,
1
29
, https://doi.org/10.1029/2009JF001514 (
2010
).
44.
D.
Pritchard
,
B. R.
Duffy
, and
S.
Wilson
, “
Shallow flows of generalised Newtonian fluids on an inclined plane
,”
J. Eng. Math.
94
,
115
133
(
2015
).
45.
C.
Di Cristo
,
M.
Iervolino
,
T.
Moramarco
, and
A.
Vacca
, “
Applicability of kinematic and diffusive models for mud-flows: A steady state analysis
,”
J. Hydrol.
559
,
585
595
(
2018
).
46.
F.
Wang
,
X.
Chen
,
J.
Chen
, and
Y.
You
, “
Experimental study on a debris-flow drainage channel with different types of energy dissipation baffle
,”
Eng. Geol.
220
,
43
233
(
2017
).
47.
M.
Greco
,
C.
Di Cristo
,
M.
Iervolino
, and
A.
Vacca
, “
Numerical simulation of mud-flows impacting structures
,”
J. Mt. Sci.
16
,
364
382
(
2019
).
48.
C.
Di Cristo
,
O.
Fecarotta
,
M.
Iervolino
, and
A.
Vacca
, “
Impact dynamics of mud flows against rigid walls
,”
J. Hydrol.
612
,
128221
(
2022
).
49.
B.
Hunt
, “
Asymptotic solution for dam break on sloping channel
,”
J. Hydraul. Eng.
110
,
1058
1071
(
1984
).
50.
C.
Di Cristo
,
M.
Iervolino
,
T.
Moramarco
, and
A.
Vacca
, “
Applicability of Kinematic model for mud-flows: An unsteady analysis
,”
J. Hydrol.
577
,
123967
(
2019
).
51.
C.
Ancey
,
S.
Cochard
, and
N.
Andreini
, “
The dam-break problem for viscous fluids in the high-capillary-number limit
,”
J. Fluid Mech.
624
,
1
22
(
2009
).
52.
A.
Ganguly
,
M.
Reza
, and
A. S.
Gupta
, “
Thin-film flow of a power-law fluid down an inclined plane
,”
J. Fluids Eng.
134
(
4
),
044502
(
2012
).
53.
A.
Hogg
and
D.
Pritchard
, “
The effects of hydraulic resistance on dam-break and other shallow inertial flows
,”
J. Fluid Mech.
501
,
179
212
(
2004
).
54.
W. S.
Saric
and
A. H.
Nayfeh
, “
Nonparallel stability of boundary layer flows
,”
Phys. Fluids
18
,
945
950
(
1975
).
55.
K. S.
Yeo
,
B. C.
Khoo
, and
W. K.
Chong
, “
The linear stability of boundary-layer flow over compliant walls: Effects of boundary-layer growth
,”
J. Fluid Mech.
280
,
199
225
(
1994
).
56.
S.
Zuccher
and
P.
Luchini
, “
Boundary-layer receptivity to external disturbances using multiple scale
,”
Meccanica
49
,
441
467
(
2014
).
57.
V.
Citro
and
P.
Luchini
, “
Multiple-scale approximation of instabilities in unsteady boundary layers
,”
Eur. J. Mech. B
50
,
1
8
(
2015
).
58.
V.
Citro
and
P.
Luchini
, “
Unsteady boundary–layer transition prediction
,” in
Proceedings of XXI AIMETA Congress 2013
(Libreria Cortina,
2013
) pp.
1
9
.
59.
M.
Gaster
, “
A note on the relation between temporally–increasing and spatially–increasing disturbances in hydrodynamic stability
,”
J. Fluid Mech.
14
,
222
224
(
1962
).
60.
S.
Gottlieb
and
C.
Shu
, “
Total variation diminishing Runge–Kutta schemes
,”
Math. Comput.
67
,
73
85
(
1998
).
61.
A.
Harten
, “
High resolution schemes for hyperbolic conservation laws
,”
J. Comput. Phys.
49
,
357
393
(
1983
).
62.
R. R.
Paz
,
M. A.
Storti
, and
L.
Garelli
, “
Local absorbent boundary condition for non-linear hyperbolic problems with unknown Riemann invariants
,”
Comput. Fluids
40
,
52
67
(
2011
).
63.
H. C.
Weed
,
F. J.
Ryerson
, and
A. J.
Piwinskii
, “
Rheological properties of Molten Kilauea Iki basalt containing suspended crystals
,” in
Mineral Matter and Ash in Coal
(
ACS Publications
,
1986
), pp.
223
233
.
64.
S.
Longo
,
V.
Di Federico
,
R.
Archetti
,
L.
Chiapponi
,
V.
Ciriello
, and
M.
Ungarish
, “
On the axisymmetric spreading of non-Newtonian power-law gravity currents of time-dependent volume: An experimental and theoretical investigation focused on the inference of rheological parameters
,”
J. Non–Newtonian Fluid Mech.
201
,
69
79
(
2013
).
65.
N. J.
Balmforth
,
J. W.
M-Bush
, and
R. V.
Craster
, “
Roll waves on flowing cornstarch suspensions
,”
Phys. Lett. A
338
,
479
484
(
2005
).
66.
I. M. R.
Sadiq
and
R.
Usha
, “
Thin Newtonian film flow down a porous inclined plane: Stability analysis
,”
Phys. Fluids
20
,
022105
(
2008
).
67.
A.
Mukhopadhyay
and
A. K.
Gaonkar
, “
Effects of the variation of viscosity on the stability of thin liquid film flows along a uniformly heated substrate under heat flux boundary condition
,”
Phys. Fluids
35
,
052109
(
2023
).
68.
L.
Trefethen
,
A.
Trefethen
,
S.
Reddy
, and
T.
Driscoll
, “
Hydrodynamic stability without eigenvalues
,”
Science
261
,
578
584
(
1993
).
69.
A.
Samantaa
, “
Non–modal stability analysis in viscous fluid flows with slippery walls
,”
Phys. Fluids
32
,
064105
(
2020
).
70.
A.
Maranzoni
and
M.
Tomirotti
, “
New formulation of the two-dimensional steep-slope shallow water equations. Part II: Numerical modeling, validation, and application
,”
Adv. Water Resources
177
,
104403
(
2023
).