水位下降誘發(fā)含水層–弱透水層1維黏彈性固結(jié)分析

徐 進，楊偉濤，陳 征，王少偉

(1.煙臺大學 土木工程學院，山東 煙臺 264005；2.武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢 430072)

地面沉降是一種多發(fā)的嚴重環(huán)境地質(zhì)問題，主要誘因是過量開采地下水引起的地下水位降深。如何對地下水位變化誘發(fā)的地面沉降災(zāi)害進行正確評估與科學防治一直備受學術(shù)界關(guān)注[1–4]。荷載作用下土體固結(jié)問題一直是土力學的傳統(tǒng)課題，已得到較深入的研究。Gray[5]較早對瞬時荷載作用下雙層地基固結(jié)問題進行了研究。欒茂田等[6]針對層狀飽和土體，運用分離變量法求解了荷載作用下的Terzaghi 1維固結(jié)方程。謝康和[[7]及楊駿[8]等進一步考慮了變荷載作用下的類似問題。馮健雪等[9–10]考慮了連續(xù)排水邊界條件下自重與線性加載時地基土的1維固結(jié)解析解。陳宗基[11]最先將流變模型引入固結(jié)分析中，劉加才等[12]基于廣義Voigt流變模型，給出了荷載作用下雙層黏彈性地基土的1維固結(jié)解析解。

與荷載作用不同，含水層中地下水開采引起的水位下降，會導致含水層和弱透水層發(fā)生層內(nèi)滲流及層間越流，伴隨這一滲流過程中的孔壓和有效應(yīng)力轉(zhuǎn)化會導致土體發(fā)生變形，最終誘發(fā)地面沉降[13]。駱冠勇等[13]給出了含水層降水引起的弱透水層1維彈性沉降固結(jié)度計算公式。Tseng等[14]考慮了土體重力的影響，給出了含水層中水位下降為定值時弱透水層的1維彈性固結(jié)解。Tao等[15]給出含水層水位瞬時下降和線性下降兩種模式下的單層土體1維固結(jié)彈性解。吳浩等[16]考慮了弱透水層黏性土的結(jié)構(gòu)性，推導了含水層降水所引起的弱透水層1維固結(jié)解。最近，Lo等[17]推導了流體通量邊界條件下非飽和土的1維彈性解。可以看出，為了反映地面沉降的真實演化特征，現(xiàn)有解析研究不斷深入，趨向于考慮水位變化下更復(fù)雜的土體變形特性。

地面沉降具有持續(xù)時間長的特點。大量現(xiàn)場觀測和試驗研究表明，控制地下水開采以后，地下水位得到穩(wěn)定甚至回升，地層變形仍會長期發(fā)展。一般認為，產(chǎn)生這種變形滯后現(xiàn)象的重要原因之一在于弱透水層黏性土具有的流變性[18–20]。Liu等[21]考慮弱透水層土體的黏滯性，并給出了水位變化下的1維固結(jié)半解析解。Li等[22]在此基礎(chǔ)上考慮了弱透水層中流體為非達西流，推導了水位變化下的單層地基土的1維固結(jié)黏彈性解。同時，近年來進一步的研究發(fā)現(xiàn)，由于顆粒間的滑移、錯動以及破碎等細觀機理，含水層砂土同樣表現(xiàn)出一定的流變性，當層厚較大時，含水層的這種蠕變變形將不可忽略[23–26]。因此，給出水位下降誘發(fā)含水層–弱透水層雙層系統(tǒng)的1維固結(jié)模型，在其中考慮含水層和弱透水層的流變性，對于地面沉降的合理評估具有一定理論意義。

為此，針對地下水位變化引起的土層固結(jié)問題，用廣義開爾文模型描述飽和弱透水層黏土與飽和含水層砂土的黏滯性，建立固結(jié)變形控制方程，給出求解條件，經(jīng)過Laplace變換等數(shù)學推導給出傳統(tǒng)矩陣傳遞法和邊界轉(zhuǎn)換法兩種計算方法，求得水位下降誘發(fā)弱透水層–含水層1維固結(jié)的半解析解。采用Stehfest數(shù)值方法實現(xiàn)Laplace逆變換計算，驗證本文計算方法的正確性。通過算例分析傳統(tǒng)矩陣傳遞法和邊界轉(zhuǎn)換法的區(qū)別，進一步探討土體黏滯性、滲透性等力學性質(zhì)與水文地質(zhì)性質(zhì)差異大小以及水位下降速率對土層變形的影響。

1 數(shù)學模型

水位變化引起的弱透水層–含水層1維固結(jié)模型如圖1所示。

圖1 水位變化下弱透水–含水層固結(jié)示意圖Fig. 1 Consolidation schematic diagram of aquitard–aquifer due to groundwater level variation

假定弱透水層頂部與含水層底部總水頭相等且為h，弱透水層頂部有足夠補給，其水頭始終保持不變，在含水層底面t時刻的水位變化為 ?h(t)。 用參數(shù)n表示土層編號（n=1,2 ），每層土厚度依次是H1和H2，且滿足總厚度為H=H1+H2；kv1、kv2分別為兩土層的滲透系數(shù)。

1.1 控制方程

根據(jù)達西定律和有效應(yīng)力原理，可得兩層土體的控制方程為[21]：

式中：kvn為土體滲透系數(shù)，其中n=1,2 ；hn為 第n層土的超孔隙水頭高度；εzn為 豎向應(yīng)變；z和t分別為土層離地面的距離和時間。

土體的黏滯性用廣義開爾文模型描述，結(jié)構(gòu)如圖2所示，由一系列虎克彈簧和牛頓黏壺組成。該模型的優(yōu)點為可以很好地描述流變性相對較弱的含水層砂土的蠕變特性，同時也可以退化成常用的Merchant模型（結(jié)構(gòu)簡圖見圖3），以描述弱透水層黏性土的流變特性。

圖2 廣義開爾文流變模型Fig. 2 Generalized Kelvin rheological model

圖3 Merchant流變模型Fig. 3 Merchant rheological model

基于廣義開爾文流變模型，豎向應(yīng)變εzn可以表示成如下積分形式：

式中，γw為單位體積水的重度， δn(0)為流變模型的柔度系數(shù) δn(t) 在t=0 時 刻的值。δn(t)表達式為：

式中， η1n=E1n/K1n， η2n=E2n/K2n，其中E0n、E1n、E2n、K1n、K2n為廣義開爾文體流變模型元件的參數(shù)（圖2）。

將式（3）代入式（2），再將式（2）代入式（1），則可以得到基于廣義開爾文流變模型的固結(jié)控制方程為：

1.2 定解條件

在兩層土體交界面（即z=H1）處，應(yīng)滿足水頭和流量連續(xù)條件：

上、下邊界條件滿足：

初始條件為：

2 問題求解

2.1 矩陣傳遞法

矩陣傳遞法是利用Laplace變換，對待求微分方程與定解條件實現(xiàn)轉(zhuǎn)換后，再根據(jù)定解條件求解方程得到變量的矩陣表達形式，最后通過矩陣運算法及Laplace逆變換求得解的一種經(jīng)典解析方法[27]。詳細求解過程如下：

首先，進行如下變量代換：

將式（10）代入控制方程（4）及求解條件（5）～（9）得：

根據(jù)初始條件（16），對控制方程（11）及邊界條件（12）～（15）關(guān)于時間t進行Laplace變換得：

對于控制方程（17）進行求解得：

將邊界條件（18）～（21）代入通解式（22）并整理成矩陣形式得：

式中：

從而可得超靜孔壓：

式中，un(z,t) 為t時刻土層各位置的超孔隙水壓。

2.2 邊界轉(zhuǎn)換法

邊界轉(zhuǎn)換法是由Chen等[28]提出的一種新的邊界處理方法，其最大的優(yōu)點在于可將復(fù)合邊界轉(zhuǎn)化為統(tǒng)一的單一邊界條件。本文的滲流連續(xù)邊界條件（5）及（6）就屬于Cauchy邊界（即邊界上同時作用第1類及第2類邊界）[29]。采用邊界轉(zhuǎn)換法將滲流連續(xù)邊界轉(zhuǎn)換為單一邊界進行求解。

根據(jù)邊界轉(zhuǎn)換法，將滲流連續(xù)邊界轉(zhuǎn)換為第1類邊界，并結(jié)合土層頂面及底面邊界可綜合表示為：

基于邊界式（26），控制方程式（17）的通解為：

將通解式（27）代入土層頂面邊界（20）、流量連續(xù)條件（19）及底面邊界（21），并整理成矩陣形式：

式中：

2.3 沉降和固結(jié)度計算公式

根據(jù)式（2）可知，任意時刻雙層土的總沉降為：

由式（29）可知，任意時刻下按沉降定義的雙層土體平均固結(jié)度可表示為：

3 兩種解法的驗證與對比

為了得到定解問題在真實物理空間的解，需要進行Laplace逆變換，本文選用Stehfest法[30]進行數(shù)值逆變換計算。利用MATLAB分別編寫了計算程序，為驗證本文計算方法和程序的正確性，進行了算例對比驗證。

3.1 單層地基下本文解法的驗證

Tao等[15]給出了水頭瞬時下降時弱透水層變形的彈性解，如下：

式中，S(t)為t時 刻下土體的沉降量，Tv為土體的固結(jié)時間因數(shù)，h0與hd分別為初始水頭高度及抽水后含水層的水頭高度，H為弱透水層厚度，N=(2n?1)（n=1,2,3,···），。

將本文解退化到這種單層情形，利用計算程序得出水頭瞬時下降時弱透水層變形的彈性解，與文獻[15]進行對比，見圖4。算例中參數(shù)取值如下：H1=H2=5m，kv1=kv2=8.64×10?4m/d，E01=E02=2MPa，E11=E12=5MPa， η11=η12近 似取無窮大（如1 ×1010d?1），初始水頭h=12m ，瞬時水位變化 ?h(t)=10m，hd=2m。由圖4可以看出：不考慮黏滯性時，本文兩種計算方法所得的解與文獻[15]所得的解析解結(jié)果均吻合良好，在變形初期（ 25d之前）邊界轉(zhuǎn)換法所得結(jié)果略優(yōu)于矩陣傳遞法。這是由于在數(shù)值逆變換中，當時間t較小時會出現(xiàn)計算結(jié)果溢出，導致計算精度降低甚至出錯（如當t取 到1 ×10?6d），邊界轉(zhuǎn)換法在求解中的未知參數(shù)少于矩陣傳遞法，所以在 25d之前邊界轉(zhuǎn)換法計算結(jié)果優(yōu)于矩陣傳遞法。

圖4 單層地基土本文解與文獻[15]所得彈性解的沉降對比Fig. 4 Comparison of Settlement for single-layer foundation obtained from solutions by the proposed methods and reference [15] method

此外，Liu等[21]基于Merchant流變模型給出了水位變換引起單層土體1維固結(jié)的半解析解。為了便于對比，將本文廣義開爾文模型退化為Merchant模型，雙層解退化成單層解，利用編寫的計算程序，分別計算含水層水頭瞬時下降，在t=100d 和t=500d時弱透水層變形的黏彈性解，將本文解與Liu等[21]的解對比，如圖5所示。算例中 η11=η12= 8 .64×10?4d?1，其余參數(shù)保持不變?？梢钥闯?，考慮黏滯性時兩種情形下本文計算結(jié)果都與Liu等[21]的解吻合良好，而且矩陣傳遞法和邊界轉(zhuǎn)換法都表現(xiàn)出良好的一致性。

圖5 單層地基土本文解與文獻[21]所得解的孔壓對比Fig. 5 Comparison of pore pressures for single-layer foundation obtained from solutions by the proposed methods and reference [21] method

3.2 雙層地基下本文解法的驗證

上述兩個退化算例驗證了兩種特殊情況，原則上，理論計算模型的合理性和適用性仍需要通過現(xiàn)場實測資料或者室內(nèi)試驗結(jié)果加以驗證?，F(xiàn)場的水文地質(zhì)條件和邊界條件復(fù)雜，無法完全確定實測結(jié)果的控制因素。相對而言，室內(nèi)模型試驗可以對試驗條件進行人為控制和簡化，去除了很多現(xiàn)場實測中的不確定因素，所得的結(jié)果更有助于驗證理論解的有效性[31]。

為此，利用本文邊界轉(zhuǎn)換法計算程序?qū)Υ迳剿防傻拇蟊壤孛娉两的Ｐ驮囼炦M行了分析[32]，試驗?zāi)Ｐ腿鐖D6所示。

