付延玲，金瑋澤，陳興賢，2，談金忠

1.河海大學地球科學與工程學院，南京 210098

2.江蘇省地質礦產勘查局，南京 210018

3.江蘇南京地質工程勘察院，南京 210041

0 引言

地面沉降是一種由多種因素引起的地面高程緩慢降低的地質現(xiàn)象，其中過量開采地下水是引發(fā)地面沉降問題的主要原因。但近年來，隨著城市建設進程的加快，高層、超高層建筑的不斷增多，建筑荷載引發(fā)的地面沉降問題也日益突出。研究者們［1－4］認為，在上海，高群體建筑引起的地面沉降占總沉降量的30%左右，且沉降速度與建筑面積增長速度之間存在線性相關關系。

高層建筑荷載一方面直接導致土體有效應力增加，發(fā)生固結壓縮，引發(fā)地面沉降；另一方面在地下水埋深較淺的地區(qū)，高層建筑荷載的作用也會對地下水滲流場產生影響，導致滲流場重新調整，滲流場的變化同樣會引起應力場的變化，引發(fā)地面沉降或隆起變形。以往學者的研究更多集中在由抽水作用引起地下水滲流場變化而引發(fā)的地面沉降變形方面［5－7］，對由高層建筑荷載引發(fā)的土體變形特征，尤其是高層建筑荷載短時間作用下引起滲流場調整變化，進而形成的土體變形的研究尚不多見。

為了準確模擬預測由高層建筑荷載引發(fā)的地面沉降與隆起變形問題，筆者以比奧固結理論為基礎，綜合考慮了含水層和隔水層的流變特性，以及水力學和土力學參數(shù)隨滲流場和應力場發(fā)生變化的動態(tài)變化問題，建立了高層建筑荷載引發(fā)地面沉降與隆起變形的三維有限元數(shù)值模型，對高層建筑荷載引起的地面沉降與隆起變形進行了模擬預測分析，探討了變形過程中土體水力參數(shù)與力學參數(shù)的變化特征。

1 理論基礎與方法

1.1 比奧固結理論

飽和土體（假定其土骨架變形為線彈性、微小變形、滲流符合達西定律、水不可壓縮或微壓縮）的三維比奧固結方程［8］如下：

式中：G為剪切模量；ν為泊松比；wx、wy、wz分別為x、y、z方向上的位移分量；u為孔隙水壓力；Kx、Ky、Kz分別為x、y、z方向上的滲透系數(shù)；γ為土的重度；γw為水的重度；▽為拉普拉斯算子，▽2＝為時間。

1.2 土體本構模型

土的本構關系是土的力學特性即應力－應變－強度－時間等關系的數(shù)學表達式。對于考慮流變特性的土體來說，其變形特征主要表現(xiàn)為變形的時間與應力水平有關，所顯示的是具有彈性、塑性和黏滯性的黏彈塑性體，若將此類土體的總應變增量dε分為彈塑性應變增量dεep、黏彈性應變增量dεve和黏塑性應變增量dεvp，則具有流變特性的土體中任意點在任意時刻的應變增量［7，9］為

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

1）彈塑性應變增量

由土體的彈塑性本構模型可得

2）黏彈性應變增量

由kelvin流變模型E1ε＋Keε＝σ，可得應力不變時復雜應力狀態(tài)下的黏彈性應變增量：

式中：ηe＝E1/Ke，E1為kelvin體黏彈性模量，Ke為黏滯系數(shù)；A為應力矩陣，

E為彈性模量。

3）黏塑性應變增量

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

式中：Qs為塑性勢函數(shù)對于摩爾－庫倫材料來說，

其中：φ為內摩擦角；c為黏聚力；

β＝為第三偏應力不變量，

σ與τ分別為主應力與剪應力。

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

數(shù)值計算與絕對穩(wěn)定的時間步與假定的破壞準則有關。

對于摩爾－庫倫材料有

塑性勢函數(shù)對應力的偏導數(shù)可以表示為

所以，應用摩爾－庫倫準則的土體黏塑性應變增量可以表示為

將式（4）、（5）、（10）代入式（3），可求得彈塑－黏彈－黏塑性體的應變增量為

