考慮材料參數空間變異性的堆石壩非侵入式隨機有限元研究

王建娥， 楊 杰， 程 琳， 馬春輝， 冉 蠡

(1. 西安理工大學 水利水電學院， 陜西 西安 710048； 2.西安理工大學 省部共建西北旱區生態水利國家重點實驗室， 陜西 西安 710048)

1 研究背景

隨著國內外水利水電資源的不斷開發，堆石壩由于可就地取材、地形適應能力強以及穩定性好等優點，成為高壩大庫的首選壩型。目前堆石壩正向300m級高度發展。堆石壩的沉降控制是此壩型設計研究的重點和難點，更是確保工程建設、運營安全的重要依據。隨著計算機技術的迅猛發展，各種計算分析方法層出不窮，其中有限元法被廣泛應用于堆石壩的數值分析中。

然而，一方面，堆石壩中的堆石材料是一種非常復雜的工程材料，其材料參數物理力學特性隨機性強且變異性大。另一方面，由于試驗測量誤差、數值計算理論以及本構模型的局限無法完全模擬實際工程中材料的受力情況等以及由測量誤差所帶來的統計不確定性等均對計算結果的精度有很大影響。因此，在常規有限元計算中假定大范圍的筑壩材料參數為唯一、確定的參數值顯然無法反映工程實際，造成計算值與實測值存在較大的誤差。因此，若能在有限元計算中考慮堆石料材料參數的空間變異性和隨機性，有助于提高堆石壩沉降控制水平，優化堆石壩設計標準、指導大壩施工和運行，從而推動壩工理論的發展。

近年來，眾多國內外專家學者致力于在數值計算中考慮材料參數空間變異性的影響，高大釗[1]、包承綱[2]、陳祖煜[3]等均認為應該考慮巖土工程中存在的不確定性因素。考慮材料參數空間變異性的理論分為隨機變量理論和隨機場理論。其中，隨機變量理論是通過離散試驗點的變異性來模擬土體參數變異性，該方法假定研究尺度內任意點處的參數相互之間完全無關，無法考慮空間不同點處局部與整體物理力學性質之間的差異性。相比之下，隨機場理論能夠較好地描述這種空間變異性。Vanmarcke[4]首先提出土體的隨機場模型，并提出用波動范圍表征土體參數的空間變異性。Griffiths等[5]通過蒙特卡羅模擬(MCS)，實現了邊坡可靠度分析的隨機有限元方法；Sett等[6]研究了考慮土石料參數空間變異性的動力反應問題；Low等[7]采用隨機方法計算土石邊坡安全系數和可靠度；在國內，蔣水華等[8]、李典慶等[9]開展了邊坡可靠度分析與隨機有限元方法結合的研究，并進行了隨機多項式展開計算邊坡安全系數的顯式表達式的研究；祁小輝等[10]采用譜展開法在考慮土體空間變異性的基礎上，進行了邊坡最危險滑動面的隨機分析。目前多數隨機場模擬方法存在嚴格的適用條件，如局部平均法等只適用于矩形或四邊形區域；當研究區域不規則、協方差函數形式復雜時，二維第二類Fredholm積分方程求解存在一定困難。相較而言，Cholesky分解方法具有編程簡單、適用于不規則邊界區域、可提高離散精度等優勢，對工程有較強適用性。

上述隨機有限元研究主要針對邊坡工程進行，而對于堆石壩材料參數的空間變異性研究還鮮有涉及。由于堆石壩材料分區較多且體型較大等原因，致使堆石壩材料參數難免存在空間差異。另一方面，國內外現行堆石壩填筑規范中，通常只對堆石料的干密度以及表征密實程度參數給出了較明確的指標，而對堆石料的級配和母巖性質等的控制則相當粗放。許多研究[11]表明，級配不同時，密實程度相同的同一種堆石料力學特性也有較大差別。因此，雖然筑壩過程中對材料進行了質量控制，但材料的物理力學性質仍具有一定的空間變異性和隨機性。