圖6 地面沉降模型試驗示意圖Fig. 6 Schematic diagram of model experiments on land subsidence

在試驗中，上覆層弱透水層始終保持自由水面，并且水頭在含水層頂板以上1.5 m處；下臥含水砂層水頭從初始水頭h0=1.5m 下 降 ?h=?1m來模擬水位降深，并在整個試驗過程中長期觀測地層的沉降量。

圖7給出了彈性和黏彈性兩種情形下的本文理論解與試驗值的對比結(jié)果。計算參數(shù)如下：上覆弱透水層的滲透系數(shù)kv1=0.00006m/d ，彈性模量Es1=0.48MPa ；下層砂土的滲透系數(shù)kv2=1.12m/d，彈性模量Es2=20MPa；黏彈性解中，弱透水層的黏滯系數(shù)η11=0.01d?1，E0=E1=0.48MPa，含水層砂土參數(shù)不變。圖7可以看出：前期（含水層抽水開始 40d以內(nèi)）試驗值與彈性解、黏彈性計算值均比較吻合，但是，隨著時間增長，彈性解越來越偏離試驗值；相比而言，黏彈性解在整體上與 2 00d的試驗結(jié)果吻合更加良好。圖7中同時給出了該雙層地基問題的黏彈性半解析數(shù)值模擬結(jié)果[4]，由圖7可以看出本文計算解與該數(shù)值結(jié)果也保持了較好的一致性，體現(xiàn)了水位下降引起的含水層與弱透水層長期變形特性。

