擋水結(jié)構(gòu)地基滲流問題的隨機(jī)有限元分析

潘健，何志虎

（華南理工大學(xué)土木與交通學(xué)院，廣東廣州510640）

擋水結(jié)構(gòu)地基滲流問題的隨機(jī)有限元分析

潘健，何志虎

（華南理工大學(xué)土木與交通學(xué)院，廣東廣州510640）

采用蒙特卡羅模擬的隨機(jī)有限元方法，并結(jié)合滲流隨機(jī)場的概念，對設(shè)置有兩道防滲墻的大壩滲流問題進(jìn)行分析.研究了滲透系數(shù)標(biāo)準(zhǔn)差和空間相關(guān)距離對壩基滲流參數(shù)的影響，對比了確定性分析和隨機(jī)有限元分析.結(jié)果表明，土體滲透系數(shù)標(biāo)準(zhǔn)差和空間相關(guān)距離對流量、揚壓力和出逸梯度有較顯著影響.

滲透系數(shù)；變異系數(shù)；隨機(jī)場；隨機(jī)有限元

0 引言

如果土的參數(shù)假定為均值常數(shù)，則巖土工程分析大多屬于確定性問題.這種確定性方法不僅用于經(jīng)典土力學(xué)計算，也廣泛應(yīng)用于有限元數(shù)值模擬分析.實際上由于土中各點的性質(zhì)都不一樣，要得到比較準(zhǔn)確的值，只有通過大量的原位測試實驗數(shù)據(jù)，但這樣做既昂貴又不切實際.隨機(jī)場模型可以用于分析巖土材料的性質(zhì)，參數(shù)也可以從有限的實驗測量數(shù)據(jù)中分析估得.

近些年來，在巖土工程分析中，土體的不確定性分析越來越得到重視[1]，且隨著隨機(jī)有限元的發(fā)展和計算機(jī)功能的強(qiáng)大，巖土工程的一些領(lǐng)域也開始運用隨機(jī)有限元分析方法.例如Beacher和Ingra[2]的應(yīng)力分析與地基沉降研究，Ishii和Suzuki[3]對邊坡穩(wěn)定性的分析，Smith和Freeze[4]對有壓滲流的研究等.

本文把對隨機(jī)場進(jìn)行離散的局部平均法與確定性有限元法相結(jié)合，用蒙特卡羅法[5]進(jìn)行滲流分析.蒙特卡羅隨機(jī)有限元法利用統(tǒng)計抽樣理論研究隨機(jī)變量，計算結(jié)果比較精確，而且不受隨機(jī)變量變異性大小的限制，因此，對于變異性相對較大的滲透系數(shù)分析，采取蒙特卡羅方法是比較合適的.

1 壩基滲流模型及其有限元網(wǎng)格

1.1 壩基滲流分析模型

研究大壩地基的滲流時，必須考慮3個重要參數(shù)：地基滲流量、梯度壓力和揚壓力，分析模型如圖1所示.

分析壩基下土體時，對于一定范圍的土體，隨著“相關(guān)距離”[6]的增加，土體的性質(zhì)會變得更加均一.本文采取局部平均法[7]，可以在給定均值、標(biāo)準(zhǔn)差和空間相關(guān)性條件下實現(xiàn)隨機(jī)滲流場的分析.滲流場的統(tǒng)計特征都反映在有限元網(wǎng)格的劃分、勢函數(shù)和流函數(shù)的邊界條件中.

圖1 大壩下的邊界條件示意圖

1.2 有限元網(wǎng)格劃分

穩(wěn)定滲流問題一般都可以用二維拉普拉斯方程求解，方程中獨立變量?是指x-y坐標(biāo)平面中任一點的測壓管水頭位置的勢能.方程如下：

式中，kx和ky分別是x，y方向的滲透系數(shù).

滲流場都假定為各向同性（kx=ky=k）.雖然這里討論的方法是簡單地通過隨機(jī)場理論擴(kuò)展到各向異性，但這種擴(kuò)展必須限定在一個特定的條件下.

方程（1）只有在取k為常數(shù)時是唯一有效的，滲透系數(shù)將在每個單元中作為常數(shù)，它的取值來自于局部一定范圍單元的幾何平均滲透系數(shù)，不同的單元k值也會有所不同，這樣就反映了滲流性的隨機(jī)性質(zhì).

本文選用一種典型的有限元網(wǎng)格劃分，如圖2所示.它包含了1 400個單元，是一個有兩道防滲墻的大壩下的二維滲流模型.上游水位和下游水位的勢能值分別為10 m和0 m.兩個防滲墻假設(shè)為零厚度，并且防滲墻兩側(cè)的節(jié)點有兩個相應(yīng)的左右側(cè)勢能值.

