- Email: [email protected]
Instability of power-law fluid flow down a porous incline
J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Instability of power-law fluid flow down a porous incline J.P. Pascal Department of Mathematics, Ryerson University, Toronto, Ont., M5B 2K3, Canada Received 25 July 2005; received in revised form 14 November 2005; accepted 19 November 2005
Abstract In this paper we investigate the generation and structure of roll waves developing on the surface of a power-law fluid layer flowing down a porous incline. The unsteady equations of motion for the power-law fluid layer are depth integrated according to the von K´arm´an momentum integral method accounting for the variation of the velocity distribution with depth. The slip boundary condition at the interface between the fluid layer and the porous plane is based on the assumption that the flow through the porous medium is governed by the modified Darcy’s law, and that the characteristic length scale of the pore space is much smaller than the depth of the fluid layer above. An analytical theory of permanent roll waves is employed to determine under what flow conditions roll waves can exist and to calculate the wavelength, wave height, and speed of the roll waves. The nonlinear stability analysis is also carried out by numerically solving the time dependent governing equations and calculating the nonlinear evolution of infinitesimal disturbances imposed on the uniform and steady flow. Conclusions are drawn regarding the effect of the permeability of the porous medium and flow conditions on the development and characteristics of roll waves arising from the instability of the uniform flow. © 2005 Elsevier B.V. All rights reserved. Keywords: Flow down a porous incline; Power-law fluid; Linear and nonlinear stability analysis; Roll waves
1. Introduction It has been observed that the steady flow of a uniform fluid layer in an inclined open channel becomes unstable under certain conditions. The instability evolves into a flow of nonconstant depth consisting of progressive bores connected by sections of gradually varying flow. This wave system is referred to as roll waves and has significant implications in various flows found in the environment. Of particular interest are situations when the fluid exhibits non-Newtonian rheological properties, as for example molten lava flows and mud flows [1]. It has been reported [2] that for many such flows, at low shear rates the shear rate–shear stress relation is of a shear-thinning nature and is appropriately modelled by the power-law model which in simple shear is expressed as [3]: ∂u ∂u n−1 τ=µ , ∂z ∂z where µ is the consistency and n is the power-law index of the fluid. For shear-thinning fluids n is between 0 and 1, and the case with n = 1 corresponding to a Newtonian fluid with µ being the dynamic viscosity. E-mail address: [email protected]. 0377-0257/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2005.11.007
Due to the relevant impact of the instability of flow down an incline, it is important to theoretically determine the critical conditions for the onset of instability and to predict the evolution of the unstable flow. By performing a linear stability analysis on the governing equations, one determines the range of the flow parameters for which infinitesimal disturbances grow exponentially in time. However, as the disturbances grow the nonlinear interactions become significant and a nonlinear analysis must be conducted in order to predict the subsequent evolution. The shallow-water equations have been widely used as the equations governing the velocity and depth of the fluid layer. Most notable is the early work by Dressler [4] who analytically obtained roll waves solutions to the shallow-water equations augmented with an empirical term to model the turbulent frictional force. Recently, Zanuttigh and Lamberti [5] obtained roll waves solutions by numerically solving the de St. Venant equations and calculating the evolution of the flow until roll waves are attained. The shallow-water equations are based on the assumption that the longitudinal velocity is uniform with depth (flat velocity profile) and are thus only appropriate for high Reynolds number flows. For laminar flows, the governing equations must account for the variation of the velocity with depth due to viscous effects, and incorporate the rheological properties of the fluid. In order to include the variation of the velocity with depth, an accurate approximation can be obtained by assuming the
110
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
velocity profile to be the same as that of the uniform and steady flow. This approximation is used in connection with von K´arm´an’s momentum integral method which consists of depth integrating the equations established under the assumption of moderate Reynolds number flow and disturbances of long wavelength (small aspect ratio). This approach was taken by Ng and Mei [1] to investigate the instability and roll waves structure of a power-law fluid flow down an impermeable incline. They carry out an analytical investigation by seeking roll waves solutions characterized as periodic shocks connected by smoothly increasing depth profiles. They obtain necessary conditions for the existence of roll waves, and present a procedure for calculating the wavelength, speed, and wave height of the roll waves. In this paper we extend the problem of instability of powerlaw fluid flow down an incline by incorporating the effect of the permeability of the inclined plane. This extension should prove useful for modelling environmental flows which mostly occur over porous surfaces. Such an extension is now being explored in the literature in connection with gravity currents. Much attention is being focussed on incorporating the effect of bottom permeability in the study of the spread of dense fluids under layers of lighter fluid, with the goal of obtaining more realistic models for gravity-driven flows over such surfaces as soil, river beds and layers of gravel. A summary of these investigations can be found in a recent paper by Thomas et al. [6]. It should also be pointed out that the layer of fluid behind the head of a gravity current resulting from the spread of dense fluid down an incline is subject to instability. As was recently discovered by Cenedese et al. [7], this uniform flow can evolve into a wave regime exhibiting a roll waves pattern. In order to determine the governing equations for our problem we apply von K´arm´an’s momentum integral method to obtain a model for the flow down an incline with a depth dependent longitudinal velocity. At the interface between the fluid layer and the porous medium we employ the slip conditions experimentally verified by Beavers and Joseph [8]. We begin our investigation by first conducting a linear stability analysis, which yields an analytical formula for the relationship between the Reynolds number, the inclination of the plane, and the permeability parameter which describes the critical conditions for the exponential growth of infinitesimal disturbances. Then by extending the analytical theory devised by Ng and Mei we study the existence of permanent wave solutions and calculate the relevant quantities describing the roll waves flow. Finally, we also carry out a nonlinear stability analysis by numerically solving the nonlinear governing equations to determine the evolution of a perturbed uniform flow. This will allow us to determine which of the analytically obtained permanent roll wave solutions do in fact arise from the evolution of infinitesimal disturbances.
u = (u, w)T and p, respectively. We assume the characteristic longitudinal length of the flow, L to be much larger than the characteristic depth, H, with the depth varying slowly in time and in the longitudinal direction. The equations of motion for the layer are obtained from the Navier–Stokes equations by neglecting terms of O(2 ) with = H/L being the aspect ratio, and are expressed as [1]:
2. The governing equations
where κ is the permeability of the porous medium, χ is a dimensionless parameter dependent on the structure of the porous medium, and up = (up , wp )T is the Darcian mean filter velocity in the porous medium. The normal velocity does not vary through the boundary-layer so the condition at the interface is
Consider the two-dimensional laminar flow of a thin layer of a power-law fluid along the surface of a porous medium incline at an angle θ with respect to the horizontal. We define an (x, z) coordinate system with the x-axis along the bottom and the z-axis normal to it. We denote the velocity and the total pressure by
∂u ∂w + = 0, ∂x ∂z ρ
(1)
∂u ∂u ∂u +u +w ∂t ∂x ∂z
∂p ∂ ∂u ∂u n−1 = − + gρ sin θ + µ , ∂x ∂z ∂z ∂z 0=−
1 ∂p − g cos θ, ρ ∂z
(2)
(3)
where g is the acceleration due to gravity, and ρ, n and µ are the mass density, the power-law index, and the consistency of the power-law fluid, respectively. If we neglect terms of O(2 ), including those arising from the effect of surface tension, then the continuity of force condition at the surface reduces to p = 0,
∂u =0 ∂z
at z = h(x, t).
(4)
The kinematic condition at the surface is given by w=
∂h ∂h +u ∂t ∂x
at z = h(x, t).
(5)
We assume that the porous medium is saturated with fluid and that the flow through it is governed by the modified Darcy’s Law. We note that even if the free flow is a high Reynolds number flow, the application of the modified Darcy’s Law to model the flow through the porous medium is validated by the assumption that the geometry of the pore space results in a low Reynolds number flow. The condition at the fluid layer-porous medium interface in this case has been investigated by Beavers and Joseph [8]. They propose that a boundary-layer is formed at the top of the porous medium and that the tangential velocity rapidly changes across this region from the Darcian velocity to the velocity at the interface. Beavers and Joseph verify experimentally that the appropriate condition for the tangential velocity at the interface is ∂u χ = √ (u − up ) ∂z κ
w = wp
at z = 0.
at z = 0,
(6)
(7)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
The conditions (6) and (7) have been used previously [9] in a model for Newtonian fluid flow down a porous incline. In the case of a general power-law fluid, condition (6) is expressed as ∂u χ = 1/(n+1) (u − up ) at z = 0, (8) ∂z κ where κ is the modified permeability of the porous medium with dimension of length to power n + 1. This extension of the Beavers–Joseph boundary condition has been recently used by Rao and Mishra [10] to study peristaltic transport of a power-law fluid in a porous tube. According to the modified Darcy’s Law which includes the nonlinear rheological properties of non-Newtonian power-law fluids, the flow in the porous medium is governed by [11]: 1/n ∂pp ∂pp 1/n−1 κ up = − − ρg sin θ , (9) µ ∂x ∂x ∂up ∂wp + = 0, ∂x ∂z ∂pp = −ρg cos θ, ∂z pp = p at z = 0,
(10) (11) (12)
where pp denotes the total pressure in the porous medium. We will consider cases when the pore space geometry is such that the flow is slower than that of the layer above to the extent that the mean filter velocity terms in the interface conditions (8) and (7) are negligible. We will now determine the order of magnitude of the permeability of the porous medium that is required for these terms to be of O(2 ). Based on the hydrostatic nature of the pressure distribution revealed by Eqs. (3) and (11) and the continuity of pressure ∂p condition (12) we conclude that we can use ∂xp = ∂p ∂x in Eq. (9). But the characteristic pressure in the fluid layer above the porous medium is n L U , P ∼µ H H where U denotes the characteristic longitudinal velocity. From Eq. (9) we determine the characteristic longitudinal velocity in the porous medium to be κ1/n U . H 1+1/n If we scale the variables in the boundary conditions at the fluid layer-porous medium interface, (8) and (7) with their characteristic values we obtain ∂u∗ χ = u∗ − χδ1/n u∗p , ∂z∗ δ and Up ∼
w∗ = δ2+1/n w∗p , 1/(n+1)
where δ = κ H is the ratio of the characteristic length scale of the pore space to the characteristic depth in the fluid layer flowing above the porous medium. If we assume that δ ≤ O(2n /χn ), then the up term in (8) is of O(2 ), and if δ ≤ O(2n/(2n+1) ) then
111
the wp term in (7) is of O(2 ). We also mention that we must also assume that δ ≤ O(3n/(n+2) ) for the Reynolds number of the porous medium flow to be of O(2 ) when the Reynolds number of the flow above is of O(1/). If we then neglect terms of O(2 ) the conditions at the bottom of the fluid layer reduce to ∂u χ (13) = 1/(n+1) u, w = 0, z = 0. ∂z κ We can now investigate the flow over the porous plane with the effect of the permeability of the surface incorporated through the slip condition in (13) and the parameters κ and χ. Integrating Eq. (3) and using the boundary condition (4) we find the expression for the total pressure to be: p = g cos θρ(h − z). Inserting these expressions into the x-momentum equation (2) we obtain ∂u ∂u ∂u +u +w ∂t ∂x ∂z ∂h µ ∂ ∂u n−1 ∂u = g sin θ − g cos θ + . (14) ∂x ρ ∂z ∂z ∂z 2.1. Depth-averaged equations A uniform and steady flow of the fluid layer is given by h = H = const. n(2n + 1)u¯ 0 (n + 1)[n + γ(2n + 1)] z 1+1/n n + 1 × 1− 1− + γ , H n
u = u0 (z) =
(15)
where γ=
κ1/(n+1) , χH
(16)
and u¯ is the depth-averaged velocity given by 1 H u¯ 0 = u0 (z) dz. H 0 Since we are dealing with a small aspect ratio flow which varies slowly in the longitudinal direction, we can apply von K´arm´an’s momentum integral method. We thus assume that the relationship between the velocity and the fluid depth for the steady flow is also true when these quantities are time dependent and nonuniform. Depth-integrating the continuity equation (1) and the momentum equation (14) over the fluid layer we obtain ∂h ∂(hu) ¯ + = 0, ∂t ∂x ∂h µ ∂ ∂(hu) ¯ + βhu¯ 2 = g sin θh − g cos θh − ∂t ∂x ∂x ρ n n u¯ 2n + 1 , × n + γ(2n + 1) h
(17)
(18)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
112
where β=
2 2n + 1 2nγ 2n2 2 γ + + . n + γ(2n + 1) 2n + 1 (3n + 2)(2n + 1) (19)
We now nondimensionalize these equations by introducing the nondimensional quantities (signified by an asterisk) according to L ∗ x = Lx∗ , u¯ = u¯ 0 u¯ ∗ , t= t , h = Hh∗ , u¯ 0 where u¯ 0 and H are the velocity and thickness of the uniform steady flow, respectively, and L=
u¯ 20 . g sin θ
Employing this transformation to nondimensionalize equations (17) and (18) and dropping the asterisks for notational convenience we obtain ∂h ∂(hu) ¯ + = 0, (20) ∂t ∂x n u¯ 1 ∂ ∂(hu) ¯ , (21) + βhu¯ 2 + αh2 = h − ∂t ∂x 2 h where
2 2n + 1 n−2 2n + 1 cot θ , α= n Re n + γ(2n + 1) 2/n 2−n n ρ 1+2/n 2/n−1 Re = H (g sin θ) , µ 2n + 1 and the expression for β is given by (19). In the case when γ = 0 ρu¯ 2−n H n
the Reynolds number can be expressed as Re = 0 µ . Eqs. (20) and (21) govern the depth and the depth-averaged velocity of the fluid layer. As parameters in these equations we have the power-law index, 0 < n ≤ 1, the uniform flow parameter, α ≥ 0, and the scaled permeability, γ ≥ 0, given by Eq. (16). As an illustration of the possible values of the parameter γ that would arise from applications we consider the case with n = 0.4 and H = 0.0024 m. This corresponds to an example presented by Ng and Mei in connection with mud flow. The particular value of H is the maximum fluid thickness that results in a laminar flow down an impermeable slope with an inclination of θ = 0.1 rad. In order to obtain possible values for the parameters related to a porous bottom we resort to the experimental conditions set up by Thomas et al. [6], who used a porous substrate with porosity Φ = 0.375 and particle diameter d = 0.00286 m. Using these values we can calculate the modified permeability associated with power-law fluid flow from the formula (see [12]) n n+1 6 Φd nΦ κ= 25 3n + 1 3(1 − Φ) which yields κ ≈ 2.37 × 10−6 m1.4 . Now, the experiments conducted by Beavers and Joseph indicate that the parameter χ ranges from 0.1 to 4. Therefore, for the current example γ ranges from 0.001 to 0.3992.
The eigenvalues of the coefficient matrix of the derivatives with respect to x of the unknowns in Eqs. (20) and (21) are (22) s1,2 = βu¯ ± β(β − 1)u¯ 2 + αh. These eigenvalues are real since it can clearly be seen that β > 0 and β−1=
n3 > 0. (3n + 2)[n + (2n + 1)γ]2
As a result the governing equations constitute a system of nonlinear hyperbolic conservation laws with a source term. The particular case of zero permeability for the surface over which the fluid layer is flowing is obtained by setting γ = 0. In this case we have β = 2(2n+1) 3n+2 and the governing equations reduce to those obtained Ng and Mei [1] for flow down an impermeable incline. 3. Linear stability analysis Suppose the scaled uniform flow is perturbed as follows: h = 1 + η(x, t),
u¯ = 1 + v(x, t),
where η, v 1. The linearized equations can be combined to produce the following equation in η: ηtt + (2β − α)ηxx + 2βηxt + (2n + 1)ηx + nηt = 0. If we consider a normal mode for the disturbance η = ηˆ ei(kx−ωt) , where k is the real wavenumber, then the growth rate is given by σ ≡ I(ω) while the phase speed is given by R(ω)/k. For the dispersion equation we obtain ω2 − (2βk − ni)ω + (β − α)k2 − (2n + 1)ik = 0. Solving this equation we obtain √ 1 ω = [2βk − ni ± a + ib], 2 where a = 4(β2 − β + α)k2 − n2 ,
b = 4(2n + 1 − nβ)k.
From this we obtain 1 1 2 I(ω) = −n ± ( a + b2 − a) 2 2 and
1 1 2 2 R(ω) = 2βk ± ( a + b + a) . 2 2
(23)
It can be shown that I(ω) > 0 if n2 (β2 − β + α) < (2n + 1 − nβ)2 . Thus the flow will be unstable if (2 n + 1)[(4 n + 5 n2 + 2 n3 + 1) γ 2 α < αs =
+ (2 n3 + 2 n + 4 n2 )γ + n2 ] [n + (2 n + 1)γ]2 n2
,
(24)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
with disturbances of any wavenumber being amplified in time. If γ = 0 then this reduces to αs =
2n + 1 , n2
(25)
which is the expression obtained by Ng and Mei [1] for the case of flow down an impermeable incline. The instability condition (24) can also be expressed as Λ < Λs =
(2 n + 1)1−n [(4 n + 5 n2 + 2 n3 + 1)γ 2 n4−n + (2 n3 + 2 n + 4 n2 ) γ + n2 ],
113
For the case of a Newtonian fluid layer flowing down a porous incline the critical conditions for the onset of instability are given by Λ < Λs = 12γ 2 + 8γ + 1. This result is in good agreement with a previous investigation θ 6 [9] where we determined that the flow is unstable if cot Re < 5 + 36 2 5 γ + O(γ ). This expression was rigorously obtained from the solution of the Orr-Sommerfeld equation under the assumption of long-wave perturbations and a Reynolds number of order unity.
(26)
where Λ ≡ cot θ/Re is a nondimensional number representing the ratio of stabilizing to destabilizing forces and can be used instead of α to describe the flow conditions of the uniform flow. The parameter Λ is independent of the permeability parameter γ, but it is a function of all the other independent parameters defining the uniform flow conditions including the fluid properties ρ, µ, and n, the thickness of the layer H, and the inclination of the plane, θ. The graphs of Λs as a function of γ for several values of n are presented in Fig. 1. From this instability condition we can deduce that increasing the permeability parameter, γ or decreasing the power-law index, n results in an increase in the instability interval [0, Λs ] for the uniform flow parameter Λ. We can thus conclude that the effect of the permeability of the inclined plane and the effect of the shear-thinning nature of the fluid is to destabilize the flow. Physically, this behavior is to be expected since we know that the effect of friction is to retard the flow and thus act as a stabilizing factor. Increasing the permeability of the porous substrate diminishes the friction of the fluid with the bottom as, due to the slip condition, the velocity gradient at the bottom of the fluid layer decreases as the modified permeability, κ is increased. Similarly, increasing the shear-thinning aspect of the flow diminishes the laminar friction in the flow which is caused by viscous forces.
4. Nonlinear analysis The linear stability analysis is insufficient in a number of important aspects which thus necessitates an analysis of the nonlinear effects. Firstly, nonlinear interactions may affect the evolution of disturbances once the amplitude is amplified to finite values. Therefore, the critical conditions for the onset of instability obtained from the linear analysis may be inaccurate. Secondly, the analysis described in the previous section indicates that in supercritical conditions disturbances of all wavelengths become unstable, however experiments [13] reveal the existence of a specific wavelength for roll waves developed spontaneously from unstable infinitesimal waves. Finally, considering the nonlinear governing equations reveals information about the structure of the roll waves. 4.1. Existence of permanent roll waves solutions In this section we analytically investigate the existence of permanent travelling wave solutions. We extend the theory developed by Ng and Mei [1] to the case of flow over a permeable surface and determine the parameter values for which the nonlinear governing equations (20) and (21) admit permanent travelling wave solutions consistent with the roll waves phenomenon. A procedure can then be employed which yields the pertinent aspects of the wave pattern including the wavelength, height of the bore, and the wave speed. We seek a travelling wave solution of the form of a smooth profile of finite length. We thus introduce the moving coordinate ξ = x − ct, moving with a constant speed c. As a result the governing equations become h(c − u) ¯ = q,
(27)
dh A(h) = , dξ B(h)
(28)
where q is a constant which represents the flow rate relative to a coordinate system moving with speed c, and c q n A(h) = h − − 2 , h h Fig. 1. Critical flow conditions for instability as a function of γ.
B(h) =
αh3 + (β − 1)c2 h2 − βq2 . h2
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
114
The observed wave pattern associated with roll waves flows consists of a smooth depth profile between successive bores which is increasing in the downhill direction. Thus, we seek solutions with dh dξ > 0, which requires that A and B have the same sign over the range of h. Examining the polynomial in the numerator in the expression for B(h) it can be seen that since β − 1 > 0, there is only one change in sign in the coefficients of the powers of h. We can thus conclude that B(h) has only one positive root and furthermore, B(h) changes sign from negative to positive across this root. We will now show that the root is a section in the flow. Consider a wavelength of the roll waves solution connecting two subsequent bores with h1 denoting the depth at the left end and h2 at the right end. The eigenvalue s1 given by (22) must be an increasing function on the interval (h1 , h2 ) since otherwise positively sloped characteristics would intersect contradicting the required smoothness of the solution. Since the bores move with speed c then s1 (h1 , u¯ 1 ) < c < s1 (h2 , u¯ 2 ). As a result there is a critical section in the flow with h = hc and u¯ = u¯ c (h1 < hc < h2 , u¯ 1 < u¯ c < u¯ 2 ) such that s1 (hc , u¯ c ) = c. This equation can also be expressed as B(hc ) = 0.
(29)
We will now show that due to the energy requirement we may assume that q is positive. In terms of the nondimensional variables the change in mechanical energy across a bore is given by h
1 h 1 1 c 3 2 2 2 ˙ = E u dz − c u dz + uh ¯ − h α , (30) 2 0 2 0 2 2
where [f ]12 ≡ f2 − f1 with f2 denoting the limit of the quantity f as the bore is approached from the left and f1 the limit of f as the bore is approached from the right. From our approximation regarding the variation of the velocity with depth given by Eq. (15) we have h h u2 dz = βhu¯ 2 , u3 dz = φhu¯ 3 , 0
0
where (2n + 1)2 [(24n3 + 46n2 + 29n + 6)γ 3 + (36n3 + 51n2 + 18n)γ 2 + (18n2 + 24n3 )γ + 6n3 ] . [n + (1 + 2n)γ]3 (2 + 3n)(4n + 3) The energy drop across the bore must be nonpositive since mechanical energy is dissipated into heat due to turbulence. Since this mechanism is not accounted for in our model, the criterion of nonpositive energy generation across the bore can be used to dismiss certain solutions as being nonphysical. Now, the shock conditions at the bore can be expressed as 1 2 1 1 1 1 2 c[h]2 = [uh] ¯ 2, c[uh] ¯ 2 = βu¯ h + αh . 2 2 φ=
Combining these two conditions and upon using Eq. (27) to eliminate u¯ we obtain 1 2 2 βq − h1 h2 c (β − 1) + α(h1 + h2 ) = 0. 2 Using this equation to simplify the expression for the energy drop gives ˙ = (−mq + ζ)(h2 − h1 ), E where m=
φα(h2 − h1 )2 φ(β − 1)c2 (h2 + h1 ) φ−β α+ + 4βh2 h1 β 2βh2 h1
and ζ=
3(φ − β)αc(h2 + h1 ) (2φβ − 3φ + β)c3 + . 4β 2β
It turns out that φ − β > 0 and 2φβ − 3φ + β > 0 for all the relevant parameter values, and as a result we have m > 0 and ζ > 0. Therefore, for roll waves solutions the value of q is positive, ˙ ≤ 0 indicates since h2 > h1 and thus the physical requirement E that ζ q≥ > 0. (31) m For positive values of q we can show that A(h) has either no positive roots or two positive roots, with A(h) < 0 between the roots and A(h) > 0 otherwise. In order for dh dξ to be finite it is necessary that A(hc ) = 0,
(32)
and hc be the only root of A(h) on the interval [h1 , h2 ]. Solving for q and c from (29) and (32) we obtain 4+2/n 2+1/n q = (β − 1)hc + β(β − 1)hc + αh3c (33) and c = h1+1/n + c
q . hc
(34)
In order for A and B to have the same sign we must have A (hc ) > 0, which can be expressed as 1+n > q. h2+1/n c n Using the expressions for c and q from above we obtain the following condition for α: (2n + 1)2 − (3n2 + 2n)β 1+2/n G(hc ) ≡ hc > α. (35) n2 For a uniform flow on λ which an infinitesimal disturbance is imposed, we have λ−1 0 h dξ ≈ 1. Then by conservation the same holds for the roll waves solution arising from this perturbed flow. We will thus assume that h = 1 is a section in the flow, i.e. h1 < 1 < h2 . So if hc < 1 then A(1) must be positive and if hc > 1 then A(1) must be negative. These conditions can be expressed as > 1 + q, hc > 1, c < 1 + q, hc < 1.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
115
Fig. 2. (a) Graphs of the functions F (hc ) and G(hc ). (b) Graphs of the function F (hc ).
Introducing the expression for c we find that both inequalities are satisfied if 2+1/n
hc
− hc > q, hc − 1
which upon substituting the expression for q gives another condition for α: −(β − 1) hc 4+2/n + βhc 2+2/n F (hc ) ≡
+ 2(β − 1)hc 2+1/n − 2 βhc 1+1/n + 1 > α. hc (hc − 1)2
(36)
It turns out that F (1) = G(1) = αs , with αs given in (24). It can be shown that for n = 1, hc = 1 is the only solution to F (hc ) = G(hc ) for all relevant values of γ. A numerical investigation indicates this to be true for general values of n. The general relationship between F (hc ) and G(hc ) is illustrated in Fig. 2a. It can be seen that for α < αs there are values of hc for which the two conditions (35) and (36) are simultaneously satisfied. Therefore, for these α values roll waves solutions exist. This is consistent with the linear analysis since α < αs corresponds to linearly unstable flows. However, these solutions may be nonphysical due to the violation of the energy requirement, or may not be the quasi-steady state of the evolution of infinitesimal disturbances, and thus do not constitute a wave pattern arising from the instability of the uniform flow. As it can be seen in Fig. 2a, G(hc ) > F (hc ) for hc > 1. So, if we let αmax = maxhc ≥1 F (hc ), then α values that satisfy α < αmax are admissible for roll waves solutions. Thus for cases with αmax > αs roll waves may form under stable flow conditions, provided that the mathematical solutions turn out to be physically relevant. If F (1) > 0 then αmax > αs as illustrated by curve 1 in Fig. 2b. It can be shown that F (1) > 0 if n[3n3 + 2n2 − 3n − 2 + n2 (2 + 3n)(1 − n)(n2 + 5n + 3)] γ > γ1 (n) ≡ . (1 − n)(2 + 3n)(1 + 2 n)(n + 1)
√ The only √ positive root of γ1 (n) is n = 1/ 2, with γ1 < 0 if 0 < n < 1/ 2. As a result we can state that in highly non-Newtonian cases αmax > αs for all relevant values of γ, including the case γ = 0 as was discovered by Ng and Mei. If F (1) ≤ 0 then F may have a maximum for hc > 1 that exceeds αs like the example illustrated by curve 2 in Fig. 2b. Or we could have F (hc ) < αs for hc > 1 as in the case of curve 3 in Fig. 2b. In this case, we can conclude that roll waves cannot be generated from a linearly stable flow. An analysis of the dependence of F (hc ) on the parameters of the problem reveal that for a given value of n there is a critical value of γ, γcrit (n) ≥ 0, such that if γ > γcrit then αmax > αs . exist for linearly stable So, if γ > γcrit then roll waves solutions √ flow conditions. Clearly, if 0 < n ≤ 1/ 2 then γcrit = 0. But √ for n > 1/ 2, γcrit is positive and, as it turns out, an increasing function of n. The dependence of γcrit on n is presented in Fig. 3.
Fig. 3. Critical value of γ for the existence of roll waves solutions in linearly stable flow as a function of n.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
116
4.2. Permanent roll waves calculations It is possible to obtain substantial information about the permanent roll waves solutions without having to determine the depth distribution, which would require a numerical scheme to solve the differential equation for h(ξ) given by (28). We can use this equation and the shock conditions to calculate the wavelength, λ, the height of the bore (wave height), ψ ≡ h2 − h1 , and the wave speed, c. Using the procedure outlined by Ng and Mei, for a specific set of values for the parameters n, γ, and α, we choose a value for hc and use (33) and (34) to obtain q and c. We then iterate over the value of h1 calculating the corresponding value of h2 from the shock conditions. Now, integrating (28) with respect to h we obtain h2 h2 dξ B(h) dh = λ= dh (37) h1 dh h1 A(h) and 1
h ≡ λ
0
λ
1 h dξ = λ
Fig. 4. The evolution of the perturbed uniform flow to roll waves solution with n = 0.6, γ = 0.2, Λ = 4, and λ = 10. h2
h1
which lead to h2 B(h) ( h − h) dh = 0. A(h) h1
hB(h) dh, A(h)
(38)
(39)
This condition is used to test the values of h1 and h2 . For the correct values of h1 and h2 the wavelength is then calculated using (37). The integrals in (37) and (39) must be evaluated numerically. A particular roll waves solution can be tested for physical validity by determining if the energy requirement is satisfied. This is accomplished by substituting the computed quantities into the condition (31). 4.3. Numerical approach The theory described above allows us to obtain the possible roll waves solutions to the governing equations. However, in order to determine if these solutions correspond to flows arising from the instability of the uniform flow, we must verify that the predicted roll waves pattern is generated by the evolution of an infinitesimal disturbance imposed on the uniform flow. We accomplish this by numerically solving Eqs. (20) and (21) and calculating the evolution of the uniform flow perturbed by certain infinitesimal disturbances. We considered periodic disturbances with various monochromatic and polychromatic profiles. The calculation is performed on the periodic domain x ∈ [0, λ], with λ being the wavelength of the disturbance. For the cases when the disturbance is amplified in time, the calculation is carried out until the wave speed and surface profile attain a steady state, providing the roll waves solution. An illustration of the evolution can be seen in Fig. 4 where we present snapshots of the surface of the fluid at various times. Eventually, snapshots at time intervals of length T ≈ 29.77 coincide. This indicates that the surface profile is unchanged and propagates with constant speed. The
solution thus describes roll waves with wavelength λ and wave speed λ/T . To solve the system of nonlinear hyperbolic conservation laws given by (20) and (21) we employ the quasi-steady wavepropagation algorithm developed by LeVeque [14]. This method is particularly effective in the context of stability analysis since it captures the balance between the source term and the flux gradient which is appropriate for cases when the solution is “quasisteady". The discretization of the equations is accomplished by a high-resolution finite volume scheme. The energy drop across the bore is calculated from (30) using the values of u¯ and h from the numerical solution, with the speed of propagation of the bore given by the formula c = [uh] ¯ 12 /[h]12 , which is derived from the shock conditions.
Fig. 5. Growth rate of disturbances as a function of the wavelength from the linear theory with n = 0.5 and Λ = 2.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Our numerical experiments indicate that the roll waves solution attained by the evolution is essentially independent of the type of disturbance if the amplitude is infinitesimal (less than 0.001). A specific solution obtained by the permanent roll waves calculation procedure can be numerically verified by selecting the periodic computation domain for the evolution of infinitesimal disturbances to be equal to the predicted wavelength of the roll waves. If roll waves result from the evolution, then the wave height and wave speed can be compared to the analytically predicted values. 5. Results Under conditions when the uniform flow is linearly unstable, i.e. Λ ≡ cot θ/Re < Λs , using the numerical approach we found that the roll waves predicted by the permanent waves theory are in fact the result of the evolution of infinitesimal distur-
117
bances imposed on the uniform flow, with excellent agreement for the wavelength, wave height and wave speed. For those linearly stable flow conditions where the analytical theory reveals the existence of roll waves solutions, the evolution of infinitesimal disturbances calculated by our numerical method indicates that they are damped in time. Thus suggesting that linearly stable flows are also nonlinearly stable. We were actually able to numerically produce roll waves that are in agreement with those obtained with the procedure described in Section 4.2 but only for disturbances with sufficiently large amplitude. For example, in the case with n = 0.4 and γ = 0.5 the uniform flow is linearly unstable if Λ < Λs ≈ 70.3 and the permanent roll waves theory indicates that roll waves solutions exist if Λ < Λmax ≈ 8404.7. With Λ = 72 > Λs one of the possible solutions is such that the wavelength is λ ≈ 169, the wave height is ψ ≈ 0.74, and the wave speed is c ≈ 6.16. We numerically generated roll waves with these parameters, but we had to use an initial disturbance with amplitude of 0.1. Infinitesimal
Fig. 6. (a) Wave height of the minimum roll wave as function of Λ for a Newtonian fluid. (b) Wave speed of the minimum roll wave as function of Λ for a Newtonian fluid. (c) Wavelength of the minimum roll wave as function of Λ for a Newtonian fluid.
118
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Fig. 7. (a) Wave height of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75. (b) Wave speed of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75. (c) Wavelength of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75.
disturbances were found to decay in time under these flow conditions. For a value of Λ for which there exist roll waves solutions with nonpositive energy drop across the bore, this quantity decreases with wavelength. So, there is a critical wavelength, λcrit such that for λ > λcrit the energy drop is nonpositive as physically required. The roll waves solution with this critical wavelength is referred to as the minimum roll wave [1]. Furthermore, as it is illustrated in Fig. 5, the linear stability analysis indicates that the growth rate is a decreasing function of wavelength, i.e. the shorter the disturbance the greater the amplification rate in time. Therefore, the minimum roll wave is the flow pattern expected to emerge from an unstable uniform flow. In Figs. 6–8 we present the various aspects of the minimum roll wave for n = 1, 0.75, and 0.4, respectively, with various values of γ. In Fig. 8 we also include some numerical results to illustrate the excellent agreement between the solutions obtained
through the permanent roll waves theory and those obtained by means of our numerical approach. The linear stability analysis yields a formula for the speed of disturbances in terms of their wavelength given by R(ω)/k, with R(ω) given by (23). In Fig. 7b we present the speed given by this formula for disturbances with wavelength equal to that of the minimum roll waves. It can be seen that the wave speed predicted by the linear theory is accurate for small values of Λ with the accuracy increasing with γ. Now, as indicated by Fig. 7c, the wavelength of the minimum roll is small for small values of Λ, and for short disturbances we have the asymptotic expression R(ω) 1 2 =β+ β −β+α+O . k k2 This formula can thus serve as an approximation for the velocity of roll waves generated from highly unstable flows.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
119
Fig. 8. (a) Wave height of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4. (b) Wave speed of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4. (c) Wavelength of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4.
From the results in Figs. 6–8 we can gain information about the effect of the permeability of the inclined plane and the correlation with the effect of the rheology of the fluid and the flow conditions of the uniform flow. The range of Λ for the results increases with the scaled permeability parameter γ, indicating that the permeability of the surface over which the fluid is flowing has the effect of increasing the range of flow conditions under which roll waves form. The results in Figs. 6a, 7a and 8a indicate that the wave height of the minimum roll wave decreases with γ, with the dependence on γ being greatest for highly unstable flow conditions (small values of Λ). There appears to be less dependence on γ for highly non-Newtonian fluids. 6. Concluding remarks We derived a model for power-law fluid flow down a porous incline and used it to study the stability of a uniform flow and
the formation of roll waves. Using an analytical permanent roll waves theory we determined flow conditions under which roll waves solutions are possible. The theory leads to a procedure for calculating quantities relevant to describing the structure of these roll waves. Our numerical method determines the evolution of the perturbed uniform flow and thus provides a complete nonlinear stability analysis. Under linearly unstable flow conditions, the numerical approach indicates that infinitesimal disturbances are amplified in time and evolve into roll waves flows. These roll wave patterns are in excellent agreement with those predicted by the permanent roll waves theory. Roll waves solutions predicted by the permanent roll waves theory under linearly stable conditions could not be generated by our numerical approach if the initial flow consist of infinitesimal disturbances imposed on the uniform flow. These disturbances were actually damped. Thus, our numerical stability analysis indicates that linearly stable flows are also nonlinearly stable.
120
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Acknowledgement Financial support for this research was provided by the Natural Sciences and Engineering Research Council of Canada. References [1] C. Ng, C.C. Mei, Roll waves on a shallow layer of mud modelled as a power-law fluid, J. Fluid Mech. 263 (1994) 151–183. [2] G. Vereet, J. Berlamont, Rheology and non-Newtonian behavior of sea and estuarine mud, in: N.P. Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, vol. 7, Gulf, Huston, TX, 1988 (Chapter 5). [3] R.I. Tanner, Engineering Rheology, Clarendon, Oxford, 1985. [4] R.F. Dressler, Mathematical solution of the problem of roll waves in inclined open channels, Commun. Pure Appl. Math. 2 (1949) 149–194. [5] B. Zanuttigh, A. Lamberti, Roll waves simulation using shallow water equations and weighted average flux method, J. Hydraul. Res. 40 (2002) 610–622. [6] L.P. Thomas, B.M. Marino, P.F. Linden, Lock-release inertial gravity currents over a thick porous layer, J. Fluid Mech. 503 (2004) 299–319.
[7] C. Cenedese, J.A. Whitehead, T.A. Ascarelli, M. Ohiwa, A dense current flowing down a sloping bottom in a rotating fluid, J. Phys. Oceanogr. 34 (2004) 188–203. [8] G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30 (1967) 197–207. [9] J.P. Pascal, Linear stability of fluid flow down a porous inclined plane, J. Phys. D: Appl. Phys. 32 (1999) 417–422. [10] A.R. Rao, M. Mishra, Peristaltic transport of a power-law fluid in a porous tube, J. Non-Newtonian Fluid Mech. 121 (2004) 163–174. [11] H. Pascal, Rheological behaviour effect of non-Newtonian fluids on steady and unsteady flow through a porous medium, Int. J. Numer. Anal. Meth. Geomechan. 7 (1983) 298–303. [12] R.H. Christopher, S. Middleman, Power-law flow through a packed tube, I & EC Fund. 4 (1965) 422–426. [13] P.Y. Julien, D.M. Hartley, Formation of roll waves in laminar sheet flow, J. Hydraul. Res. 24 (1986) 5–17. [14] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm, J. Comput. Phys. 146 (1998) 346–365.
Instability of power-law fluid flow down a porous incline J.P. Pascal Department of Mathematics, Ryerson University, Toronto, Ont., M5B 2K3, Canada Received 25 July 2005; received in revised form 14 November 2005; accepted 19 November 2005
Abstract In this paper we investigate the generation and structure of roll waves developing on the surface of a power-law fluid layer flowing down a porous incline. The unsteady equations of motion for the power-law fluid layer are depth integrated according to the von K´arm´an momentum integral method accounting for the variation of the velocity distribution with depth. The slip boundary condition at the interface between the fluid layer and the porous plane is based on the assumption that the flow through the porous medium is governed by the modified Darcy’s law, and that the characteristic length scale of the pore space is much smaller than the depth of the fluid layer above. An analytical theory of permanent roll waves is employed to determine under what flow conditions roll waves can exist and to calculate the wavelength, wave height, and speed of the roll waves. The nonlinear stability analysis is also carried out by numerically solving the time dependent governing equations and calculating the nonlinear evolution of infinitesimal disturbances imposed on the uniform and steady flow. Conclusions are drawn regarding the effect of the permeability of the porous medium and flow conditions on the development and characteristics of roll waves arising from the instability of the uniform flow. © 2005 Elsevier B.V. All rights reserved. Keywords: Flow down a porous incline; Power-law fluid; Linear and nonlinear stability analysis; Roll waves
1. Introduction It has been observed that the steady flow of a uniform fluid layer in an inclined open channel becomes unstable under certain conditions. The instability evolves into a flow of nonconstant depth consisting of progressive bores connected by sections of gradually varying flow. This wave system is referred to as roll waves and has significant implications in various flows found in the environment. Of particular interest are situations when the fluid exhibits non-Newtonian rheological properties, as for example molten lava flows and mud flows [1]. It has been reported [2] that for many such flows, at low shear rates the shear rate–shear stress relation is of a shear-thinning nature and is appropriately modelled by the power-law model which in simple shear is expressed as [3]: ∂u ∂u n−1 τ=µ , ∂z ∂z where µ is the consistency and n is the power-law index of the fluid. For shear-thinning fluids n is between 0 and 1, and the case with n = 1 corresponding to a Newtonian fluid with µ being the dynamic viscosity. E-mail address: [email protected]. 0377-0257/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2005.11.007
Due to the relevant impact of the instability of flow down an incline, it is important to theoretically determine the critical conditions for the onset of instability and to predict the evolution of the unstable flow. By performing a linear stability analysis on the governing equations, one determines the range of the flow parameters for which infinitesimal disturbances grow exponentially in time. However, as the disturbances grow the nonlinear interactions become significant and a nonlinear analysis must be conducted in order to predict the subsequent evolution. The shallow-water equations have been widely used as the equations governing the velocity and depth of the fluid layer. Most notable is the early work by Dressler [4] who analytically obtained roll waves solutions to the shallow-water equations augmented with an empirical term to model the turbulent frictional force. Recently, Zanuttigh and Lamberti [5] obtained roll waves solutions by numerically solving the de St. Venant equations and calculating the evolution of the flow until roll waves are attained. The shallow-water equations are based on the assumption that the longitudinal velocity is uniform with depth (flat velocity profile) and are thus only appropriate for high Reynolds number flows. For laminar flows, the governing equations must account for the variation of the velocity with depth due to viscous effects, and incorporate the rheological properties of the fluid. In order to include the variation of the velocity with depth, an accurate approximation can be obtained by assuming the
110
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
velocity profile to be the same as that of the uniform and steady flow. This approximation is used in connection with von K´arm´an’s momentum integral method which consists of depth integrating the equations established under the assumption of moderate Reynolds number flow and disturbances of long wavelength (small aspect ratio). This approach was taken by Ng and Mei [1] to investigate the instability and roll waves structure of a power-law fluid flow down an impermeable incline. They carry out an analytical investigation by seeking roll waves solutions characterized as periodic shocks connected by smoothly increasing depth profiles. They obtain necessary conditions for the existence of roll waves, and present a procedure for calculating the wavelength, speed, and wave height of the roll waves. In this paper we extend the problem of instability of powerlaw fluid flow down an incline by incorporating the effect of the permeability of the inclined plane. This extension should prove useful for modelling environmental flows which mostly occur over porous surfaces. Such an extension is now being explored in the literature in connection with gravity currents. Much attention is being focussed on incorporating the effect of bottom permeability in the study of the spread of dense fluids under layers of lighter fluid, with the goal of obtaining more realistic models for gravity-driven flows over such surfaces as soil, river beds and layers of gravel. A summary of these investigations can be found in a recent paper by Thomas et al. [6]. It should also be pointed out that the layer of fluid behind the head of a gravity current resulting from the spread of dense fluid down an incline is subject to instability. As was recently discovered by Cenedese et al. [7], this uniform flow can evolve into a wave regime exhibiting a roll waves pattern. In order to determine the governing equations for our problem we apply von K´arm´an’s momentum integral method to obtain a model for the flow down an incline with a depth dependent longitudinal velocity. At the interface between the fluid layer and the porous medium we employ the slip conditions experimentally verified by Beavers and Joseph [8]. We begin our investigation by first conducting a linear stability analysis, which yields an analytical formula for the relationship between the Reynolds number, the inclination of the plane, and the permeability parameter which describes the critical conditions for the exponential growth of infinitesimal disturbances. Then by extending the analytical theory devised by Ng and Mei we study the existence of permanent wave solutions and calculate the relevant quantities describing the roll waves flow. Finally, we also carry out a nonlinear stability analysis by numerically solving the nonlinear governing equations to determine the evolution of a perturbed uniform flow. This will allow us to determine which of the analytically obtained permanent roll wave solutions do in fact arise from the evolution of infinitesimal disturbances.
u = (u, w)T and p, respectively. We assume the characteristic longitudinal length of the flow, L to be much larger than the characteristic depth, H, with the depth varying slowly in time and in the longitudinal direction. The equations of motion for the layer are obtained from the Navier–Stokes equations by neglecting terms of O(2 ) with = H/L being the aspect ratio, and are expressed as [1]:
2. The governing equations
where κ is the permeability of the porous medium, χ is a dimensionless parameter dependent on the structure of the porous medium, and up = (up , wp )T is the Darcian mean filter velocity in the porous medium. The normal velocity does not vary through the boundary-layer so the condition at the interface is
Consider the two-dimensional laminar flow of a thin layer of a power-law fluid along the surface of a porous medium incline at an angle θ with respect to the horizontal. We define an (x, z) coordinate system with the x-axis along the bottom and the z-axis normal to it. We denote the velocity and the total pressure by
∂u ∂w + = 0, ∂x ∂z ρ
(1)
∂u ∂u ∂u +u +w ∂t ∂x ∂z
∂p ∂ ∂u ∂u n−1 = − + gρ sin θ + µ , ∂x ∂z ∂z ∂z 0=−
1 ∂p − g cos θ, ρ ∂z
(2)
(3)
where g is the acceleration due to gravity, and ρ, n and µ are the mass density, the power-law index, and the consistency of the power-law fluid, respectively. If we neglect terms of O(2 ), including those arising from the effect of surface tension, then the continuity of force condition at the surface reduces to p = 0,
∂u =0 ∂z
at z = h(x, t).
(4)
The kinematic condition at the surface is given by w=
∂h ∂h +u ∂t ∂x
at z = h(x, t).
(5)
We assume that the porous medium is saturated with fluid and that the flow through it is governed by the modified Darcy’s Law. We note that even if the free flow is a high Reynolds number flow, the application of the modified Darcy’s Law to model the flow through the porous medium is validated by the assumption that the geometry of the pore space results in a low Reynolds number flow. The condition at the fluid layer-porous medium interface in this case has been investigated by Beavers and Joseph [8]. They propose that a boundary-layer is formed at the top of the porous medium and that the tangential velocity rapidly changes across this region from the Darcian velocity to the velocity at the interface. Beavers and Joseph verify experimentally that the appropriate condition for the tangential velocity at the interface is ∂u χ = √ (u − up ) ∂z κ
w = wp
at z = 0.
at z = 0,
(6)
(7)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
The conditions (6) and (7) have been used previously [9] in a model for Newtonian fluid flow down a porous incline. In the case of a general power-law fluid, condition (6) is expressed as ∂u χ = 1/(n+1) (u − up ) at z = 0, (8) ∂z κ where κ is the modified permeability of the porous medium with dimension of length to power n + 1. This extension of the Beavers–Joseph boundary condition has been recently used by Rao and Mishra [10] to study peristaltic transport of a power-law fluid in a porous tube. According to the modified Darcy’s Law which includes the nonlinear rheological properties of non-Newtonian power-law fluids, the flow in the porous medium is governed by [11]: 1/n ∂pp ∂pp 1/n−1 κ up = − − ρg sin θ , (9) µ ∂x ∂x ∂up ∂wp + = 0, ∂x ∂z ∂pp = −ρg cos θ, ∂z pp = p at z = 0,
(10) (11) (12)
where pp denotes the total pressure in the porous medium. We will consider cases when the pore space geometry is such that the flow is slower than that of the layer above to the extent that the mean filter velocity terms in the interface conditions (8) and (7) are negligible. We will now determine the order of magnitude of the permeability of the porous medium that is required for these terms to be of O(2 ). Based on the hydrostatic nature of the pressure distribution revealed by Eqs. (3) and (11) and the continuity of pressure ∂p condition (12) we conclude that we can use ∂xp = ∂p ∂x in Eq. (9). But the characteristic pressure in the fluid layer above the porous medium is n L U , P ∼µ H H where U denotes the characteristic longitudinal velocity. From Eq. (9) we determine the characteristic longitudinal velocity in the porous medium to be κ1/n U . H 1+1/n If we scale the variables in the boundary conditions at the fluid layer-porous medium interface, (8) and (7) with their characteristic values we obtain ∂u∗ χ = u∗ − χδ1/n u∗p , ∂z∗ δ and Up ∼
w∗ = δ2+1/n w∗p , 1/(n+1)
where δ = κ H is the ratio of the characteristic length scale of the pore space to the characteristic depth in the fluid layer flowing above the porous medium. If we assume that δ ≤ O(2n /χn ), then the up term in (8) is of O(2 ), and if δ ≤ O(2n/(2n+1) ) then
111
the wp term in (7) is of O(2 ). We also mention that we must also assume that δ ≤ O(3n/(n+2) ) for the Reynolds number of the porous medium flow to be of O(2 ) when the Reynolds number of the flow above is of O(1/). If we then neglect terms of O(2 ) the conditions at the bottom of the fluid layer reduce to ∂u χ (13) = 1/(n+1) u, w = 0, z = 0. ∂z κ We can now investigate the flow over the porous plane with the effect of the permeability of the surface incorporated through the slip condition in (13) and the parameters κ and χ. Integrating Eq. (3) and using the boundary condition (4) we find the expression for the total pressure to be: p = g cos θρ(h − z). Inserting these expressions into the x-momentum equation (2) we obtain ∂u ∂u ∂u +u +w ∂t ∂x ∂z ∂h µ ∂ ∂u n−1 ∂u = g sin θ − g cos θ + . (14) ∂x ρ ∂z ∂z ∂z 2.1. Depth-averaged equations A uniform and steady flow of the fluid layer is given by h = H = const. n(2n + 1)u¯ 0 (n + 1)[n + γ(2n + 1)] z 1+1/n n + 1 × 1− 1− + γ , H n
u = u0 (z) =
(15)
where γ=
κ1/(n+1) , χH
(16)
and u¯ is the depth-averaged velocity given by 1 H u¯ 0 = u0 (z) dz. H 0 Since we are dealing with a small aspect ratio flow which varies slowly in the longitudinal direction, we can apply von K´arm´an’s momentum integral method. We thus assume that the relationship between the velocity and the fluid depth for the steady flow is also true when these quantities are time dependent and nonuniform. Depth-integrating the continuity equation (1) and the momentum equation (14) over the fluid layer we obtain ∂h ∂(hu) ¯ + = 0, ∂t ∂x ∂h µ ∂ ∂(hu) ¯ + βhu¯ 2 = g sin θh − g cos θh − ∂t ∂x ∂x ρ n n u¯ 2n + 1 , × n + γ(2n + 1) h
(17)
(18)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
112
where β=
2 2n + 1 2nγ 2n2 2 γ + + . n + γ(2n + 1) 2n + 1 (3n + 2)(2n + 1) (19)
We now nondimensionalize these equations by introducing the nondimensional quantities (signified by an asterisk) according to L ∗ x = Lx∗ , u¯ = u¯ 0 u¯ ∗ , t= t , h = Hh∗ , u¯ 0 where u¯ 0 and H are the velocity and thickness of the uniform steady flow, respectively, and L=
u¯ 20 . g sin θ
Employing this transformation to nondimensionalize equations (17) and (18) and dropping the asterisks for notational convenience we obtain ∂h ∂(hu) ¯ + = 0, (20) ∂t ∂x n u¯ 1 ∂ ∂(hu) ¯ , (21) + βhu¯ 2 + αh2 = h − ∂t ∂x 2 h where
2 2n + 1 n−2 2n + 1 cot θ , α= n Re n + γ(2n + 1) 2/n 2−n n ρ 1+2/n 2/n−1 Re = H (g sin θ) , µ 2n + 1 and the expression for β is given by (19). In the case when γ = 0 ρu¯ 2−n H n
the Reynolds number can be expressed as Re = 0 µ . Eqs. (20) and (21) govern the depth and the depth-averaged velocity of the fluid layer. As parameters in these equations we have the power-law index, 0 < n ≤ 1, the uniform flow parameter, α ≥ 0, and the scaled permeability, γ ≥ 0, given by Eq. (16). As an illustration of the possible values of the parameter γ that would arise from applications we consider the case with n = 0.4 and H = 0.0024 m. This corresponds to an example presented by Ng and Mei in connection with mud flow. The particular value of H is the maximum fluid thickness that results in a laminar flow down an impermeable slope with an inclination of θ = 0.1 rad. In order to obtain possible values for the parameters related to a porous bottom we resort to the experimental conditions set up by Thomas et al. [6], who used a porous substrate with porosity Φ = 0.375 and particle diameter d = 0.00286 m. Using these values we can calculate the modified permeability associated with power-law fluid flow from the formula (see [12]) n n+1 6 Φd nΦ κ= 25 3n + 1 3(1 − Φ) which yields κ ≈ 2.37 × 10−6 m1.4 . Now, the experiments conducted by Beavers and Joseph indicate that the parameter χ ranges from 0.1 to 4. Therefore, for the current example γ ranges from 0.001 to 0.3992.
The eigenvalues of the coefficient matrix of the derivatives with respect to x of the unknowns in Eqs. (20) and (21) are (22) s1,2 = βu¯ ± β(β − 1)u¯ 2 + αh. These eigenvalues are real since it can clearly be seen that β > 0 and β−1=
n3 > 0. (3n + 2)[n + (2n + 1)γ]2
As a result the governing equations constitute a system of nonlinear hyperbolic conservation laws with a source term. The particular case of zero permeability for the surface over which the fluid layer is flowing is obtained by setting γ = 0. In this case we have β = 2(2n+1) 3n+2 and the governing equations reduce to those obtained Ng and Mei [1] for flow down an impermeable incline. 3. Linear stability analysis Suppose the scaled uniform flow is perturbed as follows: h = 1 + η(x, t),
u¯ = 1 + v(x, t),
where η, v 1. The linearized equations can be combined to produce the following equation in η: ηtt + (2β − α)ηxx + 2βηxt + (2n + 1)ηx + nηt = 0. If we consider a normal mode for the disturbance η = ηˆ ei(kx−ωt) , where k is the real wavenumber, then the growth rate is given by σ ≡ I(ω) while the phase speed is given by R(ω)/k. For the dispersion equation we obtain ω2 − (2βk − ni)ω + (β − α)k2 − (2n + 1)ik = 0. Solving this equation we obtain √ 1 ω = [2βk − ni ± a + ib], 2 where a = 4(β2 − β + α)k2 − n2 ,
b = 4(2n + 1 − nβ)k.
From this we obtain 1 1 2 I(ω) = −n ± ( a + b2 − a) 2 2 and
1 1 2 2 R(ω) = 2βk ± ( a + b + a) . 2 2
(23)
It can be shown that I(ω) > 0 if n2 (β2 − β + α) < (2n + 1 − nβ)2 . Thus the flow will be unstable if (2 n + 1)[(4 n + 5 n2 + 2 n3 + 1) γ 2 α < αs =
+ (2 n3 + 2 n + 4 n2 )γ + n2 ] [n + (2 n + 1)γ]2 n2
,
(24)
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
with disturbances of any wavenumber being amplified in time. If γ = 0 then this reduces to αs =
2n + 1 , n2
(25)
which is the expression obtained by Ng and Mei [1] for the case of flow down an impermeable incline. The instability condition (24) can also be expressed as Λ < Λs =
(2 n + 1)1−n [(4 n + 5 n2 + 2 n3 + 1)γ 2 n4−n + (2 n3 + 2 n + 4 n2 ) γ + n2 ],
113
For the case of a Newtonian fluid layer flowing down a porous incline the critical conditions for the onset of instability are given by Λ < Λs = 12γ 2 + 8γ + 1. This result is in good agreement with a previous investigation θ 6 [9] where we determined that the flow is unstable if cot Re < 5 + 36 2 5 γ + O(γ ). This expression was rigorously obtained from the solution of the Orr-Sommerfeld equation under the assumption of long-wave perturbations and a Reynolds number of order unity.
(26)
where Λ ≡ cot θ/Re is a nondimensional number representing the ratio of stabilizing to destabilizing forces and can be used instead of α to describe the flow conditions of the uniform flow. The parameter Λ is independent of the permeability parameter γ, but it is a function of all the other independent parameters defining the uniform flow conditions including the fluid properties ρ, µ, and n, the thickness of the layer H, and the inclination of the plane, θ. The graphs of Λs as a function of γ for several values of n are presented in Fig. 1. From this instability condition we can deduce that increasing the permeability parameter, γ or decreasing the power-law index, n results in an increase in the instability interval [0, Λs ] for the uniform flow parameter Λ. We can thus conclude that the effect of the permeability of the inclined plane and the effect of the shear-thinning nature of the fluid is to destabilize the flow. Physically, this behavior is to be expected since we know that the effect of friction is to retard the flow and thus act as a stabilizing factor. Increasing the permeability of the porous substrate diminishes the friction of the fluid with the bottom as, due to the slip condition, the velocity gradient at the bottom of the fluid layer decreases as the modified permeability, κ is increased. Similarly, increasing the shear-thinning aspect of the flow diminishes the laminar friction in the flow which is caused by viscous forces.
4. Nonlinear analysis The linear stability analysis is insufficient in a number of important aspects which thus necessitates an analysis of the nonlinear effects. Firstly, nonlinear interactions may affect the evolution of disturbances once the amplitude is amplified to finite values. Therefore, the critical conditions for the onset of instability obtained from the linear analysis may be inaccurate. Secondly, the analysis described in the previous section indicates that in supercritical conditions disturbances of all wavelengths become unstable, however experiments [13] reveal the existence of a specific wavelength for roll waves developed spontaneously from unstable infinitesimal waves. Finally, considering the nonlinear governing equations reveals information about the structure of the roll waves. 4.1. Existence of permanent roll waves solutions In this section we analytically investigate the existence of permanent travelling wave solutions. We extend the theory developed by Ng and Mei [1] to the case of flow over a permeable surface and determine the parameter values for which the nonlinear governing equations (20) and (21) admit permanent travelling wave solutions consistent with the roll waves phenomenon. A procedure can then be employed which yields the pertinent aspects of the wave pattern including the wavelength, height of the bore, and the wave speed. We seek a travelling wave solution of the form of a smooth profile of finite length. We thus introduce the moving coordinate ξ = x − ct, moving with a constant speed c. As a result the governing equations become h(c − u) ¯ = q,
(27)
dh A(h) = , dξ B(h)
(28)
where q is a constant which represents the flow rate relative to a coordinate system moving with speed c, and c q n A(h) = h − − 2 , h h Fig. 1. Critical flow conditions for instability as a function of γ.
B(h) =
αh3 + (β − 1)c2 h2 − βq2 . h2
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
114
The observed wave pattern associated with roll waves flows consists of a smooth depth profile between successive bores which is increasing in the downhill direction. Thus, we seek solutions with dh dξ > 0, which requires that A and B have the same sign over the range of h. Examining the polynomial in the numerator in the expression for B(h) it can be seen that since β − 1 > 0, there is only one change in sign in the coefficients of the powers of h. We can thus conclude that B(h) has only one positive root and furthermore, B(h) changes sign from negative to positive across this root. We will now show that the root is a section in the flow. Consider a wavelength of the roll waves solution connecting two subsequent bores with h1 denoting the depth at the left end and h2 at the right end. The eigenvalue s1 given by (22) must be an increasing function on the interval (h1 , h2 ) since otherwise positively sloped characteristics would intersect contradicting the required smoothness of the solution. Since the bores move with speed c then s1 (h1 , u¯ 1 ) < c < s1 (h2 , u¯ 2 ). As a result there is a critical section in the flow with h = hc and u¯ = u¯ c (h1 < hc < h2 , u¯ 1 < u¯ c < u¯ 2 ) such that s1 (hc , u¯ c ) = c. This equation can also be expressed as B(hc ) = 0.
(29)
We will now show that due to the energy requirement we may assume that q is positive. In terms of the nondimensional variables the change in mechanical energy across a bore is given by h
1 h 1 1 c 3 2 2 2 ˙ = E u dz − c u dz + uh ¯ − h α , (30) 2 0 2 0 2 2
where [f ]12 ≡ f2 − f1 with f2 denoting the limit of the quantity f as the bore is approached from the left and f1 the limit of f as the bore is approached from the right. From our approximation regarding the variation of the velocity with depth given by Eq. (15) we have h h u2 dz = βhu¯ 2 , u3 dz = φhu¯ 3 , 0
0
where (2n + 1)2 [(24n3 + 46n2 + 29n + 6)γ 3 + (36n3 + 51n2 + 18n)γ 2 + (18n2 + 24n3 )γ + 6n3 ] . [n + (1 + 2n)γ]3 (2 + 3n)(4n + 3) The energy drop across the bore must be nonpositive since mechanical energy is dissipated into heat due to turbulence. Since this mechanism is not accounted for in our model, the criterion of nonpositive energy generation across the bore can be used to dismiss certain solutions as being nonphysical. Now, the shock conditions at the bore can be expressed as 1 2 1 1 1 1 2 c[h]2 = [uh] ¯ 2, c[uh] ¯ 2 = βu¯ h + αh . 2 2 φ=
Combining these two conditions and upon using Eq. (27) to eliminate u¯ we obtain 1 2 2 βq − h1 h2 c (β − 1) + α(h1 + h2 ) = 0. 2 Using this equation to simplify the expression for the energy drop gives ˙ = (−mq + ζ)(h2 − h1 ), E where m=
φα(h2 − h1 )2 φ(β − 1)c2 (h2 + h1 ) φ−β α+ + 4βh2 h1 β 2βh2 h1
and ζ=
3(φ − β)αc(h2 + h1 ) (2φβ − 3φ + β)c3 + . 4β 2β
It turns out that φ − β > 0 and 2φβ − 3φ + β > 0 for all the relevant parameter values, and as a result we have m > 0 and ζ > 0. Therefore, for roll waves solutions the value of q is positive, ˙ ≤ 0 indicates since h2 > h1 and thus the physical requirement E that ζ q≥ > 0. (31) m For positive values of q we can show that A(h) has either no positive roots or two positive roots, with A(h) < 0 between the roots and A(h) > 0 otherwise. In order for dh dξ to be finite it is necessary that A(hc ) = 0,
(32)
and hc be the only root of A(h) on the interval [h1 , h2 ]. Solving for q and c from (29) and (32) we obtain 4+2/n 2+1/n q = (β − 1)hc + β(β − 1)hc + αh3c (33) and c = h1+1/n + c
q . hc
(34)
In order for A and B to have the same sign we must have A (hc ) > 0, which can be expressed as 1+n > q. h2+1/n c n Using the expressions for c and q from above we obtain the following condition for α: (2n + 1)2 − (3n2 + 2n)β 1+2/n G(hc ) ≡ hc > α. (35) n2 For a uniform flow on λ which an infinitesimal disturbance is imposed, we have λ−1 0 h dξ ≈ 1. Then by conservation the same holds for the roll waves solution arising from this perturbed flow. We will thus assume that h = 1 is a section in the flow, i.e. h1 < 1 < h2 . So if hc < 1 then A(1) must be positive and if hc > 1 then A(1) must be negative. These conditions can be expressed as > 1 + q, hc > 1, c < 1 + q, hc < 1.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
115
Fig. 2. (a) Graphs of the functions F (hc ) and G(hc ). (b) Graphs of the function F (hc ).
Introducing the expression for c we find that both inequalities are satisfied if 2+1/n
hc
− hc > q, hc − 1
which upon substituting the expression for q gives another condition for α: −(β − 1) hc 4+2/n + βhc 2+2/n F (hc ) ≡
+ 2(β − 1)hc 2+1/n − 2 βhc 1+1/n + 1 > α. hc (hc − 1)2
(36)
It turns out that F (1) = G(1) = αs , with αs given in (24). It can be shown that for n = 1, hc = 1 is the only solution to F (hc ) = G(hc ) for all relevant values of γ. A numerical investigation indicates this to be true for general values of n. The general relationship between F (hc ) and G(hc ) is illustrated in Fig. 2a. It can be seen that for α < αs there are values of hc for which the two conditions (35) and (36) are simultaneously satisfied. Therefore, for these α values roll waves solutions exist. This is consistent with the linear analysis since α < αs corresponds to linearly unstable flows. However, these solutions may be nonphysical due to the violation of the energy requirement, or may not be the quasi-steady state of the evolution of infinitesimal disturbances, and thus do not constitute a wave pattern arising from the instability of the uniform flow. As it can be seen in Fig. 2a, G(hc ) > F (hc ) for hc > 1. So, if we let αmax = maxhc ≥1 F (hc ), then α values that satisfy α < αmax are admissible for roll waves solutions. Thus for cases with αmax > αs roll waves may form under stable flow conditions, provided that the mathematical solutions turn out to be physically relevant. If F (1) > 0 then αmax > αs as illustrated by curve 1 in Fig. 2b. It can be shown that F (1) > 0 if n[3n3 + 2n2 − 3n − 2 + n2 (2 + 3n)(1 − n)(n2 + 5n + 3)] γ > γ1 (n) ≡ . (1 − n)(2 + 3n)(1 + 2 n)(n + 1)
√ The only √ positive root of γ1 (n) is n = 1/ 2, with γ1 < 0 if 0 < n < 1/ 2. As a result we can state that in highly non-Newtonian cases αmax > αs for all relevant values of γ, including the case γ = 0 as was discovered by Ng and Mei. If F (1) ≤ 0 then F may have a maximum for hc > 1 that exceeds αs like the example illustrated by curve 2 in Fig. 2b. Or we could have F (hc ) < αs for hc > 1 as in the case of curve 3 in Fig. 2b. In this case, we can conclude that roll waves cannot be generated from a linearly stable flow. An analysis of the dependence of F (hc ) on the parameters of the problem reveal that for a given value of n there is a critical value of γ, γcrit (n) ≥ 0, such that if γ > γcrit then αmax > αs . exist for linearly stable So, if γ > γcrit then roll waves solutions √ flow conditions. Clearly, if 0 < n ≤ 1/ 2 then γcrit = 0. But √ for n > 1/ 2, γcrit is positive and, as it turns out, an increasing function of n. The dependence of γcrit on n is presented in Fig. 3.
Fig. 3. Critical value of γ for the existence of roll waves solutions in linearly stable flow as a function of n.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
116
4.2. Permanent roll waves calculations It is possible to obtain substantial information about the permanent roll waves solutions without having to determine the depth distribution, which would require a numerical scheme to solve the differential equation for h(ξ) given by (28). We can use this equation and the shock conditions to calculate the wavelength, λ, the height of the bore (wave height), ψ ≡ h2 − h1 , and the wave speed, c. Using the procedure outlined by Ng and Mei, for a specific set of values for the parameters n, γ, and α, we choose a value for hc and use (33) and (34) to obtain q and c. We then iterate over the value of h1 calculating the corresponding value of h2 from the shock conditions. Now, integrating (28) with respect to h we obtain h2 h2 dξ B(h) dh = λ= dh (37) h1 dh h1 A(h) and 1
h ≡ λ
0
λ
1 h dξ = λ
Fig. 4. The evolution of the perturbed uniform flow to roll waves solution with n = 0.6, γ = 0.2, Λ = 4, and λ = 10. h2
h1
which lead to h2 B(h) ( h − h) dh = 0. A(h) h1
hB(h) dh, A(h)
(38)
(39)
This condition is used to test the values of h1 and h2 . For the correct values of h1 and h2 the wavelength is then calculated using (37). The integrals in (37) and (39) must be evaluated numerically. A particular roll waves solution can be tested for physical validity by determining if the energy requirement is satisfied. This is accomplished by substituting the computed quantities into the condition (31). 4.3. Numerical approach The theory described above allows us to obtain the possible roll waves solutions to the governing equations. However, in order to determine if these solutions correspond to flows arising from the instability of the uniform flow, we must verify that the predicted roll waves pattern is generated by the evolution of an infinitesimal disturbance imposed on the uniform flow. We accomplish this by numerically solving Eqs. (20) and (21) and calculating the evolution of the uniform flow perturbed by certain infinitesimal disturbances. We considered periodic disturbances with various monochromatic and polychromatic profiles. The calculation is performed on the periodic domain x ∈ [0, λ], with λ being the wavelength of the disturbance. For the cases when the disturbance is amplified in time, the calculation is carried out until the wave speed and surface profile attain a steady state, providing the roll waves solution. An illustration of the evolution can be seen in Fig. 4 where we present snapshots of the surface of the fluid at various times. Eventually, snapshots at time intervals of length T ≈ 29.77 coincide. This indicates that the surface profile is unchanged and propagates with constant speed. The
solution thus describes roll waves with wavelength λ and wave speed λ/T . To solve the system of nonlinear hyperbolic conservation laws given by (20) and (21) we employ the quasi-steady wavepropagation algorithm developed by LeVeque [14]. This method is particularly effective in the context of stability analysis since it captures the balance between the source term and the flux gradient which is appropriate for cases when the solution is “quasisteady". The discretization of the equations is accomplished by a high-resolution finite volume scheme. The energy drop across the bore is calculated from (30) using the values of u¯ and h from the numerical solution, with the speed of propagation of the bore given by the formula c = [uh] ¯ 12 /[h]12 , which is derived from the shock conditions.
Fig. 5. Growth rate of disturbances as a function of the wavelength from the linear theory with n = 0.5 and Λ = 2.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Our numerical experiments indicate that the roll waves solution attained by the evolution is essentially independent of the type of disturbance if the amplitude is infinitesimal (less than 0.001). A specific solution obtained by the permanent roll waves calculation procedure can be numerically verified by selecting the periodic computation domain for the evolution of infinitesimal disturbances to be equal to the predicted wavelength of the roll waves. If roll waves result from the evolution, then the wave height and wave speed can be compared to the analytically predicted values. 5. Results Under conditions when the uniform flow is linearly unstable, i.e. Λ ≡ cot θ/Re < Λs , using the numerical approach we found that the roll waves predicted by the permanent waves theory are in fact the result of the evolution of infinitesimal distur-
117
bances imposed on the uniform flow, with excellent agreement for the wavelength, wave height and wave speed. For those linearly stable flow conditions where the analytical theory reveals the existence of roll waves solutions, the evolution of infinitesimal disturbances calculated by our numerical method indicates that they are damped in time. Thus suggesting that linearly stable flows are also nonlinearly stable. We were actually able to numerically produce roll waves that are in agreement with those obtained with the procedure described in Section 4.2 but only for disturbances with sufficiently large amplitude. For example, in the case with n = 0.4 and γ = 0.5 the uniform flow is linearly unstable if Λ < Λs ≈ 70.3 and the permanent roll waves theory indicates that roll waves solutions exist if Λ < Λmax ≈ 8404.7. With Λ = 72 > Λs one of the possible solutions is such that the wavelength is λ ≈ 169, the wave height is ψ ≈ 0.74, and the wave speed is c ≈ 6.16. We numerically generated roll waves with these parameters, but we had to use an initial disturbance with amplitude of 0.1. Infinitesimal
Fig. 6. (a) Wave height of the minimum roll wave as function of Λ for a Newtonian fluid. (b) Wave speed of the minimum roll wave as function of Λ for a Newtonian fluid. (c) Wavelength of the minimum roll wave as function of Λ for a Newtonian fluid.
118
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Fig. 7. (a) Wave height of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75. (b) Wave speed of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75. (c) Wavelength of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.75.
disturbances were found to decay in time under these flow conditions. For a value of Λ for which there exist roll waves solutions with nonpositive energy drop across the bore, this quantity decreases with wavelength. So, there is a critical wavelength, λcrit such that for λ > λcrit the energy drop is nonpositive as physically required. The roll waves solution with this critical wavelength is referred to as the minimum roll wave [1]. Furthermore, as it is illustrated in Fig. 5, the linear stability analysis indicates that the growth rate is a decreasing function of wavelength, i.e. the shorter the disturbance the greater the amplification rate in time. Therefore, the minimum roll wave is the flow pattern expected to emerge from an unstable uniform flow. In Figs. 6–8 we present the various aspects of the minimum roll wave for n = 1, 0.75, and 0.4, respectively, with various values of γ. In Fig. 8 we also include some numerical results to illustrate the excellent agreement between the solutions obtained
through the permanent roll waves theory and those obtained by means of our numerical approach. The linear stability analysis yields a formula for the speed of disturbances in terms of their wavelength given by R(ω)/k, with R(ω) given by (23). In Fig. 7b we present the speed given by this formula for disturbances with wavelength equal to that of the minimum roll waves. It can be seen that the wave speed predicted by the linear theory is accurate for small values of Λ with the accuracy increasing with γ. Now, as indicated by Fig. 7c, the wavelength of the minimum roll is small for small values of Λ, and for short disturbances we have the asymptotic expression R(ω) 1 2 =β+ β −β+α+O . k k2 This formula can thus serve as an approximation for the velocity of roll waves generated from highly unstable flows.
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
119
Fig. 8. (a) Wave height of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4. (b) Wave speed of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4. (c) Wavelength of the minimum roll wave as function of Λ for a non-Newtonian fluid having n = 0.4.
From the results in Figs. 6–8 we can gain information about the effect of the permeability of the inclined plane and the correlation with the effect of the rheology of the fluid and the flow conditions of the uniform flow. The range of Λ for the results increases with the scaled permeability parameter γ, indicating that the permeability of the surface over which the fluid is flowing has the effect of increasing the range of flow conditions under which roll waves form. The results in Figs. 6a, 7a and 8a indicate that the wave height of the minimum roll wave decreases with γ, with the dependence on γ being greatest for highly unstable flow conditions (small values of Λ). There appears to be less dependence on γ for highly non-Newtonian fluids. 6. Concluding remarks We derived a model for power-law fluid flow down a porous incline and used it to study the stability of a uniform flow and
the formation of roll waves. Using an analytical permanent roll waves theory we determined flow conditions under which roll waves solutions are possible. The theory leads to a procedure for calculating quantities relevant to describing the structure of these roll waves. Our numerical method determines the evolution of the perturbed uniform flow and thus provides a complete nonlinear stability analysis. Under linearly unstable flow conditions, the numerical approach indicates that infinitesimal disturbances are amplified in time and evolve into roll waves flows. These roll wave patterns are in excellent agreement with those predicted by the permanent roll waves theory. Roll waves solutions predicted by the permanent roll waves theory under linearly stable conditions could not be generated by our numerical approach if the initial flow consist of infinitesimal disturbances imposed on the uniform flow. These disturbances were actually damped. Thus, our numerical stability analysis indicates that linearly stable flows are also nonlinearly stable.
120
J.P. Pascal / J. Non-Newtonian Fluid Mech. 133 (2006) 109–120
Acknowledgement Financial support for this research was provided by the Natural Sciences and Engineering Research Council of Canada. References [1] C. Ng, C.C. Mei, Roll waves on a shallow layer of mud modelled as a power-law fluid, J. Fluid Mech. 263 (1994) 151–183. [2] G. Vereet, J. Berlamont, Rheology and non-Newtonian behavior of sea and estuarine mud, in: N.P. Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, vol. 7, Gulf, Huston, TX, 1988 (Chapter 5). [3] R.I. Tanner, Engineering Rheology, Clarendon, Oxford, 1985. [4] R.F. Dressler, Mathematical solution of the problem of roll waves in inclined open channels, Commun. Pure Appl. Math. 2 (1949) 149–194. [5] B. Zanuttigh, A. Lamberti, Roll waves simulation using shallow water equations and weighted average flux method, J. Hydraul. Res. 40 (2002) 610–622. [6] L.P. Thomas, B.M. Marino, P.F. Linden, Lock-release inertial gravity currents over a thick porous layer, J. Fluid Mech. 503 (2004) 299–319.
[7] C. Cenedese, J.A. Whitehead, T.A. Ascarelli, M. Ohiwa, A dense current flowing down a sloping bottom in a rotating fluid, J. Phys. Oceanogr. 34 (2004) 188–203. [8] G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30 (1967) 197–207. [9] J.P. Pascal, Linear stability of fluid flow down a porous inclined plane, J. Phys. D: Appl. Phys. 32 (1999) 417–422. [10] A.R. Rao, M. Mishra, Peristaltic transport of a power-law fluid in a porous tube, J. Non-Newtonian Fluid Mech. 121 (2004) 163–174. [11] H. Pascal, Rheological behaviour effect of non-Newtonian fluids on steady and unsteady flow through a porous medium, Int. J. Numer. Anal. Meth. Geomechan. 7 (1983) 298–303. [12] R.H. Christopher, S. Middleman, Power-law flow through a packed tube, I & EC Fund. 4 (1965) 422–426. [13] P.Y. Julien, D.M. Hartley, Formation of roll waves in laminar sheet flow, J. Hydraul. Res. 24 (1986) 5–17. [14] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm, J. Comput. Phys. 146 (1998) 346–365.