，，

(1.江蘇南京地質工程勘察院，江蘇 南京 210041; 2.江蘇省地質礦產勘查局，江蘇 南京 210018)

0 引 言

地下水開采引起地下水滲流場變化和土體自重固結壓縮引起的應變場變化相互疊加，致使區(qū)域地面沉降愈來愈嚴重。國內外學者對區(qū)域地面沉降的研究起步很早，但更多的是對基于抽水引發(fā)的地面沉降影響研究(陳希崇等，2001；駱祖江等，2008；Luo et al,2011)，對土體自重固結壓縮引起的地面沉降研究尚不多見。

為了準確評價土體自重固結壓縮對地面沉降的影響，以比奧固結理論為基礎，結合土體非線性流變理論，將比奧固結理論中的本構關系推廣到黏彈塑性，同時考慮水力學參數及土力學參數隨有效應力的動態(tài)變化過程。以江蘇南通為例，建立南通市土體自重固結壓縮、地下水開采與地面沉降三維全耦合模型，模擬預測了南通在地下水停止開采、僅在土體自重固結壓縮影響下的地面沉降發(fā)生、發(fā)展趨勢，為南通市地面沉降科學防控提供決策依據。

1 研究區(qū)概況

研究區(qū)域包括南通市區(qū)(崇川、港閘)和通州區(qū)。全區(qū)第四紀松散沉積物發(fā)育，分布范圍廣，最大厚度達300 m左右，含水層以第四紀松散層孔隙水為主。垂向上從上往下分別為潛水含水層、第I承壓含水層、第II承壓含水層、第III承壓含水層以及各含水層間的黏性土弱含水層。其中，潛水含水層由第四系全新統(tǒng)地層組成，巖性以粉砂、亞砂土為主，潛水位埋深1～3 m，厚度20～30 m；第I承壓含水層由上更新統(tǒng)沖積、沖海積松散地層組成，巖性以中細砂為主，頂板埋深50～60 m，承壓水頭埋深1～3 m；第Ⅱ承壓含水層由中更新統(tǒng)河流、河口相沉積的松散層組成，巖性以粉細砂、中粗砂及砂礫為主，頂板埋深約140 m，承壓水頭埋深3～5 m；第III承壓含水層由下更新統(tǒng)長江古河道沉積的松散層組成，巖性以中細砂、中粗砂為主，局部有含礫卵石，頂板埋深187～270 m，承壓水頭埋深30～40 m。

2 地質概念模型

研究區(qū)平面上以行政區(qū)域為界，面積為1 586 km2，垂向深度在230～350 m。垂向上，從上往下將潛水含水層、第I承壓含水層、第II承壓含水層、第III承壓含水層及各含水層之間的黏性土弱含水層剖分成獨立的7個層位參與計算，各含水層在水力學上均概化為非均質各向異性，在土力學上均概化為非均質地層。研究區(qū)南部是長江邊界，其從上往下切割潛水含水層、第I承壓含水層及之間的黏性土弱含水層，故在水力學上可概化上述含水層的長江邊界為定水頭邊界，其余在水力學上均概化為第二類流量邊界，系統(tǒng)的頂部一方面接受大氣降水及農業(yè)灌溉入滲的補給，是補給邊界；另一方面，地下水又通過其蒸發(fā)，是排泄邊界，系統(tǒng)的底部為隔水邊界。研究區(qū)所有側向邊界和底部邊界在土力學上均概化為零位移邊界。研究區(qū)地下水開采是唯一的源匯項，將其概化為大井處理。土體在自重應力作用下發(fā)生固結壓縮。

3 數學模型的建立與求解

3.1 比奧固結理論

飽和土體中假定土骨架變形為線彈性、微小變形、滲流符合達西定律、水不可壓縮或微壓縮的三維比奧固結方程如下(錢家歡等，1996)。

(1)

(2)

