坡地災害潛勢、觀測與預警技術之研究（二）-子計畫：

5 ^{成 果 報 告}

行政院國家科學委員會補助專題研究計畫

□期中進度報告

坡地災害潛勢、觀測與預警技術之研究（二）-子計畫：

多孔介質彈性力學與生態工法在邊坡穩定上之應用研究（II）

計畫類別：□ 個別型計畫

5

^{整合型計畫}

計畫編號：NSC 94－2625－Z－216－001－

執行期間： 94 年 8 月 1 日至 95 年 7 月 31 日

計畫主持人：呂志宗 共同主持人：

計畫參與人員：

成果報告類型(依經費核定清單規定繳交)：□精簡報告

5

^完整報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、

列管計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年□二年後可公開查詢

執行單位：中華大學土木與工程資訊學系

中 華 民 國 95 年 7 月 31 日

摘 要

根據本計畫所完成的研究成果，已整理並發表了 18 篇期刊暨國 際性 EI 等級的研討會論文，如本報告之第一部分所示。其中包括前 往 韓 國 清 州 參 加 「 第 10 屆亞洲科技與數學研討會（10^th Asian Technology Conference in Mathematics）」、前往中國上海參加「2006 上 海 岩 土 工 程 國 際 會 議 （GeoShanghai International Conference 2006）」、前往大陸北京參加「第 34 屆國際水文地質大會（34^th Congress of International Association of Hydrogeologists）」、及參加在本國所舉辦 之「2006 東亞放射性廢棄物管理論壇研討會（2006 East Asia Forum on Radwaste Management Conference）」等等國際間相當重要之學術會議 並發表最新之研究成果，計畫主持人所付出之努力，極有助於後續相 關研究計畫之推動與執行。本計畫亦已將一部分成果發表在「中華理 工學刊」及國內或兩岸學術研討會中，尚有其他相關研究成果，正陸 續積極整理並準備發表中。

本計畫報告之第二部分為「公路邊坡生態工法之多媒體教材建 構」，主要是以公路邊坡生態工法為案例探討對象，利用網路及多媒 體設計的技術，將生態工法之定義、精神、工法的種類以及應用生態 工法施作完成的公路邊坡案例等，製作成教材並發布在網路上，以供 各界參考。

本計畫報告之第三部分為「邊坡穩定之輔助分析軟體的應用與研 發」，是以有限差分法軟體 FLAC 3D 為架構，輔以有限元素法軟體 ANSYS，再利用 Visual Basic 程式，撰寫一應用於邊坡穩定之輔助分 析程式，用以化簡邊坡穩定數值分析之前處理步驟，使有利於FLAC 3D 軟體之使用。以此為基礎，在後續的研究中，可繼續探討邊坡土 壤中之水份等因素的影響，亦即可引用有考慮水份與土壤之交互作用 現象的多孔介質彈性力學理論，研討邊坡穩定之相關問題。

第四部分為出席國際學術會議之報告，共兩份。包括參加「第 10 屆亞洲科技與數學研討會」及「2006 上海岩土工程國際會議」之 會議報告。此一部分並非本計畫報告應具備之附件，僅供參考用。

關鍵詞：地下水、沉陷、生態工法、邊坡穩定

第 第 一 一 部 部 分 分

已 已 發 發 表 表 在 在 期 期 刊 刊 暨 暨 國 國 內 內 外 外 學 學 術 術 會 會 議 議 之 之 研 研 究 究 成 成 果 果

註： 註 ：共 共 18 1 8 篇論 篇 論文 文

I

目 錄

Lu, John C.-C. and Feng-Tsai Lin, 2006/6/6~8, “The Transient Ground Surface Displacements Due to a Point Sink/Heat Source in an Elastic Half-Space,”

Geotechnical Special Publication No. 148, GeoShanghai International Conference 2006, ASCE, Shanghai, China, pp. 210-218. (Supported by the National Science Council of NSC94-2625-Z-216-001.) (EI)

Lu, John C.-C. and Feng-Tsai Lin, 2006/10/9~13, “Analysis of Transient Ground Surface Displacements Due to a Point Sink in a Pervious/Impervious Poroelastic Half-Space,” 34^th Congress of International Association of Hydrogeologists, Beijing, China, p. 353. (Supported by the National Science Council of NSC94-2625-Z-216-001.)

Lu, John C.-C., 2006, “Long-Term Behaviors of a Buried Deep Point Heat Source in a Transversely Isotropic Thermoelastic Porous Medium,” Chung Hua Journal of Science and Engineering, Vol. 4, No. 1, pp. 11-22. (Supported by the National Science Council of NSC94-2625-Z-216-001.)

