楊 鴿，朱 晟，D. V. Griffiths

(1.杭州國(guó)家水電站大壩安全和應(yīng)急工程技術(shù)中心有限公司，浙江 杭州 311122； 2.中國(guó)電建華東勘測(cè)設(shè)計(jì)研究院有限公司，浙江 杭州 311122；3.河海大學(xué)水工結(jié)構(gòu)研究所，江蘇 南京 210098；4.科羅拉多礦業(yè)大學(xué)土木與環(huán)境工程學(xué)院，科羅拉多 戈?duì)柕?81401)

中圖分類號(hào)：TV641.41;TU435文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：0559-9342(2022)01-0077- 08

0 引 言

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

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

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

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

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

圖1 隨機(jī)有限元分析心墻壩隨機(jī)地震反應(yīng)的基本流程

局部平均細(xì)分法生成材料隨機(jī)場(chǎng)的原理見(jiàn)文獻(xiàn)[37]以及Fenton[38]和Vanmarcke[39]的研究。土石壩非線性靜動(dòng)力有限元分析的原理見(jiàn)顧淦臣等[40]的介紹。計(jì)算結(jié)果的統(tǒng)計(jì)穩(wěn)定采用主要性能指標(biāo)的滑動(dòng)平均值(Running Average，Run.ava.)和滑動(dòng)標(biāo)準(zhǔn)差(Running Standard Deviation，Run.sd.)判斷，其具體定義見(jiàn)文獻(xiàn)[37]。

2 案例概況

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

2.1 計(jì)算模型

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

圖2 有限元計(jì)算網(wǎng)格

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

2.2 堆石料的物理力學(xué)性質(zhì)及其不確定性

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

相比于力學(xué)性質(zhì)參數(shù)，巖土材料物理性質(zhì)參數(shù)的統(tǒng)計(jì)特性則通常較容易獲得。為此，筆者提出[37]采用包含材料物理指標(biāo)的力學(xué)模型描述筑壩巖土材料的力學(xué)性質(zhì)，這樣便可通過(guò)材料物理參數(shù)的統(tǒng)計(jì)特征推算出材料力學(xué)特性的統(tǒng)計(jì)特征，從而大大降低獲取材料力學(xué)參數(shù)統(tǒng)計(jì)特性的成本，使巖土結(jié)構(gòu)概率性分析中材料參數(shù)統(tǒng)計(jì)特性的獲方法更加可靠。

2.2.1 動(dòng)力本構(gòu)模型

堆石料動(dòng)力性質(zhì)方面，將采用Menq[5]提出的包含堆石料孔隙比e、不均勻系數(shù)Cu、平均粒徑D50的動(dòng)本構(gòu)模型，并將上述物理參數(shù)設(shè)置為具有空間變異性的隨機(jī)場(chǎng)，即

(1)

得到Gmax后，可通過(guò)下式計(jì)算動(dòng)剪模量G，即

(2)

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

材料的阻尼特性方面，由于Menq給出的計(jì)算公式形式過(guò)于復(fù)雜，本文采用下式計(jì)算阻尼比：

(3)

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

上述模型中，參數(shù)KG、nr、K1、K2、K3、K4需通過(guò)對(duì)材料動(dòng)力試驗(yàn)曲線擬合得到。結(jié)合室內(nèi)試驗(yàn)土料的物理指標(biāo)，對(duì)R工程堆石料動(dòng)力三軸試驗(yàn)結(jié)果進(jìn)行擬合，得到Menq模型的擬合模型參數(shù)取值見(jiàn)表1。

表1 動(dòng)力本構(gòu)模型參數(shù)

2.2.2 靜力本構(gòu)模型

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

表2 靜力本構(gòu)模型(Duncan-Chang E-B模型)參數(shù)

此外，土體干密度決定了靜力分析階段的自重荷載以及動(dòng)力分析階段的慣性力，因而也被設(shè)置為隨機(jī)場(chǎng)。