圖7 彈性和黏彈性情形下本文解與半解析數(shù)值解以及試驗值的驗證Fig. 7 Verification of the present solutions with the semianalytical numerical solutions and experimental values for elastic and viscoelastic cases

3.3 本文兩種計算方法的對比

矩陣傳遞法是從固體力學規(guī)則厚板推廣而成的經(jīng)典解析方法，可操作性強，只要列出矩陣方程，便可通過矩陣運算法則進行求解。但是，隨著土層劃分數(shù)N的增加，求解方程（20）時則有2N個未知數(shù)，這會使得計算規(guī)模翻倍，計算效率降低。邊界轉(zhuǎn)換法通過轉(zhuǎn)換邊界，使得層間接觸面水頭的連續(xù)性在解方程之前得到滿足，對于N層土而言，只需N+1個參數(shù)。因此，相比于矩陣傳遞法，邊界轉(zhuǎn)換法在計算工作量方面具有一定優(yōu)勢，而且在處理混合邊界條件時也更加方便。為進一步分析兩種方法的計算效率，通過數(shù)值算例，分別計算分析兩種解法下沉降量的長期發(fā)展過程，并將計算結(jié)果繪制于圖8。

圖8 兩種解法下的沉降值Fig. 8 Settlement values under two solutions

圖8結(jié)果表明，兩種解法得到的沉降值計算精度相近，但是計算效率相差較大，其計算效率矩陣傳遞法比邊界轉(zhuǎn)換法高74.56%。隨著土層數(shù)增加，邊界轉(zhuǎn)換法的計算優(yōu)勢將更加明顯。此外，從圖4中看出，在彈性解中，變形初期（ 25d之前）邊界轉(zhuǎn)換法的精度略高于矩陣傳遞法。