為提高堆石壩沉降控制水平，將隨機場理論引入堆石壩沉降值計算中，建立基于Cholesky分解的隨機場離散方法。隨后，建立堆石壩應力應變分析的非侵入式隨機有限元法，以考慮本構模型參數的自相關性以及互相關性。最后，嘗試探討巖土體參數的變異系數對堆石壩壩體沉降值的影響規律。

2 隨機場模擬原理

2.1 隨機場參數

對隨機場的合理離散是隨機有限元的重要內容之一，不僅結構需被離散成有限個單元，模型參數也需被離散成一定數量的隨機場單元。參數隨機場與結構單元可采用不同網格體系，通常每個隨機場單元可包含1～2個結構單元[12]，本文隨機場單元與結構單元網格劃分一致。各參數的隨機場單元參數c=(ci,i=1,2,…,ne) (ne為單元總數)的統計特性可由均值μχi和標準差σχi表征，變異系數為：

C·V·(χi)=σχi/μχi

(1)

隨機場中任意兩點的自相關系數定義為[13]：

ρ[H(xi,yi),H(xj,yj)]=

(2)

式中：H(xi,yi)、H(xj,yj)為中心點空間坐標為(xi,yi)和(xj,yj)的單元參數隨機場模擬值； Var[H(xi,yi)]和Var[H(xj,yj)]為中心點空間坐標為(xi,yi)和(xj,yj)的單元方差。

一般采用理論自相關函數來描述堆石料同一本構參數中兩個不同空間位置間的自相關性。由于高斯型自相關函數平穩性較好，本文選用高斯型自相關函數：

(3)

式中：τx為空間中兩點在水平方向上的相對距離，τx=|xi-xj|；τy為空間中兩點在垂直方向上的相對距離，τy=|yi-yj| ；δh為水平波動范圍，δv為垂直波動范圍，均需根據地質統計或經驗確定。從表達式(3)可以看出，自相關函數只與空間中兩點的相對位置有關，而與兩點的絕對位置無關。

因此，同一參數隨機場的自相關系數矩陣為：

(4)

式中：χi為某一參數隨機場在計算區域內任意點處的隨機特性值，χi=H(xi,yi),i=1,2,…,ne。

2.2 相關對數正態隨機場模擬原理

在堆石料隨機場的模擬中，常用高斯隨機場來模擬，但對于數值大于0且參數間存在一定的相關性的隨機參數，直接采用高斯隨機場模擬是不合適的。對此，應在高斯隨機場的基礎上，采用相關對數隨機場來模擬堆石料材料參數的隨機場。

不同材料參數間的互相關性，可用互相關系數矩陣R表示，m個隨機場的互相關系數矩陣表示為：

R=(ρij)m×m

(5)

(6)

(7)

3 堆石壩SFEM實現

面板堆石壩由多個分區組成，各個分區材料均不同，每一種材料的E-B模型均由9個參數確定，若對全部材料分區及材料參數進行隨機場模擬，則計算時間較長，且部分參數對計算結果影響很小。為此，首先對E-B模型參數進行敏感性分析，確定對計算結果影響明顯的材料參數，有助于簡化隨機有限元計算，節約計算時間。眾多學者分析了E-B模型參數的敏感性，本文根據Zheng Dongjian等[15]采用Morris法的敏感性分析結論，即E-B模型參數K、φ0、Kb對壩體沉降影響較大，對不同材料分區的K、φ0、Kb參數進行隨機場的模擬，其余模型參數均采用試驗值。由于面板堆石壩各分區材料參數的隨機場特征值均不同，因此隨機場模擬需按照分區分別進行，最后將模擬結果賦值給各單元，并進行有限元計算。具體步驟為：

(1)根據敏感性分析結果，將敏感性較高的參數作為待模擬的隨機場參數，敏感程度較低的參數采用試驗數據；

(2)本文選用高斯型自相關函數表征堆石料參數的自相關特征。建立K、φ0、Kb3個參數的自相關系數矩陣，并分別確定墊層區、過渡區、3BI區、3BII區和3C區的均值、變異系數；

