周躍峰,譚國(guó)煥,甄偉文,樊少鵬,3
(1.長(zhǎng)江科學(xué)院水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室,武漢 430010;2.香港大學(xué)土木工程系,香港;3.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究院樞紐設(shè)計(jì)處,武漢 430010)
非飽和滲流分析在FLAC3D中的實(shí)現(xiàn)和應(yīng)用
周躍峰1,2,譚國(guó)煥2,甄偉文2,樊少鵬2,3
(1.長(zhǎng)江科學(xué)院水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室,武漢 430010;2.香港大學(xué)土木工程系,香港;3.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究院樞紐設(shè)計(jì)處,武漢 430010)
基于FLAC3D內(nèi)置的FISH語(yǔ)言,對(duì)其滲流模型進(jìn)行了擴(kuò)展。將實(shí)驗(yàn)室測(cè)得的原狀黃土土水特征曲線,利用MATLAB實(shí)現(xiàn)了參數(shù)擬合。將擬合參數(shù)編寫為子程序,在每步滲流分析中加以調(diào)用,從而實(shí)現(xiàn)非飽和土的滲流模擬。利用該方法與有限元法分別模擬3種邊界條件下的一維入滲問題,并比較模擬結(jié)果。經(jīng)驗(yàn)證,結(jié)果吻合。其結(jié)果對(duì)研究非飽和土理論及擴(kuò)展FLAC3D在非飽和土數(shù)值計(jì)算方面的應(yīng)用均具有一定的意義。
非飽和土;滲流;土水特征曲線;FLAC3D;二次開發(fā)
滲流是指水或其他流體在巖土等孔隙或裂隙介質(zhì)中的流動(dòng)。非飽和滲流是巖土體中孔隙水的最普遍的一種滲流方式[1]。該問題廣泛存在于邊坡工程[2]、地基/路基工程[3]、隧道工程[4]中。滲流直接影響作為滲流骨架的巖土體力學(xué)性質(zhì),因而在巖土工程變形和穩(wěn)定性分析中占有重要地位。通過非飽和滲流計(jì)算能夠得到土中水的運(yùn)動(dòng)軌跡,并為流固耦合分析提供水分場(chǎng)的分布規(guī)律。
FLAC3D是美國(guó)ITASCA咨詢公司開發(fā)的三維有限差分程序[5]。該程序采用顯式拉格朗日算法和混合-離散分區(qū)技術(shù)進(jìn)行積分計(jì)算,可求解地質(zhì)材料的高度非線性、孔隙介質(zhì)的應(yīng)力-滲流耦合、熱-力耦合以及動(dòng)力學(xué)等問題。該程序依靠較快的運(yùn)算速度、強(qiáng)大的前后處理、靈活的FISH語(yǔ)言,在目前的科研及工程領(lǐng)域具有廣泛的應(yīng)用。但是其內(nèi)置的滲流模塊,在進(jìn)行非飽和滲流分析時(shí)并不夠完善。由于缺少負(fù)孔隙水壓力(或基質(zhì)吸力)同含水量(或飽和度)的關(guān)聯(lián)函數(shù),F(xiàn)LAC3D中負(fù)孔隙水壓力僅反映土體的剪脹變形[5]。因此,難以直接利用FLAC3D對(duì)非飽和土進(jìn)行準(zhǔn)確的滲流分析,一定程度上限制了FLAC3D的應(yīng)用。李毅等[6]通過為非飽和區(qū)賦以很低的滲透系數(shù),運(yùn)用FLAC3D進(jìn)行非飽和滲流模擬,提供了一個(gè)可供借鑒的方法。但該方法本質(zhì)上并未考慮負(fù)孔隙水壓力同含水量或飽和度之間的重要聯(lián)系,影響求解域中孔壓場(chǎng)的演化和進(jìn)一步分析應(yīng)力場(chǎng)的分布規(guī)律。為了解決以上不完善之處,本文介紹了非飽和滲流分析在FLAC3D中的實(shí)現(xiàn)方法。
國(guó)內(nèi)外許多學(xué)者投入了大量的精力,逐漸建立起較為完備的土的入滲理論框架,如Richard[7],Philip[8],Lumb[9],F(xiàn)redlund[10]等。三維條件下,土體的滲流控制方程可以表示為下式(1)的形式,即
式中:kwx,kwy和kwz分別表示3個(gè)方向的滲透系數(shù)(m/s);ρw表示水的密度(103kg/m3);mw為土水特性曲線的斜率;uw表示土中的孔隙水壓力(kPa)。
土的飽和-非飽和滲流問題,至少需要2個(gè)基本函數(shù):土水特征曲線和滲透系數(shù)函數(shù)。其中,土水特征曲線描述了土的含水量隨基質(zhì)吸力或負(fù)孔隙水壓力的變化;滲透系數(shù)函數(shù)則表示滲透系數(shù)隨含水量或飽和度的變化。目前,有多個(gè)經(jīng)驗(yàn)公式可用來描述土壤土水特征曲線,包括:Brooks-Corey模型、Garden模型、van Genuchten模型、Fredlund&Xing模型等[11]。其中,van Genuchten模型作為較有代表性的公式被廣泛使用,即
式中:θ,θr,θs分別表示體積含水量(滿足θ=ns;即孔隙率n和飽和度s的乘積)、殘余含水量,以及最小吸力含水量;ua和uw分別表示孔隙氣壓力和孔隙水壓力;af,nf,mf為擬合參數(shù)。
非飽和土的滲透系數(shù)函數(shù)大多描述了非飽和土滲透系數(shù)與飽和土滲透系數(shù)的關(guān)系。在數(shù)值模擬中較多采用的一個(gè)公式[4,12],可表示為
式中:ksat和kunsat分別為土的飽和滲透系數(shù)和非飽和滲透系數(shù)。
3.1 曲線擬合
土水特征曲線的測(cè)試方法較多,包括:體積壓力板法、鹽溶液法、Tempe儀法、濾紙法、GDS四維應(yīng)力路徑法等[13]。實(shí)驗(yàn)室測(cè)量的結(jié)果通常都是基質(zhì)吸力同含水量相對(duì)應(yīng)的多個(gè)離散點(diǎn)。鑒于土水特征曲線的連續(xù)性,需要對(duì)由實(shí)驗(yàn)室試驗(yàn)測(cè)得的離散點(diǎn)作進(jìn)一步分析,獲取數(shù)值模擬所需的基本參數(shù)。實(shí)驗(yàn)室中測(cè)得的離散點(diǎn)可利用Matlab R2010寫執(zhí)行腳本作曲線擬合,步驟如下:
(1)創(chuàng)建一個(gè)新的M-文件,并輸入van Genuchten公式,
function F=fun(x,xdata)
F=x(1)+(x(2)-x(1))./(1+(x(3)*xdata).^x(4)).^(1-1./x(4))。
(2)在命令窗口分別輸入試驗(yàn)測(cè)得的基質(zhì)吸力和體積含水量,
xdata=[x1 x2…xn];
ydata=[y1 y2…yn]。
(3)為迭代計(jì)算賦初值,
x0=[0.1,0.1,0.1,1]。
(4)利用函數(shù)擬合試驗(yàn)數(shù)據(jù),獲取擬合參數(shù),
[x,resnorm]=lsqcurvefit(@fun,x0,xdata,ydata)。
此時(shí),可在工作區(qū)獲得相應(yīng)的擬合參數(shù),亦可在工作區(qū)繪制出擬合圖形。
3.2 算法的修正思路
非飽和土的滲流問題是一個(gè)非線性問題,滲透系數(shù)、孔隙壓力、逸出面范圍等只能在計(jì)算過程中迭代確定。滲透系數(shù)會(huì)受到顆粒級(jí)配、孔隙率、飽和度等因素的影響。對(duì)飽和土體,滲透系數(shù)通常簡(jiǎn)化為常數(shù);而對(duì)非飽和土體,滲透系數(shù)可能因含水量的變化而變化很大。根據(jù)第2節(jié)所述,修正中需考慮負(fù)孔隙水壓力及滲透系數(shù)隨含水量的變化。需要注意的是有限差分程序中,滲透系數(shù)是單元變量,而孔隙水壓力是節(jié)點(diǎn)變量(單元孔隙水壓力為各節(jié)點(diǎn)值的映射)。
數(shù)值計(jì)算步驟如下:①設(shè)置求解域的邊界條件和初始條件(如節(jié)點(diǎn)的孔隙水壓力);②在程序中指定土水特征曲線(按上述Matlab擬合結(jié)果)和滲透系數(shù)函數(shù);③根據(jù)土水特征曲線,將節(jié)點(diǎn)的孔隙水壓力轉(zhuǎn)化為節(jié)點(diǎn)的飽和度;④計(jì)算一個(gè)滲流時(shí)步,按照質(zhì)量平衡方程,各節(jié)點(diǎn)含水量在計(jì)算中發(fā)生變化;⑤根據(jù)土水特征曲線,將新的節(jié)點(diǎn)飽和度轉(zhuǎn)化成為節(jié)點(diǎn)的孔隙水壓力;⑥按照單元中的各節(jié)點(diǎn)的孔隙水壓力和飽和度,映射得到該單元的孔隙水壓力和飽和度;⑦按照該單元的飽和度,確定相應(yīng)的非飽和滲透系數(shù);⑧計(jì)算下一個(gè)滲流時(shí)步。
為了驗(yàn)證上述方法的正確性,本研究中,利用GEO-STUDIO軟件中的有限元程序SEEP/W的計(jì)算結(jié)果與有線差分程序FLAC3D的計(jì)算結(jié)果進(jìn)行對(duì)比。GEO-STUDIO是一個(gè)二維的巖土工程仿真分析軟件,其滲流模塊發(fā)展較為成熟,應(yīng)用較廣泛,但是,它只局限于二維滲流計(jì)算,不能進(jìn)行三維分析并且在其與SIGMA/W進(jìn)行耦合計(jì)算時(shí)常出現(xiàn)收斂不理想的情況。
圖1 基于FLAC3D和SEEP/W的計(jì)算模型Fig.1 Num ericalmodel built in FLAC3D and SEEP/W
本研究中,分別利用FLAC3D和SEEP/W建立相同尺寸的一維模型。模型長(zhǎng)度10 m,包括50個(gè)計(jì)算單元(圖1)。程序驗(yàn)證中共模擬3組算例,分別比較了3種條件下2個(gè)程序的計(jì)算結(jié)果。
本文的土水特征曲線是在GDS公司生產(chǎn)的非飽和三軸儀上進(jìn)行增濕試驗(yàn),利用軸平移技術(shù)測(cè)得的原狀黃土的基質(zhì)吸力和含水量的關(guān)系。并按照第3.1節(jié)中所述方法,將實(shí)驗(yàn)室測(cè)試得到的土水特性離散點(diǎn)在Matlab程序中進(jìn)行曲線擬合(圖2)。黃土的飽和滲透系數(shù)通過現(xiàn)場(chǎng)雙環(huán)注水試驗(yàn)測(cè)得。測(cè)試過程詳見筆者的博士論文[14],在此不作贅述,相關(guān)參數(shù)概括于表1。根據(jù)擬合結(jié)果,該黃土殘余含水量為8.7%,飽和含水量為40.7%;由于采用增濕曲線,不對(duì)該黃土進(jìn)氣值進(jìn)行分析。
圖2 測(cè)量及擬合的土水特征曲線Fig.2 Measured and fitted SWCC
表1 算例中采用的參數(shù)Table 1 Parameters adopted in the simu lated cases
3組算例中,底部邊界均設(shè)為恒定孔隙水壓力為-147 kPa(-15 m水頭)。頂部邊界設(shè)置為:①算恒定的壓力邊界條件為uw=0 kPa(例1);②恒定的流量邊界條件為q=0.75 ksat(例2);③算恒定的流量邊界條件為q=2 ksat(例3)。
為了合理比較FLAC3D和SEEP/W的計(jì)算結(jié)果,數(shù)值模擬過程中還注意了以下條件:①相同的單元形狀和尺寸;②相同的土水特征曲線和滲透系數(shù)函數(shù);③相同的初始條件(各節(jié)點(diǎn)孔隙水壓力均為-147 kPa,即-15 m水頭);④相同的時(shí)步(10 s)。每組算例中均模擬了4種不同時(shí)長(zhǎng),即:20 000,40 000,60 000,80 000 s。針對(duì)不同時(shí)長(zhǎng),分別比較了2個(gè)程序的計(jì)算結(jié)果。
圖3—圖5為利用FLAC3D和SEEP/W模擬的孔隙水壓力和體積含水量(反映浸潤(rùn)范圍)在不同的模擬時(shí)間下隨深度的分布情況??梢钥闯觯涸诓煌乃憷?,分別利用修正后FLAC3D和SEEP/W進(jìn)行數(shù)值計(jì)算,在不同的計(jì)算時(shí)間下二者的模擬結(jié)果均較為相符。模擬結(jié)果表明運(yùn)用本文方法對(duì)FLAC3D的滲流模型的修正合理,且導(dǎo)入的土水特征曲線有效。2個(gè)程序計(jì)算結(jié)果的微小差異可能是由于有限差分法和有限元法計(jì)算方法精度的不同所導(dǎo)致。
圖3 uw為0條件下的體積含水量和孔隙水壓力隨深度的分布曲線Fig.3 Distributions of volumetric water content and pore-water pressure w ith depth(uw=0 kPa)
圖4 q為0.75ksat條件下的體積含水量和孔隙水壓力隨深度的分布曲線Fig.4 Distributions of volumetric water content and pore-water pressure w ith depth(q=0.75ksat)
圖5 q為2ksat條件下的體積含水量和孔隙水壓力隨深度的分布曲線Fig.5 Distributions of volumetric water content and pore-water pressure w ith depth(q=2ksat)
模擬結(jié)果表明:
(1)飽和一維入滲條件下,在浸潤(rùn)前鋒附近,土體含水量由近飽和狀態(tài)迅速降低(圖3)。數(shù)值模擬結(jié)果與Lumb[9]的現(xiàn)場(chǎng)觀測(cè)和分析解呈現(xiàn)的規(guī)律一致。
表3 流固耦合計(jì)算中所用力學(xué)參數(shù)Table 3 Mechanical parameters adopted in the fluid-mechanical coupled analysis
(2)當(dāng)入滲速度小于飽和滲透系數(shù)時(shí),其浸潤(rùn)區(qū)僅損失部分負(fù)孔隙水壓力。在浸潤(rùn)前鋒,含水量亦迅速降低(圖4)。
(3)如圖5所示,在邊界條件q=2ksat下,2個(gè)程序的模擬結(jié)果也較為相符,孔隙水壓力由正值平滑降為負(fù)值。該結(jié)果說明修正后該滲流模型從飽和到非飽和入滲過渡良好。
(4)體積含水量與孔隙水壓力在浸潤(rùn)區(qū)前緣的降低速度受所采用的土水特征曲線和滲透系數(shù)曲線影響(圖3—圖5)。其影響程度應(yīng)進(jìn)一步參照試驗(yàn)數(shù)據(jù)進(jìn)行深入研究。
在實(shí)現(xiàn)滲流計(jì)算的基礎(chǔ)上,可進(jìn)一步利用FLAC3D進(jìn)行流固耦合分析。對(duì)于非飽和土,Bishop[15]將太沙基的飽和土有效應(yīng)力公式做如下擴(kuò)展,即
式中:χ為與飽和度s有關(guān)的參數(shù)(作為簡(jiǎn)化,可取χ=s);ua為孔隙氣壓力。
利用公式(4)修正FLAC3D中的總應(yīng)力與有效應(yīng)力的關(guān)系后,將可在程序中描述非飽和土的應(yīng)力應(yīng)變關(guān)系。需注意的是,孔隙氣壓力ua與孔隙水壓力相關(guān)聯(lián),非飽和土ua=0,飽和土ua=uw。非飽和土的摩爾-庫(kù)侖強(qiáng)度公式可表示[15]為
式中c′和?′分別為有效黏聚力和有效內(nèi)摩擦角。
結(jié)合前面對(duì)FLAC3D滲流模型的擴(kuò)展,非飽和土有效應(yīng)力公式(4),強(qiáng)度公式(5),以及比奧固結(jié)理論[16],即可利用FLAC3D內(nèi)置的耦合模型進(jìn)行非飽和土流固耦合分析(弱耦合)。
運(yùn)用該方法,筆者模擬了灌溉水沿黃土裂縫入滲導(dǎo)致的土體破壞過程[14],包括裂縫中水位上升、維持、下降和消散的過程(表2)。限于篇幅,這里僅選取灌溉水完全沿裂縫入滲的典型算例簡(jiǎn)要介紹如下。計(jì)算網(wǎng)格見圖6(a),所選滲流參數(shù)及力學(xué)參數(shù)見表1、表3。當(dāng)裂縫中水位降低到裂縫底部后,其孔壓場(chǎng)、位移場(chǎng)及塑性區(qū)的最終累計(jì)量如圖6(b)和圖6(c)所示(對(duì)稱條件下僅取左部分)。可以看出:浸潤(rùn)區(qū)向深部發(fā)展,部分土體單元達(dá)到了塑性狀態(tài),且土體沿裂縫臨空面向外移動(dòng)。數(shù)值模擬結(jié)果反映了現(xiàn)場(chǎng)灌水時(shí)土體中滲流、變形破壞的趨勢(shì),取得了較好的應(yīng)用效果。
表2 裂縫參數(shù)及灌溉模擬過程Table 2 The param eters of the crack and the simulated processes of irrigation
圖6 水沿裂縫入滲過程的流固耦合分析圖Fig.6 Fluid-mechanical coupled analysis of soil failure along the crack subjected to water infiltration
本文利用FLAC3D內(nèi)置的FISH語(yǔ)言,對(duì)其滲流模型進(jìn)行了擴(kuò)展。文中所述內(nèi)容概括如下:
(1)將實(shí)驗(yàn)室測(cè)得的原狀黃土土水特征曲線,利用MATLAB實(shí)現(xiàn)了參數(shù)擬合。根據(jù)擬合結(jié)果,該黃土殘余含水量為8.7%,飽和含水量為40.7%。
(2)在FLAC3D中,將擬合參數(shù)編寫為子程序,并在滲流分析時(shí)迭代調(diào)用,可實(shí)現(xiàn)非飽和土的滲流模擬。運(yùn)用該方法模擬了不同邊界條件下的一維入滲過程。驗(yàn)證表明,該方法準(zhǔn)確。
(3)模擬結(jié)果顯示,在浸潤(rùn)前鋒附近,土體含水量由近飽和狀態(tài)迅速降低。模擬結(jié)果與前人的現(xiàn)場(chǎng)觀測(cè)和分析解呈現(xiàn)的規(guī)律一致。
(4)基于以上非飽和滲流模型,F(xiàn)LAC3D的流固耦合模塊可進(jìn)一步擴(kuò)展應(yīng)用于非飽和土流固耦合分析。本文的研究成果對(duì)研究非飽和土理論,及擴(kuò)展FLAC3D在非飽和土數(shù)值計(jì)算方面均具有一定意義。
[1] 毛昶熙.滲流計(jì)算分析與控制[M].北京:中國(guó)水利水電出版社,2003:1-2.(MAO Chang-xi.Seepage Computation and Analysis and Control[M].Beijing:China Water Power Press,2003:1-2.(in Chinese))
[2] NG CWW,SHIQ A.Numerical Investigation of the Stability of Unsaturated Soil Slopes Subjected to Transient Seepage[J].Computers and Geotechnics,1998,22(1):1-28.
[3] 王鐵行.非飽和黃土路基水分場(chǎng)的數(shù)值分析[J].巖土工程學(xué)報(bào),2008,30(1):41-45.(WANG Tie-hang.Moisture Migration in Unsaturated Loess Subgrade[J].Chinese Journal of Geotechnical Engineering,2008,30(1):41-45.(in Chinese))
[4] YOO C,KIM SB.Three-Dimensional Numerical Investigation of Multifaced Tunneling in Water-Bearing Soft Ground[J].Canadian Geotechnical Journal,2008,45(10):1467-1486.
[5] ITASCA.FLAC 3D Version 2.61 User’s Guide[K].Minneapolis,Minnesota,USA:ITASCA Inc.,2002.
[6] 李 毅,伍 嘉,李 坤.基于FLAC3D的飽和-非飽和滲流分析[J].巖土力學(xué),2012,33(2):617-622.(LIYi,WU Jia,LIKun.Saturated-Unsaturated Seepage Analysis Based on FLAC3D[J].Rock and Soil Mechanics,2012,33(2):617-622.(in Chinese))
[7] RICHARDSR A.Capillary Conduction of Liquid Through Porous Mediums[J].Physics,1931,1(5):318-333.
[8] PHILIP JR.The Theory of Infiltration 1:The Infiltration Equation and Its Solution[J].Soil Science,1957,83(5):345-358.
[9] LUMB P.Effect of Rainstorms on Slope Stability[M].Hong Kong:Hong Kong Joint Group of the Institutions of Civil,Mechanical and Electrical Engineers,1962.
[10]FREDLUND D G,RAHARDJO H.Soil Mechanics for Unsaturated Soil[M].US:Wiley,1993.
[11]LEONG EC,RAHARDJO H.Review of Soil-Water Characteristic Curve Equations[J].Journal of Geotechnical and Geoenvironmental Engineering,1997,123(12):1106-1117.
[12]ZHOU Y D,CHEUK CY,THAM LG.NumericalModelling of Soil Nails in Loose Fill Slope under Surcharge Loading[J].Computers and Geotechnics,2009,36:837-850.
[13]李志清,李 濤,胡瑞林,等.非飽和土土水特征曲線(SWCC)測(cè)試與預(yù)測(cè)[J].工程地質(zhì)學(xué)報(bào),2007,15(5):700-707.(LI Zhi-qing,LI Tao,HU Rui-lin,et al.Methods for Testing and Predicting of SWCC in Unsaturated Soil Mechanics[J].Journal of Engineering Geology,2007,15(5):700-707.(in Chinese))
[14]ZHOU Y F.Study on Landslide in Loess Slope Due to Infiltration[D].Hong Kong:The University of Hong Kong,2012.
[15]BISHOP A W.The Principle of Effective Stress[J].Teknisk Ukeblad,1959,106(39):859-863.
[16]李廣信.高等土力學(xué)[M].北京:清華大學(xué)出版社,2006.(LIGuang-xin.Advanced SoilMechanics[M].Beijing:Tsinghua University Press,2006.(in Chinese) )
(編輯:姜小蘭)
Imp lementation and Verification of Unsaturated Seepage Analysis in FLAC3D
ZHOU Yue-feng1,2,THAM L G2,YANW M2,F(xiàn)AN Shao-peng2,3
(1.Key Laboratory of Geotechnical Mechanics and Engineering of Ministry ofWater Resources,Yangtze River Scientific Research Institute,Wuhan 430010,China;2.Department of Civil Engineering,The University of Hong Kong,Hong Kong,China;3.Hydraulic Design Department,Changjiang Institute of Survey,Planning,Design and Research,Wuhan 430010,China)
The seepagemodel in FLAC3D was extended using its in-built programming language FISH.The SWCC(SoilWater Characteristic Curve)for seepage in unsaturated loess was best-fitted using MATLAB.The obtained parameterswere written as a subroutine and invoked in each step in the seepagemodel of FLAC3D to simulate the unsaturated loess.One-dimensional infiltration problem in three boundary conditionswasmodeled using both the finite differencemethod and the finite elementmethod to verify the extension in FLAC3D.The comparison results show that the abovemethod is appropriate.The findings of this paper aremeaningful in the investigation of the unsaturated soil theory and in the application of FLAC3D to simulate unsaturated soil problems.
unsaturated soil;seepage;SWCC;FLAC3D;redevelopment
TV139.1
A
1001-5485(2013)02-0057-05
10.3969/j.issn.1001-5485.2013.02.012
2012-10-11;
2012-11-22
香港研究資助局資助項(xiàng)目(HKU7140/08E)
周躍峰(1982-),男,山西侯馬人,工程師,博士,主要從事黃土滑坡機(jī)理研究與巖土工程數(shù)值計(jì)算研究,(電話)13971606626(電子信箱)zhouyuefenghku@gmail.com。