行政院國家科學委員會補助 大專學生參與專題研究計畫研究成果報告

行政院國家科學委員會補助

大專學生參與專題研究計畫研究成果報告

* ********* ********************************************** *

*

*

*

計 畫 ：

名 稱

抽水所引致壓密沉陷之數學模式的建立與解析

*

*

*

* ********* ********************************************** *

執行計畫學生： 陳馨

學生計畫編號： NSC 97-2815-C-216-003-E

研 究 期 間 ： 97 年 07 月 01 日至 98 年 02 月 28 日止，計 8 個月 指 導 教 授 ： 呂志宗

處理方式： 本計畫可公開查詢

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

中華民國 98 年 03 月 31 日

行政院國家科學委員會補助

大專學生參與專題研究計畫研究成果報告

＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

＊ 計畫 ＊

＊ ：

抽水所引致壓密沉陷之數學模式的建立與解析

＊

＊ 名稱 ＊

＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

執行計畫學生：陳 馨

學生計畫編號：NSC 97-2815-C-216-003-E

研 究 期 間 ：97 年 7 月 1 日至 98 年 2 月底止，計 8 個月 指 導 教 授 ：呂志宗

執 行 單 位：中華大學土木與工程資訊學系 中華民國 九十八 年 三 月 三十一 日

摘 要

本計畫之研究重點有三個，包括：(1)多孔介質三維壓密理論的研討與整理，

(2)抽水所引致壓密沉陷問題之數學模式的建立，(3)以積分轉換方法解析出壓密 沉陷問題的閉合解。根據文獻研討知，目前常應用之多孔介質彈力理論有三種，

其主要差異為所引用之基本力學常數並不完全相同，但均是耦合的多孔介質彈力 理論。應用本研究所研討出之結果，則不論引用何種理論模式所建立之數學模 式，均可利用所建立之基本關係式，得出數值研討所需之力學常數值，明瞭各組 基本力學常數之異同，並順利推動相關之研究工作。

本研究於建立數學模式時是以單點抽水模擬單井抽水現象，並將含水層模擬 為均質之線彈性多孔介質，其中地表之滲流邊界條件是分別模擬為透水暨不透水 條件，而抽水型態則分別考慮為以穩定速率抽水及瞬時抽水兩種情況。本研究是 引用Laplace、Hankel 與 Fourier 積分轉換方法解析所建立之數學模式，分別研討 出以上所述抽水及滲流邊界條件下之閉合解。根據本研究之探討得知以下重要研 究成果：

(1) 在穩態抽水的考量下，隨著抽水延時的增加，地層中各點之超額孔隙水的負 壓會逐漸增加，最後形成一穩態平衡。因孔隙水的負壓增加時，會形成地層 有效應力的增加，而導致地層的壓密沉陷量逐漸變大。在瞬時抽水情況下，

地層的壓密沉陷量在一開始抽水時會達到極大值，因本研究是將含水層模擬 為線彈性之多孔介質，故瞬時抽水所引致之沉陷量會逐漸消失。

(2) 在穩態抽水且地表模擬為透水條件的考量下，所引致之地表最大水平位移與

垂直沉陷量均與抽水深度無關，且地表最大水平位移發生位置與黃金比例φ

有關。由數值結果之研討得知，地表最大水平位移約為地表最大沉陷量之 30.0%。另外，在瞬時抽水條件的考慮下，其所引起的地表最大水平位移約為

所對應之地表最大沉陷量的 38.5%。由此可知，抽水所引致之地表水平位移

相當顯著，不宜忽略。

(3) 地表模擬為不透水時，抽水所引致之壓密沉陷量明顯較大。在穩態抽水條件 的考慮下，地表模擬為不透水時抽水所引致之壓密沉陷量約是地表模擬為透 水時之三倍。由此可知，地表滲流邊界條件對抽水所引起的壓密沉陷之影響 相當顯著。

關鍵詞：單點抽水、半無限域、積分轉換、閉合解。

目 錄

摘要 ……… i

目錄 ……… ii

圖目錄 ……… iii

表目錄 ……… iv

第一章 基本方程式研討 ……… 1

1-1 前言 ………...……… 1

1-2 第一種類型的多孔介質彈性力學理論 ……… 2

1-3 第二種類型的多孔介質彈性力學理論 ……… 4

1-4 第三種類型的多孔介質彈性力學理論 ……… 5

1-5 結語 ……….… 7

第二章 抽水引致沉陷之數學模式的建立與其解析 ………. 8

2-1 前言 ………...……… 8

2-2 數學模式 …...……… 9

2-3 數學模式的解析……… 11

2-3-1 對變數 t 作 Laplace 積分轉換 ……… 11

2-3-2 對變數 r 作 Hankel 積分轉換 ……… 13

2-3-3 齊性解的解析 ……….… 15

2-3-4 非齊性解的解析……….….. 16

2-3-5

(

^{z; ,ξ s}

)

定義域之特解 ……….… 19

2-3-6 實數域

(

^{r z t, ,}

)

之特解 ……….… 24

2-4 數值結果 ……….… 27

2-5 結語 ……….… 37

參考文獻 ……….………… 38

誌謝 ……….………… 45

符號說明 ……….………… 46

圖 目 錄

圖 2.1 半無限域中之單點抽水示意圖 ……… 9 圖2.2 以穩定速率抽水且地表模擬為透水情況時於 r = 0 位置之平均壓密度

U ……….... 30 圖2.3 瞬時抽水且地表模擬為透水情況時於 r = 0、h、2h、5h 位置之平均壓密 度 U ……….. 31

圖2.4 以穩定速率抽水且地表模擬為透水情況時之無因次化壓密沉陷曲

線 ………. 31

圖2.5 以穩定速率抽水且地表模擬為透水情況時之無因次化地表水平位移曲