(3)鑒于拉丁超立方抽樣(LHS)具有適用范圍廣、抽樣估值穩定、樣本具有更好的代表性和均勻性等優點，根據文獻[16]，1000次拉丁超立方抽樣誤差很小，可以滿足計算要求，因此在滿足參數分布規律的基礎上，本文采用LHS進行1000次抽樣構建獨立標準正態隨機樣本矩陣ξ；

(5)重復步驟(3)與(4)，對墊層區、過渡區、3BI區、3BII區和3C區分別進行參數隨機場的離散，得到各個分區的參數隨機場離散值，然后將離散值賦值給各對應的單元，進行有限元計算。

其次，針對侵入式隨機有限元需修改有限元軟件計算程序，且存在工作量大、適用性和靈活性弱等問題，本文在現有商業有限元軟件MSC.Marc的基礎上，基于Fortran語言實現有限元軟件的二次開發，大大提高了計算效率。NSFEM方法可以使隨機過程和有限元計算獨立進行，更有利于NSFEM方法的推廣和使用。其計算流程如圖1所示。

圖1 基于MSC.Marc二次開發的非侵入式隨機有限元(NSFEM)計算方法流程

4 工程實例

本文用公伯峽堆石壩驗證所提方法，對該壩筑壩材料的E-B本構模型參數進行隨機場模擬。公伯峽堆石壩FEM模型風格劃分如圖2所示，模型采用空間8節點等參單元，共1 430個單元、2946個節點。計算中考慮了壩體分期填筑以及分期蓄水對沉降的影響。

圖2 公伯峽二維FEM網格劃分

本文將根據布設于壩中部的電磁式沉降儀測線ES2實測數據驗證模擬結果的可靠性，測線位于壩左0+130.00 m斷面，共計24個測點，如圖3所示。

隨機有限元計算過程分為：(1)確定材料參數，進行參數隨機場模擬；(2)在所得參數隨機場基礎上，對有限元軟件進行E-B本構模型二次開發；采用UltraEdit在DOS環境下調用Marc，進行“非侵入式”有限元計算；(3)隨機有限元結果與實測值對比，確定隨機場模擬的可靠性從而得到更精確的面板堆石壩應變值。

堆石壩各分區材料的隨機場參數取值如表1和表2所示，其中均值、變異系數以及標準差取值參考文獻[15]中公伯峽面板堆石壩的均值以及變化范圍得出。由于目前隨機場理論在堆石壩中的應用鮮有涉及，大多應用于邊坡可靠度分析中，Ronold[17]、程強等[18]、吳振君等[19]得出了不同地區邊坡工程的土體參數波動范圍的統計值，但由于堆石壩各分區的堆石料均由人為粒徑控制以及碾壓等，不適用于土體參數波動范圍的統計規律。對于同一分區的筑壩材料，在經過人為篩選、控制后，材料差異較小，因此可將堆石壩各分區的水平和垂直波動范圍取為所在分區的幾何尺寸。

圖3 公伯峽堆石壩左0+130.00 m剖面材料分區及測點分布圖(單位：m)

分區KμσC.V.KbμσC.V.φ0μσC.V.波動范圍/m墊層區10503150.313003900.349.49.880.2δv=122.5, δh=3過渡區10903270.38302490.350.410.080.2δv=127, δh=33B I9502850.35501650.353.610.720.2δv=122.5, δh=387.683B II16009600.6800800.158.435.040.6δv=91, δh=149.53C7202160.3485970.256.584.751.5δv=91, δh=146.97

表2 堆石料其他材料參數取值表

在確定隨機場參數的基礎上，提取面板堆石壩二維模型的節點坐標數據，按照本文第3節提出的步驟分別進行各個材料分區參數K、φ0、Kb隨機場的模擬。將各分區的模擬值賦值給相應的單元，得到整個壩體斷面的參數隨機場，其中隨機場模擬的一次實現結果如圖4所示。

圖4 參數隨機場模擬的一次實現結果