式(1)、式(2)中，G為剪切模量；ν為泊松比；wx、wy、wz分別為x、y、z方向上的位移分量；u為孔隙水壓力；kx、ky、kz分別為主軸x、y、z方向上的滲透系數；γ為土的重度；γw為水的重度。

3.2 土體本構模型

土的本構關系是土的力學特性即應力-應變-強度-時間等關系的數學表達式。對于考慮流變特性的土體來說，其變形特征主要表現為變形的時間與應力水平有關，所顯示的是具有彈性、塑性和黏滯性的黏彈塑性體，若將此類土體的總應變增量dε分為彈塑性應變增量dεep、黏彈性應變增量dεve、黏塑性應變增量dεvp，則具有流變特性的土體中任意點在任意時刻的應變增量為(維亞諾夫，1987；錢家歡等，1996；陳曉平等，2001；姚仰平等，2009)：

dε=dεep+dεve+dεvp

(3)

式(3)中各部分的應變增量可以由下述方法確定。

3.2.1 彈塑性應變增量 由土體的彈塑性本構模型可得：

dεep=[C]dσ′

(4)

3.2.2 黏彈性應變增量 由Kelvin流變模型E1ε+Keε=σ可得應力不變時復雜應力狀態(tài)下的黏彈性應變增量為：

(5)

式(5)中，ηe=E1/Ke，E1為Kelvin體黏彈性模量；Ke為黏滯系數；[A]為應力矩陣，

3.2.3 黏塑性應變增量 利用黏塑性法確定黏塑性應變增量，在該種方法中允許材料在有限“期間”內超越破壞準則(以破壞準則函數F的值>0來表示)。在討論土體的黏塑性應變而非塑性應變時，應變的變化率與超越量有關，即有以下關系式：

(6)

式(6)中，Qs為塑性勢函數；F為破壞準則函數，對于摩爾-庫倫材料來說，

其中：φ為摩擦角，c為黏聚力。σ=(σx,σy,σz,τxy,τyz,τzx)T，

E,ν分別為彈性模量和泊松比。

如果將黏塑性應變率與一個偽時間步相乘，就可以得到累加到下一個時間步的黏塑性應變增量，于是有：

(δεvp)i=Δt(εvp)i

(7)

(Δεvp)i=(Δεvp)i-1+(δεvp)i

(8)

數值計算絕對穩(wěn)定的時間步與假定的破壞準則有關。對于摩爾-庫倫材料有：

(9)

塑性勢函數對應力的偏導數可以表示為：

(10)

(11)

將式(4)、(5)、(11)代入式(3)，可求得彈塑-黏彈-黏塑性體的應變增量為(殷宗澤等，2006；繆林昌，2007；姚仰平等，2009)：

(12)

式(12)即為土體的黏彈塑性本構方程。

3.3 比奧固結有限元方程

利用伽遼金加權余量法離散方程，考慮到土體的非線性特性，取Δt時間內的位移增量來代替位移，將式(1)、(2)離散成增量形式(李醫(yī)民等，2004)：

(13)

因為滲流取決于孔隙壓力全量的分布，而不是取決于時間內孔隙壓力增量，所以孔壓要用全量的形式表示，記時刻tn和tn+1時單元節(jié)點i的孔壓全量分別為ui(n)和ui(n+1)，且Δui=ui(n+1)-ui(n)，則式(4)可變換為：

(14)

式(14)即為三維比奧固結有限元方程。

3.4 參數的動態(tài)變化模型

3.4.1 孔隙度與滲透系數的非線性 流固耦合問題實際上是孔隙應力的消散引起土體骨架的變形，孔隙系數的變化，從而影響土體的滲透性，宏觀上表現為土體的固結變形。在比奧固結的假定條件下，根據孔隙度的相關定義和滲流力學Kozeny-Carman方程推得孔隙度n和滲透系數k的動態(tài)表達式(田杰等，2005)：

(15)

