考慮堆石料空間變異性的心墻壩地震安全性隨機有限元分析

楊 鴿，朱 晟，D. V. Griffiths

(1.杭州國家水電站大壩安全和應急工程技術中心有限公司，浙江 杭州 311122； 2.中國電建華東勘測設計研究院有限公司，浙江 杭州 311122；3.河海大學水工結構研究所，江蘇 南京 210098；4.科羅拉多礦業大學土木與環境工程學院，科羅拉多 戈爾登 81401)

中圖分類號：TV641.41;TU435文獻標識碼：A文章編號：0559-9342(2022)01-0077- 08

0 引 言

土石壩由天然巖土材料堆筑而成，雖然國內外現行的施工填筑規范中都對堆石料的干密度以及表征其密實程度的指標給出了較明確的控制標準，但對堆石料的級配和母巖性質等因素的控制則相當“粗放”[1-3]。不少研究都證明[4- 6]，密實程度相同的同一種堆石料，其級配不同時力學性質也會有較大差別；而母巖特性對堆石料力學性質的影響更是無須贅言。因此，雖然土石壩施工過程中對筑壩堆石料進行了質量控制，但壩料的物理力學性質仍很可能具有一定的不確定性。

近年來，已有很多學者致力于探討巖土材料的不確定性對結構反應影響的研究，例如Vanmarcke[7]、Baecher和Christian[8,9]、Griffiths和Fenton[10-16]、龔曉南和黃廣龍等[17-19]、Duncan和Low[20, 21]、Popescu等[22-23]、傅旭東等[24-26]以及李典慶、蔣水華、祁小輝等[27-33]的一系列研究。然而，在土石壩抗震工程領域，只有吳再光等[34]、Sanchez Lizarraga[35]以及王篤波等[36]曾探討筑壩材料的不確定性對大壩動力反應的影響。筆者曾探索采用隨機有限元法(Random Finite Element Method， RFEM)對土石壩進行隨機地震響應進行分析，并以某高混凝土面板堆石壩為例，探討了堆石料的級配、Duncan-ChangE-B模型模量系數以及初始摩擦角等物理力學性質的空間變異性對大壩地震加速度響應及永久變形的影響[37]。

本文將以某高心墻堆石壩為例，進一步探討筑壩材料的不確定性對心墻壩動力響應及抗震安全的影響。研究思路為在蒙特卡羅法的框架下，采用RFEM模擬考慮筑壩堆石料空間變異性時心墻壩的加速度響應及永久變形，通過對比隨機有限元法和常規確定性有限元法的計算結果，討論筑壩堆石料的空間變異性對心墻壩動力響應的影響，然后再對各隨機物理性質參數的敏感性進行比較分析。

1 心墻壩地震安全性隨機有限元分析方法

采用隨機有限元法分析心墻壩加速度響應及永久變形的基本流程如圖1所示。首先，采用局部平均細分法生成一組具有空間變異性的材料隨機場；然后，采用“逐單元法”將材料隨機場映射至有限元網格，使每個網格單元都找到相應的材料參數取值；之后，調用有限元靜力分析模塊，模擬壩體材料為該組隨機場時壩體的應力、變形等；最后，輸入地震動荷載，調用動力有限元分析模塊，采用Wilson-θ法逐步求解動力平衡方程，計算每個小時步內的動力加速度響應以及動應力。在此基礎上，采用沈珠江經驗公式計算殘余體應變和剪應變，利用初應變法計算殘余變形。通過上述過程即完成了一次隨機有限元分析，得到了一個隨機樣本，在蒙特卡羅法的框架下，重復上述步驟多次，直至計算得到的隨機樣本達到統計穩定后，即可通過對隨機樣本的統計分析，估計大壩地震響應的統計特性。

圖1 隨機有限元分析心墻壩隨機地震反應的基本流程

局部平均細分法生成材料隨機場的原理見文獻[37]以及Fenton[38]和Vanmarcke[39]的研究。土石壩非線性靜動力有限元分析的原理見顧淦臣等[40]的介紹。計算結果的統計穩定采用主要性能指標的滑動平均值(Running Average，Run.ava.)和滑動標準差(Running Standard Deviation，Run.sd.)判斷，其具體定義見文獻[37]。

2 案例概況