式（11）即為土體的黏彈塑性本構方程。

1.3 比奧固結有限元方程

利用伽遼金加權余量法離散方程，考慮到土體的非線性特性，取Δt時間內的位移增量來代替位移，將式（1）、（2）離散成增量形式［10］：

式中：Δδ為結點位移增量；Δu為結點孔隙壓力增量；ˉK為固體剛度矩陣；K為滲透流量矩陣；K′為應力－滲流耦合項矩陣；ΔQ為流量增量矩陣；B為自由面的積分矩陣；R為等效節(jié)點荷載，當存在高層建筑荷載時，R為包括高層建筑荷載引起的附加應力值；Rt為t時刻已經發(fā)生的位移所平衡了的那部分荷載。

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

式（13）即為三維比奧固結有限元方程。

1.4 參數(shù)的動態(tài)變化模型

1.4.1 孔隙度與滲透系數(shù)的非線性

流固耦合問題實際上是孔隙應力的消散引起土體骨架的變形，導致滲透系數(shù)變化，從而影響土體的滲透性，宏觀上表現(xiàn)為土體的固結變形。在比奧固結的假定條件下，根據(jù)孔隙度的相關定義和滲流力學Kozeny－Carman方程推得孔隙度n和滲透系數(shù)K的動態(tài)表達式［11－12］：

式中：n0為初始孔隙度；K0為初始滲透系數(shù)；εv為體應變

1.4.2 土體參數(shù)的非線性

采用鄧肯－張非線性模型，將土體的本構關系推廣到非線性，則本構關系｛Δσ｝＝D｛Δε｝中矩陣D中的彈性常數(shù)E、ν不再視為常量，而是隨著應力狀態(tài)改變而改變，其切線彈性模量Et和切線泊松比νt的表達式如下［13］：

式中：Rf為破壞比；σ1為第一主應力；σ3為第三主應力；ω為彈性模量與固結壓力曲線的斜率；logα，G為土體常規(guī)三軸壓縮實驗結果所繪曲線截距；I＝0.04，D＝3為土體實驗參數(shù)；p為大氣壓強。

1.5 定解條件

1.5.1 初始條件

1）地應力初始條件

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

式中：σx、σz為土體的初始水平向和垂向應力；z為計算點深度；K0為靜止側壓力系數(shù)，K0＝。為有效內摩擦角

2）位移初始條件

3）孔隙水壓力初始條件

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

1.5.2 邊界條件

1）孔隙水壓力邊界條件Γ1

式中：us為水頭邊界Γ1上的已知孔隙水壓力。

2）流量邊界條件Γ2

式中：qL為邊界Γ2上的已知單位面積流量。

3）自由面邊界條件Γ3

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

4）位移邊界條件Γ4

比奧固結有限元方程結合定解條件和水力學參數(shù)及土力學參數(shù)動態(tài)變化模型即可運用Fortran語言編制相應的有限元程序進行求解［5，14］。

2 應用算例

2.1 概念模型

為了研究高層建筑荷載影響下的地面沉降情況，排除其他影響因素干擾，建立一個簡單的均質含水層模型。模型長×寬×高為500m×500m×100 m，建筑物按每層自重1.2t/m2、承重每層0.3t/m2計算，地基尺寸的長×寬為70m×50m。用八節(jié)點六面體單元離散化模型，在平面上剖分為2 500個矩形網格單元，垂向上從上往下概化為：潛水含水層、第1承壓含水層，第2承壓含水層、第3承壓含水層及各含水層之間的黏性土弱含水層，共7層。每層土體劃分為一個參數(shù)分區(qū)，共7個參數(shù)分區(qū)。其中：潛水含水層巖性以粉砂、亞砂土為主，底板埋深20～30m，靜水位埋深0.9～1.2m；第1承壓含水層以粉砂為主，底板埋深50～52m，靜水位埋深3.2～3.8m；第2承壓含水層以亞砂土、細中砂為主，底板埋深79～83m，靜水位埋深7.6～8.1m；第3承壓含水層以細中砂、中粗砂為主，底板埋深96～100m，靜水位埋深11.5～13.4m。模型四周均概化為第一類已知水頭邊界，底部概化為隔水邊界。建筑物位置及模型分層如圖1所示。模型各參數(shù)分區(qū)物理力學性質參數(shù)如表1所示。其中：K0x、K0y為初始水平向滲透系數(shù)；K0z為初始垂向滲透系數(shù)；ν0為初始泊松比；E0為初始彈性模量；Sy為儲水率。