由圖4中可以看出，每個單元的材料參數取值均不同；由于本文分別對各個分區分別進行隨機場離散，因此同一個分區內同一參數在不同單元的取值服從高斯函數；由圖4也可看出，材料參數隨機場極端值較少，這是因為高斯型自相關函數穩定性較好。

采用基于MSC.Marc二次開發的隨機有限元方法進行隨機場的離散，分別利用1 000次隨機場離散結果進行有限元計算，1 000次計算結果均明顯大于常規FEM方法計算值，對1 000次模擬結果取均值，其沉降值如圖5所示。圖5為壩體斷面最大累計沉降等值線圖，采用隨機有限法計算得壩體蓄水完成時斷面最大沉降為70.12 cm，位于壩體中部偏向下游，這是由于3C區堆石質量標準略低于3BⅡ區。因此，壩體整體變形規律符合常規規律，本文提出的隨機有限元方法基本正確。圖6為常規有限元計算得到的斷面最大累計沉降等值線圖。通過比較圖5與6可知：(1)圖5中的等值線較圖6波動略明顯，這是由于非侵入式隨機有限元考慮了材料參數空間變異性，材料參數具有一定隨機性，使得不均勻沉降值增加所造成的；(2)圖5出現了零星的封閉等值線，表明考慮材料參數空間變異性之后，由于隨機性引起的不均勻沉降增加；(3)隨機有限元計算得到的壩體最大沉降值為70.12 cm，常規有限元計算值為60.05 cm，實測最大沉降值為68.52 cm，因此隨機有限元的計算結果大于常規有限元，更接近于實測值。

圖51000次隨機有限元計算的平均最大 圖6常規有限元計算的最大斷面

斷面累計沉降等值線圖(單位：m) 累計沉降等值線圖(單位：m)

將常規有限元計算值、測線ES2的實測沉降值與非侵入式隨機有限元計算值進行對比，如圖7所示。由圖7可知：(1)NSFEM計算結果明顯大于FEM，這是由于面板堆石壩在材料選用、施工填筑以及材料參數測定中有很大的隨機性，NSFEM按對數正態分布規律對材料參數進行離散，考慮了材料參數可能存在偏離均值的情況，因此所計算的沉降值大于常規有限元的結果；(2)目前的FEM方法計算值與實測沉降值有著明顯差距，尤其在壩體中部以及上部，實測沉降值明顯大于計算值，這符合現有堆石壩實測沉降值大于有限元計算值的實際情況；(3)NSFEM計算值與實測值比較接近，證明NSFEM方法通過考慮筑壩材料的隨機性和不均勻性，其計算值更符合實際工況；(4)在壩體中部以及中上部，NSFEM計算值與實測沉降值擬合較好，在壩體下部略有誤差，但整體擬合良好。

圖7 壩體斷面最大累計沉降實測值與計算值對比

NSFEM與實測值的相對誤差平均約為5%，FEM與實測值的相對誤差平均約為20%，隨機有限元模擬值尤其是壩體中上部模擬值明顯更加接近于實測值。

5 筑壩材料變異系數對壩體沉降的影響規律

在上述研究的基礎上，分析筑壩材料變異系數對壩體最大沉降值的影響，由于最大沉降值往往出現在3BII區和3C區，本文在控制其他分區其他材料參數以及參數統計特征值不變的前提下，分別對3BII區和3C區參數K、φ0、Kb的變異系數對壩體沉降的影響進行了研究。