Lu, John C.-C. and Feng-Tsai Lin, 2005/12/12~16, “Analysis of Transient Ground Surface Displacements Due to a Point Sink in a Porous Elastic Half-Space,”

Proceedings of the 10^th Conference of Advanced Technology Council in Mathematics, Cheong-Ju, Korea, pp. 135-144. (Supported by the National Science Council of NSC94-2625-Z-216-001.)

呂志宗，2005/10/13~14，「地表不透水情況下因單點抽水所引致的暫態壓密沉陷 解析」，九十四年度農業工程研討會論文摘要集，苗栗市，台灣，中華民國，

第113 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。) 呂志宗、戴雪蘭，2005/10/13~14，「溫泉自岩層裂隙滲流所引致之地層行為變化

理論解析」，九十四年度農業工程研討會論文摘要集，苗栗市，台灣，中華 民國，第114 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下 所完成。)

呂志宗、任克泰、曾柏領，2005/11/25~26，「Benford 定律與河川流量關係之探討」， 2005 年海峽兩岸科學技術研討會論文摘要集，斗六市，台灣，中華民國，

第437 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。) 呂志宗、林鳳彩，2005/11/25~26，「橫向等向性/等向性地層因溫泉抽水所引致的

壓密沉陷解析」，中國機械工程學會第二十二屆全國學術研討會論文集，中 壢市，台灣，中華民國，第 1383-1389 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。)

呂志宗、林鳳彩，2005/11/25~26，「深層溫泉開發所引致之地層力學行為變化的 暫態閉合解」，中國機械工程學會第二十二屆全國學術研討會論文集，中壢

II

市，第1417-1423 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補 助下所完成。)

呂志宗、周宜興，2005/12/16~17，「FLAC 3D 之輔助分析軟體的研發與應用」，

中華民國第二十九屆全國力學會議論文集，新竹市，台灣，中華民國，第 C003-1~C003-8 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助 下所完成。)

呂志宗、張君平，2005/12/16~17，「公路邊坡生態工法之互動式電腦輔助教學軟 體的設計」，中華民國第二十九屆全國力學會議論文集，新竹市，台灣，中 華 民 國 ， 第 C004-1~C004-8 頁 。 ( 本 文 係 在 國 科 會 研 究 計 畫 NSC94-2625-Z-216-001 補助下所完成。)

呂志宗、林鳳彩，2005/12/16~17，「深層溫泉開發所引致之地層力學行為的暫態 閉合解－考慮孔隙水與固體介質為可壓縮」，中華民國第二十九屆全國力學 會議論文集，新竹市，台灣，中華民國，第 C005-1~C005-8 頁。(本文係在 國科會研究計畫NSC94-2625-Z-216-001 補助下所完成。)

呂志宗、林鳳彩，2006/5/5，「以黃金比例觀察大自然之美」，正修科技大學 94 學年度通識教育學術研討會論文集，高雄縣，台灣，中華民國，第131-150 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。) 呂志宗、任克泰、溫志浩、廖英達、梁惠儀、楊俊丞、吳明銘，2006/10/19~20，

「以 Benford 定律探討北台灣之流量及雨量特徵」，九十五年度農業工程研 討會論文摘要集，台南市，台灣，中華民國，第 114 頁。(本文係在國科會 研究計畫NSC94-2625-Z-216-001 補助下所完成。)

呂志宗、林擎天，2006/10/19~20，「補注高溫水所引致地層彈性力學行為變化探 討」，九十五年度農業工程研討會論文摘要集，台南市，台灣，中華民國，

第115 頁。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。) Lu, John C.-C. and Feng-Tsai Lin, 2006/11/24~25, “Analytic Solutions of

Thermo-consolidation Due to a point Heat Source Buried in a Porous Elastic Half-Space,” Proceeding of the 23^rd National Conference on Mechanical Engineering, Tainan, Taiwan R.O.C., Accepted. (Supported by the National Science Council of NSC94-2625-Z-216-001.)

呂志宗、林擎天，2006/11/24~25，「溫泉在砂岩地層中之滲流機制探討」，中國機 械工程學會第二十三屆全國學術研討會，台南縣，台灣，中華民國，論文已 被接受。(本文係在國科會研究計畫 NSC94-2625-Z-216-001 補助下所完成。) Lu, John C.-C. and Feng-Tsai Lin, 2006/11/27~28, “Theoretical Solutions for a

Stratum Subjected to a Point Heat Source,” 2006 East Asia Forum on Radwaste Management Conference, Taoyuan, Taiwan R.O.C., Accepted. (Supported by the National Science Council of NSC94-2625-Z-216-001.)