圖1 模型示意圖Fig.1 Model sketch

2.2 建筑荷載影響下的地面沉降與隆起變形數(shù)值模擬

模型建筑物高度按10層計算。模型計算周期為1a，劃分為15個應力期。其中：第1到第4應力期分別為5，10，20，31d；第5到15應力期每個應力期時間分別為第2月到第12月單月時間。假定初始情況下的各參數(shù)分區(qū)孔隙度大小均相同，且初始孔隙度n0＝0.36。選取建筑物中心點所在剖面（y＝260m）模型運行15個應力期時的地面沉降量，如圖2所示，其中地面變形量正值表示隆起量，負值表示沉降量。

圖2 不同應力期末的地面沉降曲線圖Fig.2 Land subsidence curve of difference stress period

土體的固結壓縮變形較為復雜，受到土體本身性質、邊界條件、排水條件以及上部荷載方式等因素影響。由圖2可以看出，由高層建筑荷載引起的地面沉降以建筑物中心點為中心呈現(xiàn)漏斗狀：到第4個應力期末，即31d時，最大沉降達到－3.8mm；在第10個應力期末，即第7個月末時，最大沉降量為－8.5mm。由圖2b可以看出：模型運行前9個應力期時，高層建筑荷載引起的地面沉降在x＝200～300m和x＝400～500m區(qū)間段均存在隆起現(xiàn)象；在第10個應力期末，即第7個月末時，隆起現(xiàn)象消失。而由圖2a可知：在高層建筑荷載作用5，10，20d時隆起量依次增大，最大值分別為0.9mm，1.0mm，1.5mm；31d時隆起量為1.3mm。可以看出，20d累計隆起量相對于31d累計隆起量較大。這是因為由于高層建筑荷載主要影響淺部地層，在模型中，淺部土體垂向滲透系數(shù)較小，所以淺部地層的孔隙水在受到高層建筑荷載作用時，更多地以橫向流動為主；而高層建筑荷載作用時間較短，其周圍土體中的孔隙水急劇增多，導致超靜孔隙水壓力在初期急劇增大，有效應力相應減小，土體骨架膨脹變形，必然使得建筑周圍土體發(fā)生回彈，地表出現(xiàn)隆起；且隨著應力期的增加，建筑周圍土體的超靜孔隙水壓力緩慢消散，地面隆起緩慢消失。

2.3 參數(shù)變化特征分析

隨著土體的固結壓縮與回彈變形，建筑物周圍及下部土體力學及水力學參數(shù)均會伴隨著土體的變形而發(fā)生相應變化。選取模型第1161號單元和第1167號單元進行孔隙度、滲透系數(shù)、彈性模量及泊松比分析。其中：第1161號單元位于建筑物中心底部；第1167號單元中心點坐標為（430，260），在剖面y＝200m上，距離建筑物40m，且兩單元均位于模型第一層，單元中心線位于y＝260m剖面上。兩單元孔隙度及滲透系數(shù)隨應力期的變化如圖3所示。

由于1161號單元與1167號單元所處區(qū)域地面沉降變化情況不同，1161號單元位于模型第1層地面沉降中心點位置，而1167號單元所處區(qū)域在前9個應力期均出現(xiàn)隆起現(xiàn)象，所以兩單元的孔隙度變化趨勢是不同的。由圖3a可以看出：1167號單元孔隙度在20d增大到最大值，隨著應力期的增加，孔隙度緩慢減??；而1161號單元在整個應力期時間段內，孔隙度一直處于減小的趨勢。造成這種現(xiàn)象的原因可能是：在建筑荷載作用初期，建筑物下部土體的孔隙水短時間擴散到周圍淺部地層，使得孔隙水大量聚集，超靜孔隙水壓力急劇增大，土體有效應力相應減小，淺部土體發(fā)生膨脹，使得建筑物周圍土體孔隙度在20d達到最大值；但隨后超靜孔隙水壓力逐漸消散，使得孔隙度緩慢減小，在第10個應力期，即第7個月時恢復到初始孔隙度大小。建筑物中心點下部土體在受到建筑荷載作用受到擠壓后，土體孔隙水向四周擴散，使得孔隙水壓力一直減小，有效應力一直增加，所以孔隙度隨著應力期的變化始終減小。