圖8反映了3BII區和3C區參數K的變異系數，即C.V.(K)對壩體沉降的影響規律。由圖8可看出：(1)不同變異系數在壩體下部的沉降值基本重合，表明材料參數的變異性對壩體下部的沉降值影響較小。這是由于壩體下部受到巖土地基的約束，因此沉降值變化不大；(2)隨著C.V.(K)的減小，壩體最大沉降值逐漸減小。這是由于C.V.(K)減小，材料相對均勻，不均勻沉降減少；(3)隨著C.V.(K)的減小，折線之間的間距越來越小，表明壩體沉降減小的速率逐漸減緩，達到最小沉降值之后不再變化。計算表明，以高程1 952 m處的沉降值為例，C.V.(K)每減小0.1，3BII區沉降值的減小速率由5.09%減小到1.7%，3C區沉降值的減小速率由0.33%減小到0.13%；(4)圖8(a)較圖8(b)折線之間的間距更大。隨著C.V.(K)的減小，3BII區最大變化率為5.09%，3C區最大變化率為0.33%，3BII區明顯大于比3C區，表明壩體沉降對3BII區的敏感度大于3C區；(5)圖8(a)較圖8(b)最大沉降值的最小值更小，如C.V.(K)=0.05時，3BII區最大沉降值為55.91 cm，而3C區最大沉降值為68.18 cm，表明實際工程中通過提高3BII區材料參數的控制標準可以更有效控制壩體沉降值。

3BII區和3C區C.V.(φ0)對壩體沉降的影響規律與圖8相似，不同之處在于：在C.V.(φ0)與沉降的關系圖中，曲線間距更大，表明壩體沉降值對C.V.(φ0)的敏感性更大。C.V.(Kb)對壩體沉降影響規律的不同在于：(1)隨C.V.(Kb)變化，壩體沉降變化不明顯，尤其是3C區出現多條折線重合的現象，表明壩體沉降對C.V.(Kb)敏感性很?。?2)C.V.(Kb)很大時對應的壩體最大沉降值仍較小，進一步證明壩體沉降對C.V.(Kb)敏感性很小。

圖9為3BII區和3C區材料參數變異系數對最大沉降值影響規律圖。

由圖9(a)可知：(1)隨著材料參數的變異系數減小，壩體沉降值隨之減小且曲線斜率也隨之減小，減小到一定值時，壩體最大沉降值基本保持不變，進一步驗證了圖8的結論；(2)當C.V.(φ0)>0.1時，隨著C.V.(φ0)的增大，最大沉降值明顯增大，當C.V.(φ0)>0.3時，隨著C.V.(φ0)的增大，最大沉降值急劇增大；(3)當C.V.(K)>0.1時，隨著K的增大，最大沉降值明顯增大，當C.V.(K)>0.7時，隨著C.V.(K)的增大，最大沉降值快速增大；(4)參數Kb的變異系數的變化對最大沉降值基本沒有影響。

由圖9(a)與9(b)對比可知：圖9(b)相對圖9(a)3條曲線的斜率均較緩，圖9(a)中K、φ0、Kb的變異系數曲線最小值分別為66.27、73.19、72.05 cm，圖9(b)中曲線的最小值分別為70.37、74.48、70.42 cm。圖9(a)中的值較小，說明沉降值對3BII區材料參數特征值更為敏感，控制3BII區材料參數變異性更有利于控制壩體最大沉降值。

圖8 參數K的變異系數對不同高度壩體沉降的影響規律

圖9 材料參數變異性與最大沉降值的關系

6 結 論

為了考慮筑壩材料的不確定性因素對壩體沉降計算結果的影響，本文引入隨機場理論因素，取得了良好的效果。本文提出的方法適用于形狀不規則的壩體填筑料分區區域，對堆石壩設計和施工中填筑材料的控制指標的制定具有重要的參考意義。將研究內容與結論總結如下：

(1)在基于Cholesky分解的隨機場模擬方法的基礎上，結合堆石壩的材料參數特性，進行了基于有限元軟件MSC.Marc的二次開發，提出了堆石壩的非侵入式隨機有限元方法；

(2)通過工程實例驗證，表明材料參數空間變異性對堆石壩變形計算結果影響明顯，最大沉降值明顯大于常規有限元計算值，更接近實測值，更符合工程實際情況；

(3)探討了表征巖土體參數的變異系數對堆石壩沉降值的影響規律：壩體沉降值對3BII區材料參數的變異性敏感程度較3C區更大；材料參數按照敏感性從大到小排序為φ0、K、Kb。

(4)得出了最大沉降值與變異系數的曲線圖，可確定最大允許沉降對應的變異系數，進而確定各個材料分區材料參數概率分布的均值和標準差，對工程實際具有一定的參考意義。