線 …... 32 圖 2.6 瞬時抽水且地表模擬為透水情況時之無因次化壓密沉陷曲線 ...…... 32 圖 2.7 瞬時抽水且地表模擬為透水情況時之無因次化地表水平位移曲線 .. 33 圖 2.8 穩定速率抽水時透水暨不透水地表面之無因次化壓密沉陷曲線 ..… 34

圖2.9 穩定速率抽水且地表模擬為透水時之無因次化超額孔隙水壓分佈

(

^{, ,}

) [

c w ⁴

]

p r z t Qγ πkh ………....… 35

圖2.10 穩定速率抽水且地表模擬為不透水時之無因次化超額孔隙水壓分佈

(

^{, ,}

) [

c w ⁴

]

p r z t Qγ πkh ………...… 36

表 目 錄

表 1.1 各組多孔介質力學常數關係表 ……….… 6 表 2.1 飽和含水層之基本參數 ………..……….… 37

第一章 基本方程式研討

1-1 前言

Biot 之多孔介質彈性力學理論有許多不同的型式，其主要差異在於所引用

之基本力學常數有所不同，本章擬探討出不同 Biot 理論模式所引用力學常數彼

此間之關係。由研討結果知，若考慮線彈性多孔介質為飽和之均向性介質，且其

中之固體與孔隙流體均為可壓縮，則目前 之 Biot 三維多孔介質彈力理論模式常

被引用的基本力學常數有四組，分別為(µ,λ ,Q ,α )、(µ,λ ,Q , R )、(µ,λ , M ,α ) 與(µ,λ , B ,ν_u)等。本章擬分別研討出各組參數彼此間之關係，使有助於獲取多 孔介質彈性力學理論中基本力學常數之參數值，並建立問題之基本方程式。

合理之三維多孔介質彈性力學理論（Poroelasticity）最早是由 Biot[1]所提出 的，此理論模式是考慮多孔介質為均向性之線彈性介質，並引用基本力學常數 (µ,λ , Q ,α )研討多孔介質之基本方程式，其中µ與λ 分別為多孔介質之剪力係 數（Shear Modulus）與 Lame 常數；Q 與α 則為固體介質與孔隙流體間之互制 作用係數。理論模式中之基本力學常數至今仍被文獻加以引用[2, 3]。

為推廣理論模式之應用範圍，在後續的研究中，Biot[4]曾進一步探討異向性

（Anisotropy）介質情況下之多孔介質彈力理論。Biot 於文獻[4]中亦曾研討均向 性介質情況下之多孔介質彈力理論，但所引用之基本力學常數(µ,λ ,Q , R )與文 獻[1]並不完全相同，其中µ與λ 之定義與文獻[1]相同，但Q 與 R 則為新定義之 固體介質與孔隙流體間的互制作用係數。目前，仍有許多學者引用文獻[4]所建 立之理論模式，研討多孔介質之彈力問題[5, 6]。

為使 Biot 多孔介質彈性力學理論所定義之基本力學常數容易由基本力學試

驗獲得，Biot 在文獻[7]中重新整理多孔介質彈力理論中之應力與應變的組成律 (Constitutive Law)關係式，以使相關理論易於應用。此一型態之理論模式係引用 (µ,λ , M ,α )為基本力學常數，其中µ與λ 之定義仍與文獻[1]相同，M 與α 則為 新定義之固體介質與孔隙流體間的互制作用係數。此一理論模式常被應用於探討 多孔介質之彈性動力學問題[8, 9]，當然亦可應用於研討擬靜態（Quasi-static）之 多孔介質彈力問題[10]。Biot 在文獻[11]中曾對其所提出之各種力學常數仔細加 以研討。

多孔介質彈性力學理論在土壤的壓密（Consolidation）問題上有重要之應 用，為使此一理論模式中固體介質與孔隙流體間的互制係數容易由土壤力學上所 熟悉之係數獲得，Rice 與 Cleary[12]保留µ與λ 這兩個力學常數，再引進孔隙流 體為不排出（Undrained）情況下所測得之柏松比（Poisson's Ratio）ν_u與 Skempton 孔隙流體壓力參數 B[13]，改寫 Biot[1]於1941 年所建立之理論模式，所建立之 基本方程式係以(µ,λ , B ,ν_u)等為基本力學常數。許多文獻[14, 15]在探討擬靜態 之多孔介質彈力問題時，即常引用 Rice 等人[12]所建立之 Biot 多孔介質彈力理 論作解析，因其所使用之基本力學常數較容易由基本土壤力學試驗獲得。

由以上說明知，近來常被引用之 Biot 多孔介質彈力理論模式常以文獻[1, 4, 7, 12]所建立的多孔介質彈力理論為基礎，雖然理論模式所引用之基本力學常數 並不完全相同，但是均屬耦合（Coupled）多孔介質彈力理論模式，各組力學常 數彼此間之關係即為本章研討重點。應用本章所研討出之結果，可順利得出數值 分析過程所需之力學常數值，使相關研究工作順利推動。

1-2 第一種類型的多孔介質彈性力學理論

此為Biot[1]於1941 年所建立之多孔介質彈性力學理論，根據此一理論 模式所建立之基本方程式如以下所示：

(

λ µ+

)

u_{k ki}, +µu_{i kk}, −αp,_i+ =f_i 0， (A.1)

, ,

1 0

kk k k

f

k p u p q

α Q

−γ + + − = ， (A.2) 其中µ與λ分別是多孔介質之剪力係數與Lame 常數；Q 和α 則係固體介質 與孔隙流體間之互制作用係數；u 是多孔介質之位移分量； p 表超額孔隙流_i 體壓力（Excess Pore Fluid Pressure），其以壓力為正；k係多孔介質之滲透 係數（Permeability）；γ_f表孔隙流體之單位重（Unit Weight）； f 是多孔介_i 質之徹體力(Body Force)；q 為孔隙流體之補注變率。