120

The Transient Ground Surface Displacements Due to a Point Sink/Heat Source in an Elastic Half-Space

J. C.-C. Lu¹ and F.-T. Lin²

Abstract

Thermoelastic deformation due to a point heat source is the analog of poroelastic response caused by a point sink. In this paper, Biot’s three-dimensional consolidation theory is introduced to derive the analytical solutions of the transient consolidation deformation with a point sink in saturated isotropic porous elastic half-space. The transient ground surface displacement produced by a point heat source is described through analog quantities between poroelasticity and thermoelasticity. Closed-form solutions of the horizontal and vertical displacements are obtained by using Laplace and Hankel integral transforms. Attention is focused on the maximum surface horizontal displacement compared to the maximum surface settlement. Results show that the horizontal displacement is about 30% of the maximum ground surface settlement. The study concludes that horizontal displacement is significant and should be considered in prediction of the transient settlement induced by groundwater withdrawal.

Introduction

Land subsidence due to groundwater withdrawal is a well-known phenomenon (Poland 1984). The pore water pressure is reduced in the withdrawal region as water pumped from an aquifer. It leads to increase in the effective stress between the soil particles and subsidence of ground surface.

The three-dimensional consolidation theory presented by Biot (1941, 1955) is generally regarded as the fundamental theory for modeling land subsidence. Based on Biot’s theory, Booker and Carter (1986a, 1986b, 1987a, 1987b), Tarn and Lu (1991) presented solutions of subsidence by a point sink embedded in saturated elastic half-

1 Department of Civil Engineering and Engineering Informatics, Chung-Hua University, Hsinchu 30012, Taiwan, ROC; PH (886)3-518-6708; FAX (886)3-537-2188; email:

cclu@chu.edu.tw

2 Department of Naval Architecture, National Kaohsiung Marine University, Kaohsiung 81157, Taiwan, ROC; PH (866)7-361-7141 ext. 3285; FAX (866)7-362-8844; email:

ftlin@mail.nkmu.edu.tw

121

space at a constant rate. In the studies of Booker and Carter, the flow properties are considered as isotropic or cross-anisotropic whereas the elastic properties of the soil are treated as isotropic with pervious half-space boundary. Tarn and Lu found that groundwater withdrawal from an impervious half-space induces a larger amount of consolidation settlement than from a pervious one. The anisotropic permeability was proved to have significant effects on the land subsidence due to fluid extraction.

Nevertheless, transient closed-form solution of the half-space due to fluid withdrawal was not obtained in the above studies.

The governing equations of poro-mechanics are similar to the equations appeared in thermo-mechanics. Based on the analogy of poroelasticity and thermoelasticity, point heat source induced transient horizontal and vertical displacements of the ground surface are obtained. Figure 1 shows a point sink or heat source buried in a stratum at a depth h where the stratum is modeled as saturated/thermally isotropic elastic half-space. Point sink is usually introduced to simulate groundwater withdrawal and radioactive canister buried in a half-space can be treated similar to a point heat source.

Figure 1. Point sink Q (constant pumping rate)/heat source H (constant heat generation rate) induced mechanics of poroelastic/thermoelastic problem.

The pervious ground surface in poroelasticity is corresponding to a constant temperature of the ground surface in thermoelasticity. The transient horizontal and vertical displacements of the ground surface due to a point sink, corresponding to a point heat source, are obtained by using Laplace and Hankel transforms. Results are illustrated and compared to provide better understanding of the time dependent ground surface displacements due to pumping or point heat source.

Mathematical Model of Poroelasticity

Governing Equations. Figure 1 shows a point sink buried in a saturated porous stratum at a depth h. The soil mass is considered as a homogeneous isotropic porous medium.

The constitutive stress behaviors of the elastic soil skeleton are _{ij ij} G _{kk ij} p _ij

G ε δ δ

ν ε ν

σ −

+ −

= 1 2

2 2 , i,j=1,2,3, (1)

122

in which σ_ij are the total stress components; ε_ij are the strain components; ν is Poisson’s ratio and G is shear modulus of the stratum. The excess pore water pressure p is positive for compression, and δ_ij is the Kronecker delta. The strains ε_ij and displacement components u are given by the linear law _i

εij = ₂¹

(

ui_,j +uj_,i

)

, i,j=1,2,3. (2) The total stress must satisfy the equilibrium equations

σ_ij,j +b_i =0, i, j=1,2,3, (3) where b denote the body forces. Eqs. (1) and (2) are used in the equilibrium equations _i to express their forms in displacements u_i and excess pore water pressure p as follows:

0

2

1 ^{, ,}

, − =

+ − _{kki i}

jj

i G p

Gu ε

ν ^{, i=1,2,3}. (4) Consider a point sink of constant pumping rate Q that is located at point

( )

0,h . The uncoupled governing equation in axially symmetric coordinates

( )

r,z is derived from the conservation of mass and Darcy’s law as

( ) ( ) ( )

0

, +2 − =

∂ + ∂

− r z h ut

r Q t n p k p

jj w

δ π δ

γ β , (5)

where k and n are the permeability and porosity of the porous medium, respectively; β is the compressibility of pore water; γ_w is the unit weight of pore water; δ

( )

^x is the Dirac delta and u

( )

t is Heaviside unit step function. Eqs. (4) and (5) constitute the basic governing equations of the time-dependent poroelastic responses of a saturated porous medium.

Boundary Conditions and Initial Conditions. Consider the half-space surface, z = 0, which is a traction-free and pervious boundary for all time t≥0. The mathematical statements of the boundary z = 0 are:

σ_rz

(

r,0,t

)

=0, σ_zz

(

r,0,t

)

=0, and p

(

r,0,t

)

=0. (6) Assuming no initial change in displacements and seepage of the stratum, the initial conditions at time t=0 of the mathematical model due to a point sink can be treated as

u_r

(

r,z,0

)

=0, u_z

(

r,z,0

)

=0, and p

(

r,z,0

)

=0. (7)

123 Mathematical Model of Thermoelasticity

Governing Equations. The constitutive behavior of the isotropic body with a point heat source buried in a thermoelastic isotropic half-space at a depth h can be expressed by

( )

ij s ij

kk ij

ij

G G

G ϑδ

ν α δ ν

ν ε ε ν

σ 1 2

1 2 2

1 2 2

−

− + + −

= , i,j=1,2,3. (8)

Here, σ_ij are the thermal stress components and ϑ is the temperature increment measured from the reference state. The constants ν , G and α_s are the Poisson’s ratio, shear modulus, and linear thermal expansion coefficient, respectively.

The kinematic equation, Eq. (2), shows the relation between the strains and displacements, and the thermal stresses must satisfy the equilibrium relations, Eq. (3).

Using Eqs. (2) and (8), the equilibrium equations (3) can be expressed in terms of displacements u and temperature change of the thermoelastic half-space _i ϑ as follows:

( )

₀

2 1

1 2 2

1 ^{, ,}

, =

−

− +

+ − _kki ^s _i

jj i

G G

Gu ϑ

ν α ε ν

ν ^{, i=1,2,3}. (9) Consider a point heat source of constant heat generation rate H that is located at point

( )

0,h . The uncoupled governing equation, which is axially symmetry, is obtained from the conservation of energy and heat conduction law as following

( ) ( ) ( )

⁰

, −2 − =

∂ + ∂

− r z h u t

r H c t

jj

t δ δ

π ϑ ϑ

λ _ε , (10)

where λ_t is the thermal conductivity and c_ε =ρc. The constants ρ and c define the density and the specific heat of the thermoelastic medium, respectively. Eqs. (9) and (10) constitute the basic governing equations of the transient responses of a thermoelastic medium due to a point heat source.

Boundary Conditions and Initial Conditions. The half-space surface, z = 0, is considered as traction-free, and it does not have temperature change for all time t≥0. The boundary conditions on surface z = 0 are given by

σ_rz

(

r,0,t

)

=0, σ_zz

(

r,0,t

)

=0, and ϑ

(

r,0,t

)

=0. (11) Assuming there are no initial change of displacement and temperature for the thermal elastic medium, the initial conditions at time t=0 due to a point heat source can be treated as

u_r

(

r,z,0

)

=0, u_z

(

r,z,0

)

=0, and ϑ

(

r,z,0

)

=0. (12)

124

From these governing equations, the corresponding quantities of poroelasticity and thermoelasticity are shown in Table 1.

Table 1. Analogy of poroelastic and thermoelastic quantities.

Poroelasticity Thermoelasticity

p

( )

_ϑ

ν α ν 2 1

1 2

−

+ _s

G β

n

( )

( )

s

G c

α ν ν _ε +

− 1 2

2 1

w

k γ

( )

( )

^t_s

G ν α

λ ν +

− 1 2

2 1

Q − H

Analytic Solutions

Applying Laplace and Hankel integral transformations (Sneddon 1951, Erdelyi et al.

1954), the transient horizontal and vertical displacements of the ground surface z = 0 due to a point sink in axially symmetric coordinates

