地下水滲流與地面沉降耦合數(shù)值模擬模型

地下水滲流與地面沉降耦合數(shù)值模擬模型

0terzadgei的數(shù)值模型正確預測地下水開采造成的地面沉降是確定地面沉降控制規(guī)劃的關鍵。地下水開采造成的地面沉降的問題實際上是滲流場和應力場相互作用的問題。目前，國內(nèi)外關于地下水滲透和地面沉降的一般計算模型主要包括基于tergahuri有效應力原理建立的三維地下水滲透和三維垂向固結(jié)沉降的結(jié)合模型。以美國地質(zhì)調(diào)查局開發(fā)的美國地質(zhì)調(diào)查局工具商業(yè)軟件為代表（坎齊梁，羅立紅，2010；周念青等，2011；李燕等，2012）。羅祖江等（2006、2008、2009）提出了基于比奧固結(jié)理論的新理論，結(jié)合土體材料與非土體界面理論的動態(tài)變化理論，將土體結(jié)構(gòu)關系推廣到粘土彈塑化。同時，考慮到土壤水力學參數(shù)和土力學參數(shù)隨著滲流場和應力場的動態(tài)變化率的關系，建立了三維結(jié)合地下水流和地面沉降的三維模型。在這項工作中，我們使用上述兩種模型比較和分析了地下滲透和地面沉降。結(jié)果表明，無論模型的計算方法還是模型所示的地面沉降機，后者相對完整，具有更高的計算精度。1基礎和方法的理論1.1u3000含水層沉降量的初始有效性對于非均質(zhì)、空間三維非穩(wěn)定流系統(tǒng),如把坐標軸取得與各向異性的主方向一致,可用以下地下水流連續(xù)性方程及其定解條件方程式(1)來描述(駱祖江等,2007;闞京梁和羅立紅,2010):式中:kxx、kyy、kzz為各向異性主方向滲透系數(shù)(m/d);k為自由面邊界滲透系數(shù)(m/d);t為模型計算時間(d);h為點(x,y,z)在t時刻的水頭值(m);W為源匯項(1/d);μs為含水層儲水率(1/m);q為流量(m3/d);μ為飽和差(自由面上升)或給水度(自由面下降),無量綱,它表示在自由面改變單位高度下,從含水層單位截面積吸收或排出的水量;θ為自由面外法線方向與垂線的交角;Γ1、Γ2、Γ3分別為第1類邊界、第2類邊界和自由面邊界;Ω為計算區(qū)域.計算含水層沉降量的沉降模型方程為:承壓含水層的彈性變形量:承壓含水層的非彈性變形量:式中:Δb為含水層彈性壓縮量(m),正為壓縮,負為回彈;Δb*為含水層非彈性壓縮量(m);ΔH為水頭變化值(m);μfe為含水層骨架成分的彈性儲水因子,無量綱;μsu3000ke為含水層骨架成分的彈性儲水率(1/m);μfv為含水層骨架成分的非彈性儲水因子,無量綱;μsu3000kv為含水層骨架成分的非彈性儲水率(1/m);b0為可壓縮含水層的厚度(m);n為孔隙率,無量綱;nW為水位以上作為多孔介質(zhì)總體積的一部分的濕氣容量,無量綱.式中:G為剪切模量(kPa);v為泊松比,無量綱;g為重力常數(shù)(m/s2);Cc為土體的壓縮系數(shù),無量綱;σue10b0為初始有效應力(kN/m2);e0為初始孔隙比,無量綱;ρw為水的密度(kg/m3).將上述二模型通過水頭項耦合起來,即可形成地下水三維滲流與一維垂向固結(jié)的地下水滲流與地面沉降耦合數(shù)值模擬模型.將該模型采用有限差分法進行求解,并采用強隱式(SIP)聯(lián)立迭代求解法求解線性方程組.將整個求解過程采用FORTRAN語言(彭國倫,2005),編制成三維地下水滲流與一維垂向固結(jié)沉降耦合有限元計算機程序.1.2儲層參數(shù)的程序法飽和土體中假定土骨架變形為線彈性、微小變形、滲流符合達西定律、水不可壓縮或微壓縮的三維比奧固結(jié)方程如下(錢家歡和殷宗澤,1996):式中:G為剪切模量(kPa);v為泊松比,無量綱;wx、wy、wz分別為x、y、z方向上的位移分量(m);u為孔隙水壓力(kPa);kx、ky、kz分別為x、y、z方向上的滲透系數(shù)(m/d);γ為土的重度(kN/m3);γw為水的重度(kN/m3).1.3土體變形特征土的本構(gòu)關系是土的力學特性即應力-應變-強度-時間等關系的數(shù)學表達式.對于考慮流變特性的土體來說,其變形特征主要表現(xiàn)為變形的時間與應力水平有關,所顯示的是具有彈性,塑性和粘滯性的粘彈塑性體,若將此類土體的總應變增量dε分為彈塑性應變增量dεep,粘彈性應變增量dεve,粘塑性應變增量dεvp,則具有流變特性的土體中任意點在任意時刻的應變增量為(LuoandZeng,2011):1.4高層建筑應力應變的數(shù)值計算利用伽遼金加權(quán)余量法離散方程,考慮到土體的非線性特性,取Δt時間內(nèi)的位移增量來代替位移,將式(7)、(8)離散成增量形式(李醫(yī)民和周鳳燕,2004):式中:Δδ為結(jié)點位移增量;Δu為結(jié)點孔隙壓力增量;K-為固體剛度矩陣;K為滲透流量矩陣;Kue10b為應力-滲流耦合項矩陣;ΔQ為流量增量矩陣;B為自由面的積分矩陣;R為等效節(jié)點荷載,當存在高層建筑荷載時,R包括高層建筑荷載引起的附加應力值;Rt為t時刻已經(jīng)發(fā)生的位移所平衡了的那部分荷載.因為滲流取決于孔隙壓力全量的分布,而不是取決于時間內(nèi)孔隙壓力增量.所以孔壓要用全量的形式表示,記時刻tn和tn+1時單元節(jié)點i的孔壓全量分別為ui(n)和ui(n+1),且Δui=ui(n+1)-ui(n),則式(19)可變換為:上式即為三維比奧固結(jié)有限元方程.1.5土體本構(gòu)模型的建立流固耦合問題實際上是孔隙應力的消散引起土體骨架的變形,滲透系數(shù)會發(fā)生相應的變化,從而影響土體的滲透性,宏觀上表現(xiàn)為土體的固結(jié)變形.在比奧固結(jié)的假定條件下,根據(jù)孔隙度的相關定義和滲流力學Kozeny-Carman方程推得孔隙度n和滲透系數(shù)k的動態(tài)表達式(冉啟全和李士倫,1997;田杰等,2005):式中:n0為初始孔隙度,無量綱;k0為初始滲透系數(shù)(m/d);εv為體應變,無量綱;采用鄧肯-張非線性模型,將土體的本構(gòu)關系推廣到非線性,則本構(gòu)關系{Δσ}=[D]{Δε}中矩陣[D]中的彈性常數(shù)E、v不再視為常量,而是隨著應力狀態(tài)改變而改變,其切線彈性模量和切線泊松比的表達式如下(羅剛和張建民,2004):式中:Rf為破壞比,無量綱;c為粘滯力(kPa);φ為內(nèi)摩擦角(°);σ1為第1主應力(kPa);σ3為第3主應力(kPa);n為彈性模量與固結(jié)壓力曲線的斜率;logα,G為土體常規(guī)三軸壓縮實驗結(jié)果所繪曲線截距,F=0.04,D=3為土體實驗參數(shù);pa為大氣壓強(kPa).1.6固定配置條件1.6.1土體初始垂向應力采用土體的自重應力估算土體的初始應力:式中:σx、σz為土體的初始水平向和垂向應力(kPa);z為計算點深度(m);K0為靜止側(cè)壓力系數(shù),無量綱,(2)第一次重大轉(zhuǎn)移的條件(3)初始階段的間隙水壓力1.6.2邊境條件式中:us為水頭邊界Γ1上的已知孔隙水壓力(kPa).(2)流量限制條件2式中:為邊界Γ2上的已知單位面積流量(m/d).(3)自由面內(nèi)法線方向與垂線的交角式中:μ為土體給水度,無量綱;θ為自由面外法線方向與垂線的交角(°);q為通過自由面邊界Γ3的單位面積流量(m/d);Z為自由面所在的高程.(4)壓解程序式中:wx,wy,wz為位移邊界Γ4上3個方向的已知位移(m).比奧固結(jié)有限元方程結(jié)合定解條件和土體力學和水力學參數(shù)的動態(tài)變化模型即可運用Fortran語言編制相應的有限元程序進行求解(SmithandGriffiths,2003;彭國倫,2005).2應用計算示例2.1地下水含水層結(jié)構(gòu)本文以南通興益大廈深基坑場地為研究區(qū),含水層以第四系松散層孔隙水為主,按照水文地質(zhì)條件對含水層進行劃分,將場地含水層細分為7層,垂向從上往下分別為潛水含水層、第1承壓含水層、第2承壓含水層、第3承壓含水層以及各含水層間的粘性土弱含水層.其中,潛水含水層:巖性以粉砂、亞砂土為主,底板埋深一般為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.2.2模型分層及模型特征為了對地下水滲流與地面沉降進行模擬研究,選取研究區(qū)場地平面范圍為500×500m2,垂向深度為100m,建立含水層模型.用八節(jié)點六面體單元離散化模型,在平面上剖分為2500個矩形網(wǎng)格單元,垂向上剖分為:潛水含水層、第1承壓含水層,第2承壓含水層、第3承壓含水層及各含水層之間的粘性土弱含水層,共7層.取第1至第7層含水層底板標高分別為-25m,-37m,-52m,-64m,-79m,-89m,-100m.每層土體劃分為1個參數(shù)分區(qū),垂向從上往下共劃分為7個參數(shù)分區(qū).模型四周均概化為第1類已知水頭邊界,底部概化為隔水邊界.抽水井及模型分層如圖1所示.2.3地下水位下降的影響抽水井工作一段時間,使得地下水位下降后,停止抽水,隨著時間的推移地下水位會恢復至初始水位.為了探討此情況下的地面沉降特征,模型設定為2個應力期,分別為抽水水位下降期與停止抽水水位回升期.抽水水位下降期又劃分為3個時間步長,停止抽水水位回升期劃分為9個時間步長,劃分一個月為一個時間步長,共12個月.抽水井抽取第3承壓含水層地下水,對第1個應力期進行定流量抽水,抽水量為70m3/d,從第2個應力期開始停止抽水.地下水位觀測井及沉降觀測井位置與抽水井位置相同,地下水位觀測井對第3承壓含水層水位進行觀測.分別采用地下水三維滲流與一維垂向固結(jié)沉降耦合有限元計算機程序以及比奧固結(jié)三維全耦合有限元計算機程序?qū)ν怀樗^程地下水位下降、回升條件下的地面沉降進行計算,并對計算結(jié)果進行分析對比.采用地下水三維滲流與一維垂向固結(jié)沉降耦合有限元計算機程序?qū)Φ叵滤患暗孛娉两颠M行計算,模型地層各參數(shù)(水平向滲透系數(shù)Kx、Ky,垂向滲透系數(shù)Kz,給水度SS,儲水率Sy,含水層單位厚度骨架成分的彈性儲水因子與非彈性儲水因子μske、μskv)如表1所示.取模型地面沉降計算值及第3承壓含水層計算水位,如圖2所示.可以看出,從第4個月停止抽水開始,地下水位回升較為明顯,在第8個月時,地下水位已基本恢復至初始水位.整個應力期時間段的地下水位變化趨勢與地面沉降變化趨勢相同,在地下水位回升至初始水位時,地面沉降與地下水位呈現(xiàn)出同步變化的趨勢,并且當?shù)叵滤恢鸩交厣脸跏妓粫r,地面沉降也逐步回升到初始的零沉降狀態(tài).此過程也反映了地下水三維滲流與一維垂向固結(jié)沉降耦合模型計算出的土體變形為線彈性變形,沉降隨著地下水位的變化而變化,與地下水位變化趨勢完全相同.運用比奧固結(jié)三維全耦合有限元計算機程序?qū)Υ诉^程進行模擬計算,模型地層各參數(shù)(初始水平向滲透系數(shù)Kx、Ky,初始垂向滲透系數(shù)Kz,給水度SS,儲水率Sy,土體重度r,內(nèi)摩擦角φ,粘聚力c,初始泊松比ν0,初始彈性模量E0)如表2所示.同樣取地面沉降值及第3承壓含水層計算水位如圖3所示.地下水位在第9個月恢復至初始水位,整個地下水位下降趨勢也和地下水三維滲流與一維垂向固結(jié)沉降耦合有限元計算機程序計算結(jié)果不同.地下水三維滲流與一維垂向固結(jié)沉降耦合有限元計算機程序所計算地下水位在第4個應力期停止抽水后迅速下降,而比奧固結(jié)三維全耦合有限元計算機程序則變化較為平緩.出現(xiàn)這種現(xiàn)象的原因是,地下水三維滲流與一維垂向固結(jié)沉降耦合有限元計算機程序在整個計算過程中,土體力學及水力學參數(shù)均未隨土體的固結(jié)變形而發(fā)生變化,而比奧固結(jié)有限元程序考慮了土體力學及水力學參數(shù)隨土體的固結(jié)變化而變化.抽水使得土體固結(jié)變形,此時土體壓密,其孔隙度減小,滲透系數(shù)減小,而在停止抽水后,地下水位回升會使土體孔隙度增大,此時土體的滲透系數(shù)又會有所增大,所以第4個月到第6個月,地下水位回升速度相對較慢,而第7個月到第9個月地下水位回升相對較快,隨著應力期的增加,地下水位必然恢復至初始狀態(tài).由于比奧固結(jié)三維全耦合有限元計算機程序在土體的本構(gòu)模型中考慮了土體的粘彈塑性,因此,土體在地下水位下降發(fā)生固結(jié)變形后,將產(chǎn)生部分永久的殘余沉降量,即土體在地下水位恢復至初始狀態(tài)時,仍存在地面沉降.并且地面沉降及回彈伴隨著地下水位的變化存在滯后性,由圖3可以看出,當?shù)?個月停止抽水時,其地面沉降仍在繼續(xù),在第5個月地面沉降值才有所減小,地面沉降滯后于地下水位變化約1個月左右.抽水導致土體中的孔隙水大量排出,地下水滲流場發(fā)生變化,孔隙水壓力消散,土體的有效應力增加,導致土體的固結(jié),宏觀上反映為地面沉降.反之,隨著地面沉降的發(fā)生,土體固結(jié),其微觀上表現(xiàn)為土體的壓密,土體孔隙減小,土體的滲透性相應降低,從而影響地下水滲流.地下水三維滲流與一維垂向固結(jié)的地下水滲流與地面沉降耦合數(shù)值模擬模型將滲流與變形單獨考慮,只通過計算地下水位將地下水滲流模型與地面沉降模型耦合起來,是部分耦合模型,并未將土體變形與孔隙壓力消散同時考慮,沒有從機理上實現(xiàn)地下水滲流與地面沉降的耦合.而且隨著土體固結(jié)壓縮,不考慮土體力學參數(shù)及水力參數(shù)隨土體中應力場改變的動態(tài)變化問題,與實際情況不符,這些不足對地下水滲流引起的地面沉降計算的準確性存在較大的影響,比奧固結(jié)地下水滲流與地面沉降三維全耦合數(shù)值模擬模型考慮土體變形與孔隙壓力消散的耦合作用,從機理上實現(xiàn)了滲流場與應力場的耦合,并且考慮土體力學及水力學參數(shù)的動態(tài)變化問題,引入土體的粘彈塑性本構(gòu)模型,大大提高了地面沉降模擬的準確性,使得計算結(jié)果更加符合實際情況.參數(shù)的準確性深刻地影響著模型計算的準確性.地下水開采引發(fā)地面沉降問題實際是一個滲流場和應力場相互影響、相互作用的問題.隨著地面沉降的發(fā)生,土體水力學參數(shù)和土力學參數(shù)均會發(fā)生相應變化.因此,如何正確模擬參數(shù)的變化過程至關重要.由于模型計算結(jié)果為孔隙水壓力,為了更好的反應計算過程中土體力學參數(shù)和水力學參數(shù)隨地下水位的變化關系,將孔隙壓力進行換算,得出相應的地下水位值.相應的地下水位其中u為孔隙水壓力,ue4d1為水的重度.選取模型第15559號單元進行孔隙度和滲透系數(shù)變化分析.該單元位于模型第7層,單元中心點與抽水井位置相同.由于抽水影響,導致土體孔隙水壓力變化,從而引起有效應力變化,土體發(fā)生固結(jié)變形,引起孔隙度的變化.由圖4可以看出在整個應力期時間段內(nèi),孔隙度的變化趨勢與地面沉降變化趨勢相同.孔隙度在前4個月均呈現(xiàn)減小趨勢,從第5個月開始,孔隙度緩慢增大,第9個月到第12個月孔隙度變化較小,基本趨于穩(wěn)定.滲透系數(shù)隨著孔隙度的變化而發(fā)生變化.孔隙度減小,使土體出水能力變?nèi)?從而導致滲透系數(shù)減小;反之,孔隙度增大,土體出水能力增強,滲透系數(shù)減小.由圖5和圖6可以看出所選單元的滲透系數(shù)變化趨勢與孔隙度變化趨勢相同,呈現(xiàn)出先減小后增大,而后趨于穩(wěn)定的趨勢.同樣,選取模型第15559號單元的彈性模量及泊松比進行分析.彈性模量及泊松比隨應力期變化如圖7和圖8所示.由于地下水開采對土體的固結(jié)變形壓縮的影響,在前4個應力期彈性模量增大,而泊松比呈現(xiàn)減小趨勢.隨著水位的回升,土體發(fā)生部分回彈變形,在第5個月到第9個月,彈性模量減小,泊松比增大,但是在第9個月