或

此為Biot[7]於1956 年所建立之多孔介質彈性力學理論，根據此一理論 模式所建立之基本方程式如以下所示：

(

λ µ+

)

u_{k ki}, +µu_{i kk}, −αp,_i+ =f_i 0， (B.1)

, ,

1 0

kk k k

f

k p u p q

α M

−γ + + − = ， (B.2) 其中µ與λ分別是多孔介質之剪力係數與Lame 常數；其中之互制力學常數

M 即為 Biot[1]在1941 年所定義之力學常數 Q ，而力學常數α亦與 Biot[1]

於 1941 年所定義之力學常數α 相同。u 是多孔介質之位移分量； p 表超額_i 孔隙流體壓力（Excess Pore Fluid Pressure），其係以壓力為正；k是多孔介 質之滲透係數（Permeability）；γ_f表孔隙流體之單位重（Unit Weight）；f 是_i 多孔介質之徹體力(Body Force)； q 為孔隙流體之補注變率。

式(A.1)-(A.2)與式(B.1)-(B.2)之型式極為相似，其由來與多孔介質之組成律

（Constitutive Law）等有關，說明如下。Biot[1]於 1941 年 所 建 立 之 理 論 模

式 中 ， 於研討多孔介質之組成律關係式時，即有考慮固體介質與孔隙流體間之 互制作用關係。若考慮多孔介質為均向性之飽和線彈性體，則其組成律關係式可 表為：

( )

, , ,

ij uk k ij ui j uj i p ij

τ =λ δ +µ + −α δ ， (1.1a)

,

1 uk k p

θ α= +Q ， (1.1b) 其中τ_ij為作用於多孔介質之總應力（Total Stress）；θ 係單位多孔介質於單位體積 內所增加之孔隙流體體積；u 是多孔介質之位移分量； p 表超額孔隙流體壓力_i

（Excess Pore Fluid Pressure），其係以壓力為正；δ_ij為 Kronecker delta 函數；µ 與λ分別是多孔介質之剪力係數與Lame 常數；Q 和α 則是固體介質與孔隙流體 間之互制作用係數。

若考慮孔隙流體的流動遵守 Darcy 定律，則：

,

i i

f

v k p

= −γ ， (1.2) 其中v 為孔隙流體之流速；_i k係多孔介質之滲透係數（Permeability）；γ_f 表孔隙 流體之單位重（Unit Weight）。將Darcy 定律式(1.2)與組成律關係式(1.1a)與(1.1b) 代入以下有考慮多孔介質徹體力（Body Force） f 與孔隙流體補注變率 q 之力平_i 衡方程式(1.3a)與流量連續方程式(1.3b)中：

, 0

ij j fi

τ + = ， (1.3a)

,

vk k q

θ + = ， (1.3b) 則理論模式之基本方程式可表為：

(

λ µ+

)

u_{k ki}, +µu_{i kk}, −αp,_i+ =f_i 0， (1.4a)

, ,

1 0

kk k k

f

k p u p q

α Q

−γ + + − = 。 (1.4b) 式(1.4a)-(1.4b)即引用(µ ,λ , Q ,α )為基本力學常數，以試驗獲取這些常數 時，必須考慮固體介質與孔隙流體間之互制力學行為變化。其他係數n、k、γ_f 本單元並未將其列為基本力學常數，因其與固體介質和孔隙流體間之力學交互作 用現象無關。

另一組與式(1.4a)-(1.4b)極為相似之基本方程式則是由 Biot[7]於1956 年所建 立的，係引用基本力學常數(µ,λ , M ,α )建立理論模式之基本方程式，其中之互 制力學常數M 即為 Biot[1]在 1941 年所定義之力學常數 Q ，而力學常數α 亦與 Biot[1]於 1941 年所定義之力學常數α 相同。然而 Biot 於文獻[7]中之研討的主

要目的是擬介紹多孔介質之彈性動力學模式，並且仍引用文獻[4]中之固體介質 與孔隙流體具相對位移的觀念建立基本方程式。為方便作比較，本節仍將相關方 程式以式(1.1a)-(1.1b)的方式加以表達。基於此，線彈性飽和均向性多孔介質的 組成律關係式可表為[7]：

( )

, , ,

ij uk k ij ui j uj i p ij

τ =λ δ +µ + −α δ ， (1.5a)

,

1 uk k p

θ α= +M ， (1.5b) 其中單位多孔介質體積內所增加之孔隙流體體積θ在文獻[7]中係以ζ 表示；其他 各項物理量或參數之符號定義仍與式(1.1a)-(1.1b)相同。

所建立之基本方程式亦可由將式(1.5a)、式(1.5b)與 Darcy 定律式(1.2)代入力 平衡方程式(1.3a)與流量連續方程式(1.3b)中獲得，其結果如以下所示：

(

λ µ+

)

u_{k ki}, +µu_{i kk}, −αp,_i+ =f_i 0， (1.6a)

, ,

1 0

kk k k

f

k p u p q

α M

−γ + + − = 。 (1.6b)

式(1.6a)-(1.6b)即是引用(µ,λ,M ,α )為基本力學常數所建立之基本方程式。

1-3 第二種類型的多孔介質彈性力學理論

此為Biot[4]於1955 年所建立之多孔介質彈性力學理論，根據此一理論 模式所建立之基本方程式如以下所示：

( )

_{k ki}, _{i kk}, Q R ,_{i i} 0

u u n p f

λ µ+ +µ − R⁺ + = ， (C.1)

2

,_{kk k k}, 0

f

k Q R n

p n u p q

R R

γ

− + + + − = ， (C.2)