( )

r, are obtained as follows: z

( ) ( )

( )







− +

= − ₃₂

2

2 2

2 , 1

0

, h r

ctr Gk

t Q r

u_r ^w

π γ ν

( )





 



 



 

 



− 



 



 



 



− +

+

∫

0^{ct ct}− ^{hr r h r r d}

2 1 2 0 2 2

3 I 8

I 8 8

exp 2

16 τ

τ τ

τ τ

τ , (13a)

( ) ( )

( )

















 +

+

= −

ct r h r

h cth Gk

t Q r

u_z ^w

erf 2 2

2 , 1

0 ,

2 2 2

2 3

π 2

γ ν

















 +

+ +



 



− +

− +

ct r h r

h h ct

r h ct

r h

h

erfc 2 4 2

exp

2 2

2 2 2

2 2

2 π ^{, (13b)}

where c=k nβγ_w. The long-term ground surface horizontal and vertical displacements can be obtained by letting t→∞ in Eqs. (13a) and (13b):

( ) ( )

) 4 (

2 , 1

0

, h2 r2 h2 r2 h

hr Gk

r Q

u_r ^w

+ + +

− −

=

∞ π

γ

ν , (14a)

( ) ( )

2

4 2

2 , 1

0

, h r

h Gk

r Q

u_z ^w

+

= −

∞ π

γ

ν . (14b)

125

The maximum surface settlement u_z

(

0,0,∞

) (

= 1−2ν

)

Qγ_w 4πGk can be found from Eq.

(14b). Solutions of buried point heat source induced ground surface displacements can be easily derived through the parameters described in Table 1.

Numerical Results

For the point sink induced land subsidence, the particular interest is the settlement of stratum at each stage of the consolidation process. The average consolidation ratio U is defined as

n compressio of

end at settlement

time at settlement t

U= . (15)

This ratio U can be expressed as bellow,











 +

+



 



− +

− +











 +

= +

ct r h ct

r h ct

r ct h

r h r

h U ct

erfc 2 exp 4

2 erf 2

2 ^{2 2 2 2}

2 2 2

2 2

2 π . (16)

Figure 2 shows the average consolidation ratio U at r=0. Note that U initially decreases rapidly then the rate of settlement slows down. As U asymptotically approaches 1, its consolidation is theoretically never achieved.

Figure 2. Graphical interpretation of average consolidation ratio U at r=0. The profiles of normalized vertical and horizontal displacements at the ground surface z=0 are shown in Figures 3 and 4, respectively. Figure 4 illustrates that the ground surface has significant horizontal displacement. The maximum surface horizontal displacement is about 30% of the maximum vertical settlement. The displacements of the thermoelastic half-space due to a point heat source display a similar phenomenon as shown in Figures 2-4.

126

Figure 3. Normalized vertical displacement profile uz at the ground surface z=0.

Figure 4. Normalized horizontal displacement profile ur at the ground surface z=0. Conclusions

Closed-form solutions of the transient consolidation due to pumping from pervious elastic half-space were obtained using Laplace and Hankel transformations. Vertical settlement and ground surface horizontal displacement were investigated.

Based on numerical results, the maximum surface horizontal displacement is found to be about 30% of the maximum surface settlement. From the average consolidation ratio U at r=0, the trend revealed by this model agrees with previous models by Sivaram and Swamee (1977) that U initially decreases rapidly then the rate of settlement slows down. It concludes that horizontal displacement should be properly considered for better prediction of the transient settlement induced by groundwater withdrawal.

The corresponding quantities of poroelasticity and thermoelasticity are discussed through their governing equations. The ground surface horizontal and vertical

127

displacements of the half-space due to a buried point heat source can be derived through their corresponding analogy.

Acknowledgments

This work is supported by the National Science Council of Taiwan, ROC, through grant NSC94-2625-Z-216-001.

References

Biot, M. A. (1941). “General theory of three-dimensional consolidation.” J. Appl. Phys., 12(2), 155-164.

Biot, M. A. (1955). “Theory of elasticity and consolidation for a porous anisotropic solid.” J. Appl. Phys., 26(2), 182-185.

Booker, J. R., and Carter, J. P. (1986a). “Analysis of a point sink embedded in a porous elastic half space.” Int. J. Numer. Anal. Methods Geomech., 10(2), 137-150.

Booker, J. R., and Carter, J. P. (1986b). “Long term subsidence due to fluid extraction from a saturated, anisotropic, elastic soil mass.” Q. J. Mech. Appl. Math., 39(1), 85-97.

Booker, J. R., and Carter, J. P. (1987a). “Elastic consolidation around a point sink embedded in a half-space with anisotropic permeability.” Int. J. Numer. Anal. Methods Geomech., 11(1), 61-77.

Booker, J. R., and Carter, J. P. (1987b). “Withdrawal of a compressible pore fluid from a point sink in an isotropic elastic half space with anisotropic permeability.” Int. J. Solids Struct., 23(3), 369-385.

Erdelyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F. G. (1954). Tables of integral transforms, McGraw-Hill, New York.

Poland, J. F. (1984). Guidebook to studies of land subsidence due to ground-water withdrawal, The United Nations Educational Scientific and Cultural Organization, Paris, 3-16.

Sivaram, B., and Swamee, P. K. (1977). “A computational method for consolidation coefficient.” Soils Found., 17(2), 48-52.

Sneddon, I.N. (1951). Fourier transforms, McGraw-Hill, New York, 48-70.

Tarn, J.-Q., and Lu, C.-C. (1991). “Analysis of subsidence due to a point sink in an anisotropic porous elastic half space.” Int. J. Numer. Anal. Methods Geomech., 15(8), 573-592.

1

Analysis of Transient Ground Surface Displacements Due to a Point Sink in a Pervious/Impervious

Poroelastic Half-Space

John C.-C. Lu¹ and Feng-Tsai Lin²

1Department of Civil Engineering and Engineering Informatics Chung Hua University

Hsinchu 30012, Taiwan, ROC Email: cclu@chu.edu.tw

2Department of Naval Architecture National Kaohsiung Marine University

Kaohsiung 81157, Taiwan, ROC Email: ftlin@mail.nkmu.edu.tw

ABSTRACT

Based on Biot’s three-dimensional consolidation theory of porous media, analytical solutions of the transient consolidation deformation due to a point sink in saturated isotropic porous elastic half-space are presented.

Using the Laplace and Hankel integral transform techniques, closed-form solutions of the horizontal and vertical displacements of the ground surface are obtained. In the analysis, cases of pervious and impervious half-space boundary are studied. The consolidation effected by the consolidation parameters are illustrated and discussed.

The numerical results show that the maximum surface horizontal displacement is about 30% of the maximum surface settlement for the pervious ground surface. The results indicate pumping from an impervious half-space leads to a much larger amount of land subsidence than from a pervious one. The study concludes that horizontal displacement and boundary condition are significant and should be considered in the prediction of the transient settlement induced by groundwater withdrawal.

Keywords: Point sink, Closed-form solution, Porous medium, Groundwater withdrawal

1. Introduction

Land subsidence due to groundwater withdrawal is a well known phenomenon (Poland, 1984). As water pumps from an aquifer, the pore water pressure is reduced in the withdrawal region. It leads to increase in the effective stress between the soil particles and subsidence of ground surface.

The three-dimensional consolidation theory presented by Biot (1941, 1955) is generally regarded as the fundamental theory for modeling land subsidence. Based on Biot’s theory, Booker and Carter (1986a, 1986b, 1987a, 1987b), Tarn and Lu (1991) presented solutions of

subsidence by a point sink embedded in saturated elastic half-space at a constant rate. In the studies of Booker and Carter (1986a, 1986b, 1987a, 1987b), the flow properties are considered as isotropic or cross-anisotropic whereas the elastic properties of the soil are treated as isotropic with pervious half-space boundary. Tarn and Lu (1991) found that groundwater withdrawal from an impervious half-space induces a larger amount of consolidation settlement than from a pervious one. The anisotropic permeability was proved to have significant effects on the land subsidence due to fluid extraction. Lu and Lin (2005, 2006) displayed transient closed-form solutions of

2 the pervious half-space due to fluid withdrawal.

Nevertheless, transient closed-form solutions of the impervious half-space due to pumping were not obtained in the above studies.

In this paper, the soil mass is modeled as an isotropic saturated elastic half-space. Cases of pervious and impervious ground surface boundary are investigated.

Using the Laplace and Hankel transform techniques, transient horizontal and vertical displacements of the ground surface due to a point sink are obtained. Results are illustrated and compared to provide better understanding of the time dependent consolidation settlement due to pumping.

2. Mathematical Models

Basic Equations

Figure 1 presents a point sink buried in a saturated porous stratum at a depth h. The soil mass is considered as a homogeneous isotropic porous medium with a vertical axis of symmetry. The constitutive behavior of the elastic soil skeleton for linear axially symmetric deformation in the cylindrical coordinates (r ,,θ z) are expressed by

( ) _p

z u G r u G r u

G _{r r z}

rr −

∂

∂ + −

+ −

∂

∂

−

= −

ν ν ν

ν ν

τ ν

2 1

2 2

1 2 2

1 1

2 , (1a)

( ) _p

z u G r u G r u

G _{r r z} −

∂

∂ + −

− + −

∂

∂

= −

ν ν ν

ν ν

τ_θθ ν

2 1

2 2

1 1 2 2

1

2 , (1b)

( ) _p

z u G

r u G r u

G _{r r z}

zz −

∂

∂

− + − + −

∂

∂

= −

ν ν ν

ν ν

τ ν

2 1

1 2 2

1 2 2

1

2 , (1 c)



 





∂ +∂

∂

= ∂

r u z G u^{r z}

τrz , (1d) where τ , _rr τ , etc., are the total stress components; the _θθ excess pore fluid pressure p is positive for compression;

ur, u_z are the displacements in the radial and axial directions, respectively; ν and G are the Poisson’s ratio and shear modulus of the stratum, respectively. The

shear stress components τ_rθ and τ vanish by locating _θz the vertical z-axis through the point sink.

Figure 1. Point sink induced land subsidence problem.

The total stresses must satisfy the following equilibrium relations

=0

∂ + +∂ + −

∂

∂

rz r rr

rr b

z r

r

τ τ τ

τ _θθ

, (2a)

=0

∂ + +∂

∂ +

∂

zz z rz

rz b

z r r

τ τ

τ , (2b)

in which b_r and b_z denote the body forces. By using equations (1a)-(1d), the equilibrium equations for axially symmetric problem without body forces b_i can be expressed in terms of displacements u_i and excess pore water pressure p as follows:

2 0

1 ²

2 =

∂

−∂

∂ −

∂ + −

∇ r

p r Gu r u G

G _r ε ^r

ν , (3a) 2 0

1

2 =

∂

−∂

∂

∂ + −

∇ z

p z u G

G _z ε

ν , (3b) where ∇² =∂² ∂r²+1r∂ ∂r+∂² ∂z² is the Laplacian operator and the volume strain of the porous medium is

z u r u r

u_r ∂ + _r +∂ _z ∂

∂

ε = .

Assuming the problem is decoupled, i.e., the flow field is sought independently from that of the displacement field. The third relation between u_r, u_z and p obtained from the conservation of mass can be expressed as:

3

( )

[ ]

+ =0

∂ + ∂

−

⋅

∇ q

t n p

n v_w v_s β , (4) where n is the porosity of the porous medium; v_w and v_s are the velocities of pore water and solid matrix, respectively;

β is the compressibility of pore water; q is the rate of water extracted from the porous medium per unit volume.

Assuming that the pore water is governed by Darcy’s law, we have

( ) 



 





∂ +∂

∂

− ∂

=

− _{r z}

w s

w z

p r p

nv v k i i

γ , (5) in which k denotes the permeability of the soil mass;

γ is the unit weight of pore water. w

Let us consider a point sink of constant strength Q located at point ( )0,h . Substituting (5) into (4) yields

t n p z

p r p r r

p k

w ∂

+ ∂



 





∂ +∂

∂ + ∂

∂

− ∂²₂ ²₂ β γ

1

( ) ( ) ( ) 0

2 − =

+ r z hut r

Q δ δ

π , (6) in which δ( )x and u( )t are Dirac delta and Heaviside unit step function, respectively. Eqs. (3a), (3b) and (6) constitute the basic governing equations of the axially symmetric time-dependent poroelastic responses of a saturated porous medium.

Boundary Conditions

Consider the half-space surface, z=0, as a traction-free boundary for all time t≥0. From Eqs. (1c) and (1d), the mechanical boundary conditions at z=0 are expressed in terms of u_r and u_z by

( ) ₀

2 1

1 2 2

1

2 =

∂

∂

− + −



 



 +

∂

∂

− z

u G

r u r u

G _{r r z}

ν ν ν

ν , (7a)

=0



 





∂ +∂

∂

∂ r u z

G u^{r z} . (7b)

An additional condition is provided by considering the half-space as either pervious or impervious. The

mathematical statements of the two types of flow conditions at the boundary z = 0 are given by

For pervious half-space:

=0

p , (7c) For impervious half-space:

=0

∂

∂ z

p . (7d)

The boundary conditions at z→∞ due to the effect of the point sink must vanish at any time.

Initial Conditions

Assuming no initial changes in displacements and seepage of the stratum, the initial conditions at time t=0 of the mathematical model can be treated as

=0

ur , u_z =0, p=0. (8)

3. Analytic Solutions

Laplace and Hankel Transforms Solutions

The governing partial differential equations (3a), (3b) and (6) can be reduced to ordinary differential equations by performing appropriate Laplace and Hankel transforms (Sneddon, 1951) with respect to the time variable t and the radial coordinate r:

( ) ~ 1 ~ 0 1

~ 2 2 ²

2

2  − − + =



 



 − p

G dz

u u d

dz

d _z

r η ξ ξ

ηξ , (9a)

( ) ~ 2 ~ 1 ~ 0 1

2 ₂ ²

2  − =



 



 −

+

− dz

p d u G dz

d dz

u d

z

r η ξ

ξ

η , (9b)

( ) 0 2

~

~

2 2

2  + + − =



 



 −

− z h

s p Q s n dz p

d k

w

π δ β

γ ξ , (9c)

where ξ and s are Hankel and Laplace transform parameters. The parameter ^η=(1−^ν) (1−2^ν) and the symbols u~_r, u~_z, p~ are defined as

( ; , )⁼

∫

0^∞ { ( , , )} ( )1

~ z s rLu r zt J rdr

u_r ξ _r ξ , (10a)

4

( ; , )⁼

∫

0^∞ { ( , , )} ( )0

~ z s rLu r z t J rdr

u_z ξ _z ξ , (10b)

( ; , )⁼

∫

0^∞ { ( , , )} ( )0

~ z s rL p r zt J r dr

p ξ ξ , (10c) in which J_n( )x represents the first kind of Bessel function of order n and the Laplace transformations with respect to u_r, u_z and p are denoted by

( )

{u r,z,t}⁼

∫

0^∞u (r,z,t) ( )exp⁻stdt

L _{r r} , (11a)

( )

{u r,z,t}⁼

∫

0^∞u (r,z,t) ( )exp⁻st dt

L _{z z} , (11b)

( )

{pr,z,t}⁼

∫

0^∞p(r,z,t) ( )exp⁻st dt

L . (11c) The general solutions of equations (9a)-(9c) are obtained as

( )z C z ( )z C

u_r expξ expξ

~ = 1 + 2

( )z C z ( )z

C −ξ + −ξ

+ ₃exp ₄ exp











− +

+











 +

+ z

c C s

cz

C₅exp ξ² s ₆exp ξ²

(

z h

)

s c Gk Q _w

−

−

− ξ

πη

γ exp

8 ²











− + − +

+ z h

c s

c s s

c Gk

Q _{w 2}

2

2 exp

8 ξ

ξ ξ πη

γ , (12a)

( )z C z ( )z C

C

u_z ξ ξ

ξ η

η 1 exp exp

1 2

1

~ 2

2 2

1  −



 





− + +

−

=

( )z C z ( )z C

C ξ ξ

ξ η

η  − + −



 





− + +

+ 1 exp exp

1 2

1 2

4 4

3











 +

+

− z

c C s

c

s 2

5

2 exp

1 ξ ξ

ξ











− + +

+ z

c C s

c

s 2

6

2 exp

1 ξ ξ

ξ

(

^{z h}

)

s c Gk Q _w

−

−ξ πη

γ exp

8 ²

m











− + −

± z h

c s s

c Gk

Q _{w 2}

2 exp

8 ξ

πη

γ , (12b)











 +

−

= z

c C s

c G s

p 2 1 ₅exp ²

~ ξ

η ξ











− +

− z

c C s

c

G1s ₆exp ²

2 ξ

η ξ











− + −

+

− z h

c s

c s s

k

Q _w 2

2

1 exp 1

4 ξ

π ξ

γ , (12c)

where the parameterc=k nβγ_w ; and C_i(i=1,2,_L,6)

are functions of the transformed variables ξ and s which must be determined from the transformed boundary conditions. The upper and lower signs in equation (12b) are for the conditions of (z− h)≥0 and (z− h)<0, respectively.

Transformed Boundary Conditions

Taking Hankel and Laplace transforms for Eqs. (7a)-(7d) yield the mechanical boundary conditions at z = 0 of the transformed domains as follows:

~ 0

~^r − u_z = dz

u

d ξ , (13a)

( 1) ~ 0

~^z + − u_r = dz

u

d η ξ

η . (13b) The transformed flow boundary conditions at z = 0 are given by

For pervious half-space:

~ =0

p . (13c) For impervious half-space:

~ 0 dz = p

d , (13d)

where u~_r, u~_z and p~ follows the definitions of Eqs.

(10a)-(10c).

The constants C_i(i=1,2,_L,6) of the general solutions can be determined by the transformed half-space