滲透系數(shù)隨著孔隙度的變化而發(fā)生變化。孔隙度增大，土體出水能力增強，導致滲透系數(shù)變大，而孔隙度減小，使土體出水能力變弱，從而導致滲透系數(shù)減小。由圖3也可以看出，1161號單元和1167號單元的滲透系數(shù)與孔隙度變化趨勢相同。

表1 地層參數(shù)Table 1 Stratum parameters

圖3 1161號單元和1167號單元孔隙度、水平向滲透系數(shù)和垂向滲透系數(shù)隨應力期變化曲線圖Fig.3 Variation curve among porosity and horizontal hydraulic conductivity and vertical hydraulic conductivity and stress periods of 1161element and 1167element

伴隨著土體的回彈及壓縮，彈性模量、泊松比必然發(fā)生變化。由圖4可以看出，1167號單元彈性模量、泊松比的變化與孔隙度的變化相對應，在20d左右土體彈性模量達到最小值，泊松比達到最大值，隨著應力期增加彈性模量緩慢增大，泊松比緩慢減小。而1161號單元彈性模量及泊松比的變化均較為平緩，彈性模量不斷增大，泊松比不斷減小。通過對土體力學及水力學參數(shù)的變化趨勢分析可以發(fā)現(xiàn)，在前9個應力期內各參數(shù)的變化相比后幾個應力期要快，這是因為隨著應力期的增長，土體固結變形量逐漸減小，導致土體力學及水力參數(shù)的變化相對于前期而言越來越不明顯。

3 結論

1）以比奧固結理論為基礎，結合土體的非線性流變理論，將土體的本構關系推廣到黏彈塑性，同時考慮土體水力學參數(shù)和土力學參數(shù)隨有效應力的動態(tài)變化問題，建立起的地下水滲流與土體變形三維全耦合模型更加符合實際，進一步提高了模型計算的精度。

2）由高層建筑荷載引起的地面沉降呈現(xiàn)漏斗狀，以建筑物中心點為漏斗中心。建筑荷載作用初期會引起周圍土體隆起，建筑物周圍淺部土層土體孔隙度、滲透系數(shù)及泊松比也隨著增加，到達峰值后，隨著時間的增加又緩慢減??；而彈性模量逐漸減小，到達最小值后，隨著時間的增加緩慢增大。建筑荷載下部土層，孔隙度、滲透系數(shù)及泊松比隨時間的增加均緩慢減小，而彈性模量隨時間的增加緩慢增大。

（References）：

［1］龔士良.上海城市建設對地面沉降的影響［J］.中國地質災害與防治學報，1998，9（2）：108－111.Gong Shiliang.Urban Construction’s Impact on Land Subsidence in Shanghai［J］.The Chinese Journal of Geological Hazard and Control，1998，9（2）：108－111.

［2］許燁霜，馬磊，沈水龍.上海市城市化進程引起的地面沉降因素分析［J］.巖土力學，2011，32（增刊1）：578－582.Xu Yeshuang，Ma Lei，Shen Shuilong.Influential Factors on Development of Land Subsidence with Process of Urbanization in Shanghai［J］.Rock and Soil Mechanics，2011，32（Sup.1）：578－582.

［3］嚴學新，沈國平.上海城區(qū)建筑密度與地面沉降關系分析［J］.水文地質工程地質，2002，29（6）：21－25.Yan Xuexin，Shen Guoping.Relationship Between Building Density and Land Subsidence in Shanghai Urban Zone［J］.Hydrogeology and Engineering Geology，2002，29（6）：21－25.