圖2 有限元網(wǎng)格劃分示意圖（單元尺寸：0.2 m×0.2 m）

有限元程序用拉普拉斯方程求解.單元傳導(dǎo)矩陣往往被假定為一個整體矩陣，假設(shè)后的整體傳導(dǎo)關(guān)系變?yōu)椋?/p>

式中，K為整體剛度矩陣，Φ為整體勢能矢量，Q為整體流量矢量.一旦整體傳導(dǎo)方程求解后再推導(dǎo)出節(jié)點勢能值，那么相關(guān)的輸出量流量、上揚壓力和梯度壓力就能很容易地推導(dǎo)計算出來.

隨機(jī)滲流場特征由3個參數(shù)確定，分別是均值μk、標(biāo)準(zhǔn)差σk和相關(guān)距離θk.為了得到更為合理穩(wěn)定的統(tǒng)計特征，每個參數(shù)組合將會用1 000份樣本點作統(tǒng)計分析，這樣可以計算出比較精確的結(jié)果.

1.3 滲流勢能分布

滲透系數(shù)的現(xiàn)場測試結(jié)果表明，它近似服從一個對數(shù)正態(tài)分布[8-9].因此，本文的模型也采取同樣的分布.滲流場通過如下公式轉(zhuǎn)化：

式中，ki是分配到第i個單元的滲透系數(shù)，gi是在第i個單元中標(biāo)準(zhǔn)高斯隨機(jī)場的局部均值，μlnk和σlnk分別是ki基于對數(shù)的均值和標(biāo)準(zhǔn)差.

2 確定性分析

確定性分析時，所有單元的滲透系數(shù)都假定為一個恒定值且等于10-5m/s.取這個值是因為它是隨后隨機(jī)分析得到的平均值.勢能和流線的反推可以推出圖3所示的流網(wǎng)圖.其中流槽數(shù)nf=5，等勢線間隔數(shù)nd=20.

圖3 確定性分析流網(wǎng)圖（nf=5，nd=20）

所有的輸出量都采取無量綱計算.就流量來說，整體流向量Q是通過計算整體傳導(dǎo)矩陣的矩陣方程（2）得到的.這些值最后都計入到流量Q（m3/s）中，由此推導(dǎo)出一個無量綱，定義為：

式中，μk是滲透性的均值（各向同性），H是上下游的總水頭差高度.大壩基礎(chǔ)的揚壓力U的計算要結(jié)合兩防滲墻間沿著壩基的壓力分布情況，它可以通過一個數(shù)值積分和節(jié)點的勢能值很容易地計算出來.無量綱揚壓力定義為：

式中，γw為水的重度，L是兩防滲墻的間距，是大壩淹沒在水中時作為浮力出現(xiàn)的揚壓力.下游防滲墻如圖4所示.出逸坡降ie是最接近大壩下游尾端的溢出口處初始值，采用四點反向差分法準(zhǔn)則計算：

圖4 下游防滲墻示意圖

式中，?i對應(yīng)于圖4所示出口點下垂直線上的4個節(jié)點處測壓管水頭，b是節(jié)點間恒定的垂直距離.下游處初始勢能值設(shè)為0，即?0=0.0 m.

圖1、圖2中，上述常數(shù)值分別為：H=10 m，μk＝10-5m/s，γw=9.81 kN/m3，L=6 m，b=0.2 m.

考慮到大多數(shù)土體的臨界水力梯度近似相等或一致，上面ie的取值在實際情況中是不能被接受的.ie值大小和壓頭差H成正比，H在此實例中被簡化設(shè)定為10 m.

3 隨機(jī)有限元分析

3.1 勢能場的統(tǒng)計特征

圖5給出了由1 000個取θk=1.0 m的樣本得到的勢能場標(biāo)準(zhǔn)差和均值輪廓線.這里采用了由Smith和Freeze通過一個一維和一個二維承壓流問題的模擬實驗得到的一系列數(shù)值結(jié)果[10].圖5（a）中的等勢線非常類似于圖3流網(wǎng)確定性分析得到的等勢線，圖5（b）中勢能的標(biāo)準(zhǔn)差顯示了關(guān)于勢能值最大不確定性的區(qū)域.

在隨后所有的二維隨機(jī)分析中，土被假定為各向同性，隨機(jī)場把每個有限元單元的滲透性幾何均值取為常數(shù)10-5m/s.參數(shù)分析時會考慮到土滲流場的不同標(biāo)準(zhǔn)差（σk）和相關(guān)距離（θk）的影響.每組參數(shù)都根據(jù)1 000個樣本點統(tǒng)計分析，最后計算得到了統(tǒng)計