R工程初步設計擋水建筑物為礫石土直心墻堆石壩，最大壩高315 m，上游壩坡為1∶2.1，下游壩坡坡比為1∶2.0。心墻順河向最大寬度為162 m，上、下游坡比均為0.23；心墻與上、下游壩殼堆石之間設反濾層、過渡層。

2.1 計算模型

采用如圖2所示的有限元計算模型，共有4 595個實體單元，另在心墻與反濾I區之間設置124個接觸面單元。靜力分析模擬壩體的填筑和蓄水過程，壩體自重荷載的施加考慮施工填筑的順序分34級施加；水荷載則分8級施加至正常蓄水位。動力荷載方面，根據工程場地地震安全性評價報告，場地100年超越概率2%的基巖地震加速度峰值為0.434g，典型場地地震動時程如圖3a所示，對應場地加速度反應譜如圖3b所示。

圖2 有限元計算網格

圖3 地震動輸入(100年超越概率為2%)

2.2 堆石料的物理力學性質及其不確定性

筑壩堆石料的物理力學性質都具有一定程度的不確定性以及空間變異性。為了解堆石料力學性質(特別是變形特性)的不確定性，需要對大量的力學試驗結果進行統計分析，成本極高。受此影響，目前對巖土材料不確定性的研究多集中在其物理性質上，只有少數研究涉及到了材料的強度特性和變形特性；而在考察材料不確定性對結構反應的影響的研究中，也一般只考慮材料物理性質或強度特性的不確定性，極少研究涉及到材料變形特性不確定性[38,41,42]。特別值得注意的是，在幾乎所有考察巖土材料力學特性不確定性對結構反應影響的研究中，描述力學特性不確定性的統計參數都是作者根據有限的已發表的文獻中的數據“擬定”而來，并非根據對具體研究對象的試驗測量結果獲得。

相比于力學性質參數，巖土材料物理性質參數的統計特性則通常較容易獲得。為此，筆者提出[37]采用包含材料物理指標的力學模型描述筑壩巖土材料的力學性質，這樣便可通過材料物理參數的統計特征推算出材料力學特性的統計特征，從而大大降低獲取材料力學參數統計特性的成本，使巖土結構概率性分析中材料參數統計特性的獲方法更加可靠。

2.2.1 動力本構模型

堆石料動力性質方面，將采用Menq[5]提出的包含堆石料孔隙比e、不均勻系數Cu、平均粒徑D50的動本構模型，并將上述物理參數設置為具有空間變異性的隨機場，即

(1)

得到Gmax后，可通過下式計算動剪模量G，即

(2)

式中，γ為動剪應變；γr為參考應變，γr=K1×Cu-nr×(σ0/Pa)0.5，其中，nr、K1均為擬合參數；a為指數，a=K2+0.1×lg(σ0/Pa)，其中，K2為參數。

材料的阻尼特性方面，由于Menq給出的計算公式形式過于復雜，本文采用下式計算阻尼比：

(3)

式中，λr為參考阻尼比，λr=K3(σ0/Pa)+K4，其中，K3、K4均為擬合參數。

上述模型中，參數KG、nr、K1、K2、K3、K4需通過對材料動力試驗曲線擬合得到。結合室內試驗土料的物理指標，對R工程堆石料動力三軸試驗結果進行擬合，得到Menq模型的擬合模型參數取值見表1。

表1 動力本構模型參數

2.2.2 靜力本構模型

材料的靜力本構模型方面，由于目前暫時尚未有描述堆石料物理指標與靜應力-應變關系的模型，本文將繼續沿用在心墻壩有限元分析中常用的Duncan-ChangE-B模型，并參考以往經驗將對結構反應影響較大的參數(模量系數K和初始摩擦角φ0)設置為隨機場。對R工程堆石料靜三軸試驗結果進行擬合，得到Duncan-ChangE-B模型參數見表2。

表2 靜力本構模型(Duncan-Chang E-B模型)參數

此外，土體干密度決定了靜力分析階段的自重荷載以及動力分析階段的慣性力，因而也被設置為隨機場。

2.2.3 隨機物理力學參數的統計特性