4 影響土體沉降因素分析

在上述邊界轉(zhuǎn)換法計算程序基礎(chǔ)上，通過數(shù)值算例進一步探討了土體的黏滯性、滲透性、層狀土性質(zhì)的差異性以及水頭下降速率對水位變化引起土層變形過程的影響。

算例具體描述如下：某雙層土地基，土層總厚H=12m，H1= 9m，H2=3m ，初始總水頭高度h=14m，t時刻承壓層瞬時降水高度 ?h(t)=12m，土體豎向滲透系數(shù)kv1=9×10?4m/d ，kv2= 3 ×10?4m/d。兩層土均采用Merchant模型：第1層土模型參數(shù)為E01=2000kPa，E11= 5 000kPa ，η11=10?4d?1；第2層土模型參數(shù)為E02=4800kPa，E12=4800kPa ，η12= 1 0?3d?1。

4.1 土體黏滯性的影響

黏滯系數(shù)η1n是 描述土體流變性的主要參數(shù)，其大小反映了土體黏滯性的程度。為了充分探討?zhàn)詫ν馏w長期變形的影響，參考文獻[21]給出了關(guān)于黏滯系數(shù)η1n一 般取值范圍，本文所采用的黏滯系數(shù)η1n取值見表1。算例中其余參數(shù)不變。當η1n接近1010倍時，近似認為η1n趨 于∞，此時流變模型等效于彈性模量為E0E1/(E0+E1)的串聯(lián)彈性模型。用計算程序進行了數(shù)值算例，將計算得到的固結(jié)度–時間關(guān)系繪于圖9。