其中µ與λ分別是多孔介質之剪力係數與Lame 常數；Q 與 R 分別為 Biot[4]

新定義之固體介質與孔隙流體間的互制作用係數，亦即這裏的力學常數Q 與

式(1.1b)所定義之力學常數Q 並不相同。u 是多孔介質之位移分量；_i p 表超 額孔隙流體壓力（Excess Pore Fluid Pressure），其以壓力為正；n與k係多孔 介質之孔隙率（Porosity）與滲透係數（Permeability）；γ_f表孔隙流體之單位 重（Unit Weight）； f 是多孔介質之徹體力(Body Force)； q 為孔隙流體之補_i 注變率。

式(C.1)與式(C.2)之由來與多孔介質之組成律（Constitutive Law）有關，說明 如下。Biot[4]在此引進孔隙率（Porosity）n的觀念說明作用於多孔介質之孔隙

流體壓力變化，再以孔隙流體與固體介質會產生相對位移的現象分析多孔介質之 位移變化，並將其理論模式推廣至異向性介質情況。本單元將僅研討均向性介質 之理論模式，並將相關方程式整理成類似前一節所述之型式，以方便作比較。

基於以上之說明，考慮飽和多孔介質為線彈性體，則 Biot 所建立之多孔介 質組成律關係式如以下所示[4]：

( )

, , ,

ij k k ij i j j i ij

u u u nQ R p

τ =λ δ +µ + − R⁺ δ ， (1.7a)

2 , k k

Q R n

n u p

R R

θ = ⁺ + ， (1.7b)

式中τ_ij、θ 、u 、 p 、_i δ_ij、µ、λ 之定義仍和式(1.1a)與式(1.1b)相同；Q 與 R 分

別為 Biot[4]新定義之固體介質與孔隙流體間的互制作用係數，亦即這裏的力學

常數Q 與式(1.1b)所定義之力學常數 Q 並不相同。上式中所引用之基本力學常數 (µ,λ, Q , R )與文獻[4]中曾出現之力學常數A 、N、S的關係為：

µ =N， (1.8a) Q2

S A

λ = = − R 。 (1.8b) 將式(1.7a)、式(1.7b)與 Darcy 定律式(1.2)代入力平衡方程式(1.3a)與流量連續方 程式(1.3b)中，則理論模式之基本方程式可表為：

( )

_{k ki}, _{i kk}, Q R ,_{i i} 0

u u n p f

λ µ+ +µ − R⁺ + = ， (1.9a)

2

,_{kk k k}, 0

f

k Q R n

p n u p q

R R

γ

− + + + − = ， (1.9b)

上式係引用基本力學常數(µ,λ, Q , R )所建立之基本方程式。以試驗取得其他係 數n、k、γ_f 時，不必考慮固體介質與孔隙流體間之交互作用現象。

1-4 第三種類型的多孔介質彈性力學理論

近年來，由於 Rice 與 Cleary[12]所引用之多孔介質力學常數容易由基本力

學試驗獲得，因此，其所建立之 Biot 多孔介質彈性力學理論已漸被廣泛引用。

其理論模式中，引進流體不排出情況下所測得之柏松比ν_u與Skemptom 孔隙流體

壓力參數 B[13]，改寫 Biot[1]於1941 年所建立之理論模式，所建立之基本方程

式係以(µ,λ, B ,ν_u)等為基本力學常數。基於此，飽和均向性之線彈性多孔介質 的組成律關係式可表為[12]：

( ) ₍ ⁽ ₎₍ ⁾ ₎

, , ,

3 1 2 1

u

ij k k ij i j j i ij

u

u u u p

B

τ λ δ µ ν ν δ

ν ν

= + + − −

− + ， (1.10a)

( )

( )( ) ( )( )

( )( )

, 2 2

3 9 1 2

1 2 1 2 1 2 1

u u u

k k

u u

u p

B B

ν ν ν ν ν

θ ν ν µ ν ν

− − −

= +

− + − + ， (1.10b) 其中θ在文獻[12]中仍以ζ 表示；ν 為考慮流體係排出情況下（Drained）所測得 之多孔介質柏松比，其與力學常數µ、λ的關係為：

( )

2 ν λ

= λ µ

+ 。 (1.11) 將式(1.10a)、式(1.10b)與 Darcy 定律式(1.2)代入力平衡方程式(1.3a)與流量 連續方程式(1.3b)中，則多孔介質彈性力學理論之基本方程式可表為：

( ) ( )

( )( )

, , ,

3 0

1 2 1

u

k ki i kk i i

u

u u p f

B λ µ µ ν ν

ν ν

+ + − − + =

− + ， (1.12a)

( )

( )( ) ( )( )

( )( )

, , 2 2

3 9 1 2

1 2 1 2 1 2 1 0

u u u

kk k k

f u u

k p u p q

B B

ν ν ν ν ν

γ ν ν µ ν ν

− − −

− + + − =

− + − + ， (1.12b)

其中(µ,λ , B ,ν_u)等力學常數最易經由試驗取得，許多文獻均有列示不同種類多 孔介質之(µ,λ , B ,ν_u)值。因此，以(µ,λ , B ,ν_u)為基本力學常數所建立之多孔介 質彈力理論之基本方程式漸漸廣被引用。為清楚起見，本單元擬利用如下所建立

之表1.1 列示各組力學常數彼此間之關係，使相關研究人員容易獲得所需之基本

力學常數，並熟悉各組基本力學常數之異同。

表1.1 各組多孔介質力學常數關係表

(µ,λ, Q ,α)^[1] (µ,λ, Q , R )^[4] (µ,λ, M ,α)^[7] (µ,λ, B ,ν_u)^[12]

(µ,λ, Q ,α)^[1] − Q^[4]=nQ(α−n)^[1]

[4] 2 [1]