隨機有限元分析中將考慮堆石料孔隙比e、不均勻系數Cu、平均粒徑D50、干密度、E-B模型的模量系數以及初始摩擦角的空間變異性。由于示例工程尚未建成，并無法通過對壩體材料的檢驗確定上述隨機變量的統計特性，因而本文將參考現有填筑標準、相似工程現場檢測資料以及相關研究成果綜合估計上述隨機參數的統計特性，統計參數的具體擬定過程見文獻[37]，擬定結果見表3。

表3 堆石料物理參數

隨機場內不同位置對應的變量間的相關性強弱用相關距離描述。目前僅有少量針對天然土層的相關距離的研究成果[43-47]，而并沒有針對心墻壩工程的研究成果。為此，本研究參考天然土層的相關距離取值，將各個參數隨機場的水平和垂直相關距離分別設定為150 m和50 m。圖4為隨機參數為干密度的典型參數隨機場。

圖4 不確定性筑壩堆石料參數分布示意(隨機參數為干密度)

3 心墻壩地震安全性分析結果

3.1 計算結果的統計穩定性

在根據RFEM計算結果估計總體統計規律之前，先要保證計算結果達到“統計穩定”。為此，選擇壩體最大震后沉降以及壩頂水平峰值加速度作為性能指標，繪出他們的滑動平均值和滑動標準差隨RFEM模擬次數的變化，如圖5所示，壩體最大震后沉降在模擬次數達到300時已經非常穩定，壩頂水平峰值加速度在模擬次數達到400以后也已非常穩定。因而，可以認為本文中對600次RFEM模擬結果的分析能夠代表總體的統計特性。

圖5 隨機有限元計算結果的穩定性

3.2 加速度響應

圖6統計了考慮材料隨機性時壩頂典型節點的水平向峰值加速度的頻率分布情況，其中實線為600次RFEM計算結果的均值，虛線為不考慮材料隨機性時的計算結果。表4為水平向加速度峰值計算結果的統計特征值。由表4可知，考慮材料隨機性時，壩頂峰值加速度計算結果的均值(0.998 22g)略大于確定性計算結果(0.971 91g)，對應放大倍數為2.30。計算結果的離散性方面，水平向壩頂峰值加速度的極差高達約0.6g，但變異系數僅為7.44%，而90%、95%以及98%保證率對應的峰值加速度取值分別為0.109 431g、0.114 279g以及0.118 187g，分別為確定性分析結果的1.126倍、1.176倍和1.216倍。

圖6 壩頂水平向峰值加速度頻率分布

表4 壩頂水平向峰值加速度統計參數

筑壩材料具有不確定性時，除了加速度響應的極值，加速度響應的頻譜特性也具有不確定性。圖7所示為壩頂典型節點加速度時程的反應譜，其中曲線為確定性有限元計算的結果，其余為4次典型RFEM的計算結果。可以看到，雖然大多數響應加速度反應譜的形狀與確定性計算結果相似，但不同方案反應譜的高峰區所對應的頻率范圍或者峰值仍有所不同。例如，順河向加速度方面，多數方案在0.8 Hz附近的反應譜取值已經較小，但方案3、4在0.8 Hz附近卻出現了一個尖峰；方案2在1.5 Hz附近出現峰值，方案3在0.8 Hz附近出現峰值，而其他方案的峰值則都在0.2～0.3 Hz附近。類似的情況在豎向加速度反應譜圖中也可觀察到。加速度響應反應譜高峰區對應的頻率反映了壩體的自振特性，因而上述高峰區位置不同的現象反映了筑壩材料的不確定性引起了壩體自振特性的變化。

圖7 壩頂加速度反應譜(典型案例)

3.3 永久變形

圖8統計了考慮材料隨機性時壩體最大震后沉降頻率分布情況，其中實線為600次RFEM計算結果的均值，虛線為不考慮材料隨機性時的計算結果，為327.3 cm。表5為最大沉降計算結果的統計特征值。由表5可知，考慮材料隨機性時最大沉降計算結果的均值明顯大于確定性計算結果；而確定性有限元計算結果的保證率約為25%，即在材料存在隨機性的情況下，最大沉降大于確定性有限元預測值的概率約有75%。計算結果的離散性方面，最大沉降的變異系數為17.2%。若要設計值的保證率達到90%、95%和98%，最大沉降的取值分別為455.2、490.7、527.7 cm，分別為確定性分析結果的1.39倍、1.50倍和1.61倍。

圖8 震后壩體最大沉降頻率分布

表5 震后壩體最大沉降統計特性