3.4.2 土體參數的非線性 采用鄧肯-張非線性模型，將土體的本構關系推廣到非線性，則本構關系{Δσ}=[D]{Δε}中矩陣[D]中的彈性常數E、ν不再視為常量，而是隨著應力狀態(tài)改變而改變，其切線彈性模量和切線泊松比的表達式如下(羅剛等，2004)：

(16)

(17)

式(16)、式(17)中，Rf為破壞比；c為黏聚力；φ為內摩擦角；σ1為第一主應力；σ3為第三主應力；n為彈性模量與固結壓力曲線的斜率(lgα)；G為土體常規(guī)三軸壓縮實驗結果所繪曲線截距；F=0.04，D=3為土體實驗參數；pa為大氣壓強。

3.5 定解條件

3.5.1 初始條件 初始條件是指在初始時刻(t=0)所研究對象各個求解變量的空間分布情況。初始條件越接近真實值，迭代收斂所用的時間越短，模型計算結果越符合實際情況。

(1) 地應力初始條件：

采用土體的自重應力估算土體的初始應力：

(18)

(2) 位移初始條件：

w(x,y,z,t)|t=0=0

(19)

(3) 孔隙水壓力初始條件：

u(x,y,z,t)|t=0=u0(x,y,z)

(20)

式(20)中，u0(x,y,z)為研究區(qū)域內已知初始孔隙水壓力。

3.5.2 邊界條件 (1) 孔隙水壓力邊界條件Γ1：

u(x,y,z)|Γ1=us

(21)

式(21)中，us為水頭邊界Γ1上的已知孔隙水壓力。

(2) 流量邊界條件Γ2：

(22)

(3) 自由面邊界條件Γ3：

(23)

式(23)中，μ為土體給水度；θ為自由面外法線方向與垂線的交角；q為通過自由面邊界Γ3的單位面積流量；Z為自由面所在的高程。

(4) 位移邊界條件Γ4：

(24)

式(24)中，wx,wy,wz為位移邊界Γ4上3個方向的已知位移。

上述比奧固結有限元數學模型即可運用Fortran語言編制相應的有限元程序進行求解(錢家歡等，1996；毛昶熙等，1999)。

4 模型識別與驗證

對整個研究區(qū)域用八節(jié)點六面體單元進行離散化，在平面上剖分成6 489個矩形網格單元，剖面上按地層巖性(含水沙層和黏性土相對隔水層)剖分成7層。單元總數為45 423個，節(jié)點總數為13 514個。研究區(qū)三維網格剖分詳見圖1。

圖1 三維網格剖分圖

利用前述開發(fā)的有限元計算機模型，選取2009年12月31日—2010年12月31日作為模型識別、驗證的時段，每個月作為1個應力期，共分12個應力期，每個應力期為1個時間步長。各含水層的初始孔隙水壓力由初始水位值換算得出，初始位移值及邊界上的位移均為0。對觀測井的水位進行擬合、反演，獲得了各含水層的水文地質參數，整個模型參數分區(qū)共58個，以第III承壓含水層組為例，參數分區(qū)如圖2所示，其參數分區(qū)參數值見表1。

圖2 第III承壓含水層組參數分區(qū)

圖3 地下水位擬合圖(南通市煙濾嘴廠)

圖4 地下水位擬合圖(通州二甲水廠)

圖5 2010年地面沉降量計算值與觀測值對比(mm)

圖3、圖4、圖5舉例說明了模型的地下水位值和地面沉降值擬合精度較高，從識別的結果看，各參數分區(qū)參數值的級別大小均符合常規(guī)，識別過程中所涉及的物理量都是在實際調查的基礎上確定的，因此模型的識別具有一定的精度和可信度。

5 模型預測

