- Email: [email protected]
A model for two-layer moderate Reynolds number flow down an incline
International Journal of Non-Linear Mechanics 36 (2001) 977}985
A model for two-layer moderate Reynolds number #ow down an incline J.P. Pascal* Applied Mathematics Institute, Department of Mathematical Sciences, University of Alberta, Edmonton, Alberta, Canada T6G 2G1 Received 15 May 2000; accepted 12 June 2000
Abstract In this paper we investigate the gravity current resulting from the spread of a #uid of non-Newtonian power law rheology down an incline of constant slope under a less dense shallow Newtonian #uid layer. The theoretical model employed in this investigation accounts for the e!ect of both inertial and viscous forces, and is suited for moderately high Reynolds number #ows down gentle slopes. The governing equations are obtained from the unsteady equations of motion by applying von KaH rmaH n's momentum integral method. Results are obtained by a well-established numerical scheme for systems of non-linear hyperbolic conservation laws. For a particular case analytical results are obtained by employing an asymptotic matching approach. A good agreement is obtained between the numerical and analytical results. The e!ects of the thickness of the ambient layer, Reynolds number, and rheology on the gravity are discussed. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Gravity current; Two-layer model; Free surface #ow; Non-Newtonian power law #uid
1. Introduction Gravity-driven #ows consisting of a #uid #owing into another #uid of di!erent density are a common occurrence in many natural and industrial situations [1]. If the di!erence in the densities of the two #uids is &small', the resulting #ow is referred to as a gravity current. In this case the #ow is more complex due to the coupling of the dynamics of the #ow of the two #uids. It turns out that many gravity currents occurring in natural phenomena and man-made situations are high Reynolds
* Fax: #1-780-492-6826. E-mail address: [email protected] (J.P. Pascal).
number #ows. For this reason most theoretical investigations concentrate on the balance between gravitational and inertial forces, neglecting the effect of viscous forces. Such an inviscid-#uid theory investigation was recently carried out by D'Alessio et al. [2] to study the spread of #uid along a horizontal bottom under a shallow layer of less-dense #uid with a free surface. They utilized a two-layer #uid model based on shallow-water theory to consider the coupling of the dynamics of the two layers. The inviscid-#uid theory however, proves to be inadequate for gravity currents down an incline of constant slope. This was pointed out by Montgomery and Moodie [3] who proposed using a drag condition at the front of the gravity current to account for the e!ect of the viscous forces and
0020-7462/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 6 2 ( 0 0 ) 0 0 0 6 3 - 9
978
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
balance the gravitational acceleration down the slope. Viscous forces must also be considered when studying mud #ows, which fall in the range of moderately high Reynolds number #ows. Fluid mud consists of a cohesive suspension containing clay and other mineral and organic particles in water. These particles are advected by rivers to estuarine and costal regions where they settle on the bottom, thus generating a layer of mud which slowly spreads along the coast [4,5]. Numerous experiments have been carried out with the aim of determining the rheological properties of #uid mud in laminar #ow. It has been found that in high particle concentrations #uid mud exhibits strong non-Newtonian properties [4,5]. In particular, it has been reported that 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 [5,6] which is expressed as "2(D)D, where is the deviatoric stress tensor and D is the shear rate tensor, with the viscosity being given by (D)"[2tr(D)]L\ , where is the consistency and n is the power-law index of the #uid. For shear thinning #uids n is between 0 and 1, and the case with n"1 corresponding to a Newtonian #uid with being the dynamic viscosity. In the literature one "nds several theoretical investigations of mud #ows down an incline. Hunt [7] uses a matched asymptotic perturbation technique to obtain "rst-order solutions to the depth averaged single-layer equations for the #ow resulting from a point source. This approach is also utilized by Huang and Garcia [8] to investigate #ows of Herschel}Bulkley #uids which are #uids with an yield stress and a power-law relation between shear rate and shear stress. Ng and Mei [5] utilize the power-law model in their investigation of the hydrodynamic stability of a layer of mud #owing down an incline. They use the depth-averaged single-layer equations to investigate the formation and evolution of roll waves on the surface of the mud layer.
In the investigations mentioned above the e!ect of the ambient medium on the #ow was neglected. This is consistent with the applications under consideration which consist of mud #ows under the atmosphere. However, in the case of submarine mud #ows the #ow constitutes a gravity current spreading under a less-dense ambient #uid. If the ambient layer is shallow we expect a signi"cant coupling of the dynamics of the two #uids. In this paper we derive and solve a two-layer model that incorporates the dynamics of the gravity current and of the ambient layer and includes the e!ect of inertial and viscous forces. The model assumes the gravity current #uid to be of powerlaw rheology and the ambient #uid to be a much less-viscous Newtonian #uid. Thus, the model is applicable to mud #ows in water. We solve the depth-averaged equations of motion numerically and for a particular case compare the results with the analytical ones obtained by the matchedasymptotic expansion used by Hunt [7] and Huang and Garcia [8].
2. The governing equations Consider a two-dimensional laminar gravity current of density intruding into a shallow ambient #uid of lesser density along an incline at an angle with respect to the horizontal. As is illustrated in Fig. 1 we de"ne 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"(u, w)2 and p, respectively, and use the subscript a to denote quantities pertaining to the ambient layer. We expect the characteristic longitudinal length of the #ow, ¸ to be much larger than the
Fig. 1. A sketch of the two-layer model used in this paper.
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
979
characteristic depth, H, with the depth varying slowly in time and in the longitudinal direction. As a result we neglect terms in the equations of motion that are of O(), where "H/¸ is the aspect ratio. The equations of motion are then expressed as (see [5])
obtained from the continuity of force requirement. Assuming that the medium above the ambient layer does not exert a force on it we have the conditions
u w # "0, x z
where z"h (x, t) is the position of the free surface of the ambient layer. At the interface between the two #uid layer we have
(1)
u u p u #u #w "! #g sin x z x t
(2)
(4)
1 p 0"! !g cos , z
(5) (6)
where g is the acceleration due to gravity, is the consistency of the power law #uid, and is the dynamic viscosity of the Newtonian ambient #uid. In obtaining these approximations we assumed that 1 1 ¸g sin , &O() and &O(1), Re Re ; where ; is the characteristic longitudinal velocity of the #ow and the Reynolds numbers are de"ned as ;\LHL Re"
(8)
u "0 at z"h(x, t), z
(9)
where z"h(x, t) is the position of the interface. In obtaining the condition (9) we made the assumption that Re &O() Re
u u p u #u #w "! #g sin x z x t u # , z
(7)
and
(3)
u w # "0, z x
u "0 at z"h (x, t), z
p "p at z"h(x, t)
u L\u , # z z z
1 p !g cos , 0"! z
p "0,
;H and Re " .
Since is small these assumptions are appropriate for moderately high Reynolds number #ows down gentle slopes. Boundary conditions at the top of the ambient layer and at the interface between the two #uids are
which corresponds to cases when the gravity current is much more viscous than the ambient #uid. At the interface between the #uid we also have the condition of continuity of velocity which is stated as u "u, w "w at z"h(x, t).
(10)
Finally, at the bottom of the gravity current the appropriate boundary condition is u"w"0 at z"0.
(11)
The kinematic conditions at the surface of the ambient #uid and at the interface between the two #uids can be expressed as h h w " #u t x h h w" #u t x
at z"h (x, t),
at z"h(x, t).
Integrating Eqs. (3) and (6) and using the boundary conditions (7) and (8) we "nd the expressions for the total pressures to be p "g cos (h !z), p"!g cos z#g cos h #g cos (! )h.
980
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Inserting these expressions into the x-momentum equations (2) and (5) we obtain u u h u #u #w "g sin !gcos x z x t #(g!g)cos
h x
u L\u # , z z z
(12)
u u u #u #w x z t h u "g sin !g cos # , x z
(13)
(14)
h !z h !h
3 2n#1 # u ! u , 2 2(n#1)
q (h !h)# "0, x t
(17)
q q 1 # # g cos h t x h 2
2n#1 L q L ! , n h
(18)
q q(h !h) 2 qq # ! t h x 5 5 h
6q 1 # # g cos h 5(h !h) 2
!3 (15)
where u and u are the depth-averaged velocities given by
1 F u(z) dz u " h and
(16)
h "g sin (h !h)#g cos h x
and 3 2n#1 u (z)" u !u 2 n#1
h q # "0, t x
h "g sin h#(g!g)cos h x
where g"g(! )/ is the reduced gravity. For a uniform steady #ow these equations together with the boundary conditions (7)}(11) can be solved to give 2n#1 h!z > L u 1! u(z)" n#1 h
Depth-integrating the continuity equations (1) and (4) and the momentum equations (12) and (13) in the two layers we obtain
F 1 u " u (z) dz. h !h F Since we are dealing with a small aspect ratio #ow which varies slowly in the longitudinal direction, we can apply von KaH rmaH n's momentum integral method. We thus assume that the relationship between the velocity and the #uid depth for the steady #ow is also true when these quantities are timedependent and non-uniform.
q h! q(h !h) , h(h !h)
(19)
where q"huN ,
q "(h !h)uN ,
2(2n#1) , " 3n#2
2n#1 " . n#1
We assume that the surface of the ambient layer is initially horizontal, and express its position as z"h (x, t),H#tan x# (x, t), where H is the undisturbed depth at x"0 and
(x, t) represents the displacement of the surface from its undisturbed con"guration. Since we are considering cases when the Reynolds number in the ambient layer is much larger than that in the gravity current, the displacement of this surface due to the gravity current will be proportional to the relative density di!erence between the two #uids which is given by g/g"(! )/.
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
We now nondimensionalize the system of equations (16)}(19) using a scaling that focuses on the non-linear internal gravity wave processes. To this end we introduce the nondimensional quantities (signi"ed by a tilde) according to x"¸x ,
z"Hz ,
¸ t" tI , (h, h )"H(hI , hI ), ;
" H , (q, q )";H(q , q ), where ; ;"(Hg cos ) , ¸" g sin
and "g/g.
Employing this transformation to nondimensionalize equations (16)}(19) and dropping the tildes for notational convenience we obtain h q # "0, t x
(20)
h q q # # "0, t x x
(21)
q q 1 # # h t x h 2
(22)
q(h !h) 2 qq q # ! x 5 h 5 h t
h h 1 q h! q(h !h) "h !h# ! , x R h(h !h)
(23)
where
.
3.1. Numerical method To obtain numerical solutions to the equations (20)}(23) we consider the problem of the instantaneous release of a "xed volume of power-law #uid initially at rest into a sloped channel containing a less-dense Newtonian #uid. We assume that the power-law #uid is initially contained in the region 0)x)x and has a horizontal surface, as depic ted by the dash line in Fig. 1. Then the appropriate initial conditions are
h # x for 0)x)x , 0 for x'x , q(x, 0)"q (x, 0)" (x, 0)"0, x*0, h(x, 0)"
where x "1 was used in all the computations. The parameter h is the fraction of the total depth initially occupied by the heavy #uid at x"0. If we assume that a solid boundary is located at the left end of the channel then we have the boundary conditions
L n H tan 1 Re, R " Re , " " 2n#1 3 ¸
The nondimensional position of the free surface above the ambient layer is z"h (x, t)"1# x# (x, t).
t*0.
Eqs. (20)}(23) constitute a system of four "rst-order quasilinear partial di!erential equations governing the unknowns q , q, h, and . In the particular case of "1 the eigenvalues of the coe$cient matrix of the derivatives of the unknowns can be expressed as 1
" (1$5 ), 2
6q 1 # # h 5(h !h) 2
R"
3. Method of solution
q(0, t)"q (0, t)"0,
h 1 q L "h#(1!1/ )h ! , x R h
981
q
"h! , h q 6q
" ! #h !h 5h 5(h !h) which are clearly real. As a result the system of equations (20)}(23) is hyperbolic if "1. In the general case, although the eigenvalues of the coe$cient matrix can be stated explicitly, the expressions are too complex to allow for an analysis that would provide a condition for hyperbolicity in terms of the variables of the problem. However, we expect the system of equations to be hyperbolic
982
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
in which case the initial discontinuity at x"x will persist and under certain circumstances the solution will develop discontinuities in "nite time. We therefore must employ a numerical procedure that correctly captures discontinuities and e!ectively dampens spurious numerical oscillations while maintaining a high order of accuracy. An e$cient strategy for dealing with these di$culties that was proposed by Harten [9] is to modify a second-order accurate scheme by adding arti"cial viscosity to dampen numerical oscillations. Because adding the necessary arti"cial viscosity reduces the accuracy to "rst order, the arti"cial viscosity term is solution dependent in such a way as to add signi"cant arti"cial viscosity only around discontinuities where the oscillations occur. The resulting scheme then remains second order accurate where the solution is smooth and is "rst order accurate only near discontinuities. To obtain numerical solutions to our problem we use a scheme constructed by adding a selfadjusting arti"cial viscosity term to MacCormack's scheme [9]. MacCormack's is a conservative second-order scheme appropriate for systems of hyperbolic conservation laws with source terms [10]. We may thus assert that our numerical results are second order accurate on smooth regions, "rst order accurate near discontinuities, and converge to the physical weak solution of the problem. For all our numerical work we used a space increment of x"0.02. The corresponding values for the time increment that were required for stability ranged from t"0.002 to 0.005. 3.2. Analytical solution In the special case "1 the equations governing the #ow of the two layers decouple and the two equations governing the gravity current reduce to the single-layer model. These two equations can be solved analytically if beside terms of O() we also neglect terms of O(). This solution can be used as a check of the numerical results. An analytical solution for the single-layer case in connection with the Herschel}Bulkley #uid model was obtained by Huang and Garcia [8]. This model reduces to the power law one if the yield stress parameter is set to zero. Therefore, below we outline their approach
and present the solution for the power-law #uid case. To identify the terms of O() in Eqs. (20)}(23) for the case "1 we return to the dimensional form and normalize the variables by applying the transformation ¸ h"HhH, t" tH, q"H;qH, ;
x"¸xH, where
;"(Hg cos ) . We thus obtain h q # "0, t x
(24)
q q 1 # # h t x h 2
"tan h!
2n#1 L 1 q L , Re h n
(25)
where the stars have been dropped for notational convenience. Neglecting terms of O() and solving by the method of characteristics along with the initial condition
h(x, 0)"
x for 0)x)x , 0 for x'x
we obtain x"(Re tan ) Lh> Lt#h,
0)x)x (t), (26) where x (t) is the location of the front of the gravity current which is given by
(2n#1)x #h , x" 2(n#1)h with h being the depth of the gravity current at the front. The solution given by (26) is the "rst-order outer approximation in a singular perturbation expansion. However, this solution is not "rst-order accurate near the shock at the front where h/x becomes large and therefore the equations have been scaled incorrectly. An inner solution that is valid near the front can be obtained by introducing
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
983
the variable x!x " in which case h /&O(1). The inner solution h , is given implicitly by
h (Re tan ) h#1 log(h !h ) tan "h # 45 2
h (Re tan ) # h!1 [log(h #h ) 2 45 (Re tan ) h h(log h !1). !log 4]! ! 45 2
Fig. 2. Comparison between the numerical and analytical solutions with n"1, R"1, and "1, at t"10.
A composite solution h , can be obtained by substracting the common matching term h , from the sum of the inner and outer solutions, and is thus expressed as h "h#h !h . The composite solution can be rescaled so that a comparison can be made with the numerical results. This comparison is illustrated in Fig. 2.
4. Results In this section we present and interpret various numerical results for the two-layer model derived in this investigation. Of particular interest is the e!ect of the thickness of the ambient layer above the gravity current on the evolution of the shape of the current, and how this e!ect correlates with other aspects of the #ow such as the Reynolds number, the rheology of the gravity current #uid and the density di!erence of the two #uids. Our results indicate that a rear shock forms on the interface between the two #uids if the ambient layer above the gravity current is thin compared to the gravity current. The relative di!erence in the depths of the two layers is indicated by the value of the parameter h , which is a measure of the initial depth ratio at x"0. The rear shock propagates downstream with a greater speed than the front, which results in the eventual collapse of the head of the gravity current. This behavior is illustrated in
Fig. 3. Time evolution of the gravity current with n"1, h "0.8, R"10, R "1000, and "0.2.
Fig. 3 which corresponds to the case h "0.8. As the value of h is decreased the strength of the rear shock decreases and for lower values a shock does not form at all, as it can be seen in Fig. 4 where h "0.3. This is due to the lower reverse #ow speed in the ambient layer as h is decreased. The formation of the rear shock and its dependence on the depth ratio of the two layers can be determined by employing an inviscid two-layer model as it was done by D'Alessio et al. [2]. The current model however also takes into consideration the e!ect of viscous forces and allows us to determine how the Reynolds number of the #ow a!ects the formation and propagation of the rear
984
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Fig. 4. Time evolution of the gravity current with n"1, h "0.3, R"10, R "1000, and "0.2.
Fig. 5. Time evolution of the gravity current with n"1, h "0.8, R"5, R "1000, and "0.2.
shock. Our results suggest that as the Reynolds number is decreased the strength of the rear shock decreases and is eventually replaced by a wave of depression for small Reynolds numbers akin to that observed when the fraction depth of the gravity current is small. This behavior can be seen in Figs. 3, 5, and 6 which display the evolution of the gravity current with R"10, 5, 1, respectively. The parameter R is related to the Reynolds number of the #ow by R"
L n Re. 2n#1
For mud #ow of power law rheology [5] Re&O(10), and
L n 1 ) (1 2n#1 3 for 0(n)1. As a result, we considered values of R ranging from O(1) to O(10). An especially relevant aspect of a non-Newtonian #uid #ow problem is the e!ect of the rheology of the #uid on the #ow behavior. For power-law #uids the rheology is strongly related to the parameter n, the power-law index. Our results indicated that increasing the shear-thinning nature of the #uid, i.e. decreasing the value of n, results in a decrease in the rate of spreading of the gravity current after the collapse of the head of the current. This
Fig. 6. Time evolution of the gravity current with n"1, h "0.8, R"1, R "1000, and "0.2.
e!ect can be observed by comparing the evolution of the gravity current for n"1 which is shown in Fig. 5 to that for n"0.4 which is shown in Fig. 7. In the stage following the collapse of the head, the gravity current experiences a decrease in the rate of shear. As a result a #uid with a lower value of n would be thicker in this stage and thus spreads at a slower rate. The e!ect of the relative density di!erence of the two #uids is determined by varying the parameter . The results reveal that, as can be seen in Fig. 8, the rate of spreading of the denser #uid is an increasing function of . This observation is consistent with the underlying physical mechanism, since
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Fig. 7. Time evolution of the gravity current with n"0.4, h "0.8, R"5, R "1000, and "0.2.
985
#uids and is thus appropriate for various depths of the ambient layer including the case when this layer is shallow. This aspect makes the model relevant in a wide range of situations with realistic applications. The model derived in this paper also accounts for both inertial and viscous forces and is more useful than the two-layer #uid model based on the inviscid #uid theory. This is made evident by the fact that our results indicated that varying the Reynolds number has a qualitative e!ect on the shape of the gravity current. Also retaining the viscous terms in the equations of motion prevents the creation of runaway kinetic energy which results when the inviscid equations of motion are used to describe #ows down a constant slope.
References
Fig. 8. Shape of the gravity current for various values of with n"1, h "0.8, R"5, and R "1000 at t"4.
increasing the density di!erence increases the buoyancy force which drives the spreading process.
5. Concluding remarks In this paper we have presented a model for the unsteady spread of a non-Newtonian power-law #uid down an incline under a less-dense Newtonian #uid layer. The model is a two-layer #uid model that couples the dynamics of the #ow of the two
[1] J.E. Simpson, Gravity Currents: In the Environment and the Laboratory, Cambridge University Press, Cambridge, 1997. [2] S.J.D. D'Alessio, T.B. Moodie, J.P. Pascal, G.E. Swaters, Gravity currents produced by sudden release of a "xed volume of heavy #uid, Stud. Appl. Math. 96 (1996) 359}385. [3] P.J. Montgomery, T.B. Moodie, Two-layer gravity currents with topography, Stud. Appl. Math. 102 (1999) 221}266. [4] G. Vereet, J. Berlamont, Rheology and non-Newtonian behavior of sea and estuarine mud, in: N.P. Cheremisino!, (Ed.), Encyclopedia of Fluid Mechanics, Vol. 7, Gulf, Huston, TX, 1988 (Chapter 5). [5] C. Ng, C.C. Mei, Roll waves on a shallow layer of mud modelled as a power-law #uid, J. Fluid Mech. 263 (1994) 151}183. [6] R.I. Tanner, Engineering Rheology, Clarendon, Oxford, 1985. [7] B. Hunt, Newtonian #uid mechanics treatment of debris #ows and avalanches, J. Hydraul. Engng ASCE 120 (1994) 1350}1363. [8] X. Huang, M.H. Garcia, A Herschel}Bulkley model for mud #ow down a slope, J. Fluid Mech. 374 (1998) 305}333. [9] A. Harten, The arti"cial compression method for computation of shocks and contact discontinuities III. Selfadjusting hybrid schemes, Math. Comp. 32 (1978) 363}389. [10] R.J. LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with sti! source terms, J. Comp. Phys. 86 (1990) 187}210.
A model for two-layer moderate Reynolds number #ow down an incline J.P. Pascal* Applied Mathematics Institute, Department of Mathematical Sciences, University of Alberta, Edmonton, Alberta, Canada T6G 2G1 Received 15 May 2000; accepted 12 June 2000
Abstract In this paper we investigate the gravity current resulting from the spread of a #uid of non-Newtonian power law rheology down an incline of constant slope under a less dense shallow Newtonian #uid layer. The theoretical model employed in this investigation accounts for the e!ect of both inertial and viscous forces, and is suited for moderately high Reynolds number #ows down gentle slopes. The governing equations are obtained from the unsteady equations of motion by applying von KaH rmaH n's momentum integral method. Results are obtained by a well-established numerical scheme for systems of non-linear hyperbolic conservation laws. For a particular case analytical results are obtained by employing an asymptotic matching approach. A good agreement is obtained between the numerical and analytical results. The e!ects of the thickness of the ambient layer, Reynolds number, and rheology on the gravity are discussed. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Gravity current; Two-layer model; Free surface #ow; Non-Newtonian power law #uid
1. Introduction Gravity-driven #ows consisting of a #uid #owing into another #uid of di!erent density are a common occurrence in many natural and industrial situations [1]. If the di!erence in the densities of the two #uids is &small', the resulting #ow is referred to as a gravity current. In this case the #ow is more complex due to the coupling of the dynamics of the #ow of the two #uids. It turns out that many gravity currents occurring in natural phenomena and man-made situations are high Reynolds
* Fax: #1-780-492-6826. E-mail address: [email protected] (J.P. Pascal).
number #ows. For this reason most theoretical investigations concentrate on the balance between gravitational and inertial forces, neglecting the effect of viscous forces. Such an inviscid-#uid theory investigation was recently carried out by D'Alessio et al. [2] to study the spread of #uid along a horizontal bottom under a shallow layer of less-dense #uid with a free surface. They utilized a two-layer #uid model based on shallow-water theory to consider the coupling of the dynamics of the two layers. The inviscid-#uid theory however, proves to be inadequate for gravity currents down an incline of constant slope. This was pointed out by Montgomery and Moodie [3] who proposed using a drag condition at the front of the gravity current to account for the e!ect of the viscous forces and
0020-7462/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 6 2 ( 0 0 ) 0 0 0 6 3 - 9
978
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
balance the gravitational acceleration down the slope. Viscous forces must also be considered when studying mud #ows, which fall in the range of moderately high Reynolds number #ows. Fluid mud consists of a cohesive suspension containing clay and other mineral and organic particles in water. These particles are advected by rivers to estuarine and costal regions where they settle on the bottom, thus generating a layer of mud which slowly spreads along the coast [4,5]. Numerous experiments have been carried out with the aim of determining the rheological properties of #uid mud in laminar #ow. It has been found that in high particle concentrations #uid mud exhibits strong non-Newtonian properties [4,5]. In particular, it has been reported that 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 [5,6] which is expressed as "2(D)D, where is the deviatoric stress tensor and D is the shear rate tensor, with the viscosity being given by (D)"[2tr(D)]L\ , where is the consistency and n is the power-law index of the #uid. For shear thinning #uids n is between 0 and 1, and the case with n"1 corresponding to a Newtonian #uid with being the dynamic viscosity. In the literature one "nds several theoretical investigations of mud #ows down an incline. Hunt [7] uses a matched asymptotic perturbation technique to obtain "rst-order solutions to the depth averaged single-layer equations for the #ow resulting from a point source. This approach is also utilized by Huang and Garcia [8] to investigate #ows of Herschel}Bulkley #uids which are #uids with an yield stress and a power-law relation between shear rate and shear stress. Ng and Mei [5] utilize the power-law model in their investigation of the hydrodynamic stability of a layer of mud #owing down an incline. They use the depth-averaged single-layer equations to investigate the formation and evolution of roll waves on the surface of the mud layer.
In the investigations mentioned above the e!ect of the ambient medium on the #ow was neglected. This is consistent with the applications under consideration which consist of mud #ows under the atmosphere. However, in the case of submarine mud #ows the #ow constitutes a gravity current spreading under a less-dense ambient #uid. If the ambient layer is shallow we expect a signi"cant coupling of the dynamics of the two #uids. In this paper we derive and solve a two-layer model that incorporates the dynamics of the gravity current and of the ambient layer and includes the e!ect of inertial and viscous forces. The model assumes the gravity current #uid to be of powerlaw rheology and the ambient #uid to be a much less-viscous Newtonian #uid. Thus, the model is applicable to mud #ows in water. We solve the depth-averaged equations of motion numerically and for a particular case compare the results with the analytical ones obtained by the matchedasymptotic expansion used by Hunt [7] and Huang and Garcia [8].
2. The governing equations Consider a two-dimensional laminar gravity current of density intruding into a shallow ambient #uid of lesser density along an incline at an angle with respect to the horizontal. As is illustrated in Fig. 1 we de"ne 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"(u, w)2 and p, respectively, and use the subscript a to denote quantities pertaining to the ambient layer. We expect the characteristic longitudinal length of the #ow, ¸ to be much larger than the
Fig. 1. A sketch of the two-layer model used in this paper.
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
979
characteristic depth, H, with the depth varying slowly in time and in the longitudinal direction. As a result we neglect terms in the equations of motion that are of O(), where "H/¸ is the aspect ratio. The equations of motion are then expressed as (see [5])
obtained from the continuity of force requirement. Assuming that the medium above the ambient layer does not exert a force on it we have the conditions
u w # "0, x z
where z"h (x, t) is the position of the free surface of the ambient layer. At the interface between the two #uid layer we have
(1)
u u p u #u #w "! #g sin x z x t
(2)
(4)
1 p 0"! !g cos , z
(5) (6)
where g is the acceleration due to gravity, is the consistency of the power law #uid, and is the dynamic viscosity of the Newtonian ambient #uid. In obtaining these approximations we assumed that 1 1 ¸g sin , &O() and &O(1), Re Re ; where ; is the characteristic longitudinal velocity of the #ow and the Reynolds numbers are de"ned as ;\LHL Re"
(8)
u "0 at z"h(x, t), z
(9)
where z"h(x, t) is the position of the interface. In obtaining the condition (9) we made the assumption that Re &O() Re
u u p u #u #w "! #g sin x z x t u # , z
(7)
and
(3)
u w # "0, z x
u "0 at z"h (x, t), z
p "p at z"h(x, t)
u L\u , # z z z
1 p !g cos , 0"! z
p "0,
;H and Re " .
Since is small these assumptions are appropriate for moderately high Reynolds number #ows down gentle slopes. Boundary conditions at the top of the ambient layer and at the interface between the two #uids are
which corresponds to cases when the gravity current is much more viscous than the ambient #uid. At the interface between the #uid we also have the condition of continuity of velocity which is stated as u "u, w "w at z"h(x, t).
(10)
Finally, at the bottom of the gravity current the appropriate boundary condition is u"w"0 at z"0.
(11)
The kinematic conditions at the surface of the ambient #uid and at the interface between the two #uids can be expressed as h h w " #u t x h h w" #u t x
at z"h (x, t),
at z"h(x, t).
Integrating Eqs. (3) and (6) and using the boundary conditions (7) and (8) we "nd the expressions for the total pressures to be p "g cos (h !z), p"!g cos z#g cos h #g cos (! )h.
980
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Inserting these expressions into the x-momentum equations (2) and (5) we obtain u u h u #u #w "g sin !gcos x z x t #(g!g)cos
h x
u L\u # , z z z
(12)
u u u #u #w x z t h u "g sin !g cos # , x z
(13)
(14)
h !z h !h
3 2n#1 # u ! u , 2 2(n#1)
q (h !h)# "0, x t
(17)
q q 1 # # g cos h t x h 2
2n#1 L q L ! , n h
(18)
q q(h !h) 2 qq # ! t h x 5 5 h
6q 1 # # g cos h 5(h !h) 2
!3 (15)
where u and u are the depth-averaged velocities given by
1 F u(z) dz u " h and
(16)
h "g sin (h !h)#g cos h x
and 3 2n#1 u (z)" u !u 2 n#1
h q # "0, t x
h "g sin h#(g!g)cos h x
where g"g(! )/ is the reduced gravity. For a uniform steady #ow these equations together with the boundary conditions (7)}(11) can be solved to give 2n#1 h!z > L u 1! u(z)" n#1 h
Depth-integrating the continuity equations (1) and (4) and the momentum equations (12) and (13) in the two layers we obtain
F 1 u " u (z) dz. h !h F Since we are dealing with a small aspect ratio #ow which varies slowly in the longitudinal direction, we can apply von KaH rmaH n's momentum integral method. We thus assume that the relationship between the velocity and the #uid depth for the steady #ow is also true when these quantities are timedependent and non-uniform.
q h! q(h !h) , h(h !h)
(19)
where q"huN ,
q "(h !h)uN ,
2(2n#1) , " 3n#2
2n#1 " . n#1
We assume that the surface of the ambient layer is initially horizontal, and express its position as z"h (x, t),H#tan x# (x, t), where H is the undisturbed depth at x"0 and
(x, t) represents the displacement of the surface from its undisturbed con"guration. Since we are considering cases when the Reynolds number in the ambient layer is much larger than that in the gravity current, the displacement of this surface due to the gravity current will be proportional to the relative density di!erence between the two #uids which is given by g/g"(! )/.
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
We now nondimensionalize the system of equations (16)}(19) using a scaling that focuses on the non-linear internal gravity wave processes. To this end we introduce the nondimensional quantities (signi"ed by a tilde) according to x"¸x ,
z"Hz ,
¸ t" tI , (h, h )"H(hI , hI ), ;
" H , (q, q )";H(q , q ), where ; ;"(Hg cos ) , ¸" g sin
and "g/g.
Employing this transformation to nondimensionalize equations (16)}(19) and dropping the tildes for notational convenience we obtain h q # "0, t x
(20)
h q q # # "0, t x x
(21)
q q 1 # # h t x h 2
(22)
q(h !h) 2 qq q # ! x 5 h 5 h t
h h 1 q h! q(h !h) "h !h# ! , x R h(h !h)
(23)
where
.
3.1. Numerical method To obtain numerical solutions to the equations (20)}(23) we consider the problem of the instantaneous release of a "xed volume of power-law #uid initially at rest into a sloped channel containing a less-dense Newtonian #uid. We assume that the power-law #uid is initially contained in the region 0)x)x and has a horizontal surface, as depic ted by the dash line in Fig. 1. Then the appropriate initial conditions are
h # x for 0)x)x , 0 for x'x , q(x, 0)"q (x, 0)" (x, 0)"0, x*0, h(x, 0)"
where x "1 was used in all the computations. The parameter h is the fraction of the total depth initially occupied by the heavy #uid at x"0. If we assume that a solid boundary is located at the left end of the channel then we have the boundary conditions
L n H tan 1 Re, R " Re , " " 2n#1 3 ¸
The nondimensional position of the free surface above the ambient layer is z"h (x, t)"1# x# (x, t).
t*0.
Eqs. (20)}(23) constitute a system of four "rst-order quasilinear partial di!erential equations governing the unknowns q , q, h, and . In the particular case of "1 the eigenvalues of the coe$cient matrix of the derivatives of the unknowns can be expressed as 1
" (1$5 ), 2
6q 1 # # h 5(h !h) 2
R"
3. Method of solution
q(0, t)"q (0, t)"0,
h 1 q L "h#(1!1/ )h ! , x R h
981
q
"h! , h q 6q
" ! #h !h 5h 5(h !h) which are clearly real. As a result the system of equations (20)}(23) is hyperbolic if "1. In the general case, although the eigenvalues of the coe$cient matrix can be stated explicitly, the expressions are too complex to allow for an analysis that would provide a condition for hyperbolicity in terms of the variables of the problem. However, we expect the system of equations to be hyperbolic
982
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
in which case the initial discontinuity at x"x will persist and under certain circumstances the solution will develop discontinuities in "nite time. We therefore must employ a numerical procedure that correctly captures discontinuities and e!ectively dampens spurious numerical oscillations while maintaining a high order of accuracy. An e$cient strategy for dealing with these di$culties that was proposed by Harten [9] is to modify a second-order accurate scheme by adding arti"cial viscosity to dampen numerical oscillations. Because adding the necessary arti"cial viscosity reduces the accuracy to "rst order, the arti"cial viscosity term is solution dependent in such a way as to add signi"cant arti"cial viscosity only around discontinuities where the oscillations occur. The resulting scheme then remains second order accurate where the solution is smooth and is "rst order accurate only near discontinuities. To obtain numerical solutions to our problem we use a scheme constructed by adding a selfadjusting arti"cial viscosity term to MacCormack's scheme [9]. MacCormack's is a conservative second-order scheme appropriate for systems of hyperbolic conservation laws with source terms [10]. We may thus assert that our numerical results are second order accurate on smooth regions, "rst order accurate near discontinuities, and converge to the physical weak solution of the problem. For all our numerical work we used a space increment of x"0.02. The corresponding values for the time increment that were required for stability ranged from t"0.002 to 0.005. 3.2. Analytical solution In the special case "1 the equations governing the #ow of the two layers decouple and the two equations governing the gravity current reduce to the single-layer model. These two equations can be solved analytically if beside terms of O() we also neglect terms of O(). This solution can be used as a check of the numerical results. An analytical solution for the single-layer case in connection with the Herschel}Bulkley #uid model was obtained by Huang and Garcia [8]. This model reduces to the power law one if the yield stress parameter is set to zero. Therefore, below we outline their approach
and present the solution for the power-law #uid case. To identify the terms of O() in Eqs. (20)}(23) for the case "1 we return to the dimensional form and normalize the variables by applying the transformation ¸ h"HhH, t" tH, q"H;qH, ;
x"¸xH, where
;"(Hg cos ) . We thus obtain h q # "0, t x
(24)
q q 1 # # h t x h 2
"tan h!
2n#1 L 1 q L , Re h n
(25)
where the stars have been dropped for notational convenience. Neglecting terms of O() and solving by the method of characteristics along with the initial condition
h(x, 0)"
x for 0)x)x , 0 for x'x
we obtain x"(Re tan ) Lh> Lt#h,
0)x)x (t), (26) where x (t) is the location of the front of the gravity current which is given by
(2n#1)x #h , x" 2(n#1)h with h being the depth of the gravity current at the front. The solution given by (26) is the "rst-order outer approximation in a singular perturbation expansion. However, this solution is not "rst-order accurate near the shock at the front where h/x becomes large and therefore the equations have been scaled incorrectly. An inner solution that is valid near the front can be obtained by introducing
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
983
the variable x!x " in which case h /&O(1). The inner solution h , is given implicitly by
h (Re tan ) h#1 log(h !h ) tan "h # 45 2
h (Re tan ) # h!1 [log(h #h ) 2 45 (Re tan ) h h(log h !1). !log 4]! ! 45 2
Fig. 2. Comparison between the numerical and analytical solutions with n"1, R"1, and "1, at t"10.
A composite solution h , can be obtained by substracting the common matching term h , from the sum of the inner and outer solutions, and is thus expressed as h "h#h !h . The composite solution can be rescaled so that a comparison can be made with the numerical results. This comparison is illustrated in Fig. 2.
4. Results In this section we present and interpret various numerical results for the two-layer model derived in this investigation. Of particular interest is the e!ect of the thickness of the ambient layer above the gravity current on the evolution of the shape of the current, and how this e!ect correlates with other aspects of the #ow such as the Reynolds number, the rheology of the gravity current #uid and the density di!erence of the two #uids. Our results indicate that a rear shock forms on the interface between the two #uids if the ambient layer above the gravity current is thin compared to the gravity current. The relative di!erence in the depths of the two layers is indicated by the value of the parameter h , which is a measure of the initial depth ratio at x"0. The rear shock propagates downstream with a greater speed than the front, which results in the eventual collapse of the head of the gravity current. This behavior is illustrated in
Fig. 3. Time evolution of the gravity current with n"1, h "0.8, R"10, R "1000, and "0.2.
Fig. 3 which corresponds to the case h "0.8. As the value of h is decreased the strength of the rear shock decreases and for lower values a shock does not form at all, as it can be seen in Fig. 4 where h "0.3. This is due to the lower reverse #ow speed in the ambient layer as h is decreased. The formation of the rear shock and its dependence on the depth ratio of the two layers can be determined by employing an inviscid two-layer model as it was done by D'Alessio et al. [2]. The current model however also takes into consideration the e!ect of viscous forces and allows us to determine how the Reynolds number of the #ow a!ects the formation and propagation of the rear
984
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Fig. 4. Time evolution of the gravity current with n"1, h "0.3, R"10, R "1000, and "0.2.
Fig. 5. Time evolution of the gravity current with n"1, h "0.8, R"5, R "1000, and "0.2.
shock. Our results suggest that as the Reynolds number is decreased the strength of the rear shock decreases and is eventually replaced by a wave of depression for small Reynolds numbers akin to that observed when the fraction depth of the gravity current is small. This behavior can be seen in Figs. 3, 5, and 6 which display the evolution of the gravity current with R"10, 5, 1, respectively. The parameter R is related to the Reynolds number of the #ow by R"
L n Re. 2n#1
For mud #ow of power law rheology [5] Re&O(10), and
L n 1 ) (1 2n#1 3 for 0(n)1. As a result, we considered values of R ranging from O(1) to O(10). An especially relevant aspect of a non-Newtonian #uid #ow problem is the e!ect of the rheology of the #uid on the #ow behavior. For power-law #uids the rheology is strongly related to the parameter n, the power-law index. Our results indicated that increasing the shear-thinning nature of the #uid, i.e. decreasing the value of n, results in a decrease in the rate of spreading of the gravity current after the collapse of the head of the current. This
Fig. 6. Time evolution of the gravity current with n"1, h "0.8, R"1, R "1000, and "0.2.
e!ect can be observed by comparing the evolution of the gravity current for n"1 which is shown in Fig. 5 to that for n"0.4 which is shown in Fig. 7. In the stage following the collapse of the head, the gravity current experiences a decrease in the rate of shear. As a result a #uid with a lower value of n would be thicker in this stage and thus spreads at a slower rate. The e!ect of the relative density di!erence of the two #uids is determined by varying the parameter . The results reveal that, as can be seen in Fig. 8, the rate of spreading of the denser #uid is an increasing function of . This observation is consistent with the underlying physical mechanism, since
J.P. Pascal / International Journal of Non-Linear Mechanics 36 (2001) 977}985
Fig. 7. Time evolution of the gravity current with n"0.4, h "0.8, R"5, R "1000, and "0.2.
985
#uids and is thus appropriate for various depths of the ambient layer including the case when this layer is shallow. This aspect makes the model relevant in a wide range of situations with realistic applications. The model derived in this paper also accounts for both inertial and viscous forces and is more useful than the two-layer #uid model based on the inviscid #uid theory. This is made evident by the fact that our results indicated that varying the Reynolds number has a qualitative e!ect on the shape of the gravity current. Also retaining the viscous terms in the equations of motion prevents the creation of runaway kinetic energy which results when the inviscid equations of motion are used to describe #ows down a constant slope.
References
Fig. 8. Shape of the gravity current for various values of with n"1, h "0.8, R"5, and R "1000 at t"4.
increasing the density di!erence increases the buoyancy force which drives the spreading process.
5. Concluding remarks In this paper we have presented a model for the unsteady spread of a non-Newtonian power-law #uid down an incline under a less-dense Newtonian #uid layer. The model is a two-layer #uid model that couples the dynamics of the #ow of the two
[1] J.E. Simpson, Gravity Currents: In the Environment and the Laboratory, Cambridge University Press, Cambridge, 1997. [2] S.J.D. D'Alessio, T.B. Moodie, J.P. Pascal, G.E. Swaters, Gravity currents produced by sudden release of a "xed volume of heavy #uid, Stud. Appl. Math. 96 (1996) 359}385. [3] P.J. Montgomery, T.B. Moodie, Two-layer gravity currents with topography, Stud. Appl. Math. 102 (1999) 221}266. [4] G. Vereet, J. Berlamont, Rheology and non-Newtonian behavior of sea and estuarine mud, in: N.P. Cheremisino!, (Ed.), Encyclopedia of Fluid Mechanics, Vol. 7, Gulf, Huston, TX, 1988 (Chapter 5). [5] C. Ng, C.C. Mei, Roll waves on a shallow layer of mud modelled as a power-law #uid, J. Fluid Mech. 263 (1994) 151}183. [6] R.I. Tanner, Engineering Rheology, Clarendon, Oxford, 1985. [7] B. Hunt, Newtonian #uid mechanics treatment of debris #ows and avalanches, J. Hydraul. Engng ASCE 120 (1994) 1350}1363. [8] X. Huang, M.H. Garcia, A Herschel}Bulkley model for mud #ow down a slope, J. Fluid Mech. 374 (1998) 305}333. [9] A. Harten, The arti"cial compression method for computation of shocks and contact discontinuities III. Selfadjusting hybrid schemes, Math. Comp. 32 (1978) 363}389. [10] R.J. LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with sti! source terms, J. Comp. Phys. 86 (1990) 187}210.