- Email: [email protected]
Two-dimensional numerical model of two-layer shallow water equations for confluence simulation
Advances in Water Resources 29 (2006) 1608–1617 www.elsevier.com/locate/advwatres
Two-dimensional numerical model of two-layer shallow water equations for confluence simulation Su-Chin Chen a, Szu-Hsien Peng b
b,*
a Department of Soil and Water Conservation, National Chung-Hsing University, Taichung, Taiwan Department of Space Design, Chienkuo Technology University, No. 1, Chieh Shou N. Rd., Changhua, Taiwan
Received 25 April 2005; received in revised form 19 October 2005; accepted 2 December 2005 Available online 7 February 2006
Abstract This study presents a finite-volume explicit method to solve 2D two-layer shallow water equations. This numerical model is intended to describe two-layer shallow flows in which the superposed layers differ in velocity, density and rheology in a two-dimensional domain. The rheological behavior of mudflow or debris flow is called the Bingham fluid. Thus, the shear stress on rigid bed can be derived from the constitutive equation. The computational approach adopts the HLL scheme, a novel approach for the purpose of computing a Godunov flux and solving the Riemann problem approximately proposed by Harten, Lax and van Leer, as a basic building block, treats the bottom slope by lateralizing the momentum flux, and refines the scheme using the Strang splitting to manage the frictional source term. This study successfully performed 2D two-layer shallow water computations on a rigid bed. The proposed numerical model can describe the variety of depths and velocities of substances including water and mud, when the hyperconcentrated tributary flows into the main river. The analytical results in this study will be valuable for further advanced research and for designing or planning hydraulic engineering structures. 2005 Elsevier Ltd. All rights reserved. Keywords: Two-layer shallow water equations; Bingham fluid; Finite-volume method; HLL scheme
1. Introduction Taiwan is characterized by mountainous terrain over plain areas, steep topography, fragile geology, dense distribution of torrents and gullies and frequent earthquakes, combined with heavy rainfall. Floods carrying sediment frequently cause collapses, landslides and debris flows, possibly leading to the formation of landslide dams due to the enormous sediment carried by the hyperconcentrated tributary as it enters the main river. This confluent behavior of the tributary with hyperconcentrated flow (including debris flow and mud flow) into the main river is a frequent occurrence in mountainous regions. However, studies on this
*
Corresponding author. Tel.: +886 4 7111111x2259; fax: +886 4 7111126. E-mail address: [email protected] (S.-H. Peng). 0309-1708/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.12.001
subject are very limited in Taiwan, due to the complexity of the confluence mechanism and lack of field data. To study two-dimensional unsteady gradually varied flows, Fennema and Chaudhry [1,2] used implicit and explicit finite-difference methods to simulate transient free-surface flows, e.g. flows in channels with variable cross section and/or alignment and on the upstream and downstream sides of structures occupying part of the channel width. Zhao et al. [3,4] presented RBFVM-2D, a twodimensional unsteady flow model based on the finite-volume method combined with unstructured triangle and quadrilateral grids in a river basin system. An advantage of this model is that it treats the calculation of the mass and momentum flux across each side of elements as a Riemann problem, which is solved with the Osher scheme. Recently, many studies have been developed which employ the finite-volume method with second-order accuracy to solve two-dimensional shallow water flows and have good
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
results [5–7]. However, the above models can only deal with a clear-water flood. O’Brien et al. [8] proposed FLO-2D, a two-dimensional finite difference model, which simulates clear-water flood hazards, mudflows and debris flows on alluvial fans and urban floodplains, and can route channel flow using variable area cross sections, predict channel overbank discharge and simulate floodplain flow over complex topography. Liu and Huang [9] developed a numerical model utilizing the nonlinear constitutive law proposed by Julien and Lan [10] and fluid conservation laws, and used it in the field to simulate the debris flow that occurred in Shenmu village. Additionally, Tsai and Shieh [11] studied the morphological similarity of debris-flow fans based on experimental and numerical results. Previous studies provide significant background for this study, but few studies of confluence problems due to different layers in velocity, density and rheology are still very limited. Therefore, this study proposes a 2D two-layer shallow water computation model based on a 1D model [12] and successfully simulates the confluence phenomenon qualitatively as compared with experimental results. In the first, the present paper introduced the governing equations and numerical method. This study also presents a new expression of bottom friction for Bingham fluidity in a two-layer shallow water model. Secondly, the numerical model simulated 2D partial breach problems for clear water and mudflow in top layer and bottom layer, respectively. Finally, using this two-layer shallow water numerical model computed the confluence phenomenon induced by mudflow tributary entering the main channel.
2. Governing equations The flows in Fig. 1 feature a shallow water layer on top of a shallow mud layer, in turn flowing over a rigid bottom. The two flowing layers are assumed to have distinct but constant densities q1, q2, with independent horizontal velocities u1, u2 which are uniform over the corresponding depths h1 and h2. Layer boundaries have elevations z0, z1 = z0 + h1 and z2 = z0 + h1 + h2, and are treated as sharp interfaces across which no mass transfer takes place. Momentum exchange occurs through shear stresses s0 and s1, which are applied along the interfaces and depend on the rheology of the shearing fluid. The governing equations of 2D two-layer transient freesurface flows can be developed from one-dimensional shallow water equations. Application of mass conservation and balance of momentum to the different layers yields the following governing equations [13]: oh1 o o þ ðh1 u1 Þ þ ðh1 v1 Þ ¼ 0; ox ox ot oh2 o o þ ðh2 u2 Þ þ ðh2 v2 Þ ¼ 0; ox ox ot
ð1Þ ð2Þ
1609
Fig. 1. Idealized two-layer flow structure assumed in the present work: (a) definition sketch; (b) occurrence of ‘‘contact points’’ at locations where one of the layers approaches zero depth.
o o 1 2 o 2 ðh1 u1 Þ þ h1 u1 þ gh1 þ ðh1 u1 v1 Þ ot ox 2 oy o q s1x s0x z 0 þ 2 h2 ¼ þ gh1 ; ox q1 q1 o o o 1 ðh1 v1 Þ þ ðh1 u1 v1 Þ þ h1 v21 þ gh21 ot ox oy 2 o q s1y s0y z 0 þ 2 h2 ¼ ; þ gh1 oy q1 q1 o o 1 o ðh2 u2 Þ þ h2 u22 þ gh22 þ ðh2 u2 v2 Þ ot ox 2 oy o s1x ðz0 þ h1 Þ ¼ ; ox q2 o o o 1 ðh2 v2 Þ þ ðh2 u2 v2 Þ þ h2 v22 þ gh22 ot ox oy 2 þ gh2
þ gh2
o s1y ðz0 þ h1 Þ ¼ ; oy q2
ð3Þ
ð4Þ
ð5Þ
ð6Þ
where horizontal distances x, y and time t are the three independent variables, and where g is the gravitational acceleration. The usual hydrostatic pressure assumption is adopted throughout both layers and along the interfaces. Subscripts 1 and 2 represent the bottom and top layer respectively. The water layer (on the top) has velocity components (u2, v2) in x and y directions of the horizontal plane, while the slurry layer (on the bottom) has corresponding velocity components (u1, v1). The parameter sab denotes the b = x or y components of shear stresses s1
1610
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
and s0 exerted on the interface and bed. The shear stresses s1b are achieved by specifying the shear stress functions: qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð7Þ s1x ¼ q2 C f u2r þ v2r ur ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffi s1y ¼ q2 C f u2r þ v2r vr ; ð8Þ where ur = u2 u1 and vr = v2 v1 are the two components of the velocity jump across the interface between the upper and lower layers. These two functions are Che´zy-type relations featuring friction coefficient pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Cf. However, the relationship between the V 1 ¼ u21 þ v21 and the shear stresses s0 are derived from the constitutive relation of the Bingham fluid model for two-layer flows (see Appendix A) 2 h1 s0 smin V1 ¼ ð2s0 þ smin 3s1 Þ; ð9Þ 6K s0 s1 where sminqisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the yield stress, K is the dynamic viscosity, 2 and s1 ¼ s1x þ s21y . To determine the incipient motion of a static mud layer, substitute u1 = 0 and v1 = 0 into Eqs. (3) and (4), and lead to o/ o/ þ s1y s21 2q1 gh1 s1x ox oy " 2 # 2 o/ o/ 2 þ ðq1 gh1 Þ þ ð10Þ > s2min ; ox oy where / ¼ z0 þ h1 þ qq21 h2 . Conversely, the lower mud layer freezes if its velocity V1 reaches zero. Because of the two-layer governing equations, four waves propagate upstream and downstream at speeds k1, k2, k3 and k4 in both x and y directions respectively as follows: ð11Þ S min 6 k1 6 k2 6 k3 6 k4 6 S max ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where S min ¼ffi minðu1 ;u2 Þ gðh1 þ h2 Þ, S max ¼ maxðu1 ; u2 Þþ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p gðh þ h2 Þffi in the x direction and S min ¼ minðv 1 1 ; v2 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi gðh1 þ h2 Þ; S max ¼ maxðv1 ; v2 Þ þ gðh1 þ h2 Þ in the y direction. Harten, Lax and van Leer [14] ignore the details of the wave structure, and retain only the waves in the extreme left and right directions propagating at speeds SL and SR. Namely, applying the concept of HLL scheme to two-layer model, it is the most important to determine the speeds SL and SR in the extreme left and right directions without calculating all propagate waves. 3. Computational scheme The proposed computational scheme follows that of Young and Capart [15]. In two horizontal dimensions, the governing equations can be cast in the following vector form oU oFx oFy oU oU þ þ Gy ¼ H; þ þ Gx ot ox oy ox oy
ð12Þ T
where U ¼ ½ h1 h2 h1 u1 h1 v1 h2 u2 h2 v2 is the vector of conservation variables, Fx(U) and Fy(U) are flux vectors;
Gx(U) and Gy(U) are the Jacobian matrices of the nonconservative products, and H(U) is the vector of source terms. A finite-volume approach is adopted to discretise Eq. (12) on an orthogonal grid composed of square cells. Operator splitting is used to decompose each time update into separate substeps. First, the hyperbolic lefthand-side of it is split into its two directional components [15]: oU oFx oU þ ¼ 0; þ Gx ot ox ox oU oFy oU þ ¼ 0: þ Gy ot oy oy
ð13Þ ð14Þ
Each component is integrated independently using the lateralized version of the Harten–Lax–Van Leer (LHLL) scheme developed by Fraccarollo et al. [16], an explicit scheme which is subject to the Courant time step constraint and has excellent shock-capturing properties. The above statements are summarized from the documents [13,15]. Having dealt with the hyperbolic operator, the remaining source vector is decomposed into separate components [15,17] oU ¼ Hf 1 ; ot oU ¼ Hf 0 ; ot
ð15Þ ð16Þ
linked respectively with the friction along the interface between the upper and lower layers, and the friction along the bottom bed. Each component is integrated independently using an implicit scheme, which guarantees unconditional stability. The components of the x- and y-directions in Eq. (15) are given respectively by o s1x ð u1 Þ ¼ ; ot q 1 h1 o s1y ðv1 Þ ¼ ; ot q 1 h1 o s1x ðu2 Þ ¼ ; ot q 2 h2 o s1y ðv2 Þ ¼ ; ot q 2 h2
ð17Þ ð18Þ ð19Þ ð20Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where s1x ¼ q2 C f u2r þ v2r ur , s1y ¼ q2 C f u2r þ v2r vr , ur = u2 u1, and vr = v2 v1. With reference to Saurel and Abgrall [17], subtracting Eqs. (17) and (18) from Eqs. (19) and (20), respectively, and leads to: o 1 1 ðu2 u1 Þ ¼ s1x þ ; ð21Þ ot q2 h2 q1 h1 o 1 1 ðv2 v1 Þ ¼ s1y þ : ð22Þ ot q2 h2 q1 h1 By contrast, adding Eqs. (17) and (18) into Eqs. (19) and (20), respectively: o ðq h1 u1 þ q2 h2 u2 Þ ¼ 0; ot 1 o ðq h1 v1 þ q2 h2 v2 Þ ¼ 0; ot 1
ð23Þ ð24Þ
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1611
implying that q1h1u1 + q2h2u2 and q1h1v1 + q2h2v2 are invariant parameters. Eqs. (21) and (22) are then combined to obtain o q h 1 þ q2 h 2 V r ¼ js1 j 1 ; ð25Þ ot q1 h1 q2 h2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where V r ¼ ðu2 u1 Þ2 þ ðv2 v1 Þ2 , and the negative sign implies that the shear stress js1j acts on the opposite direction to Vr. Then, an implicit backward Euler scheme is adopted to substitute js1 j ¼ q2 C f V 2r for Eq. (25), and the solution can be derived from solving an algebraic equation. Therefore, the velocity components u0r and v0r of the computed V 0r can be obtained by u0r ¼ V 0r Vurr and v0r ¼ V 0r Vvrr . Thus, the velocity components u2, v2, u1 and v1 including upper and lower layers can be solved by considering Eqs. (23) and (24), u0r and v0r . The final calculation is Eq. (16) for the source term of the bottom friction in the lower layer. Considering only s0x and s0y, then o s0x u1 ¼ ; ot q1 h1 o s0y v1 ¼ : ot q1 h1
ð26Þ ð27Þ
Similarly, the above equations are combined to obtain: o 1 V 1 ¼ s0 ; ð28Þ ot q1 h1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where s0 ¼ s20x þ s20y . Applying the implicit backward Euler scheme and substituting Eq. (9) for the above equation generates the resultant V 01 using the Newton–Raphson method. Therefore, the velocity components u01 and v01 can be derived from u01 ¼ V 01 Vu11 and v01 ¼ V 01 Vv11 . Meanwhile, the velocity components u2 and v2 are assumed to be invariant in this step.
Fig. 2. Results computed by this study at t = 7.2 s: (a) water surface profile; (b) velocity vectors.
4. Numerical modeling results 4.1. 2D partial breach modeling for top layer The computational domain consists of a 200 m long and 200 m wide channel. The nonsymmetrical breach or sluice gates are 75 m wide, and the structure or dam is 10 m thick in the flow direction. This problem assumes that the dam fails instantaneously and that the sluice gates have opened instantly [18]. The grid has 40 · 40 elements, generating an individual mesh size of 5 m · 5 m. A horizontal channel with the friction loss Cf = 0.12, and initial conditions with a tailwater/reservoir ht/hr = 0.5, were used in the simulation (the parameters in bottom layer were all zeros). Fig. 2 plots the calculated results including the water surface profile and velocity vectors. To validate the proposed model, the RBFVM-2D model was used to simulate the same case. Figs. 3 and 4 compare the computed results. Fig. 3 compares the computed trans-
Fig. 3. Comparison of transverse profiles at different locations (thick line computed by this study; dashed line computed by RBFVM-2D model).
verse water-surface profiles at x = 80 m (upstream of the breach), x = 100 m (inside the breach) and x = 120 m (downstream of the breach). Fig. 4 illustrates the time
1612
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
4.2. Deposition of mud flow on 2D plane
Fig. 4. Comparison of variation of computed flow depth with time (thick line computed by this study; dashed line computed by RBFVM-2D model).
To test the numerical computation of the bottom layer, the computational region and initial state (see Fig. 6) was defined according to Pastor et al. [19]. The original case is the flow of liquefied tailings from Gypsum Tailings Impoundment in East Texas, which occurred in 1966 [20]. Properties of the tailings were taken from [20]: the yield stress smin = 103 Pa, viscosity K = 50 Pa s and tailing density q1 = 1400 kg/m3. Similarly, the parameters in the top layer were all set to zero. Fig. 7 shows the numerical results of the mud surface. Fig. 8 shows the mud surface profiles along the vertical plane passing through the centerline of the breach. Although the depths are smaller than those documented in Pastor et al. [19], the results agree well with the observed inundation distance (300 m) and freezing time (60–120 s) values reported by Jeyapalan et al. [20]. Fig. 9 illustrates the difference between these two sets of computed results, in which the average error of maximum mud depth was less than 5.5%.
variation of flow depth at (x = 100 m, y = 75 m), i.e. inside the breach, and at (x = 120 m, y = 75 m), i.e. downstream of the breach. The difference between the computed results in the two models is determined by the friction coefficient. In the RBFVM-2D model, the Manning’s coefficient is set to n = 0.15, but Cf = 0.12 in the proposed model. When n = 0.15, Cf ranges from 0.14 to 0.11 when the hydraulic radius is in the range 4.1–7.9 m and approximate water depths are in the range 5–10 m. The difference between RBFVM-2D model and the proposed model is insignificant while the water depths lie between 6.5 m and 8 m, because the friction coefficient Cf = 0.12 is very close to n = 0.15. Fig. 5 compares the longitudinal water-surface profile at t = 7.1 s in this study with values from other studies clearly demonstrating the results computed by this study agree well with other numerical model predictions.
Fig. 5. Comparison of computed water-surface profiles for dam breach.
Fig. 6. Deposition modeling of mudflow: (a) computational region; (b) initial state.
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1613
Fig. 7. Results of mud surface profile computed by this study: (a) t = 30 s; (b) t = 60 s; (c) t = 90 s; (d) t = 120 s.
Fig. 8. Longitudinal profiles of free surface along centerline at t = 30, 60, 90, and 120 s (the surface computed by Pastor et al. is shown as a dashed line, and the surface computed by the current study as a solid line).
4.3. Numerical simulation of confluence The 2D numerical model developed from the two-layer shallow water equations was used to simulate the confluence phenomenon. The computational region derived from
Fig. 9. Comparison between hmax values from this study and Pastor et al. [19].
the experiments [21] included a main channel (1 m · 0.2 m) and a branch (0.8 m · 0.1 m). The main channel is horizontal, and the branch with a slope of 10–30% joins to the main channel at x = 0.4–0.5 m. The grid size was set to 0.01 m · 0.01 m, and the initial conditions were: (1) water depth h2 = 0.038–0.05 m, velocity u2 = 0.0088–1.0 m/s,
1614
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
and v2 = 0 in the main channel, and (2) mud surface elevation z1 = 0.1655–0.2865 m before a gate at the middle in the branch (see Table 1). Parameters were set as q1 = 1900 kg/m3, smin = 42.5 Pa, K = 9.55 Pa s and Cf = 0.01. Figs. 10 and 11 illustrate the computational results (cases A1 and A2) in Table 1. The Courant–Friedrichs–Levy condition was set to Cr = 0.3 to ensure convergence of the numerical computation. Initially, the parameters smin,
Table 1 Initial conditions of numerical simulation on confluence No.
Branch
Main channel
Slope (%)
z1
h2
u2
Fr
A1 A2 B1 B2 B3
30 30 10 10 10
0.2865 0.2865 0.1655 0.1655 0.1655
0.038 0.043 0.038 0.04 0.05
0.0088 0.0506 0.0088 0.5 1.0
0.014 0.078 0.014 0.798 1.428
Fig. 11. Experimental data (top) and numerical results (bottom) concerning submerged alluvial fan: (a) contour; (b) topography (where h2 = 0.043 m, u2 = 0.0506 m/s, z1 = 0.2865 m, t = 5 s).
Fig. 10. Experimental data (top) and numerical results (bottom) concerning the submerged alluvial fan: (a) contour; (b) topography (where h2 = 0.038 m, u2 = 0.0088 m/s, z1 = 0.2865 m, t = 5 s).
K and Cf (see Fig. 10) were calibrated to simulate the submerged fan (case A1). Then the case A2 was verified based on the same parameters, and the numerical results were compared with experimental data (see Fig. 11). Agreement between the computational results and experimental data is acceptable in morphological features. The formation of submerged alluvial fan induced by tributary inflow under the water surface of the main channel was computed and described in this numerical simulation. The upstream inflows of the main channel are subcritical (cases A1, A2 and B1 in Table 1). Fig. 12 shows the deposition of alluvial fan after confluence behavior of tributary and main flow (cases B1, B2 and B3 in which B3 represents supercritical). When the mudflow in the branch enters the main channel, the surge is formed by the intrusion of the tributary, and the alluvial fan develops gradually under the water surface. Additionally, the computed velocity vectors (see Fig. 13) of the water and mud layers concerning confluence are close to the observations of experiments in qualitative analysis. Initially, two vortices formed and
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1615
Fig. 12. Simulation of confluence at t = 5 s: (a) case B1; (b) case B2; (c) case B3 (main flow direction from left to right; mud released by gate in tributary).
Fig. 13. Velocity vectors of confluence (case B1) in top layer (a) and (b), bottom layer (c) and (d), and tributary region (e) and (f) with different time. (a) t = 0.7 s, (b) t = 3.5 s, (c) t = 0.7 s, (d) t = 3.5 s, (e) t = 0.7 s and (f) t = 3.5 s.
1616
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
nomenon using two-layer shallow water equations, and the numerical results can capture some major characteristics of confluence phenomenon with observation. Nevertheless, a number of challenges remain. From a numerical point of view, extension to second order accuracy would be desirable to better resolve fragile features such as weak shocks and contact discontinuities. Meanwhile, more precise and complete experiments need still to be continued. Acknowledgements The work is supported by the National Science Council of the Republic of China under grant no. NSC 93-2811-B005-022. Special thanks are given to Dr. H. Capart, Associate Professor of Civil Engineering at National Taiwan University, for making several valuable suggestions. Fig. 14. Side view of alluvial fans under water surface for case B1.
developed gradually, and then the vortices dispersed and moved downstream with the main channel flow [22]. Simulations were performed to model the confluence phenomena with different main channel flow conditions. Fig. 14 indicates the simulation result of Case B1, and Cases B1, B2 and B3 show that the Froude numbers of upstream inflow in main channel were about 0.01, 0.8 and 1.4, respectively. Due to the increasing flow intensity of the main flow, the water surface on the confluence rose and the alluvial fan scale decreased. The side view of the alluvial fans (see Fig. 14), clearly show that the formation of alluvial fan induced by mud and debris tributary inflow was subject to the flow intensity of the main flow as well as the tributary conditions. Additionally, the main channel flow was compressed by the alluvial fan and was concentrated at the opposite side of this fan, while the flow velocity in the concentrated region (rose simultaneously in the numerical modeling of case B1 (see Fig. 13b)). Similarly, cases B2 and B3 generated the same results even though their Froude numbers of upstream inflow in the main channel increased.
Appendix A Consider the superposed layers which differ in velocity, density and rheology flow down on an inclined slope h (see Fig. 15). The rheology of the bottom layer is described with the Bingham model and the top layer is regarded as the clear-water. Thus, the velocity of the top layer is assumed to be uniform according to a shallow water description widely used in pure water hydrodynamics. The flowing bottom layer is subdivided into (1) an upper plug-flow region (H < y < h1), moving rigidly at velocity u = up, and (2) a lower shear-flow region, where u(y) varies from zero to up as y ranges from 0 to H. If both the hypothetical shear stress s1 and the resisting stress s0 are distributed in straight lines as shown in Fig. 15, then the shear stress s(y) in the shear-flow region is expressed as following (modified from [23]): sðyÞ ¼ g sin h½cv ðqs q2 Þðh1 yÞ þ q2 ðh1 þ h2 yÞ; ðA:1Þ where g is the acceleration due to gravity, cv is the sediment concentration by volume, and qs is the density of the y
5. Conclusions Field and laboratory observations indicate that valley topography and crossflow interactions exert a very significant influence on confluence events. Horizontal two-dimensional simulations could contribute to a better understanding of these effects. This study successfully described a two-dimensional model for two-layer shallow water computation. Applying this model to the confluence formed by a mud tributary flowing into the main river, the surfaces and velocity vectors of each layer, including the mud and water layer, can be obtained by numerical simulation. Previous investigations often did not simulate numerically the confluence of mudflow and clear-water. In this study, the main contribution is to simulate the phe-
y h2
u2 ρ2
τ1 up
ρ1
u(y)
F low
τ1
h1
τ(y)
H
θ
Flow
τ0
Fig. 15. Shear stress and velocity distribution of Bingham fluid for twolayer flows.
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
particle. Substitute q1 = cvqs + (1 cv)q2 for Eq. (A.1) to yield sðyÞ ¼ g sin h½q2 h2 þ q1 ðh1 yÞ:
ðA:2Þ
Assuming that the bottom layer behaves as a Bingham fluid, its rheology is described by the relationship s ¼ smin þ K
du ; dy
ðA:3Þ
whereby beyond a certain yield stress smin the shear stress increases linearly with the strain rate, with a corresponding dynamic viscosity K. In the shear-flow region, integrating Eq. (A.3) and substituting u(y = 0) = 0 yields 1 1 2 ðs1 smin þ q1 gh1 sin hÞy q1 g sin hy ; uðyÞ ¼ K 2 ðA:4Þ where s1 = q2gh2 sin h. Eq. (A.4) can fulfill the boundary condition s(y = h1) = s1 automatically. The depth of the shear-flow region H is the locus y where s = smin, given by smin s1 s0 smin H ¼ h1 ¼ ; ðA:5Þ q1 g sin h q1 g sin h where s0 ¼ s1 þ q1 gh1 sin h. The velocity up of the rigid-like plug is obtained by substituting y = H in Eq. (A.4): 2
up ¼
ðs0 smin Þ : 2Kq1 g sin h
ðA:6Þ
The mean velocity u1 in the bottom layer can thus be obtained from Z H Z h1 1 u1 ¼ uðyÞ dy þ up dy ; ðA:7Þ h1 0 H in which results (A.4)–(A.6) can be substituted to yield 2 h1 s0 smin u1 ¼ ð2s0 þ smin 3s1 Þ: ðA:8Þ 6K s0 s1 Let n ¼ sjsmin and 1 ¼ jsjs10 jj, then 0j 2 h1 s 0 1 n s1 2 þ n 31 u1 ¼ ; 6K 1 1 s1
ðA:9Þ
where both n and 1 lie in the interval [0, 1]. If only a single layer exists, i.e. 1 = 0, then u1 ¼
h1 s 0 ð1 nÞ2 ð2 þ nÞ: 6K
ðA:10Þ
The result is the same as that given in Pastor et al. [19]. References [1] Fennema RJ, Chaudhry MH. Implicit methods for two-dimensional unsteady free-surface flows. J Hydr Res 1989;27(3):321–32.
1617
[2] Fennema RJ, Chaudhry MH. Explicit methods for 2-D transient freesurface flows. J Hydr Engrg 1990;116(8):1013–34. [3] Zhao DH, Shen HW, Tabios III GQ, Lai JS, Tan WY. Finite-volume two-dimensional unsteady-flow model for river basins. J Hydr Engrg 1994;120(7):863–83. [4] Zhao DH, Tabios III GQ, Shen, HW. RBFVM-2D Model, River basin two-dimensional flow model using finite volume method: program documentation. University of California, Berkeley, California, USA, 1995. [5] Tseng MH, Chu CR. Two-dimensional shallow water flows simulation using TVD-MacCormack schemes. J Hydr Res 2000;38(2): 123–31. [6] Gottardi G, Venutelli M. Central scheme for two-dimensional dambreak flow simulation. Adv Water Resour 2004;27:259–68. [7] Lin GF, Lai JS, Guo WD. Finite-volume component-wise TVD schemes for 2D shallow water equations. Adv Water Resour 2003;26:861–73. [8] O’Brien JS, Julien PY, Fullerton WT. Two-dimensional water flood and mudflow simulation. J Hydr Engrg 1993;119(2):244–61. [9] Liu KF, Huang MC. The numerical simulation of debris flows and its application in Shen-Mu village. J Chin Soil Water Conserv 2002;33:215–21 [in Chinese]. [10] Julien PY, Lan Y. Rheology of hyperconcentrations. J Hydr Engrg 1991;117(3):346–53. [11] Tsai YF, Shieh CL. Experimental and numerical studies on the morphological similarity of debris-flow fans. J Chinese Inst Engrg 1997;20(6):629–42. [12] Chen SC, Peng SH, Capart H. Two-layer shallow water computation of mud flow intrusion into quiescent water. J Hydraul Res, accepted for publication. [13] Capart H, Young DL. Two-layer shallow water computations of torrential geomorphic flows. In: Bousmar, Zech, editors, Proceedings of the international conference on fluvial hydraulics (River Flow 2002), Louvain, Belgium, 2002. p. 1003–12. [14] Harten A, Lax PD, van Leer B. On upstream difference and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev 1983;25:35–61. [15] Young DL, Capart H. Development and validation of a geomorphic flood routing framework. Research Report, Council of Agriculture Executive Yuan, Taiwan, ROC, 2001. [16] Fraccarollo L, Capart H, Zech Y. A Godunov method for the computation of erosional shallow water transients. Int J Numer Meth Fluids 2003;41:951–76. [17] Saurel R, Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows. J Comput Phys 1999;150: 425–67. [18] Chaudhry MH. Open-channel flow. Englewood Cliffs: Prentice-Hall; 1993. [19] Pastor M, Quecedo M, Gonza´lez E, Herreros MI, Ferna´dez Merodo JA, Mira P. Simple approximation to bottom friction for Bingham fluid depth integrated models. J Hydr Engrg 2004;130(2):149–55. [20] Jeyapalan JK, Duncan JM, Seed HB. Investigation of flow failures of tailings dams. J Geotech Engrg 1983;109(2):172–89. [21] Peng SH. Geomorphic change induced by the flow of mud and debris surges into river channels. PhD thesis, Department of Soil and Water Conservation, National Chung-Hsing University, Taiwan, ROC, 2004. [22] Chen SC, Peng SH, Capart H. Morphology of alluvial fans formed by hyperconcentrated tributaries. In: Greco, Carravetta, Della Morte, editors. Proceedings of the 2nd international conference on fluvial hydraulics (River Flow 2004), Naples, Italy, 2004. p. 1095–102. [23] Takahashi T. Debris flow. Rotterdam: Balkema; 1991.
Two-dimensional numerical model of two-layer shallow water equations for confluence simulation Su-Chin Chen a, Szu-Hsien Peng b
b,*
a Department of Soil and Water Conservation, National Chung-Hsing University, Taichung, Taiwan Department of Space Design, Chienkuo Technology University, No. 1, Chieh Shou N. Rd., Changhua, Taiwan
Received 25 April 2005; received in revised form 19 October 2005; accepted 2 December 2005 Available online 7 February 2006
Abstract This study presents a finite-volume explicit method to solve 2D two-layer shallow water equations. This numerical model is intended to describe two-layer shallow flows in which the superposed layers differ in velocity, density and rheology in a two-dimensional domain. The rheological behavior of mudflow or debris flow is called the Bingham fluid. Thus, the shear stress on rigid bed can be derived from the constitutive equation. The computational approach adopts the HLL scheme, a novel approach for the purpose of computing a Godunov flux and solving the Riemann problem approximately proposed by Harten, Lax and van Leer, as a basic building block, treats the bottom slope by lateralizing the momentum flux, and refines the scheme using the Strang splitting to manage the frictional source term. This study successfully performed 2D two-layer shallow water computations on a rigid bed. The proposed numerical model can describe the variety of depths and velocities of substances including water and mud, when the hyperconcentrated tributary flows into the main river. The analytical results in this study will be valuable for further advanced research and for designing or planning hydraulic engineering structures. 2005 Elsevier Ltd. All rights reserved. Keywords: Two-layer shallow water equations; Bingham fluid; Finite-volume method; HLL scheme
1. Introduction Taiwan is characterized by mountainous terrain over plain areas, steep topography, fragile geology, dense distribution of torrents and gullies and frequent earthquakes, combined with heavy rainfall. Floods carrying sediment frequently cause collapses, landslides and debris flows, possibly leading to the formation of landslide dams due to the enormous sediment carried by the hyperconcentrated tributary as it enters the main river. This confluent behavior of the tributary with hyperconcentrated flow (including debris flow and mud flow) into the main river is a frequent occurrence in mountainous regions. However, studies on this
*
Corresponding author. Tel.: +886 4 7111111x2259; fax: +886 4 7111126. E-mail address: [email protected] (S.-H. Peng). 0309-1708/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.12.001
subject are very limited in Taiwan, due to the complexity of the confluence mechanism and lack of field data. To study two-dimensional unsteady gradually varied flows, Fennema and Chaudhry [1,2] used implicit and explicit finite-difference methods to simulate transient free-surface flows, e.g. flows in channels with variable cross section and/or alignment and on the upstream and downstream sides of structures occupying part of the channel width. Zhao et al. [3,4] presented RBFVM-2D, a twodimensional unsteady flow model based on the finite-volume method combined with unstructured triangle and quadrilateral grids in a river basin system. An advantage of this model is that it treats the calculation of the mass and momentum flux across each side of elements as a Riemann problem, which is solved with the Osher scheme. Recently, many studies have been developed which employ the finite-volume method with second-order accuracy to solve two-dimensional shallow water flows and have good
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
results [5–7]. However, the above models can only deal with a clear-water flood. O’Brien et al. [8] proposed FLO-2D, a two-dimensional finite difference model, which simulates clear-water flood hazards, mudflows and debris flows on alluvial fans and urban floodplains, and can route channel flow using variable area cross sections, predict channel overbank discharge and simulate floodplain flow over complex topography. Liu and Huang [9] developed a numerical model utilizing the nonlinear constitutive law proposed by Julien and Lan [10] and fluid conservation laws, and used it in the field to simulate the debris flow that occurred in Shenmu village. Additionally, Tsai and Shieh [11] studied the morphological similarity of debris-flow fans based on experimental and numerical results. Previous studies provide significant background for this study, but few studies of confluence problems due to different layers in velocity, density and rheology are still very limited. Therefore, this study proposes a 2D two-layer shallow water computation model based on a 1D model [12] and successfully simulates the confluence phenomenon qualitatively as compared with experimental results. In the first, the present paper introduced the governing equations and numerical method. This study also presents a new expression of bottom friction for Bingham fluidity in a two-layer shallow water model. Secondly, the numerical model simulated 2D partial breach problems for clear water and mudflow in top layer and bottom layer, respectively. Finally, using this two-layer shallow water numerical model computed the confluence phenomenon induced by mudflow tributary entering the main channel.
2. Governing equations The flows in Fig. 1 feature a shallow water layer on top of a shallow mud layer, in turn flowing over a rigid bottom. The two flowing layers are assumed to have distinct but constant densities q1, q2, with independent horizontal velocities u1, u2 which are uniform over the corresponding depths h1 and h2. Layer boundaries have elevations z0, z1 = z0 + h1 and z2 = z0 + h1 + h2, and are treated as sharp interfaces across which no mass transfer takes place. Momentum exchange occurs through shear stresses s0 and s1, which are applied along the interfaces and depend on the rheology of the shearing fluid. The governing equations of 2D two-layer transient freesurface flows can be developed from one-dimensional shallow water equations. Application of mass conservation and balance of momentum to the different layers yields the following governing equations [13]: oh1 o o þ ðh1 u1 Þ þ ðh1 v1 Þ ¼ 0; ox ox ot oh2 o o þ ðh2 u2 Þ þ ðh2 v2 Þ ¼ 0; ox ox ot
ð1Þ ð2Þ
1609
Fig. 1. Idealized two-layer flow structure assumed in the present work: (a) definition sketch; (b) occurrence of ‘‘contact points’’ at locations where one of the layers approaches zero depth.
o o 1 2 o 2 ðh1 u1 Þ þ h1 u1 þ gh1 þ ðh1 u1 v1 Þ ot ox 2 oy o q s1x s0x z 0 þ 2 h2 ¼ þ gh1 ; ox q1 q1 o o o 1 ðh1 v1 Þ þ ðh1 u1 v1 Þ þ h1 v21 þ gh21 ot ox oy 2 o q s1y s0y z 0 þ 2 h2 ¼ ; þ gh1 oy q1 q1 o o 1 o ðh2 u2 Þ þ h2 u22 þ gh22 þ ðh2 u2 v2 Þ ot ox 2 oy o s1x ðz0 þ h1 Þ ¼ ; ox q2 o o o 1 ðh2 v2 Þ þ ðh2 u2 v2 Þ þ h2 v22 þ gh22 ot ox oy 2 þ gh2
þ gh2
o s1y ðz0 þ h1 Þ ¼ ; oy q2
ð3Þ
ð4Þ
ð5Þ
ð6Þ
where horizontal distances x, y and time t are the three independent variables, and where g is the gravitational acceleration. The usual hydrostatic pressure assumption is adopted throughout both layers and along the interfaces. Subscripts 1 and 2 represent the bottom and top layer respectively. The water layer (on the top) has velocity components (u2, v2) in x and y directions of the horizontal plane, while the slurry layer (on the bottom) has corresponding velocity components (u1, v1). The parameter sab denotes the b = x or y components of shear stresses s1
1610
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
and s0 exerted on the interface and bed. The shear stresses s1b are achieved by specifying the shear stress functions: qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð7Þ s1x ¼ q2 C f u2r þ v2r ur ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffi s1y ¼ q2 C f u2r þ v2r vr ; ð8Þ where ur = u2 u1 and vr = v2 v1 are the two components of the velocity jump across the interface between the upper and lower layers. These two functions are Che´zy-type relations featuring friction coefficient pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Cf. However, the relationship between the V 1 ¼ u21 þ v21 and the shear stresses s0 are derived from the constitutive relation of the Bingham fluid model for two-layer flows (see Appendix A) 2 h1 s0 smin V1 ¼ ð2s0 þ smin 3s1 Þ; ð9Þ 6K s0 s1 where sminqisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the yield stress, K is the dynamic viscosity, 2 and s1 ¼ s1x þ s21y . To determine the incipient motion of a static mud layer, substitute u1 = 0 and v1 = 0 into Eqs. (3) and (4), and lead to o/ o/ þ s1y s21 2q1 gh1 s1x ox oy " 2 # 2 o/ o/ 2 þ ðq1 gh1 Þ þ ð10Þ > s2min ; ox oy where / ¼ z0 þ h1 þ qq21 h2 . Conversely, the lower mud layer freezes if its velocity V1 reaches zero. Because of the two-layer governing equations, four waves propagate upstream and downstream at speeds k1, k2, k3 and k4 in both x and y directions respectively as follows: ð11Þ S min 6 k1 6 k2 6 k3 6 k4 6 S max ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where S min ¼ffi minðu1 ;u2 Þ gðh1 þ h2 Þ, S max ¼ maxðu1 ; u2 Þþ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p gðh þ h2 Þffi in the x direction and S min ¼ minðv 1 1 ; v2 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi gðh1 þ h2 Þ; S max ¼ maxðv1 ; v2 Þ þ gðh1 þ h2 Þ in the y direction. Harten, Lax and van Leer [14] ignore the details of the wave structure, and retain only the waves in the extreme left and right directions propagating at speeds SL and SR. Namely, applying the concept of HLL scheme to two-layer model, it is the most important to determine the speeds SL and SR in the extreme left and right directions without calculating all propagate waves. 3. Computational scheme The proposed computational scheme follows that of Young and Capart [15]. In two horizontal dimensions, the governing equations can be cast in the following vector form oU oFx oFy oU oU þ þ Gy ¼ H; þ þ Gx ot ox oy ox oy
ð12Þ T
where U ¼ ½ h1 h2 h1 u1 h1 v1 h2 u2 h2 v2 is the vector of conservation variables, Fx(U) and Fy(U) are flux vectors;
Gx(U) and Gy(U) are the Jacobian matrices of the nonconservative products, and H(U) is the vector of source terms. A finite-volume approach is adopted to discretise Eq. (12) on an orthogonal grid composed of square cells. Operator splitting is used to decompose each time update into separate substeps. First, the hyperbolic lefthand-side of it is split into its two directional components [15]: oU oFx oU þ ¼ 0; þ Gx ot ox ox oU oFy oU þ ¼ 0: þ Gy ot oy oy
ð13Þ ð14Þ
Each component is integrated independently using the lateralized version of the Harten–Lax–Van Leer (LHLL) scheme developed by Fraccarollo et al. [16], an explicit scheme which is subject to the Courant time step constraint and has excellent shock-capturing properties. The above statements are summarized from the documents [13,15]. Having dealt with the hyperbolic operator, the remaining source vector is decomposed into separate components [15,17] oU ¼ Hf 1 ; ot oU ¼ Hf 0 ; ot
ð15Þ ð16Þ
linked respectively with the friction along the interface between the upper and lower layers, and the friction along the bottom bed. Each component is integrated independently using an implicit scheme, which guarantees unconditional stability. The components of the x- and y-directions in Eq. (15) are given respectively by o s1x ð u1 Þ ¼ ; ot q 1 h1 o s1y ðv1 Þ ¼ ; ot q 1 h1 o s1x ðu2 Þ ¼ ; ot q 2 h2 o s1y ðv2 Þ ¼ ; ot q 2 h2
ð17Þ ð18Þ ð19Þ ð20Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where s1x ¼ q2 C f u2r þ v2r ur , s1y ¼ q2 C f u2r þ v2r vr , ur = u2 u1, and vr = v2 v1. With reference to Saurel and Abgrall [17], subtracting Eqs. (17) and (18) from Eqs. (19) and (20), respectively, and leads to: o 1 1 ðu2 u1 Þ ¼ s1x þ ; ð21Þ ot q2 h2 q1 h1 o 1 1 ðv2 v1 Þ ¼ s1y þ : ð22Þ ot q2 h2 q1 h1 By contrast, adding Eqs. (17) and (18) into Eqs. (19) and (20), respectively: o ðq h1 u1 þ q2 h2 u2 Þ ¼ 0; ot 1 o ðq h1 v1 þ q2 h2 v2 Þ ¼ 0; ot 1
ð23Þ ð24Þ
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1611
implying that q1h1u1 + q2h2u2 and q1h1v1 + q2h2v2 are invariant parameters. Eqs. (21) and (22) are then combined to obtain o q h 1 þ q2 h 2 V r ¼ js1 j 1 ; ð25Þ ot q1 h1 q2 h2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where V r ¼ ðu2 u1 Þ2 þ ðv2 v1 Þ2 , and the negative sign implies that the shear stress js1j acts on the opposite direction to Vr. Then, an implicit backward Euler scheme is adopted to substitute js1 j ¼ q2 C f V 2r for Eq. (25), and the solution can be derived from solving an algebraic equation. Therefore, the velocity components u0r and v0r of the computed V 0r can be obtained by u0r ¼ V 0r Vurr and v0r ¼ V 0r Vvrr . Thus, the velocity components u2, v2, u1 and v1 including upper and lower layers can be solved by considering Eqs. (23) and (24), u0r and v0r . The final calculation is Eq. (16) for the source term of the bottom friction in the lower layer. Considering only s0x and s0y, then o s0x u1 ¼ ; ot q1 h1 o s0y v1 ¼ : ot q1 h1
ð26Þ ð27Þ
Similarly, the above equations are combined to obtain: o 1 V 1 ¼ s0 ; ð28Þ ot q1 h1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where s0 ¼ s20x þ s20y . Applying the implicit backward Euler scheme and substituting Eq. (9) for the above equation generates the resultant V 01 using the Newton–Raphson method. Therefore, the velocity components u01 and v01 can be derived from u01 ¼ V 01 Vu11 and v01 ¼ V 01 Vv11 . Meanwhile, the velocity components u2 and v2 are assumed to be invariant in this step.
Fig. 2. Results computed by this study at t = 7.2 s: (a) water surface profile; (b) velocity vectors.
4. Numerical modeling results 4.1. 2D partial breach modeling for top layer The computational domain consists of a 200 m long and 200 m wide channel. The nonsymmetrical breach or sluice gates are 75 m wide, and the structure or dam is 10 m thick in the flow direction. This problem assumes that the dam fails instantaneously and that the sluice gates have opened instantly [18]. The grid has 40 · 40 elements, generating an individual mesh size of 5 m · 5 m. A horizontal channel with the friction loss Cf = 0.12, and initial conditions with a tailwater/reservoir ht/hr = 0.5, were used in the simulation (the parameters in bottom layer were all zeros). Fig. 2 plots the calculated results including the water surface profile and velocity vectors. To validate the proposed model, the RBFVM-2D model was used to simulate the same case. Figs. 3 and 4 compare the computed results. Fig. 3 compares the computed trans-
Fig. 3. Comparison of transverse profiles at different locations (thick line computed by this study; dashed line computed by RBFVM-2D model).
verse water-surface profiles at x = 80 m (upstream of the breach), x = 100 m (inside the breach) and x = 120 m (downstream of the breach). Fig. 4 illustrates the time
1612
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
4.2. Deposition of mud flow on 2D plane
Fig. 4. Comparison of variation of computed flow depth with time (thick line computed by this study; dashed line computed by RBFVM-2D model).
To test the numerical computation of the bottom layer, the computational region and initial state (see Fig. 6) was defined according to Pastor et al. [19]. The original case is the flow of liquefied tailings from Gypsum Tailings Impoundment in East Texas, which occurred in 1966 [20]. Properties of the tailings were taken from [20]: the yield stress smin = 103 Pa, viscosity K = 50 Pa s and tailing density q1 = 1400 kg/m3. Similarly, the parameters in the top layer were all set to zero. Fig. 7 shows the numerical results of the mud surface. Fig. 8 shows the mud surface profiles along the vertical plane passing through the centerline of the breach. Although the depths are smaller than those documented in Pastor et al. [19], the results agree well with the observed inundation distance (300 m) and freezing time (60–120 s) values reported by Jeyapalan et al. [20]. Fig. 9 illustrates the difference between these two sets of computed results, in which the average error of maximum mud depth was less than 5.5%.
variation of flow depth at (x = 100 m, y = 75 m), i.e. inside the breach, and at (x = 120 m, y = 75 m), i.e. downstream of the breach. The difference between the computed results in the two models is determined by the friction coefficient. In the RBFVM-2D model, the Manning’s coefficient is set to n = 0.15, but Cf = 0.12 in the proposed model. When n = 0.15, Cf ranges from 0.14 to 0.11 when the hydraulic radius is in the range 4.1–7.9 m and approximate water depths are in the range 5–10 m. The difference between RBFVM-2D model and the proposed model is insignificant while the water depths lie between 6.5 m and 8 m, because the friction coefficient Cf = 0.12 is very close to n = 0.15. Fig. 5 compares the longitudinal water-surface profile at t = 7.1 s in this study with values from other studies clearly demonstrating the results computed by this study agree well with other numerical model predictions.
Fig. 5. Comparison of computed water-surface profiles for dam breach.
Fig. 6. Deposition modeling of mudflow: (a) computational region; (b) initial state.
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1613
Fig. 7. Results of mud surface profile computed by this study: (a) t = 30 s; (b) t = 60 s; (c) t = 90 s; (d) t = 120 s.
Fig. 8. Longitudinal profiles of free surface along centerline at t = 30, 60, 90, and 120 s (the surface computed by Pastor et al. is shown as a dashed line, and the surface computed by the current study as a solid line).
4.3. Numerical simulation of confluence The 2D numerical model developed from the two-layer shallow water equations was used to simulate the confluence phenomenon. The computational region derived from
Fig. 9. Comparison between hmax values from this study and Pastor et al. [19].
the experiments [21] included a main channel (1 m · 0.2 m) and a branch (0.8 m · 0.1 m). The main channel is horizontal, and the branch with a slope of 10–30% joins to the main channel at x = 0.4–0.5 m. The grid size was set to 0.01 m · 0.01 m, and the initial conditions were: (1) water depth h2 = 0.038–0.05 m, velocity u2 = 0.0088–1.0 m/s,
1614
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
and v2 = 0 in the main channel, and (2) mud surface elevation z1 = 0.1655–0.2865 m before a gate at the middle in the branch (see Table 1). Parameters were set as q1 = 1900 kg/m3, smin = 42.5 Pa, K = 9.55 Pa s and Cf = 0.01. Figs. 10 and 11 illustrate the computational results (cases A1 and A2) in Table 1. The Courant–Friedrichs–Levy condition was set to Cr = 0.3 to ensure convergence of the numerical computation. Initially, the parameters smin,
Table 1 Initial conditions of numerical simulation on confluence No.
Branch
Main channel
Slope (%)
z1
h2
u2
Fr
A1 A2 B1 B2 B3
30 30 10 10 10
0.2865 0.2865 0.1655 0.1655 0.1655
0.038 0.043 0.038 0.04 0.05
0.0088 0.0506 0.0088 0.5 1.0
0.014 0.078 0.014 0.798 1.428
Fig. 11. Experimental data (top) and numerical results (bottom) concerning submerged alluvial fan: (a) contour; (b) topography (where h2 = 0.043 m, u2 = 0.0506 m/s, z1 = 0.2865 m, t = 5 s).
Fig. 10. Experimental data (top) and numerical results (bottom) concerning the submerged alluvial fan: (a) contour; (b) topography (where h2 = 0.038 m, u2 = 0.0088 m/s, z1 = 0.2865 m, t = 5 s).
K and Cf (see Fig. 10) were calibrated to simulate the submerged fan (case A1). Then the case A2 was verified based on the same parameters, and the numerical results were compared with experimental data (see Fig. 11). Agreement between the computational results and experimental data is acceptable in morphological features. The formation of submerged alluvial fan induced by tributary inflow under the water surface of the main channel was computed and described in this numerical simulation. The upstream inflows of the main channel are subcritical (cases A1, A2 and B1 in Table 1). Fig. 12 shows the deposition of alluvial fan after confluence behavior of tributary and main flow (cases B1, B2 and B3 in which B3 represents supercritical). When the mudflow in the branch enters the main channel, the surge is formed by the intrusion of the tributary, and the alluvial fan develops gradually under the water surface. Additionally, the computed velocity vectors (see Fig. 13) of the water and mud layers concerning confluence are close to the observations of experiments in qualitative analysis. Initially, two vortices formed and
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
1615
Fig. 12. Simulation of confluence at t = 5 s: (a) case B1; (b) case B2; (c) case B3 (main flow direction from left to right; mud released by gate in tributary).
Fig. 13. Velocity vectors of confluence (case B1) in top layer (a) and (b), bottom layer (c) and (d), and tributary region (e) and (f) with different time. (a) t = 0.7 s, (b) t = 3.5 s, (c) t = 0.7 s, (d) t = 3.5 s, (e) t = 0.7 s and (f) t = 3.5 s.
1616
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
nomenon using two-layer shallow water equations, and the numerical results can capture some major characteristics of confluence phenomenon with observation. Nevertheless, a number of challenges remain. From a numerical point of view, extension to second order accuracy would be desirable to better resolve fragile features such as weak shocks and contact discontinuities. Meanwhile, more precise and complete experiments need still to be continued. Acknowledgements The work is supported by the National Science Council of the Republic of China under grant no. NSC 93-2811-B005-022. Special thanks are given to Dr. H. Capart, Associate Professor of Civil Engineering at National Taiwan University, for making several valuable suggestions. Fig. 14. Side view of alluvial fans under water surface for case B1.
developed gradually, and then the vortices dispersed and moved downstream with the main channel flow [22]. Simulations were performed to model the confluence phenomena with different main channel flow conditions. Fig. 14 indicates the simulation result of Case B1, and Cases B1, B2 and B3 show that the Froude numbers of upstream inflow in main channel were about 0.01, 0.8 and 1.4, respectively. Due to the increasing flow intensity of the main flow, the water surface on the confluence rose and the alluvial fan scale decreased. The side view of the alluvial fans (see Fig. 14), clearly show that the formation of alluvial fan induced by mud and debris tributary inflow was subject to the flow intensity of the main flow as well as the tributary conditions. Additionally, the main channel flow was compressed by the alluvial fan and was concentrated at the opposite side of this fan, while the flow velocity in the concentrated region (rose simultaneously in the numerical modeling of case B1 (see Fig. 13b)). Similarly, cases B2 and B3 generated the same results even though their Froude numbers of upstream inflow in the main channel increased.
Appendix A Consider the superposed layers which differ in velocity, density and rheology flow down on an inclined slope h (see Fig. 15). The rheology of the bottom layer is described with the Bingham model and the top layer is regarded as the clear-water. Thus, the velocity of the top layer is assumed to be uniform according to a shallow water description widely used in pure water hydrodynamics. The flowing bottom layer is subdivided into (1) an upper plug-flow region (H < y < h1), moving rigidly at velocity u = up, and (2) a lower shear-flow region, where u(y) varies from zero to up as y ranges from 0 to H. If both the hypothetical shear stress s1 and the resisting stress s0 are distributed in straight lines as shown in Fig. 15, then the shear stress s(y) in the shear-flow region is expressed as following (modified from [23]): sðyÞ ¼ g sin h½cv ðqs q2 Þðh1 yÞ þ q2 ðh1 þ h2 yÞ; ðA:1Þ where g is the acceleration due to gravity, cv is the sediment concentration by volume, and qs is the density of the y
5. Conclusions Field and laboratory observations indicate that valley topography and crossflow interactions exert a very significant influence on confluence events. Horizontal two-dimensional simulations could contribute to a better understanding of these effects. This study successfully described a two-dimensional model for two-layer shallow water computation. Applying this model to the confluence formed by a mud tributary flowing into the main river, the surfaces and velocity vectors of each layer, including the mud and water layer, can be obtained by numerical simulation. Previous investigations often did not simulate numerically the confluence of mudflow and clear-water. In this study, the main contribution is to simulate the phe-
y h2
u2 ρ2
τ1 up
ρ1
u(y)
F low
τ1
h1
τ(y)
H
θ
Flow
τ0
Fig. 15. Shear stress and velocity distribution of Bingham fluid for twolayer flows.
S.-C. Chen, S.-H. Peng / Advances in Water Resources 29 (2006) 1608–1617
particle. Substitute q1 = cvqs + (1 cv)q2 for Eq. (A.1) to yield sðyÞ ¼ g sin h½q2 h2 þ q1 ðh1 yÞ:
ðA:2Þ
Assuming that the bottom layer behaves as a Bingham fluid, its rheology is described by the relationship s ¼ smin þ K
du ; dy
ðA:3Þ
whereby beyond a certain yield stress smin the shear stress increases linearly with the strain rate, with a corresponding dynamic viscosity K. In the shear-flow region, integrating Eq. (A.3) and substituting u(y = 0) = 0 yields 1 1 2 ðs1 smin þ q1 gh1 sin hÞy q1 g sin hy ; uðyÞ ¼ K 2 ðA:4Þ where s1 = q2gh2 sin h. Eq. (A.4) can fulfill the boundary condition s(y = h1) = s1 automatically. The depth of the shear-flow region H is the locus y where s = smin, given by smin s1 s0 smin H ¼ h1 ¼ ; ðA:5Þ q1 g sin h q1 g sin h where s0 ¼ s1 þ q1 gh1 sin h. The velocity up of the rigid-like plug is obtained by substituting y = H in Eq. (A.4): 2
up ¼
ðs0 smin Þ : 2Kq1 g sin h
ðA:6Þ
The mean velocity u1 in the bottom layer can thus be obtained from Z H Z h1 1 u1 ¼ uðyÞ dy þ up dy ; ðA:7Þ h1 0 H in which results (A.4)–(A.6) can be substituted to yield 2 h1 s0 smin u1 ¼ ð2s0 þ smin 3s1 Þ: ðA:8Þ 6K s0 s1 Let n ¼ sjsmin and 1 ¼ jsjs10 jj, then 0j 2 h1 s 0 1 n s1 2 þ n 31 u1 ¼ ; 6K 1 1 s1
ðA:9Þ
where both n and 1 lie in the interval [0, 1]. If only a single layer exists, i.e. 1 = 0, then u1 ¼
h1 s 0 ð1 nÞ2 ð2 þ nÞ: 6K
ðA:10Þ
The result is the same as that given in Pastor et al. [19]. References [1] Fennema RJ, Chaudhry MH. Implicit methods for two-dimensional unsteady free-surface flows. J Hydr Res 1989;27(3):321–32.
1617
[2] Fennema RJ, Chaudhry MH. Explicit methods for 2-D transient freesurface flows. J Hydr Engrg 1990;116(8):1013–34. [3] Zhao DH, Shen HW, Tabios III GQ, Lai JS, Tan WY. Finite-volume two-dimensional unsteady-flow model for river basins. J Hydr Engrg 1994;120(7):863–83. [4] Zhao DH, Tabios III GQ, Shen, HW. RBFVM-2D Model, River basin two-dimensional flow model using finite volume method: program documentation. University of California, Berkeley, California, USA, 1995. [5] Tseng MH, Chu CR. Two-dimensional shallow water flows simulation using TVD-MacCormack schemes. J Hydr Res 2000;38(2): 123–31. [6] Gottardi G, Venutelli M. Central scheme for two-dimensional dambreak flow simulation. Adv Water Resour 2004;27:259–68. [7] Lin GF, Lai JS, Guo WD. Finite-volume component-wise TVD schemes for 2D shallow water equations. Adv Water Resour 2003;26:861–73. [8] O’Brien JS, Julien PY, Fullerton WT. Two-dimensional water flood and mudflow simulation. J Hydr Engrg 1993;119(2):244–61. [9] Liu KF, Huang MC. The numerical simulation of debris flows and its application in Shen-Mu village. J Chin Soil Water Conserv 2002;33:215–21 [in Chinese]. [10] Julien PY, Lan Y. Rheology of hyperconcentrations. J Hydr Engrg 1991;117(3):346–53. [11] Tsai YF, Shieh CL. Experimental and numerical studies on the morphological similarity of debris-flow fans. J Chinese Inst Engrg 1997;20(6):629–42. [12] Chen SC, Peng SH, Capart H. Two-layer shallow water computation of mud flow intrusion into quiescent water. J Hydraul Res, accepted for publication. [13] Capart H, Young DL. Two-layer shallow water computations of torrential geomorphic flows. In: Bousmar, Zech, editors, Proceedings of the international conference on fluvial hydraulics (River Flow 2002), Louvain, Belgium, 2002. p. 1003–12. [14] Harten A, Lax PD, van Leer B. On upstream difference and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev 1983;25:35–61. [15] Young DL, Capart H. Development and validation of a geomorphic flood routing framework. Research Report, Council of Agriculture Executive Yuan, Taiwan, ROC, 2001. [16] Fraccarollo L, Capart H, Zech Y. A Godunov method for the computation of erosional shallow water transients. Int J Numer Meth Fluids 2003;41:951–76. [17] Saurel R, Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows. J Comput Phys 1999;150: 425–67. [18] Chaudhry MH. Open-channel flow. Englewood Cliffs: Prentice-Hall; 1993. [19] Pastor M, Quecedo M, Gonza´lez E, Herreros MI, Ferna´dez Merodo JA, Mira P. Simple approximation to bottom friction for Bingham fluid depth integrated models. J Hydr Engrg 2004;130(2):149–55. [20] Jeyapalan JK, Duncan JM, Seed HB. Investigation of flow failures of tailings dams. J Geotech Engrg 1983;109(2):172–89. [21] Peng SH. Geomorphic change induced by the flow of mud and debris surges into river channels. PhD thesis, Department of Soil and Water Conservation, National Chung-Hsing University, Taiwan, ROC, 2004. [22] Chen SC, Peng SH, Capart H. Morphology of alluvial fans formed by hyperconcentrated tributaries. In: Greco, Carravetta, Della Morte, editors. Proceedings of the 2nd international conference on fluvial hydraulics (River Flow 2004), Naples, Italy, 2004. p. 1095–102. [23] Takahashi T. Debris flow. Rotterdam: Balkema; 1991.