利用上述經過識別和驗證的土體自重固結壓縮、地下水開采與土體變形三維全耦合數值模型，模擬預測了研究區(qū)在地下水停采、僅考慮土體自重固結壓縮的影響下，自2010年12月31日—2025年12月31日各含水層地下水流場變化特征和固結壓縮沉降的發(fā)展趨勢。模型計算得出：南通市各層最大固結壓縮量均位于通州區(qū)，累計由土體自重固結壓縮導致的最大地面沉降量為2.42 mm，最大地面沉降速率為0.16 mm/a。圖6、圖7分別舉例說明了2025年12月31日第III承壓含水層的地下水預測流場變化和預測固結壓縮量分布特征，圖8舉例說明了2025年12月31日的預測地面沉降量分布特征。

從計算結果可以看出，南通市由土體自重固結壓縮引發(fā)的地面沉降，在平面上遠離市區(qū)的地方地面沉降值逐漸增大，符合遠離市區(qū)的土體固結程度偏低，彈性模量較低的規(guī)律。

表1 第Ⅲ承壓含水層參數分區(qū)表

圖6 2025年12月31日第III承壓含水層預測流場(m)

圖7 2025年12月31日第III承壓含水層壓縮量預測等值線圖(mm)

圖8 2025年12月31日預測地面沉降量等值線圖(mm)

6 結 論

(1) 以比奧固結理論為基礎，并結合土體流變理論，將土體水力學參數以及土力學參數隨有效應力的動態(tài)變化引入到模型中，建立的土體自重固結壓縮、地下水開采與地面沉降三維全耦合模型，可準確刻畫土體自重應力、地下水開采與地面沉降之間的相互作用機制，提高模型計算的置信度。

(2) 南通市在地下水停采、僅在土體自重固結壓縮影響下，自2010年12月31日—2025年12月31日的最大地面沉降量為2.42 mm，最大地面沉降速率為0.16 mm/a，土體自重固結壓縮對南通市地面沉降的影響極為有限。

(3) 南通市由土體自重固結壓縮引發(fā)的地面沉降，在平面上遠離市區(qū)的地方，地面沉降值逐漸增大，符合土體固結程度偏低，彈性模量較低的規(guī)律。

陳崇希，裴順平.2001.地下水開采-地面沉降模型研究[J].水文地質工程地質，28(2)：5-8.

陳曉平,白世偉.2001.軟黏土地基黏彈塑性比奧固結的數值分析[J].巖土工程學報,23(4):481-484.

李醫(yī)民,周鳳燕.2004.一類三維偏微分方程邊值問題的解法[J].江蘇大學學報:自然科學版,25(4):328-331.

羅剛,張建民.2004.鄧肯-張模型和沈珠江雙屈服面模型的改進[J].巖土力學,25(6):887-890.

駱祖江，劉金寶，李朗.2008.第四紀松散沉積層地下水疏降與地面沉降三維全耦合數值模擬[J].巖土工程學報，30(2):193-198.

毛昶熙,段祥寶,李祖貽,等.1999.滲流數值計算與程序應用[M].南京:河海大學出版社.

繆林昌.2007.非飽和土的本構模型研究[J].巖土力學,28(5):855-860.

錢家歡,殷宗澤.1996.土工原理與計算[M].北京:水利水電出版社.

SMITHI M, GRIFFITHS D V. 2003.有限元方法編程[M].王菘,周堅鑫,王來,等譯.北京:電子工業(yè)出版社.

田杰,劉先貴,尚根華.2005.基于流固耦合理論的套損力學機理分析[J].水動力學研究與進展:A輯,20(2):221-225.

維亞諾夫C C.1987.土力學的流變學原理[M].杜余培，譯.北京:科學出版社.

殷宗澤,周建,趙仲輝,等.2006.非飽和土本構關系及變形計算[J].巖土工程學報,28(2):137-146.

姚仰平,侯偉.2009.土的基本力學特性及其彈塑性描述[J].巖土力學,30(10):2881-2902.

LUO ZUJIANG, ZENG FENG. 2011. Finite element numerical simulation of land subsidence and groundwater exploitation based on visco-elastic-plastic biot’s consolidation theory [J]. Journal of Hydrodynamics, 23(5): 615-624.