表1 黏滯系數(shù)的取值Tab. 1 Values of viscosity coefficients

圖9 黏滯系數(shù)對固結(jié)度的影響Fig. 9 Influence of viscosity cofficients on consolidation process

由圖9可以看出：相同時間下，黏滯系數(shù)小的固結(jié)度低，并且，隨著η1的減小，土體的蠕變性越明顯；相同水位變化下土層達到最終沉降所需的時間越長，即完成固結(jié)所需的時間越長，特別是在后期（t=100d之后）這種黏滯性的影響效應(yīng)越明顯。

4.2 滲透性的影響

滲透系數(shù)是描述孔隙水在土體內(nèi)部流動能力的主要參數(shù)，為了探討其對土體變形的影響，對兩層土的滲透系數(shù)kvn同時逐級擴大，取值見表2，其余參數(shù)保持不變。用計算程序進行了數(shù)值算例，將計算得到的滲透系數(shù)影響下，固結(jié)–時間關(guān)系曲線繪于圖10。從圖10中可以看出：在固結(jié)前期（1 00d之前），影響固結(jié)的主要因素為滲透系數(shù)，滲透系數(shù)越大，固結(jié)完成的越快；但在后期該影響趨弱，影響的主要因素是土體的黏滯性。

表2 滲透系數(shù)取值Tab. 2 Values of permeability coefficients

圖10 滲透系數(shù)對固結(jié)度的影響Fig. 10 Influence of permeability coefficients on consolidation process

4.3 土層性質(zhì)差異的影響

土體為非均質(zhì)各向異性體。為了研究土層間力學性質(zhì)差異的影響，進行了數(shù)值算例分析，參數(shù)取值為：土層厚度均取6 m，初始總水頭高度為h=14m，t時刻承壓層瞬時降水高度 ?h(t)=12m，土體豎向滲透系數(shù)kv1=kv2=6×10?4m/d，E01=E02= 2 000kPa，E11=E12=5000kPa，土層的黏滯系數(shù)取值如表3所示。

表3 層狀土黏滯性的差異Tab. 3 Difference of layered soils in viscosity property