2.2.3 隨機(jī)物理力學(xué)參數(shù)的統(tǒng)計(jì)特性

隨機(jī)有限元分析中將考慮堆石料孔隙比e、不均勻系數(shù)Cu、平均粒徑D50、干密度、E-B模型的模量系數(shù)以及初始摩擦角的空間變異性。由于示例工程尚未建成，并無(wú)法通過(guò)對(duì)壩體材料的檢驗(yàn)確定上述隨機(jī)變量的統(tǒng)計(jì)特性，因而本文將參考現(xiàn)有填筑標(biāo)準(zhǔn)、相似工程現(xiàn)場(chǎng)檢測(cè)資料以及相關(guān)研究成果綜合估計(jì)上述隨機(jī)參數(shù)的統(tǒng)計(jì)特性，統(tǒng)計(jì)參數(shù)的具體擬定過(guò)程見(jiàn)文獻(xiàn)[37]，擬定結(jié)果見(jiàn)表3。

表3 堆石料物理參數(shù)

隨機(jī)場(chǎng)內(nèi)不同位置對(duì)應(yīng)的變量間的相關(guān)性強(qiáng)弱用相關(guān)距離描述。目前僅有少量針對(duì)天然土層的相關(guān)距離的研究成果[43-47]，而并沒(méi)有針對(duì)心墻壩工程的研究成果。為此，本研究參考天然土層的相關(guān)距離取值，將各個(gè)參數(shù)隨機(jī)場(chǎng)的水平和垂直相關(guān)距離分別設(shè)定為150 m和50 m。圖4為隨機(jī)參數(shù)為干密度的典型參數(shù)隨機(jī)場(chǎng)。

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

3 心墻壩地震安全性分析結(jié)果

3.1 計(jì)算結(jié)果的統(tǒng)計(jì)穩(wěn)定性

在根據(jù)RFEM計(jì)算結(jié)果估計(jì)總體統(tǒng)計(jì)規(guī)律之前，先要保證計(jì)算結(jié)果達(dá)到“統(tǒng)計(jì)穩(wěn)定”。為此，選擇壩體最大震后沉降以及壩頂水平峰值加速度作為性能指標(biāo)，繪出他們的滑動(dòng)平均值和滑動(dòng)標(biāo)準(zhǔn)差隨RFEM模擬次數(shù)的變化，如圖5所示，壩體最大震后沉降在模擬次數(shù)達(dá)到300時(shí)已經(jīng)非常穩(wěn)定，壩頂水平峰值加速度在模擬次數(shù)達(dá)到400以后也已非常穩(wěn)定。因而，可以認(rèn)為本文中對(duì)600次RFEM模擬結(jié)果的分析能夠代表總體的統(tǒng)計(jì)特性。

圖5 隨機(jī)有限元計(jì)算結(jié)果的穩(wěn)定性

3.2 加速度響應(yīng)

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

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

表4 壩頂水平向峰值加速度統(tǒng)計(jì)參數(shù)

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

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

3.3 永久變形

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

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

表5 震后壩體最大沉降統(tǒng)計(jì)特性

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

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

4 材料物理參數(shù)隨機(jī)性的敏感性分析

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

(4)

敏感率R的取值越大，表明壩體響應(yīng)對(duì)該材料參數(shù)的敏感性越高；反之亦然。

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

表6 考慮各個(gè)材料性質(zhì)不確定性的計(jì)算結(jié)果統(tǒng)計(jì)(永久變形)

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

表7 考慮各個(gè)材料性質(zhì)不確定性的計(jì)算結(jié)果統(tǒng)計(jì)(水平峰值加速度)

5 結(jié) 論

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

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

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

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

綜上所述，忽略材料的不確定性可能導(dǎo)致低估大壩的地震反應(yīng)。目前大壩填筑中嚴(yán)格控制孔隙率和干密度是非常合理的，而施工過(guò)程中也宜加強(qiáng)對(duì)堆石料級(jí)配的控制，避免其出現(xiàn)過(guò)大的隨機(jī)性。