R =n Q M=Q ( )

2

2 2 u

Q Q ν λ α

λ µ α

= + + +

2

3

3 2 3

B Q

Q α

λ µ α

= + +

(µ,λ, Q , R )^[4]

[4]

[1] nQ R

α = R⁺

[4]

[1]

2

Q R

=n

−

nQ R α= R⁺

2

M R

=n

( )

( )

2

2

2

u

Q R R Q R

R ν λ

λ µ + +

=  + 

+ +

 

 

 

( )

( )²

3

3 2 3

B Q R n Q R

λ µ R

= +

 + + + 

 

 

 

(µ,λ, M ,α)^[7] Q M= Q nM= (α−n)

R n M= 2 − ( )

2 2

M

2 M

u

ν λ α λ µ α

= + + +

2

3 M

3 2 3 M

B α

λ µ α

= + +

(µ,λ, B ,ν_u)^[12]

( )

(1 2³ )(1 )

u

B u

ν ν

α ν ν

= −

− +

( )( )

( )( )

2 2

2 1 2 1

9 1 2

u

u u

Q µB ν ν

ν ν ν

− +

= − −

( )

( )

2 1

3 1 2

u u

Q µnB ν

ν

= +

−

( )( )

( )( )

2 2 2

2 1 2 1

9 1 2

u

u u

µn B ν ν

ν ν ν

− +

− − −

( )( )

( )( )

2 2 2

2 1 2 1

9 1 2

u

u u

R µn B ν ν

ν ν ν

− +

= − −

( )

(1 2³ )(1 )

u

B u

ν ν

α ν ν

= −

− +

( )( )

( )( )

2 2

2 1 2 1

9 1 2

u

u u

M µB ν ν

ν ν ν

− +

= − −

−

1-5 結語

本單元旨在研討 Biot 多孔介質彈力理論之基本力學常數的關係。由文獻研

討知，目前常應用之多孔介質彈力理論有三種，其主要差異為所引用之基本力學 常數並不完全相同，但均是耦合的多孔介質彈力理論。應用本文所研討出之結

果，則不論引用何種理論模式所建立之數學模式，均可利用表1.1 之關係式，得

出數值研討所需之力學常數值，明瞭各組基本力學常數之異同，並順利推動相關 研究工作。

第二章 抽水引致沉陷之數學模式的建立與其解析

2-1 前言

地下水超抽所引起的地層下陷問題，直到如今仍為國人所深刻關切[16, 17]，也是世界各國在經濟發展過程中常見的問題[18-27]。為解決此一問題，以 往國內外已有許多學者專家進行許多相關之研究，近年來國科會亦曾補助許多與 地層下陷議題相關之計畫案[28-41]進行相關之研究。此外，水利署等相關單位亦 持續進行地層下陷區之沉陷監測，並有設置宣導網站[42-44]，進行地層下陷之防 治與宣導，且已具有一定之成效，值得肯定。然而，經研討過相關之文獻後發現，

目前國內較少針對抽水所引致地層下陷問題的解析解進行相關之研究，但是很多 以數值模擬為基礎之相關研究，常需以解析解為基礎，用以研判其數值模擬之分 析結果是否正確無誤。基於此，本計畫在指導老師的指導下，進行相關之研究，

並已獲致重要之研究成果。

Terzaghi[45]首先引用有效應力觀念（Effective Stress Concept）說明土壤的壓 密過程，在其理論模式中，須先解析出平衡孔隙水壓力，然後再間接計算出壓密 沉陷量。然而 Biot[1, 4]所考慮之孔隙水的平衡過程則與土壤固體直接相關，所 得出之結果亦被證實較為合理可靠。Biot[1, 4]所建立之壓密理論模式中，若於質 量守衡方程式中同時考慮孔隙水與固體介質的影響，則所建立之壓密模式稱為耦 合壓密模式；若質量守衡方程式中僅考慮孔隙水的影響，則所建立之壓密模式稱 為非耦合壓密模式。本研究於探討單點抽水所引致的地層超額孔隙水壓分佈時，

擬引用 Biot 之非耦合壓密模式建立基本方程式，再以積分轉換方法解析所建立

之數學模式的解。Booker 與 Carter[46-49]、Tarn 與 Lu[50]、Lu 與 Lin[51-55]、

Lin 與 Lu[56]、Chen[57, 58]、Worsak 與 Chau[59]等人均曾研討過單點抽水所引 致之地層位移與超額孔隙水壓。但以上所述相關文獻之研究過程中，均未曾同時 探討穩定速率抽水與瞬時抽水之暫態壓密沉陷問題的閉合解，故本計畫擬探討其 沉陷行為之差異性；此外，本計畫亦擬探討穩定速率抽水與瞬時抽水所引致的地 層孔隙水壓分佈。

根據各界對抽水所引致地層下陷問題之學理分析發現，抽水時除會引起地層 下陷外，亦會引起地層之水平位移，國內外亦已有許多相關之文獻[60-65]的監測 結果印證此一結論。然而，地層水平位移常被忽略卻是事實，根據本計畫案之指 導老師的研究發現[51-56]，當地表模擬為透水情況時，抽水所引致之最大地表水 平位移並不宜忽略。以學理分析為基礎之相關研究成果，亦已有許多文獻[66-68]

證明抽水所引致之地表水平位移相當顯著。基於以上考量，本計畫所引用之分析 模式即是非單向度之壓密模式；另外，由文獻研討得知，以上各類問題之相關研 究成果，均較少完整討論地表邊界分別模擬為透水暨不透水、地層滲流條件為橫 向等向性時之暫態壓密解析解，這些研究重點均已呈現在本計畫案所建立之數學 模式中，並引用積分轉換等方法加以解析，以研討出單點抽水所引致壓密沉陷之 暫態閉合解。