對RFEM計算得到的最大震后永久變形進行Kolmogorov-Smirnov檢驗以確定其是否滿足正態分布，經檢驗，順河向和豎直向最大永久變形的雙側漸近顯著性都小于0.05%。因此，可以認為隨機最大永久變形不滿足正態分布。繪制最大永久變形與對數正態分布的P-P圖以檢驗其是否滿足對數正態分布，如圖9所示，計算點與45°直線重合極好，由此可知隨機最大永久變形滿足對數正態分布。

圖9 震后壩體最大沉降P-P示意

4 材料物理參數隨機性的敏感性分析

對材料物理參數(干密度、孔隙率、不均勻系數及平均粒徑)的不確定性做敏感性分析，為堆石料填筑質量控制標準的設定提供一定參考。為此進行4組蒙特卡羅隨機有限元模擬，每組分別考慮一種物理參數的隨機性，各參數的變異系數取值見表3。此外，用敏感率R表征地震響應對各個物理參數隨機性的敏感性，即

(4)

敏感率R的取值越大，表明壩體響應對該材料參數的敏感性越高；反之亦然。

表6為考慮各材料參數隨機性時最大豎向永久變形的頻率分布情況。孔隙率和干密度本身的變異系數為5%，由表6可知，孔隙率和干密度引起最大沉降的變異系數分別為6.56%和2.89%，為其本身變異系數的1.3倍和0.577；因而，在目前大壩填筑中嚴格控制孔隙率和干密度是非常合理的。相比之下，不均勻系數和平均粒徑的變異系數為40%時引起的最大豎向永久變形的變異系數分別為11.07%和6.78%；因此，現行土石壩填筑規范并不要求對這兩項參數進行非常嚴格控制也是合理的。但是不均勻系數和平均粒徑的隨機性過大也還是會對結構響應造成影響，由表6可知，考慮不均勻系數隨機性時結構響應的變異系數最大，其次是考慮平均粒徑隨機性時，隨后才是考慮孔隙率隨機性時。因此，不能完全不對不均勻系數和平均粒徑的隨機性進行控制。

表6 考慮各個材料性質不確定性的計算結果統計(永久變形)

表7為考慮各材料參數隨機性時壩頂水平峰值加速度響應的分布情況。由表6、7可知，壩頂水平峰值加速度對應各個材料參數的敏感率R普遍略低于最大豎向永久變形對應的各R值，即壩頂水平峰值加速度對材料參數隨機性的敏感度略低。此外，與最大永久變形類似，峰值加速度響應對孔隙率的敏感性相對較高，其次是干密度，而對不均勻系數和平均粒徑的敏感性相對較低。但是，并不代表就不需要對筑壩材料的不均勻系數和平均粒徑進行控制——考慮不均勻系數的隨機性時響應的變異系數仍然最大，而考慮平均粒徑隨機性時對應的變異系數也大于考慮干密度隨機性時。

表7 考慮各個材料性質不確定性的計算結果統計(水平峰值加速度)

5 結 論

本文以某高心墻堆石壩為例，采用基于局部平均細分法的隨機有限元法模擬考慮筑壩堆石料空間變異性時心墻壩的地震響應以及永久變形，據此探討了堆石料的空間變異性對心墻壩地震安全性的影響，以及心墻壩地震響應對各隨機物理參數的敏感性。主要結論為

(1)考慮筑壩材料空間變異性時，大壩峰值加速度響應的取值和頻譜特性都發生一定程度的離散，為使設計方案的保證率達到90%以上，需將確定性分析得到的峰值加速度放大1.1倍～1.2倍。

(2)大壩震后最大沉降的離散程度相比更大，且為使設計方案的保證率達到90%，需將確定性分析結果放大1.4倍～1.6倍。

(3)心墻壩地震響應對堆石料孔隙率和干密度隨機性的敏感性相對較高，對不均勻系數和平均粒徑的隨機性的敏感性較低。但是，不均勻系數和平均粒徑的隨機性仍然會對大壩的地震反應帶來變異性，因此有必要對其隨機性進行控制。

綜上所述，忽略材料的不確定性可能導致低估大壩的地震反應。目前大壩填筑中嚴格控制孔隙率和干密度是非常合理的，而施工過程中也宜加強對堆石料級配的控制，避免其出現過大的隨機性。