圖5 等勢線和勢能標(biāo)準(zhǔn)差示意圖（σk/μk＝1，θk=1.0 m）

應(yīng)該指出的是，上游和下游邊界上的勢能值是確定性的，所以在這些邊界上的勢能標(biāo)準(zhǔn)差等于0.最大的標(biāo)準(zhǔn)差出現(xiàn)在流型的中間，流型反映了兩防滲墻間大壩下的那塊區(qū)域.該區(qū)域的標(biāo)準(zhǔn)差實際上為一常數(shù).圖5中顯示的勢能場的統(tǒng)計特征與后面要考慮的揚壓力的統(tǒng)計特征有著密切關(guān)系.

3.2 參數(shù)分析

基于圖2的網(wǎng)格單元劃分的參數(shù)研究顯示了滲透系數(shù)的標(biāo)準(zhǔn)差（σk）和相關(guān)距離（θk）對輸出的流量、和ie的影響.在所有情況下，滲透性的均值μk保持為一個恒定值10-5m/s.

為了便于分析σk的影響，滲透性的變化采用無量綱系數(shù)而且考慮如下取值：σk/μk=0.125，0.25，0.50，1，2，4，8，16；相關(guān)距離考慮如下取值：θk=0，1，2，4，8和∞（m）.

分析所有取值后的結(jié)果，以σk/μk的常用對數(shù)形式反映在圖6~8中.、和ie的均值和標(biāo)準(zhǔn)差分別記作

圖6 滲透性變異系數(shù)對流量的影響

3.2.1 流量分析

由圖6（a）可知，在θk＜8 m時，隨著σk/μk的增加，μQ會顯著下降.隨著θk接近無限大，μQ也逐漸接近一個恒定值0.226.當(dāng)θk=8 m，隨著σk/μk增加，μQ也會減少，但減幅不如θk＜8 m時大.盡管如此，對于典型的相關(guān)距離，平均流量的影響是微小的.作為土體量變化的一個函數(shù)，流量的增加對設(shè)計者來說是一個很重要的考慮因素.

從圖6（b）可以看出σk/μk變化時，σQ的變化情況.特別值得注意的是，在θk≤8 m的情況下，σQ達(dá)到最大值時，σk/μk的范圍是1.0~2.0.同樣，當(dāng)θk=∞時，隨著σk/μk變大σQ持續(xù)變大.而σQ最大變化發(fā)生在這樣一種典型條件下：相關(guān)距離θk在1~4之間而且滲透性變異系數(shù)σk/μk大約為1或者2（即log（σk/μk）為0或0.3左右）.

3.2.2 揚壓力分析

由圖7（b）可知，當(dāng)σk/μk和θk增加時，σU不斷增加，當(dāng)θk→∞，σU→0時，滲流場變得完全一致.θk變大時曲線也明顯變陡.

對于給定的σk/μk和θk值，σU的實際值能輕易地從勢能值的標(biāo)準(zhǔn)差中推斷出來.圖5（b）列舉了對于特定值σk/μk=1.0，θk=1.0 m整個流域的勢能值標(biāo)準(zhǔn)差的輪廓.圖5（b）中，大壩下方的勢能標(biāo)準(zhǔn)差近似為常數(shù)0.8.在無量綱除以H=10 m后，這些值就很接近圖7（b）中對應(yīng)的值.

圖7（b）中給出的揚壓力標(biāo)準(zhǔn)差最大值也不是很大，這就意味著揚壓力可以在一個合理的置信度下被估計出來，因為揚壓力是通過大壩下大量節(jié)點的勢能值計算出來的.這個取平均的過程將會使局部勢能值波動更小，結(jié)果導(dǎo)致了方差的減小.

圖7 滲透性變異系數(shù)對揚壓力的影響

3.2.3 出逸梯度分析

圖8（a）、（b）給出了考慮參數(shù)變化時，μie和σie的變化情況.ie對σk/μk變化的敏感性顯示得非常清楚.在圖8（a）中，當(dāng)σk/μk值在0.0~1.0間波動時，μie基本接近一個確定值0.668，但是超出這個范圍，μie就明顯不穩(wěn)定，和恒定值發(fā)生偏離.值得注意的是，對于θk≤1，隨著σk/μk增加，μie曲線變化趨勢為向下，而且一直在恒定值下方.而對于較大的θk，μie值會隨著σk/μk的增加而增大，且一直高于恒定值.當(dāng)用于預(yù)測出逸梯度的差分長度為0.6 m時候，在相同σk/μk范圍內(nèi)，μie在θk取2、4、8 m時的變化比θk取0和1時大得多.