本單元旨在考慮如圖2.1 所示之單點抽水所引致的地表壓密沉陷，其係將飽

和含水地層模擬為均質之半無限域多孔介質，並以單點抽水模擬單井抽水的現

象，因讓座標z 軸通過抽水點，故問題可進一步化簡為軸對稱情況。此外，其地 表力學邊界條件是考慮為無正向應力與剪應力變化的情況，地表滲流邊界條件則

模擬為透水暨不透水兩種情況；而在無限遠處之邊界上（z→ ∞），則考慮各種

物理變化量均不受單點抽水的影響。本計畫是引用 Laplace、Hankel 與 Fourier 積分轉換方法解析所建立之數學模式，所研討出之解為閉合解（Closed-form Solution）。

圖2.1 半無限域中之單點抽水示意圖 2-2 數學模式

如圖2.1 所示，飽和含水地層係模擬為均質之線彈性多孔介質，但其水平方

向的滲流速度與垂直方向的滲流速度則考慮為不同，並以單點抽水模擬單井抽水

的現象。若考慮座標z 軸通過抽水點，且抽水深度為 h，則抽水問題的非耦合控

制方程式可以軸對稱圓柱座標(r, z)表為[51]：

2

2 0

1 2

r r

u

G p

G u G

r r r

ε ν

∂ ∂

∇ + − − =

− ∂ ∂ ， (2.1a)

2 0

z 1 2

G p

G u z z

ε ν

∂ ∂

∇ + − =

− ∂ ∂ ， (2.1b)

2 2

2 2

1 0

r z

w w

k p p k p p

n q

r r r z β t

γ γ

∂ ∂ ∂ ∂

 

− ∂ + ∂ − ∂ + ∂ + = ， (2.1c)

其中k 、_r k 、n、_z G與ν 為飽和含水地層之水平滲透係數、垂直滲透係數、孔隙 率（Porosity）、剪力模數（Shear Modulus）與柏松比（Poisson’s Ratio）；γ_w與β 是孔隙水之單位重（Unit Weight）與體積壓縮係數（Compressibility）；u 與_r u 為_z 飽和含水層之水平位移與垂直位移；p 則是超額孔隙水壓（Excess Pore Water Pressure）；q 表抽水作用源。

本研究所考慮的抽水作用源q 包含以下兩種情況。

  穩定速率抽水情況

亦即所考慮的抽水作用源q 為：

(

^{, ,}

) ( ) ( ) ( )

2 Qc

q q r z t r z h u t rδ δ

= = π − ， (2.2a) 其中Q 是於單位時間內以穩定速率抽水時所抽出的地下水體積，因此_c Q 的單位_c 是體積除以時間；^δ

( )

^{x 為}Dirac delta 函數；^{u t}

( )

則是單位階梯函數（Unit Step Function 或 Heaviside Step Function）。

  瞬時抽水情況

亦即所考慮的抽水作用源q 為：

(

^{, ,}

)

⁰

( ) ( ) ( )

2

q q r z t Q r z h t

rδ δ δ

= = π − ， (2.2b) 其中Q 是當₀ t→ 時，自飽和含水層瞬間抽出的地下水之體積，也就是說0⁺ Q 的₀ 單位是體積單位，Q 與₀ Q 的單位並不相同。因本單元是考慮飽和含水地層為線_c 彈性之多孔介質，故瞬間抽水所引致之沉陷長期而言會逐漸消失。

如圖2.1 所示之抽水問題中，其地表邊界（z = 0）是考慮為無任何應力變化，

亦即：

( ) (

^{, 0,}

) (

^{, 0,}

)

, 0, ^{r z} 0

rz

u r t u r t

r t G

z r

σ′ = ^^∂ ∂ +^∂ ∂ ^= ， (2.3a)

(

_,0,

)

₂ z

(

^,0,

) ( )

_{2 1} r

(

^,0,

)

r

(

^,0,

)

₀

zz

u r t u r t u r t

r t G G

z r r

σ′ = η ^∂ ∂ + η− ^^∂ ∂ + ^= ， (2.3b)

其中參數^{η = −}

(

^{1 ν}

) (

^{1 2− ν}

)

^；σij′（張力定義為正的量）表作用於飽和含水層之有 效應力（Effective Stress），且σ_ij′ = + ，τ_ij p τ_ij（張力定義為正的量）是作用於飽 和含水層之總應力（Total Stress），超額孔隙水壓p 則將壓力定義為正的量。此 外，z = 0 之地表邊界的滲流性質可分別模擬為透水暨不透水兩種情況。

  若考慮地表邊界完全透水

則地表之滲流邊界條件可表為：

(

^,0,

)

⁰

p r t = 。 (2.3c)   若考慮地表邊界為完全不透水

則由Darcy 定律得知，地表之滲流邊界條件可表為：

(

^,0,

)

₀

p r t z

∂ =

∂ 。 (2.3d) 式(2.3a)-(2.3d)組成半無限域之單點抽水問題於地表（z = 0）力學與滲流的邊界條 件。

此外，可考慮如圖2.1 所示半無限域中之無限深遠處（z→ ∞）的邊界不受 點抽水的影響，亦即，問題之另一組邊界條件可表為：

( )

lim _r , , 0

z u r z t

→∞ = ， (2.4a)

( )

lim _z , , 0

z u r z t

→∞ = ， (2.4b)

( )

lim , , 0

z p r z t

→∞ = 。 (2.4c) 關於問題之初始條件則是考慮所有的物理量在一開始的時候都沒有變化，因 此，問題之初始條件可表為：

( )

lim0 _r , , 0

t u r z t

→ = ， (2.5a)

( )

lim0 _z , , 0

t u r z t

→ = ， (2.5b)

( )

lim0 , , 0

t p r z t

