- Email: [email protected]
The Role Fore Air Flow in Soil Slope Stability Analysis
640
2009,21(5):640-646 DOI: 10.1016/S1001-6058(08)60195-X
THE ROLE FORE AIR FLOW IN SOIL SLOPE STABILITY ANALYSIS* ZHANG Xiao-yue, ZHU Yue-ming, FANG Chun-hui College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China, E-mail: [email protected]
(Received December 3, 2008, Revised January 18, 2009)
Abstract: In order to investigate the effect of the pore air flow in the soil slope stability analysis, a water-air two-phase flow model, based on the multi-phase flow theory, is proposed and with the model, the water-phase and air-phase seepages of the soil slope in the stable seepage and rainfall situations are simulated. The soil slope safety coefficients are obtained according to the simulated pore air pressure and pore water pressure. The calculation results show that in the stable seepage situation, the influence of the pore air pressure on the soil slope stability can be neglected, while in the rainfall situation, the pore air pressure generated in the unsaturated zone will reduce the safety coefficient, and the larger the distance between the slip surface and the underground water level is, the greater the influence of the pore air pressure on the soil slope stability will be. Key words: unsaturated soil, water-air two-phase flow, rainfall, slope stability
1. Introduction The safety of the natural and artificial soil slopes is an important issue. The unsaturated soil is a four-phase mixture, composed of solid-phase, air-phase, liquid-phase and water-air interface called shrinkable film[1], so one has to deal with not only water-phase flow but also air-phase flow in the soil[2,3]. Because the pore air pressure is hard to measure, it is not taken into consideration in the general slope stability analysis. In resent years, the multi-component multi-phase flow theory has been made a great progress. For example, Petros Gaganis advanced a multi-phase model[4], Kees proposed a time and spatial discretization method for the two-phase flow model[5], Oostrom analyzed the relationship between saturation and relative seepage coefficient in porous * Project supported by the China-Austria scientific and technological cooperative project of National Science and Technology Ministry (Grant No. CN 01/2007), the Ministry of Water Resources’ 948 plan, namely technical innovation and transformation projects (Grant No. CT200515). Biography: ZHANG Xiao-yue (1984-), Female, Ph. D. Candidate
media[6], Laroche studied the relationship between capillary pressure and saturation in a multiphase model[7], and all of these make it possible to analyze the air-phase flow in soil slopes and consider the influence of pore air pressure on slope stability. This article proposes a water-air two-phase flow model to simulate the water-air seepage in a homogeneous soil slope under stable seepage and rainfall conditions. The anti-sliding safety coefficients of the slope are calculated from the results obtained by the two-phase flow simulation, and the effect of the pore air pressure on the slope stability is analyzed.
Fig.1 Phase states and their conversions
641
transformed into the following form: 2. Water-air two-phase flow model 2.1 Definitions and assumptions The seepage in an unsaturated soil is a water-air two-phase flow, and the air-phase and water-phase both involve two components: water and air[8]. In this article, the subscripts g and w represent the air-phase and water-phase, respectively, while the superscripts a and w represent the air and water components, respectively. The air-phase and water-phase are transformed to each other by evaporation and dissolution (Fig.1). Some assumptions about the model are made. The soil skeleton is undeformable, the water phase is incompressible, the temperature keeps constant, hysteresis effect in the constitutive relationship is not considered, the chemical and biological reactions are neglected, the ideal gas law is valid for the air component, the partial pressure of the water component in the air phase is equal to the saturation pressure. 2.2 The governing equation and its solution With the mass transfer between the water-phase and the air-phase being taken into consideration, the governing equation of the two-phase flow model takes the form of the mass conservation of the component N [9]:
(2)
where the vector x represents the primary variables. Expanding Eq.(2) in Taylor’s series, neglecting high-order terms, we have at the iteration step m +1 and time step k + 1 :
§ wF · F ( x k +1,m +1 ) | F ( x k +1, m ) + ¨ < ¸ © wx ¹ k +1,m ( x k +1,m +1 x k +1,m )
(3)
As F ( x k +1, m +1 ) should be equal to zero, Eq.(3) can be rewritten as:
K ( x k +1,m )u = F ( x k +1,m )
(4)
where K = wF/wx is the coefficient matrix, u = x k +1, m +1 x k +1, m is the vector of the primary variable increments, F ( x k +1, m ) is the residual term at time step k + 1 and iteration step m . The elements K ij of the matrix K are obtained by finite difference method[11]:
§ · w ¨ ¦ SD U mol ,D xDN ¸ ª krD U mol ,D xDN D © ¹ ¦ div « K I < wt PD D ¬
(gradpD UD g )@ QN = 0
F ( x) = 0
K ij = (1)
where N = w, a represent the water component and air component, respectively, D = w, g represent the water-phase and air-phase, respectively, I , S ,
U mol , xDN and t are the porosity, saturation, molar density, mole fraction of component N in phase D and the time, respectively, K = k sw P w / U w g represents the tensor of absolute permeability of the soil, k sw and kr are the saturated and relative permeability coefficients of the soil, respectively, and p , P , g and QN are the pressure, the absolute viscosity, the gravity acceleration and the source term, respectively. In the governing equation, there are primary and secondary parameters, and the secondary parameters are functions of the primary parameters, and these parameters will be discussed later on. The nonlinear governing equation can be solved by Newton-Raphson method[10]. Equation (1) is first
wFi k +1,m | ª Fi (..., x j 1 , x j + 'x j , x j +1 ,...) wx kj +1,m ¬
Fi (..., x j 1 , x j 'x j , x j +1 ,...) ¼º / 2'x j
(5)
where 'x j is a small increment, j = 1, 2,! , n , n is equal to the number of primary variables times the element number. The more detailed calculation steps of the Newton-Raphson method can be described as follows. At time step k + 1 and iteration step m , the initial value x k +1,0 is the solution at time step k , while F ( x k +1, m ) ! H , solve
K ( x k +1, m )u = F ( x k +1,m ) ,
2
and let x k +1, m +1 = x k +1, m + K u , m = m +1 , do the next iteration until F ( x k +1, m ) H , then do calculation at 2
the next time step. In the above description, H is the allowable error, K 1 is the step coefficient of the iteration. 2.3 Primary variables The phase states in a given control volume are variables in space and time, and are directly related with the primary variables. According to the Gibbsian
642
phase rule[10], the number of degrees of freedom in a multi-component multi-phase system is F = C P + 2 + ( P 1) = C +1 , where C and P are the component number and phase number, respectively. As the temperature is kept as a constant, there should be two independent primary variables. With only air phase, the primary variables are xgw and p g , with only water phase, the primary variables are xwa and p g , with both water and air phases, the
is the Henry coefficient, T is the absolute temperature. Using the idea gas law, the mole density and the density of the air phase can be expressed as
U mol , g =
(9)
U g = U mol , g ¦ xgk M N
(10)
N
primary variables are S w and p g . The algorithm for phase transformation in the program is shown in Fig.2.
p n = g Vg RT
where R = 8.314J /(mol K) is the universal gas constant, n , Vg , T and M N are the quantity of the gas, the volume occupied by n , the absolute temperature and the molecular weight of component N , respectively. The difference between the water-phase pressure and the air-phase pressure, which, respectively, act on the two sides of the water-air interface, is usually called capillary pressure[13] pc = pw pg d 0 . The capillary pressure and the relative permeability coefficient are both the functions of the saturation degree, and they can be evaluated by using the approximate formulas proposed by Van Genuchten[14]:
pc ( S w ) = P0 ( S e 1/ m 1)1/ n Fig.2 Algorithm for phase transformation in the program
2.4 Secondary variables The mole fractions satisfy the following relation:
¦N x
N D
xgw =
=1
(6)
w sat
p pg
(7)
w is the saturation pressure of the water where psat vapor. According to Henry’s law and Dalton’s law for partial pressures[12], we have
w xwa = K H ( pg psat )
(8)
where
K H = [0.8942 +1.47 exp(0.04394T )] u 10
10
(11)
m k rw ( S w ) = Se1/ 2 ª«1 1 Se1/ m º» ¬ ¼
2
(12)
krg ( S w ) = (1 Se )1/ 3 1 Se1/ m
2m
(13)
where m and n are the empirical parameters, m = 1 1/ n , P0 is the air entry value of the soil, S e = ( S w S rw ) /( S sw S rw ) is the effective saturation, Srw and S sw are the residual and maximum water saturation, respectively. 2.5 Conditions for a unique solution Thin layer elements, called boundary condition elements, are added to the model boundary. In order to keep the boundary condition fixed in the calculation process, the volumes of the thin layer elements are assumed to take a maximum value. The initial condition of the stable seepage flow calculation is as follows: the soil elements in the slope are water-saturated, with primary variables of xwa and
pg ,
xwa = 1010 and
pg
equal to the standard
atmospheric pressure. The boundary condition elements above the upstream and downstream water
643
levels are air-saturated, with primary variables of xga and p g , xga = 0.9999 and p g equal to the standard atmospheric pressure. The boundary condition elements below the upstream and downstream water levels are water-saturated, with primary variables of xwa and p g , xwa = 1010 and p g equal to the water pressure of the centroid of the boundary condition element. The initial condition of the unstable seepage flow calculation is what obtained by the stable flow calculation. There are two kinds of boundary conditions: Dirichlet boundary conditions, as the air pressure and water pressure boundary conditions, which specify initial values to the boundary condition elements, Neumann boundary conditions, as the water and air flux boundary conditions, which specify source terms to the boundary condition elements. In the rainfall condition, the source term can be obtained by:
Q w (t ) = U w Aq (t )
(14)
where A is the effective area of the soil surface in the vertical direction of the rainwater infiltration, q (t ) is the rainfall intensity.
3. Stability analysis method With the circular slip surface, according to the simplified Bishop method, it is assumed that X R X L = 0 , with X R and X L denoting the inter-slice vertical shears. With the pore pressures at the grid nodes obtained from the water-air two-phase flow calculation, the pore pressures at the midpoints of the soil slice hemlines are determined by interpolation based on the inverse isoparametric mapping in FEM[15], then the slope anti-sliding stability coefficient F , as the ratio of the shear strength to the shear stress of the slip surface, is obtained by the moment balance. According to the vertical force balance, the total normal force N i acting on the soil slice hemline i is:
Ni =
1 mD ,i
oblique length of the soil slice hemline i , c c is the effective cohesion, I c is the effective angle of internal friction, I b represents the increasing rate of shear strength with the capillary pressure, and if the soil slice hemline is in the saturated zone, I b is equal to I c , D i is the angle between the horizontal plane and the tangent of the soil slice hemline through its midpoint, Wi is the gravity of the soil slice i ,
pw,i and pg , j are the pore water and pore air pressures at the midpoint of the soil slice hemline, respectively. According to the moment balance, we have
F=
1 ¦Wi sin Di
¦ ^ccE i
Ei sin D i F
º (tan I c tan I b ) » ¼
where mD ,i = cos D i + (sin D i tan I c) / F , E i
§ tan I b · º tan I b p g ,i E i ¨ 1 ¸ » tan I c` tan I c tan I c ¹ ¼ ©
(16)
The soil slope safety coefficients can be determined by Eqs.(15) and (16). Let the negative pore water pressure pw,i and pore air pressure pg ,i be zero in the calculation, then Eqs.(15) and (16) are turned into the calculation formulas in the saturated soil slope stability analysis. 4. Examples A homogeneous soil slope is taken as an example, with upstream water level of 11 m, downstream water level of 16 m, slope angle of 27o, and an impervious bottom. Here, two slip surfaces are studied. The slip surface 1 is deep with most of the soil slice hemlines in the saturated zone and in shape of an arc with center at (16.0, 25.0) m, radius of 16.16 m, and the slip surface 2 is shallow with most of the soil slice hemlines in the unsaturated zone and also in shape of an arc with center at (20.0, 27.0) m, radius of 13.74 m. The cross section of the slope is shown in Fig.3.
(15) is the
+ ª¬ N i pw,i Ei <
i
ccE i sin D i E sin D i ª + pw , i i < «Wi F F ¬
tan I b + pg ,i
i
Fig.3 Cross section of the slope
644
The parameters of the soil used in the simulation are shown in Table 1. Table1 Material parameters Parameter
Symbol
Value
Saturated seepage coefficient
k sw (m/s)
10-6
Air entry value
P0 (kPa)
20
Residual water saturation
S rw
0.15
Maximum water saturation
S sw
1.0
Temperature
T (ć )
15.0
Empirical parameter
m
0.457
Porosity
I
0.421
Effective cohesion
cc (kPa)
10.1
Effective angle of internal friction
Ic
42.6
Increasing rate angle
Ib
35.0
Density of soil
U (kN/m3)
18.94
Density of saturated soil
Usat (kN/m3)
19.63
the stable seepage situation and after 5 h of rainfall are shown in Fig.6 and Fig.7, respectively. By comparing the calculation results obtained by the water-air two-phase flow method and the FEM, the water-air two-phase flow model can be verified.
Fig.4 Calculation results in the stable seepage situation
Using the water-air two-phase flow model, the pore pressure and the saturation of the soil slope in the stable seepage situation and the rainfall situation with rainfall intensity of 2.0×10-5 m/s and rainfall time of 5 h, are obtained (Fig.4 and Fig.5), in which the pore air pressure and the pore water pressure are the corresponding water head hg = ( pg p0 ) / U w g and hw = ( pw p0 ) / U w g , respectively, and p0 is the standard atmospheric pressure. From Figs.4 and 5, it can be seen that in the stable seepage situation, the pore air pressure is about zero in the unsaturated zone, after 5 h of rainfall, there appears a transient saturated zone in the soil surface layer, and an increased water saturation in the unsaturated zone, which reduces the negative pore water pressure. In the rain water infiltration process, not all the pore air is dissipated, the pore air hinders the water infiltration and increases the pore air pressure significantly. Considering only the water phase flow by the FEM, the pore water pressure isolines for the slope in
Fig.5 Calculation results in the rainfall situation
According to the pore air pressure and pore water pressure obtained by the two-phase flow program and by formulas (15) and (16), the anti-sliding safety coefficients of the slip Surfaces 1 and 2 of the following four cases are obtained(Fig.8 and Fig.9): (1) only the part pw ! 0 is considered, as in the method
645
used in the saturated soil, (2) both the part pw ! 0 and p g are considered, (3) only pw is considered, (4) both pw and p g are considered.
Fig.6
slope safety can be neglected, and the increase of the safety coefficient is mainly due to the negative pore water pressure, in the rainfall situation, the influence of the pore air pressure on the slope safety increases, but is also less than that due to the negative pore water pressure. Figure 10 shows that the larger the distance between the slip surface and the underground water level is, the greater the influence of the pore pressure in the unsaturated zone on the slope stability will be.
Pore water pressure in the stable seepage situation considering only the water-phase flow
Fig.7 Pore water pressure in the rainfall situation considering only the water-phase flow
Fig.8 Safety coefficients of sliding Surface 1
Fig.9 Safety coefficients of sliding Surface 2
Comparing Cases 2 and 4 with Case 1, the percentage increase of safety coefficient can be obtained(Fig.10). Figures 8 and 9 show that when the pore air pressure and the pore water pressure are both considered(Case 4), the safety coefficient is significantly larger than that obtained by the method used in saturated soil(Case 1), in the stable seepage situation, the influence of the pore air pressure on the
1—F of slip Surface 1 in stable seepage situation, 2—F of slip Surface 1 at the time the rainfall stops, 3—F of slip Surface 2 in stable seepage situation, 4—F of slip Surface 2 at the time the rainfall stops Fig.10 Percentage increase of safety coefficient as compared with Case 1
5. Conclusions In this article, a water-air two-phase flow model is proposed to simulate the water-air two-phase seepage and to calculate the pore pressure in the soil, and a method is suggested to calculate the safety coefficients of soil slopes by the pore pressure. The pore fluid flow of a soil slope in stable seepage and rainfall situations is simulated and the safety coefficients of the slip surfaces are determined, and some conclusions can be drawn as follows: (1) In the stable seepage situation, the pore air pressure in the unsaturated zone is about zero, so the pore air pressure can be neglected in the slope safety coefficient calculation. In the rainfall situation, the negative pore water pressure in the unsaturated zone is reduced, and the pore air pressure increases because the pore air provides some resistance to the rainwater infiltration, thus the pore air and pore water pressure both affect the slope stability. (2) After considering the shear strength provided by the pore air pressure and the negative pore water pressure in the unsaturated zone, the slope safety coefficient will become significantly larger than that obtained by the method used in the saturated soil, and the pore air pressure reduces the safety coefficient while the negative pore water pressure increases it. The pore pressure in the unsaturated zone has especially great influence on the slope safety while the
646
slip surface is mainly in the unsaturated zone. (3) In the soil pores, there are both air and water, and in the general slope stability analysis, only the water flow is considered, however, by the water-air two-phase flow model, the water-phase and air-phase are both taken into consideration, so the method proposed in the paper to solve the slope safety coefficient is more accurate than the conventional method. References [1]
[2] [3]
[4]
[5] [6]
SUN Dong-mei, ZHU Yue-ming and ZHANG Ming-jin. Study on numerical for water-air two-phase flow in unsaturated soil[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(4): 560-565(in Chinese). FU Jun-feng, JIN Sheng. The seepage computation based on saturation description[J]. Chinese Journal of Hydrodynamics, 2008, 23(6): 668-674(in Chinese). XU Jiong, WANG Tong and YANG Bo et al. The measurement and analysis of motion behavior of bubbles in calm water[J]. Chinese Journal of Hydrodynamics, 2008, 23(6): 709-714(in Chinese). GAGANIS P., KARAPANAGIOTI H. K. and BURGANOS V. N. Modeling multicomponent NAPL transport in the unsaturated zone with the constituent averaging technique[J]. Advances in Water Resources, 2002, 25(7): 723-732. KEES C. E., MILLER C. T. Higher order time integration methods for two-phase flow[J]. Advances in Water Resources, 2002, 25(2): 159-177. OOSTROM M., LENHARD R. J. Comparison of relative permeability-saturation-pressure parametric models for infiltration and redistribution of a light nonaqueous-phase liquid in sandy porous media[J]. Advances in Water Resources, 1998, 21(2): 145-157.
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14] [15]
LAROCHE C., VIZIKA O. Two-phase flow properties prediction from small-scale date using pore-network modeling[J]. Transport in Porous Media, 2005, 61(1): 77-91. HU Li-ming, XING Wei-wei and ZHOU Xiao-wen. Laboratory testing and numerical simulation of multiphase flow in unsaturated soils[J]. Engineering Mechanics, 2008, 25(11): 162-166(in Chinese). SCHROTH M. H., ISTOK J. and SELKER J. S. et al. Multifluid flow in bedded porous media: Laboratory experiments and numerical simulation[J]. Advances in Water Resources, 1998, 22(2): 169-183. CLASS H., HELMIG R. and BASTIAN P. Numerical simulation of non-isothermal multiphase multicomponent processes in porous media[J]. Advances in Water Resources, 2002, 25(5): 533-550. MA Cui-lin, ZHU Ming and WANG Jian-hua. Application of finite difference method FLAC in analyzing slope stability[J]. Chinese Mine Engineering, 2008, 37(5): 19-22(in Chinese). BETTINA A. A multiscaling problem for a model of partially saturated soils[C]. Proceedings of the 1st International Conference on Long Term Effects and Seepage Behavior of Dams. Nanjing: Hohai University Press, 2008, 203-213. LU Dong-qiang. Unsteady waves due to an impulsive oseenlet beneath the capillary surface of a viscous fluid[J]. Journal of Hydrodynamics, 2008, 20(1): 23-29. KOLDITZ O., De JONGE J. Non-isothermal two-phase flow in low-permeable porous media[J]. Computational Mechanics, 2004, 33(5): 345-364. ZHU Yue-ming, CHEN Jian-yu and GONG Dao-yong et al. Refined solution to seepage in arch dam foundation with FEM[J]. Chinese Journal of Geotechnical Engineering, 2003, 25(3): 326-330(in Chinese).
2009,21(5):640-646 DOI: 10.1016/S1001-6058(08)60195-X
THE ROLE FORE AIR FLOW IN SOIL SLOPE STABILITY ANALYSIS* ZHANG Xiao-yue, ZHU Yue-ming, FANG Chun-hui College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China, E-mail: [email protected]
(Received December 3, 2008, Revised January 18, 2009)
Abstract: In order to investigate the effect of the pore air flow in the soil slope stability analysis, a water-air two-phase flow model, based on the multi-phase flow theory, is proposed and with the model, the water-phase and air-phase seepages of the soil slope in the stable seepage and rainfall situations are simulated. The soil slope safety coefficients are obtained according to the simulated pore air pressure and pore water pressure. The calculation results show that in the stable seepage situation, the influence of the pore air pressure on the soil slope stability can be neglected, while in the rainfall situation, the pore air pressure generated in the unsaturated zone will reduce the safety coefficient, and the larger the distance between the slip surface and the underground water level is, the greater the influence of the pore air pressure on the soil slope stability will be. Key words: unsaturated soil, water-air two-phase flow, rainfall, slope stability
1. Introduction The safety of the natural and artificial soil slopes is an important issue. The unsaturated soil is a four-phase mixture, composed of solid-phase, air-phase, liquid-phase and water-air interface called shrinkable film[1], so one has to deal with not only water-phase flow but also air-phase flow in the soil[2,3]. Because the pore air pressure is hard to measure, it is not taken into consideration in the general slope stability analysis. In resent years, the multi-component multi-phase flow theory has been made a great progress. For example, Petros Gaganis advanced a multi-phase model[4], Kees proposed a time and spatial discretization method for the two-phase flow model[5], Oostrom analyzed the relationship between saturation and relative seepage coefficient in porous * Project supported by the China-Austria scientific and technological cooperative project of National Science and Technology Ministry (Grant No. CN 01/2007), the Ministry of Water Resources’ 948 plan, namely technical innovation and transformation projects (Grant No. CT200515). Biography: ZHANG Xiao-yue (1984-), Female, Ph. D. Candidate
media[6], Laroche studied the relationship between capillary pressure and saturation in a multiphase model[7], and all of these make it possible to analyze the air-phase flow in soil slopes and consider the influence of pore air pressure on slope stability. This article proposes a water-air two-phase flow model to simulate the water-air seepage in a homogeneous soil slope under stable seepage and rainfall conditions. The anti-sliding safety coefficients of the slope are calculated from the results obtained by the two-phase flow simulation, and the effect of the pore air pressure on the slope stability is analyzed.
Fig.1 Phase states and their conversions
641
transformed into the following form: 2. Water-air two-phase flow model 2.1 Definitions and assumptions The seepage in an unsaturated soil is a water-air two-phase flow, and the air-phase and water-phase both involve two components: water and air[8]. In this article, the subscripts g and w represent the air-phase and water-phase, respectively, while the superscripts a and w represent the air and water components, respectively. The air-phase and water-phase are transformed to each other by evaporation and dissolution (Fig.1). Some assumptions about the model are made. The soil skeleton is undeformable, the water phase is incompressible, the temperature keeps constant, hysteresis effect in the constitutive relationship is not considered, the chemical and biological reactions are neglected, the ideal gas law is valid for the air component, the partial pressure of the water component in the air phase is equal to the saturation pressure. 2.2 The governing equation and its solution With the mass transfer between the water-phase and the air-phase being taken into consideration, the governing equation of the two-phase flow model takes the form of the mass conservation of the component N [9]:
(2)
where the vector x represents the primary variables. Expanding Eq.(2) in Taylor’s series, neglecting high-order terms, we have at the iteration step m +1 and time step k + 1 :
§ wF · F ( x k +1,m +1 ) | F ( x k +1, m ) + ¨ < ¸ © wx ¹ k +1,m ( x k +1,m +1 x k +1,m )
(3)
As F ( x k +1, m +1 ) should be equal to zero, Eq.(3) can be rewritten as:
K ( x k +1,m )u = F ( x k +1,m )
(4)
where K = wF/wx is the coefficient matrix, u = x k +1, m +1 x k +1, m is the vector of the primary variable increments, F ( x k +1, m ) is the residual term at time step k + 1 and iteration step m . The elements K ij of the matrix K are obtained by finite difference method[11]:
§ · w ¨ ¦ SD U mol ,D xDN ¸ ª krD U mol ,D xDN D © ¹ ¦ div « K I < wt PD D ¬
(gradpD UD g )@ QN = 0
F ( x) = 0
K ij = (1)
where N = w, a represent the water component and air component, respectively, D = w, g represent the water-phase and air-phase, respectively, I , S ,
U mol , xDN and t are the porosity, saturation, molar density, mole fraction of component N in phase D and the time, respectively, K = k sw P w / U w g represents the tensor of absolute permeability of the soil, k sw and kr are the saturated and relative permeability coefficients of the soil, respectively, and p , P , g and QN are the pressure, the absolute viscosity, the gravity acceleration and the source term, respectively. In the governing equation, there are primary and secondary parameters, and the secondary parameters are functions of the primary parameters, and these parameters will be discussed later on. The nonlinear governing equation can be solved by Newton-Raphson method[10]. Equation (1) is first
wFi k +1,m | ª Fi (..., x j 1 , x j + 'x j , x j +1 ,...) wx kj +1,m ¬
Fi (..., x j 1 , x j 'x j , x j +1 ,...) ¼º / 2'x j
(5)
where 'x j is a small increment, j = 1, 2,! , n , n is equal to the number of primary variables times the element number. The more detailed calculation steps of the Newton-Raphson method can be described as follows. At time step k + 1 and iteration step m , the initial value x k +1,0 is the solution at time step k , while F ( x k +1, m ) ! H , solve
K ( x k +1, m )u = F ( x k +1,m ) ,
2
and let x k +1, m +1 = x k +1, m + K u , m = m +1 , do the next iteration until F ( x k +1, m ) H , then do calculation at 2
the next time step. In the above description, H is the allowable error, K 1 is the step coefficient of the iteration. 2.3 Primary variables The phase states in a given control volume are variables in space and time, and are directly related with the primary variables. According to the Gibbsian
642
phase rule[10], the number of degrees of freedom in a multi-component multi-phase system is F = C P + 2 + ( P 1) = C +1 , where C and P are the component number and phase number, respectively. As the temperature is kept as a constant, there should be two independent primary variables. With only air phase, the primary variables are xgw and p g , with only water phase, the primary variables are xwa and p g , with both water and air phases, the
is the Henry coefficient, T is the absolute temperature. Using the idea gas law, the mole density and the density of the air phase can be expressed as
U mol , g =
(9)
U g = U mol , g ¦ xgk M N
(10)
N
primary variables are S w and p g . The algorithm for phase transformation in the program is shown in Fig.2.
p n = g Vg RT
where R = 8.314J /(mol K) is the universal gas constant, n , Vg , T and M N are the quantity of the gas, the volume occupied by n , the absolute temperature and the molecular weight of component N , respectively. The difference between the water-phase pressure and the air-phase pressure, which, respectively, act on the two sides of the water-air interface, is usually called capillary pressure[13] pc = pw pg d 0 . The capillary pressure and the relative permeability coefficient are both the functions of the saturation degree, and they can be evaluated by using the approximate formulas proposed by Van Genuchten[14]:
pc ( S w ) = P0 ( S e 1/ m 1)1/ n Fig.2 Algorithm for phase transformation in the program
2.4 Secondary variables The mole fractions satisfy the following relation:
¦N x
N D
xgw =
=1
(6)
w sat
p pg
(7)
w is the saturation pressure of the water where psat vapor. According to Henry’s law and Dalton’s law for partial pressures[12], we have
w xwa = K H ( pg psat )
(8)
where
K H = [0.8942 +1.47 exp(0.04394T )] u 10
10
(11)
m k rw ( S w ) = Se1/ 2 ª«1 1 Se1/ m º» ¬ ¼
2
(12)
krg ( S w ) = (1 Se )1/ 3 1 Se1/ m
2m
(13)
where m and n are the empirical parameters, m = 1 1/ n , P0 is the air entry value of the soil, S e = ( S w S rw ) /( S sw S rw ) is the effective saturation, Srw and S sw are the residual and maximum water saturation, respectively. 2.5 Conditions for a unique solution Thin layer elements, called boundary condition elements, are added to the model boundary. In order to keep the boundary condition fixed in the calculation process, the volumes of the thin layer elements are assumed to take a maximum value. The initial condition of the stable seepage flow calculation is as follows: the soil elements in the slope are water-saturated, with primary variables of xwa and
pg ,
xwa = 1010 and
pg
equal to the standard
atmospheric pressure. The boundary condition elements above the upstream and downstream water
643
levels are air-saturated, with primary variables of xga and p g , xga = 0.9999 and p g equal to the standard atmospheric pressure. The boundary condition elements below the upstream and downstream water levels are water-saturated, with primary variables of xwa and p g , xwa = 1010 and p g equal to the water pressure of the centroid of the boundary condition element. The initial condition of the unstable seepage flow calculation is what obtained by the stable flow calculation. There are two kinds of boundary conditions: Dirichlet boundary conditions, as the air pressure and water pressure boundary conditions, which specify initial values to the boundary condition elements, Neumann boundary conditions, as the water and air flux boundary conditions, which specify source terms to the boundary condition elements. In the rainfall condition, the source term can be obtained by:
Q w (t ) = U w Aq (t )
(14)
where A is the effective area of the soil surface in the vertical direction of the rainwater infiltration, q (t ) is the rainfall intensity.
3. Stability analysis method With the circular slip surface, according to the simplified Bishop method, it is assumed that X R X L = 0 , with X R and X L denoting the inter-slice vertical shears. With the pore pressures at the grid nodes obtained from the water-air two-phase flow calculation, the pore pressures at the midpoints of the soil slice hemlines are determined by interpolation based on the inverse isoparametric mapping in FEM[15], then the slope anti-sliding stability coefficient F , as the ratio of the shear strength to the shear stress of the slip surface, is obtained by the moment balance. According to the vertical force balance, the total normal force N i acting on the soil slice hemline i is:
Ni =
1 mD ,i
oblique length of the soil slice hemline i , c c is the effective cohesion, I c is the effective angle of internal friction, I b represents the increasing rate of shear strength with the capillary pressure, and if the soil slice hemline is in the saturated zone, I b is equal to I c , D i is the angle between the horizontal plane and the tangent of the soil slice hemline through its midpoint, Wi is the gravity of the soil slice i ,
pw,i and pg , j are the pore water and pore air pressures at the midpoint of the soil slice hemline, respectively. According to the moment balance, we have
F=
1 ¦Wi sin Di
¦ ^ccE i
Ei sin D i F
º (tan I c tan I b ) » ¼
where mD ,i = cos D i + (sin D i tan I c) / F , E i
§ tan I b · º tan I b p g ,i E i ¨ 1 ¸ » tan I c` tan I c tan I c ¹ ¼ ©
(16)
The soil slope safety coefficients can be determined by Eqs.(15) and (16). Let the negative pore water pressure pw,i and pore air pressure pg ,i be zero in the calculation, then Eqs.(15) and (16) are turned into the calculation formulas in the saturated soil slope stability analysis. 4. Examples A homogeneous soil slope is taken as an example, with upstream water level of 11 m, downstream water level of 16 m, slope angle of 27o, and an impervious bottom. Here, two slip surfaces are studied. The slip surface 1 is deep with most of the soil slice hemlines in the saturated zone and in shape of an arc with center at (16.0, 25.0) m, radius of 16.16 m, and the slip surface 2 is shallow with most of the soil slice hemlines in the unsaturated zone and also in shape of an arc with center at (20.0, 27.0) m, radius of 13.74 m. The cross section of the slope is shown in Fig.3.
(15) is the
+ ª¬ N i pw,i Ei <
i
ccE i sin D i E sin D i ª + pw , i i < «Wi F F ¬
tan I b + pg ,i
i
Fig.3 Cross section of the slope
644
The parameters of the soil used in the simulation are shown in Table 1. Table1 Material parameters Parameter
Symbol
Value
Saturated seepage coefficient
k sw (m/s)
10-6
Air entry value
P0 (kPa)
20
Residual water saturation
S rw
0.15
Maximum water saturation
S sw
1.0
Temperature
T (ć )
15.0
Empirical parameter
m
0.457
Porosity
I
0.421
Effective cohesion
cc (kPa)
10.1
Effective angle of internal friction
Ic
42.6
Increasing rate angle
Ib
35.0
Density of soil
U (kN/m3)
18.94
Density of saturated soil
Usat (kN/m3)
19.63
the stable seepage situation and after 5 h of rainfall are shown in Fig.6 and Fig.7, respectively. By comparing the calculation results obtained by the water-air two-phase flow method and the FEM, the water-air two-phase flow model can be verified.
Fig.4 Calculation results in the stable seepage situation
Using the water-air two-phase flow model, the pore pressure and the saturation of the soil slope in the stable seepage situation and the rainfall situation with rainfall intensity of 2.0×10-5 m/s and rainfall time of 5 h, are obtained (Fig.4 and Fig.5), in which the pore air pressure and the pore water pressure are the corresponding water head hg = ( pg p0 ) / U w g and hw = ( pw p0 ) / U w g , respectively, and p0 is the standard atmospheric pressure. From Figs.4 and 5, it can be seen that in the stable seepage situation, the pore air pressure is about zero in the unsaturated zone, after 5 h of rainfall, there appears a transient saturated zone in the soil surface layer, and an increased water saturation in the unsaturated zone, which reduces the negative pore water pressure. In the rain water infiltration process, not all the pore air is dissipated, the pore air hinders the water infiltration and increases the pore air pressure significantly. Considering only the water phase flow by the FEM, the pore water pressure isolines for the slope in
Fig.5 Calculation results in the rainfall situation
According to the pore air pressure and pore water pressure obtained by the two-phase flow program and by formulas (15) and (16), the anti-sliding safety coefficients of the slip Surfaces 1 and 2 of the following four cases are obtained(Fig.8 and Fig.9): (1) only the part pw ! 0 is considered, as in the method
645
used in the saturated soil, (2) both the part pw ! 0 and p g are considered, (3) only pw is considered, (4) both pw and p g are considered.
Fig.6
slope safety can be neglected, and the increase of the safety coefficient is mainly due to the negative pore water pressure, in the rainfall situation, the influence of the pore air pressure on the slope safety increases, but is also less than that due to the negative pore water pressure. Figure 10 shows that the larger the distance between the slip surface and the underground water level is, the greater the influence of the pore pressure in the unsaturated zone on the slope stability will be.
Pore water pressure in the stable seepage situation considering only the water-phase flow
Fig.7 Pore water pressure in the rainfall situation considering only the water-phase flow
Fig.8 Safety coefficients of sliding Surface 1
Fig.9 Safety coefficients of sliding Surface 2
Comparing Cases 2 and 4 with Case 1, the percentage increase of safety coefficient can be obtained(Fig.10). Figures 8 and 9 show that when the pore air pressure and the pore water pressure are both considered(Case 4), the safety coefficient is significantly larger than that obtained by the method used in saturated soil(Case 1), in the stable seepage situation, the influence of the pore air pressure on the
1—F of slip Surface 1 in stable seepage situation, 2—F of slip Surface 1 at the time the rainfall stops, 3—F of slip Surface 2 in stable seepage situation, 4—F of slip Surface 2 at the time the rainfall stops Fig.10 Percentage increase of safety coefficient as compared with Case 1
5. Conclusions In this article, a water-air two-phase flow model is proposed to simulate the water-air two-phase seepage and to calculate the pore pressure in the soil, and a method is suggested to calculate the safety coefficients of soil slopes by the pore pressure. The pore fluid flow of a soil slope in stable seepage and rainfall situations is simulated and the safety coefficients of the slip surfaces are determined, and some conclusions can be drawn as follows: (1) In the stable seepage situation, the pore air pressure in the unsaturated zone is about zero, so the pore air pressure can be neglected in the slope safety coefficient calculation. In the rainfall situation, the negative pore water pressure in the unsaturated zone is reduced, and the pore air pressure increases because the pore air provides some resistance to the rainwater infiltration, thus the pore air and pore water pressure both affect the slope stability. (2) After considering the shear strength provided by the pore air pressure and the negative pore water pressure in the unsaturated zone, the slope safety coefficient will become significantly larger than that obtained by the method used in the saturated soil, and the pore air pressure reduces the safety coefficient while the negative pore water pressure increases it. The pore pressure in the unsaturated zone has especially great influence on the slope safety while the
646
slip surface is mainly in the unsaturated zone. (3) In the soil pores, there are both air and water, and in the general slope stability analysis, only the water flow is considered, however, by the water-air two-phase flow model, the water-phase and air-phase are both taken into consideration, so the method proposed in the paper to solve the slope safety coefficient is more accurate than the conventional method. References [1]
[2] [3]
[4]
[5] [6]
SUN Dong-mei, ZHU Yue-ming and ZHANG Ming-jin. Study on numerical for water-air two-phase flow in unsaturated soil[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(4): 560-565(in Chinese). FU Jun-feng, JIN Sheng. The seepage computation based on saturation description[J]. Chinese Journal of Hydrodynamics, 2008, 23(6): 668-674(in Chinese). XU Jiong, WANG Tong and YANG Bo et al. The measurement and analysis of motion behavior of bubbles in calm water[J]. Chinese Journal of Hydrodynamics, 2008, 23(6): 709-714(in Chinese). GAGANIS P., KARAPANAGIOTI H. K. and BURGANOS V. N. Modeling multicomponent NAPL transport in the unsaturated zone with the constituent averaging technique[J]. Advances in Water Resources, 2002, 25(7): 723-732. KEES C. E., MILLER C. T. Higher order time integration methods for two-phase flow[J]. Advances in Water Resources, 2002, 25(2): 159-177. OOSTROM M., LENHARD R. J. Comparison of relative permeability-saturation-pressure parametric models for infiltration and redistribution of a light nonaqueous-phase liquid in sandy porous media[J]. Advances in Water Resources, 1998, 21(2): 145-157.
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14] [15]
LAROCHE C., VIZIKA O. Two-phase flow properties prediction from small-scale date using pore-network modeling[J]. Transport in Porous Media, 2005, 61(1): 77-91. HU Li-ming, XING Wei-wei and ZHOU Xiao-wen. Laboratory testing and numerical simulation of multiphase flow in unsaturated soils[J]. Engineering Mechanics, 2008, 25(11): 162-166(in Chinese). SCHROTH M. H., ISTOK J. and SELKER J. S. et al. Multifluid flow in bedded porous media: Laboratory experiments and numerical simulation[J]. Advances in Water Resources, 1998, 22(2): 169-183. CLASS H., HELMIG R. and BASTIAN P. Numerical simulation of non-isothermal multiphase multicomponent processes in porous media[J]. Advances in Water Resources, 2002, 25(5): 533-550. MA Cui-lin, ZHU Ming and WANG Jian-hua. Application of finite difference method FLAC in analyzing slope stability[J]. Chinese Mine Engineering, 2008, 37(5): 19-22(in Chinese). BETTINA A. A multiscaling problem for a model of partially saturated soils[C]. Proceedings of the 1st International Conference on Long Term Effects and Seepage Behavior of Dams. Nanjing: Hohai University Press, 2008, 203-213. LU Dong-qiang. Unsteady waves due to an impulsive oseenlet beneath the capillary surface of a viscous fluid[J]. Journal of Hydrodynamics, 2008, 20(1): 23-29. KOLDITZ O., De JONGE J. Non-isothermal two-phase flow in low-permeable porous media[J]. Computational Mechanics, 2004, 33(5): 345-364. ZHU Yue-ming, CHEN Jian-yu and GONG Dao-yong et al. Refined solution to seepage in arch dam foundation with FEM[J]. Chinese Journal of Geotechnical Engineering, 2003, 25(3): 326-330(in Chinese).