，，

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

0 引 言

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

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

1 研究區(qū)概況

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

2 地質(zhì)概念模型

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

3 數(shù)學(xué)模型的建立與求解

3.1 比奧固結(jié)理論

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

(1)

(2)

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

3.2 土體本構(gòu)模型

土的本構(gòu)關(guān)系是土的力學(xué)特性即應(yīng)力-應(yīng)變-強(qiáng)度-時(shí)間等關(guān)系的數(shù)學(xué)表達(dá)式。對(duì)于考慮流變特性的土體來(lái)說(shuō)，其變形特征主要表現(xiàn)為變形的時(shí)間與應(yīng)力水平有關(guān)，所顯示的是具有彈性、塑性和黏滯性的黏彈塑性體，若將此類(lèi)土體的總應(yīng)變?cè)隽縟ε分為彈塑性應(yīng)變?cè)隽縟εep、黏彈性應(yīng)變?cè)隽縟εve、黏塑性應(yīng)變?cè)隽縟εvp，則具有流變特性的土體中任意點(diǎn)在任意時(shí)刻的應(yīng)變?cè)隽繛?維亞諾夫，1987；錢(qián)家歡等，1996；陳曉平等，2001；姚仰平等，2009)：

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

(3)

式(3)中各部分的應(yīng)變?cè)隽靠梢杂上率龇椒ù_定。

3.2.1 彈塑性應(yīng)變?cè)隽?由土體的彈塑性本構(gòu)模型可得：

dεep=[C]dσ′

(4)

3.2.2 黏彈性應(yīng)變?cè)隽?由Kelvin流變模型E1ε+Keε=σ可得應(yīng)力不變時(shí)復(fù)雜應(yīng)力狀態(tài)下的黏彈性應(yīng)變?cè)隽繛椋?/p>

(5)

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

3.2.3 黏塑性應(yīng)變?cè)隽?利用黏塑性法確定黏塑性應(yīng)變?cè)隽浚谠摲N方法中允許材料在有限“期間”內(nèi)超越破壞準(zhǔn)則(以破壞準(zhǔn)則函數(shù)F的值>0來(lái)表示)。在討論土體的黏塑性應(yīng)變而非塑性應(yīng)變時(shí)，應(yīng)變的變化率與超越量有關(guān)，即有以下關(guān)系式：

(6)

式(6)中，Qs為塑性勢(shì)函數(shù)；F為破壞準(zhǔn)則函數(shù)，對(duì)于摩爾-庫(kù)倫材料來(lái)說(shuō)，

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

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

如果將黏塑性應(yīng)變率與一個(gè)偽時(shí)間步相乘，就可以得到累加到下一個(gè)時(shí)間步的黏塑性應(yīng)變?cè)隽?，于是有?/p>

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

(7)

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

(8)

數(shù)值計(jì)算絕對(duì)穩(wěn)定的時(shí)間步與假定的破壞準(zhǔn)則有關(guān)。對(duì)于摩爾-庫(kù)倫材料有：

(9)

塑性勢(shì)函數(shù)對(duì)應(yīng)力的偏導(dǎo)數(shù)可以表示為：

(10)

(11)

將式(4)、(5)、(11)代入式(3)，可求得彈塑-黏彈-黏塑性體的應(yīng)變?cè)隽繛?殷宗澤等，2006；繆林昌，2007；姚仰平等，2009)：

(12)

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

3.3 比奧固結(jié)有限元方程

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

(13)

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

(14)

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

3.4 參數(shù)的動(dòng)態(tài)變化模型

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

(15)

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

(16)

(17)

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

3.5 定解條件

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

(1) 地應(yīng)力初始條件：

采用土體的自重應(yīng)力估算土體的初始應(yīng)力：

(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ū)域內(nèi)已知初始孔隙水壓力。

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

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

(21)

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

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

(22)

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

(23)

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

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

(24)

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

上述比奧固結(jié)有限元數(shù)學(xué)模型即可運(yùn)用Fortran語(yǔ)言編制相應(yīng)的有限元程序進(jìn)行求解(錢(qián)家歡等，1996；毛昶熙等，1999)。

4 模型識(shí)別與驗(yàn)證

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

圖1 三維網(wǎng)格剖分圖

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

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

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

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

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

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

5 模型預(yù)測(cè)

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

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

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

圖6 2025年12月31日第III承壓含水層預(yù)測(cè)流場(chǎng)(m)

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

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

6 結(jié) 論

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

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

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

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

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

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

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

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

毛昶熙,段祥寶,李祖貽,等.1999.滲流數(shù)值計(jì)算與程序應(yīng)用[M].南京:河海大學(xué)出版社.

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

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

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

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

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

殷宗澤,周建,趙仲輝,等.2006.非飽和土本構(gòu)關(guān)系及變形計(jì)算[J].巖土工程學(xué)報(bào),28(2):137-146.

姚仰平,侯偉.2009.土的基本力學(xué)特性及其彈塑性描述[J].巖土力學(xué),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.