→ = 。 (2.5c) 式(2.1a)-(2.1c)、(2.3a)-(2.3d)、(2.4a)-(2.4c)、(2.5a)-(2.5c)組成單點抽水問題 之數學模式。本研究於解析數學模式之解時，擬採用Laplace、Hankel 與 Fourier 積分轉換方法等，說明如下。

2-3 數學模式的解析

2-3-1 對變數 t 作 Laplace 積分轉換

式(2.1a)-(2.1c)是與自變數r 、 z 、 t 有關之聯立偏微分方程式，此為線性

（Linear）的非齊次（Non-homogeneous）微分方程式。解析問題時，首先引用 式(2.5a)-(2.5c)之初始條件，對式(2.1a)-(2.1c)作 Laplace 積分轉換，基於此，可得 [69]：

( )

²

2 2

2 2 2

1 1 1

2 _r 2 1 U^z P 0

r r r r z U r z G r

η η

  ∂ ∂  ∂  ∂ ∂

+ − + + − − =

 ∂ ∂  ∂  ∂ ∂ ∂

  ， (2.6a)

( )

²2 ²2

1 1 1

2 1 _r 2 _z P 0

U U

z r r r r r z G z

η− ∂ ∂^{∂ ∂} + ^ +^∂^∂ + ∂^∂ + η∂^{∂ } − ^∂∂ = ， (2.6b)

2 2

2 2

1 0

r z

w w

k P P k P

n sP Q

r r r z β

γ γ

∂ ∂ ∂

 

− ∂ + ∂ − ∂ + + = ， (2.6c)

其中符號s是Laplace 積分轉換參數；符號^{η = −}

(

^{1 ν}

) (

^{1 2− ν}

)

^；U 、r U 、_z P 、Q

分別定義如下：

(

^{, ;}

)

₀

(

^{, ,}

)

^st

r r

U r z s =

∫

^∞u r z t e dt⁻ ， (2.7a)

(

^{, ;}

)

₀

(

^{, ,}

)

^st

z z

U r z s =

∫

^∞u r z t e dt⁻ ， (2.7b)

(

^{, ;}

)

₀

(

^{, ,}

)

^st

P r z s =

∫

^∞p r z t e dt⁻ ， (2.7c)

(

^{, ;}

)

₀

(

^{, ,}

)

^st

Q r z s =

∫

^∞q r z t e dt⁻ 。 (2.7d)

上式中符號U 、_r U 、 P 之 Laplace 反轉換則定義為： _z

(

^{, ,}

)

¹

(

^{, ;}

)

2

i st

r i r

u r z t U r z s e ds i

α

π α + ∞

=

∫

− ∞ ， (2.7a*)

(

^{, ,}

)

¹

(

^{, ;}

)

2

i st

z i z

u r z t U r z s e ds i

α

π α + ∞

=

∫

− ∞ ， (2.7b*)

(

^{, ,}

)

¹

(

^{, ;}

)

2

i st

p r z t i P r z s e ds i

α

π α + ∞

=

∫

− ∞ 。 (2.7c*) 式(2.6c)中之抽水作用源 Q 包含以下兩種情況。

  穩定速率抽水情況

( ) ( )

2 Qc

Q r z h

rsδ δ

= π − 。 (2.8a)   瞬時抽水情況

( ) ( )

0

2

Q Q r z h

rδ δ

= π − 。 (2.8b) 另外，地表z=0之邊界條件與無限深遠處z→ ∞之邊界條件可分別表為式 (2.9a)-(2.9d)與(2.10a)-(2.10c)的型式。茲對式(2.3a)-(2.3d)作 Laplace 積分轉換，可 得地表z=0之力學邊界條件如以下所示：

(

^,0;

) (

^,0;

)

₀

r z

U r s U r s

z r

∂ ∂

+ =

∂ ∂ ， (2.9a)

(

^,0;

) ( ) (

^,0;

) (

^,0;

)

1 0

z r r

U r s U r s U r s

z r r

η^∂ ∂ + η− ^^∂ ∂ + ^= 。 (2.9b)

  若考慮地表邊界完全透水

則地表之滲流邊界條件可表為：

(

^,0;

)

⁰

P r s = 。 (2.9c)

  若考慮地表邊界完全不透水

則地表之滲流邊界條件可表為：

(

^,0;

)

₀

P r s z

∂ =

∂ 。 (2.9d) 對式(2.4a)-(2.4c)作 Laplace 積分轉換，可得無限深遠處（z→ ∞）之邊界條 件如以下所示：

( )

lim _r , ; 0

z U r z s

→∞ = ， (2.10a)

( )

lim _z , ; 0

z U r z s

→∞ = ， (2.10b)

( )

lim , ; 0

z P r z s

→∞ = 。 (2.10c) 式(2.9a)-(2.9d)與式(2.10a)-(2.10c)中之符號U 、_r U 、 P 的定義如式(2.7a)-(2.7c)所_z 示。

2-3-2 對變數 r 作 Hankel 積分轉換

經Laplace 積分轉換後，式(2.6a)-(2.6c)、(2.9a)-(2.9d)與(2.10a)-(2.10c)是與自 變數r 、 z 有關之數學模式，解析過程中可繼續對其中之自變數 r 進行 Hankel 積 分轉換。式(2.6a)-(2.6c)分別進行 1 階、0 階與 0 階之 Hankel 積分轉換後可得[70]：

( )

2

2 2

2 _r 2 1 dU^z 1 0

d U P

dz ηξ η ξ dz Gξ

 −  − − + =

 

  ， (2.11a)

(

^{2 1}

)

^{r 2 2}₂ ² z ^{1 0}

dU d dP

dz dz U G dz

η− ξ +^ η −ξ ^ − =

  ， (2.11b)

2 2

2 0

r z

w w

k k d

P n sP Q