為了探討土層間水文地質(zhì)性質(zhì)差異對固結(jié)的影響，在保持其他參數(shù)不變的前提下，土層的黏滯系數(shù)取η11=η12=6×10?4d?1，土層的豎向滲透系數(shù)的取值見表4。土層黏滯系數(shù)和滲透性差異對固結(jié)度的影響如圖11～12所示。

表4 層狀土滲透性的差異Tab. 4 Difference of layered soils in permeability

圖11 土體黏滯性差異對固結(jié)度的影響Fig. 11 Influence of soil viscosity difference on consolidation process

從圖11可以看出，土層黏滯性的差異對固結(jié)前期（ 400d之前）影響較小，只有在后期這種差異性才逐漸表現(xiàn)出來。區(qū)別于黏滯性差異，圖12表明：土層的滲透性差異對于固結(jié)的影響主要表現(xiàn)在前期（400 d之前）；滲透性差異越大，前期固結(jié)度越高，并且在100d之后固結(jié)速率開始逐漸變慢；在固結(jié)后期，滲透性差異的影響趨弱。

圖12 土體滲透性差異對固結(jié)度的影響Fig. 12 Influence of soil permeability difference on consolidation process

4.4 水位下降速率的影響

以往研究大多假設(shè)抽水作用下含水層的水位下降是瞬時完成的，而實際情形中，地下水水位的降深不可能瞬時發(fā)生，應(yīng)該存在一個漸變的過程[13,16]。Liu等[21]以指數(shù)函數(shù)?h(t)=H(1?e?bt)描述抽水作用下水位變化過程，其中，參數(shù)b綜合反映了井的抽水效率和開采含水層的水文地質(zhì)條件，且b＞0，單位為d?1。當參數(shù)b足 夠大（如b＞10000）時，相當于水頭下降瞬時完成這種理想情形。

圖13給出了參數(shù)b取不同值時的固結(jié)度變化曲線，其余參數(shù)取值同上述算例描述。由圖13可以看出：當參數(shù)b越小，即抽水水位下降速率越慢時，固結(jié)度越??；前期沉降發(fā)生時間與參數(shù)b有 關(guān)，即參數(shù)b越小，地面沉降越遲發(fā)生，最終沉降量不受水位下降速率的影響。

圖13 參數(shù)b對固結(jié)度的影響Fig. 13 Influence of parameter b on consolidation process

5 結(jié) 論

1）利用廣義開爾文流變模型和1維固結(jié)理論建立了水位變化下弱透水層–含水層越流系統(tǒng)土體變形的控制方程。通過數(shù)學推導給出了傳統(tǒng)矩陣傳遞法和邊界轉(zhuǎn)換法兩種解法，得到了該固結(jié)方程的半解析解。編制了計算程序，通過與已有的單層解、雙層室內(nèi)試驗數(shù)據(jù)、數(shù)值解進行對比，驗證了本文兩種解法的正確性。

2）對比兩種解法發(fā)現(xiàn)，固結(jié)初期，同種數(shù)值逆變下邊界轉(zhuǎn)換法計算精度略高于傳統(tǒng)矩陣傳遞法。邊界轉(zhuǎn)換法可以將復(fù)雜的混合邊界處理為統(tǒng)一的單一邊界條件，更適合于雙層地基土固結(jié)模型中層間連續(xù)性條件的處理，當土層數(shù)較多時，邊界轉(zhuǎn)換法更加適宜。

3）基于邊界轉(zhuǎn)換法編寫的計算程序進行了數(shù)值算例。結(jié)果表明：土體的蠕變性越明顯，完成固結(jié)所需時間越長，土層變形越滯后，而且土層力學性質(zhì)差異對后期固結(jié)影響較大；相反地，水文地質(zhì)參數(shù)及其差異主要影響中前期固結(jié)過程，土層滲透性越強，土體固結(jié)度越高。值得注意的是，抽水作用下水位下降速率對地面沉降過程也有較大影響，可能會導致土層變形發(fā)生的延遲。