［4］唐益群，宋壽鵬，陳斌，等.不同建筑容積率下密集建筑群區(qū)地面沉降規(guī)律研究［J］.巖石力學與工程學報，2010，29（增刊1）：3425－3431.Tang Yiqun，Song Shoupeng，Chen Bin，et al.Study of Land Subsidence Rule of Dense Buildings Under Different Floor Area Ratios［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（Sup.1）：3425－3431.

［5］陳興賢，駱祖江，安曉宇，等.深基坑降水三維變參數(shù)非穩(wěn)定滲流與地面沉降耦合模型［J］.吉林大學學報：地球科學版，2013，43（5）：1572－1578.Chen Xingxian，Luo Zujiang，An Xiaoyu，et al.Coupling Model of Groundwatrer Three Dimensional Variable Parametrics Non－Steady Seepage and Land－Subsidence［J］.Journal of Jilin University：Earth Science Edition，2013，43（5）；1572－1578.

［6］Luo Zujiang，Zeng Feng.Finite Element Numerical Simulation of Land Subsidence and Groundwater Exploitation Based on Visco－Elastic－Plastic Biot’s Consolidation Theory［J］.Journal of Hydrodynamics，2011，23（5）：615－624.

［7］錢家歡，殷宗澤.土工原理與計算［M］.北京：水力水電出版社，1996.Qian Jiahuan，Yin Zongze.Principle and Canculation of Geotechnics［M］.Beijing：China Waterpower Press，1996.

［8］陳曉平，白世偉.軟黏土地基黏彈塑性比奧固結的數(shù)值分析［J］.巖土工程學報，2001，23（4）：481－484.Chen Xiaoping，Bai Shiwei.The Numerical Analysis of Visco－Elastic－Plastic Biot’s Consolidation for Soft Clay Foundation［J］.Chinese Journal of Geotechnical Engineering，2001，23（4）：481－484.

［9］李醫(yī)民，周鳳燕.一類三維偏微分方程邊值問題的解法［J］.江蘇大學學報：自然科學版，2004，25（4）：328－331.Li Yimin，Zhou Fengyan.Solution to a Class of Boundary Value Problem of Three Dimensional Partial Differential Equation［J］.Journal of Jiangsu University：Natural Science Edition，2004，25（4）：328－331.

［10］冉啟全，李士倫.流固耦合油藏數(shù)值模擬中物性參數(shù)動態(tài)模型研究［J］.石油勘探與開發(fā)，1997，24（3）：61－65.Ran Qiquan，Li Shlun.Study on Dynamic Models of Reservoir Parameters in the Coupled Simulation of Multiphase Flow and Reservoir Deformation［J］.Petroleum Exploration and Development，1997，24（3）：61－65.

［11］田杰，劉先貴，尚根華.基于流固耦合理論的套損力學機理分析［J］.水動力學研究與進展：A 輯，2005，20（2）：221－225.Tian Jie，Liu Xiangui，Shang Genhua.Casing Damage Mechanism Based on Theory of Fluid－Solid Coupling Flow Through Underground Rock［J］.Journal of Hydrodynamics：Series A，2005，20（2）：221－225.

［12］羅剛，張建民.鄧肯－張模型和沈珠江雙屈服面模型的改進［J］.巖土力學，2004，25（6）：887－890.Luo Gang，Zhang Jianmin.Improvement of Duncan－Chang Nonlinear Model and Shen Zhujiang’s Elastoplastic Model for Granular Soils［J］.Rock and Soil Mechanics，2004，25（6）：887－890.

［13］彭國倫.Fortran 95程序設計［M］.北京：中國電力出版社，2005.Peng Guolun.Fortran 95［M］.Beijing：China Electric Power Press，2005.

［14］Smith I M，Griffiths.D V.有限元方法編程［M］.3版.王菘，周堅鑫，王來，等，譯.北京：電子工業(yè)出版社，2003.Smith I M，Griffiths D V.Programming the Finite Element Method［M］.3rd ed.Translated by Wang Song，Zhou Jianxin， Wang Lai，et al.Beijing：Electronic Industry Press，2003.