ξ dz β

γ γ

 

− + + =

 

  ， (2.11c)

其中符號ξ是Hankel 積分轉換參數；U 、_r U 、 P 、 Q 分別定義為： _z

(

; ,

)

0

(

, ;

) ( )

1

r r

U z ξ s =

∫

^∞rU r z s J ξr dr， (2.12a)

(

; ,

)

0

(

, ;

) ( )

0

z z

U z ξ s =

∫

^∞rU r z s J ξr dr， (2.12b)

(

; ,

)

0

(

, ;

) ( )

0

P z ξ s =

∫

^∞rP r z s J ξr dr ， (2.12c)

(

; ,

)

0

(

, ;

) ( )

0

Q z ξ s =

∫

^∞rQ r z s J ξr dr。 (2.12d) U 、r U 、 P 之 Hankel 反轉換則分別定義為： _z

(

, ;

)

0

(

; ,

) ( )

1

r r

U r z s =

∫

^∞ξU z ξ s J ξr dξ， (2.12a*)

(

, ;

)

0

(

; ,

) ( )

0

z z

U r z s =

∫

^∞ξU z ξ s J ξr dξ ， (2.12b*)

(

, ;

)

0

(

; ,

) ( )

0

P r z s =

∫

^∞ξP z ξ s J ξr dξ 。 (2.12c*)

式(2.11c)中之抽水作用源 Q 包含以下兩種情況。

  穩定速率抽水情況

( )

2 Qc

Q z h

sδ

= π − 。 (2.13a)   瞬時抽水情況

( )

0

2

Q Q δ z h

= π − 。 (2.13b) 另外，再對式(2.9a)-(2.9d)所示之地表z=0的邊界條件中之變數r 分別作 1 階、0 階、0 階、0 階之 Hankel 積分轉換後，可得地表力學邊界條件(2.14a)-(2.14b) 與滲流邊界條件(2.14c)-(2.14d)如以下所示：

(

^{0; ,}

) (

^{0; ,}

)

⁰

r

z

dU s

U s

dz

ξ −ξ ξ = ， (2.14a)

(

^{0; ,}

) ( ) ( )

^{1 0; , 0}

z

r

dU s

U s

dz

η ξ + η− ξ ξ = 。 (2.14b)

  若考慮地表邊界完全透水

則地表之滲流邊界條件可表為：

(

^{0; ,}

)

⁰

P ξ s = 。 (2.14c)   若考慮地表邊界完全不透水

則地表之滲流邊界條件可表為：

(

^{0; ,}

)

dP s 0 dz

ξ = 。 (2.14d)

茲再對式(2.10a)-(2.10c)所示之無限深遠處（z→ ∞）邊界條件中的變數r 分 別作1 階、0 階、0 階之 Hankel 積分轉換後可得：

( )

lim _r ; , 0

z U z ξ s

→∞ = ， (2.15a)

( )

lim _z ; , 0

z U z ξ s

→∞ = ， (2.15b)

( )

lim ; , 0

z P z ξ s

→∞ = 。 (2.15c) 式(2.14a)-(2.14d)與式(2.15a)-(2.15c)中之符號U 、_r U 、P 的定義如式(2.12a)-(2.12c)_z 所示。

2-3-3 齊性解的解析

式(2.11a)-(2.11c) 為 聯 立 之 非 齊 性 常 微 分 方 程 式 ， 其 通 解 包 含 齊 性 解

（Homogeneous Solution）與非齊性解（Non-homogeneous Solution），這兩部分 的解均可依據傳統的微分方程解析技巧推求其解。基於此，式(2.11a)-(2.11c)的齊 性解可假設為：

1

( )

2

3 rh

zh

h

U x

U x exp z

P x

λ

   

   =

   

    

 

， (2.16)

其中

{

1 2 3

}

x x x T是特徵向量（Eigenvector）；λ 為特徵根（Eigenvalue）。將式(2.16) 代入式(2.11a)-(2.11c)之齊性微分方程式部分（不考慮Q 的影響），可得問題之特 徵方程式（Characteristic Equation）如以下所示：

(

^{2 2}

)

^{2 r 2 z 2 0}

w w

k k

λ ξ ξ n sβ λ

γ γ

 

−  + − =

  。 (2.17) 解析式(2.17)可得問題之六個特徵根λ 分別為：

λ = ± 、ξ ± 、ξ ^{r 2}

z

k s

k ξ c

± + ， (2.18)

其中符號c k n= _z βγ_w。根據所推導出之特徵根，可得出問題之齊性解U 、_rh U 、_zh P 如以下所示： h

(

; ,

)

1 2 3 4 5 ² 6 ²

r r

z z

k s k s

z z

k c k c

z z z z

Urh z ξ s = A e^ξ +A ze^ξ +A e^−ξ +A ze^−ξ +A e ^{ξ +} +A e^{− ξ +} ， (2.19a)

(

; ,

)

1 2 3 4 5 ² 6 ²

r r

z z

k s k s

z z

k c k c

z z z z

Uzh z ξ s =B e^ξ +B ze^ξ +B e^−ξ +B ze^−ξ +B e ^{ξ +} +B e^{− ξ +} ， (2.19b)

(

; ,

)

1 2 3 4 5 ² 6 ²

r r

z z

k s k s

z z

k c k c

z z z z

P zh ξ s =C e^ξ +C ze^ξ +C e^−ξ +C ze^−ξ +C e ^{ξ +} +C e^{− ξ +} ， (2.19c) 其中係數A ii

(

=^{1, ,6}

)

、B ii

(

=^{1, ,6}

)

、C ii

(

=^{1, ,6}

)

僅六個是線性獨立的，其 係數間之彼此關係可由式(2.19a)-(2.19c)需滿足式(2.11a)-(2.11c)之聯立齊性微分