圖8 滲透性變異系數(shù)對出逸梯度的影響

4 結(jié)語

對在隨機(jī)土地基上設(shè)置兩道防滲墻的擋水結(jié)構(gòu)滲流問題進(jìn)行分析，考慮了均值、標(biāo)準(zhǔn)差和相關(guān)距離對輸出變量的影響，并把這些值映射到包含有1 400個單元的有限元網(wǎng)格中，對每一組參數(shù)的確定都采取1 000個樣本點進(jìn)行分析.將有限單元的平均滲透系數(shù)設(shè)為常數(shù)并對變異系數(shù)和相關(guān)距離進(jìn)行研究，得到了流量、揚壓力和出逸梯度3個輸出變量，其中流量和揚壓力采取了無量綱化.結(jié)果表明：

1）流量的均值對于相關(guān)距離的變化相對不敏感，但是隨其增加而持續(xù)下降，這對擋水結(jié)構(gòu)的設(shè)計有一定的參考意義.流量的標(biāo)準(zhǔn)差會隨著相關(guān)距離的增加先增大后減小，但其峰值的具體位置尚待研究.

2）揚壓力均值對參數(shù)變化非常不敏感，最大的變化也只有大約10%.然而，揚壓力的標(biāo)準(zhǔn)差均隨著相關(guān)距離和變異系數(shù)的增大而不斷增大.

3）出逸梯度均值對參數(shù)變化甚為敏感，基于溢出點處與長度相關(guān)的測壓水頭的一階導(dǎo)數(shù)分析，出逸梯度均值對于隨機(jī)場生成的勢能值的局部變化非常敏感.通過四點數(shù)值微分公式引入局部平均法，并算出了一些局部均值，但均值波動性和標(biāo)準(zhǔn)差還需要進(jìn)一步研究.

[1] 松尾稔［日］.地基工程學(xué)可靠性設(shè)計的理論和實際[M].北京：人民交通出版社，1990：43-62.

[2] Baecher G B，Ingra T S.Stochastic FEM in settlement predictions[J].Geotech Eng ASCE，1981，107（GT4）：449-463.

[3] Ishii K，Suzuki M.Stochastic finite element method for slope stability analysis[J].Structural Safety，1987（4）：111-129.

[4] Smtih L，F(xiàn)reeze R A.Stochastic analysis of steady state groundwater flow inabounded domain.1.one dimensional simulation[J].Water Resour Res，1979，15（3）：521-528.

[5] 雷桂媛.關(guān)于蒙特卡羅及擬蒙特卡羅方法的若干研究[D].杭州：浙江大學(xué)博士論文，2003：2-14.

[6] 程強(qiáng)，羅書學(xué)，彭雄志.相關(guān)距離與土性參數(shù)的關(guān)系及計算方法[J].西南交通大學(xué)學(xué)報，2000，35（5）：55-59.

[7] 朱位秋，任永堅.隨機(jī)場的局部平均與隨機(jī)有限元[J].航空學(xué)報，1986，7（7）：604-611.

[8] 施小清，吳吉春，袁永生.滲透系數(shù)空間變異性研究[J].水科學(xué)進(jìn)展，2005，16（2）：210-215.

[9] 束龍倉，李偉.北塘水庫庫底地層滲透系數(shù)的隨機(jī)特分析[J].吉林大學(xué)學(xué)報（地球科學(xué)版），2007，37（2）：216-220.

[10] Smith L，F(xiàn)reeze R A.Stochastic analysis of steady state groundwater flow in abounded domain，2.two-dimensional simulations[J].Water Resour Res，1979，15（6）：1 543-1 559.

Seepage Beneath Water Retaining Structures Based on Stochastic Finite Element Analysis

PAN Jian，HE Zhi-hu

(College of Civil Engineering and Transportation，South China University of Technology，Guangzhou 510640，Guangdong，China)

Random field concepts for the generation of soil permeability properties with finite element methods were used to perform Monte Carlo simulations of the seepage problem.Analyses have been performed for the case of a dam with two cut-off walls.The effect of the standard deviation and correlation structure of the permeability on the output statistics relating to seepage quantities，exit gradients and uplift pressures are studied.In all cases，comparisons are made with results that would be achieved on a deterministicbasis.Flowratesandotherquantitiesofinterestareshowntobe significantly affected by both the standard deviation and the correlation structure of soil permeability.

permeability；variance；random fields；stochastic finite element method

TV 640

A

1001-4217（2010）03-0066-08

2010-04-14

潘?。?963-），男，廣東廣州人，博士，副教授.研究方向：巖土工程地下結(jié)構(gòu)設(shè)計方法和風(fēng)險評估.E-mail：cvpan@scut.edu.cn