擋水結構地基滲流問題的隨機有限元分析

潘健，何志虎

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

擋水結構地基滲流問題的隨機有限元分析

潘健，何志虎

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

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

滲透系數；變異系數；隨機場；隨機有限元

0 引言

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

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

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

1 壩基滲流模型及其有限元網格

1.1 壩基滲流分析模型

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

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

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

1.2 有限元網格劃分

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

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

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

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

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

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

有限元程序用拉普拉斯方程求解.單元傳導矩陣往往被假定為一個整體矩陣，假設后的整體傳導關系變為：

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

隨機滲流場特征由3個參數確定，分別是均值μk、標準差σk和相關距離θk.為了得到更為合理穩定的統計特征，每個參數組合將會用1 000份樣本點作統計分析，這樣可以計算出比較精確的結果.

1.3 滲流勢能分布

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

式中，ki是分配到第i個單元的滲透系數，gi是在第i個單元中標準高斯隨機場的局部均值，μlnk和σlnk分別是ki基于對數的均值和標準差.

2 確定性分析

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

圖3 確定性分析流網圖（nf=5，nd=20）

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

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

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

圖4 下游防滲墻示意圖

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

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

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

3 隨機有限元分析

3.1 勢能場的統計特征

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

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

圖5 等勢線和勢能標準差示意圖（σk/μk＝1，θk=1.0 m）

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

3.2 參數分析

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

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

分析所有取值后的結果，以σk/μk的常用對數形式反映在圖6~8中.、和ie的均值和標準差分別記作

圖6 滲透性變異系數對流量的影響

3.2.1 流量分析

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

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

3.2.2 揚壓力分析

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

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

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

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

3.2.3 出逸梯度分析

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

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

4 結語

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

1）流量的均值對于相關距離的變化相對不敏感，但是隨其增加而持續下降，這對擋水結構的設計有一定的參考意義.流量的標準差會隨著相關距離的增加先增大后減小，但其峰值的具體位置尚待研究.

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

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

[1] 松尾稔［日］.地基工程學可靠性設計的理論和實際[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，Freeze 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] 雷桂媛.關于蒙特卡羅及擬蒙特卡羅方法的若干研究[D].杭州：浙江大學博士論文，2003：2-14.

[6] 程強，羅書學，彭雄志.相關距離與土性參數的關系及計算方法[J].西南交通大學學報，2000，35（5）：55-59.

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

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

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

[10] Smith L，Freeze 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

潘健（1963-），男，廣東廣州人，博士，副教授.研究方向：巖土工程地下結構設計方法和風險評估.E-mail：cvpan@scut.edu.cn