周方明，潘紅梅，占良紅王 鋒

(1.宿州學(xué)院資源與土木工程學(xué)院，安徽 宿州 234000；2.南昌大學(xué)建筑工程學(xué)院，江西 南昌 330031；3.江西水利職業(yè)學(xué)院，江西 南昌 330013)

0 引 言

我國現(xiàn)役混凝土壩無論是數(shù)量還是壩高均居世界首位，這些重力壩工程健康服役直接關(guān)系到整個水工程的安危，已成為關(guān)乎國計民生的重大公共安全問題[1]。因此，對重力壩的可靠性評估研究是目前大壩安全控制領(lǐng)域高度關(guān)注的重要科學(xué)問題。

重力壩特殊的結(jié)構(gòu)性質(zhì)及其復(fù)雜的工作環(huán)境，決定了對其進(jìn)行可靠性分析的復(fù)雜性。自1926年德國學(xué)者邁耶最早將概率論引入到分析結(jié)構(gòu)可靠度以來，結(jié)構(gòu)可靠度理論的運(yùn)用取得了重大進(jìn)展[2]。我國首次把結(jié)構(gòu)可靠性分析方法引入大壩可靠度計算的是河海大學(xué)“水工可靠度”科研組。武清璽等[3]研究提出了基于有限元法的重力壩可靠性分析，并進(jìn)一步計算了重力壩抗滑穩(wěn)定的可靠性；吳世偉等[4]基于概率可靠度理論對重力壩最大可能破壞模式進(jìn)行了探討，提出了地基、壩體、地基與壩體交界面3種主要的破壞模式。將概率可靠性理論用于評估重力壩的安全狀態(tài)時，需突出解決好兩大問題：一是在重力壩可靠性分析過程中其參數(shù)通常選用設(shè)計標(biāo)準(zhǔn)值，未考慮重力壩服役期間在內(nèi)外環(huán)境及自身老化持續(xù)影響下，其參數(shù)因發(fā)生演化而改變；二是對于重力壩這種極為復(fù)雜工程，其極限狀態(tài)方程常表現(xiàn)為高度非線性且難以用顯式的方程給以刻畫，采用材料力學(xué)分析的手段構(gòu)建極限狀態(tài)方程與重力壩真實(shí)響應(yīng)存在較大誤差，這些因素勢必影響可靠性真實(shí)計算結(jié)果。

工程中常用響應(yīng)面法(RSM)以二項(xiàng)式模型來近似擬合其極限狀態(tài)方程，江勝華等[5]運(yùn)用加權(quán)響應(yīng)面法，計算了重力壩各斷層抗滑穩(wěn)定的可靠指標(biāo)，并與安全系數(shù)進(jìn)行了比較，表明其安全系數(shù)與可靠指標(biāo)并非存在完全的對等關(guān)系。然已有研究表明，在高非線性問題中，RSM法需大量的計算工作，且對于高非線性可靠性問題常出現(xiàn)無法收斂和精度不足的情況[6]。Gaspar等[7]指出，Kriging模型較其它方法具有更優(yōu)的預(yù)測能力及靈活性，與多項(xiàng)式回歸模型相比，Kriging模型求解的失效概率更為準(zhǔn)確。文獻(xiàn)[8]基于Kriging模型精確擬合了某航空結(jié)構(gòu)的極限狀態(tài)方程，并對其可靠性實(shí)現(xiàn)了高效計算。

據(jù)此，考慮到重力壩參數(shù)隨服役而發(fā)生演化改變，本文基于重力壩參數(shù)的衰減函數(shù)獲取重力壩服役后其真實(shí)概率分布信息。在此基礎(chǔ)上，借鑒Kriging模型強(qiáng)預(yù)測及擬合能力，探研重力壩極限狀態(tài)方程精確擬合方法，并根據(jù)可靠指標(biāo)幾何定義，擬建立重力壩壩體單元強(qiáng)度破壞和沿壩基面抗滑失穩(wěn)的概率可靠度優(yōu)化模型，以實(shí)現(xiàn)深入分析重力壩可靠性。

1 重力壩Kriging模型的建立及檢驗(yàn)

1.1 Kriging模型的基本原理

Kriging模型是一種統(tǒng)計近似模型[9]，對于任意輸入變量x∈Rn與響應(yīng)值y有如下關(guān)系，即：

(1)

式中：f(x)T=[f1(x)，f2(x)，…，fp(x)]為回歸模型，常采用二次多項(xiàng)式函數(shù)形式；β=[β1，β2，…，βp]T為回歸模型系數(shù)；z(x)為隨機(jī)誤差，其均值為0、方差為σ2。

假定隨機(jī)抽樣m個計算樣本點(diǎn)組成的矩陣為S=[s1，s2，…，sm]T，si∈Rn(i=1，2，…，m)，其每個樣本點(diǎn)組成的響應(yīng)值矩陣為Y=[y1，y2，…，ym]T。通過精確模擬樣本點(diǎn)矩陣S和響應(yīng)值矩陣Y之間對應(yīng)規(guī)律，可借助MATLAB強(qiáng)大計算優(yōu)勢完成Kriging模型的建立。

1.2 重力壩Kriging模型

以2.1節(jié)所述Kriging模型理論為基礎(chǔ)，通過合理選取重力壩計算參數(shù)和功能函數(shù)(響應(yīng)函數(shù))形式，構(gòu)建重力壩的Kriging模型。本文重點(diǎn)考慮重力壩單元強(qiáng)度破壞和沿壩基面抗滑失穩(wěn)兩種破壞形式。

本文參考文獻(xiàn)[10]，根據(jù)單元的第一、第二和第三主應(yīng)力(σ1、σ2和σ3)建立單元破壞的響應(yīng)值函數(shù)，即：

(2)

式中：g(x)為重力壩單元強(qiáng)度的功能函數(shù)(響應(yīng)值)；α=ft/fc；ft和fc分別為重力壩單元的抗拉強(qiáng)度和抗壓強(qiáng)度。

根據(jù)壩體與壩基抗滑面上所有單元應(yīng)力，建立失穩(wěn)破壞的響應(yīng)值函數(shù)，即：

(3)

式中：n為抗滑面上單元總數(shù)；f和c為分別為材料的摩擦系數(shù)和粘聚力；σyi和τxyi分別為單元i的正應(yīng)力和剪應(yīng)力；bi為單元i沿滑動面的邊長。

1.3 Kriging模型的檢驗(yàn)

Kriging模型近似擬合重力壩極限狀態(tài)方程，必然與真實(shí)響應(yīng)面存在一定誤差[11]。本文用復(fù)相關(guān)系數(shù)來檢驗(yàn)Kriging模型擬合效果，即：

(4)

(5)

由定義可知R2介于0到1之間，其中R2越趨于1表明擬合效果越佳。其檢驗(yàn)步驟如下：

Step1：根據(jù)重力壩參數(shù)統(tǒng)計分布信息，抽樣獲取足夠樣本點(diǎn)，得到總樣本點(diǎn)集U，再從U中隨機(jī)抽取適量樣本點(diǎn)組成子樣本集C1，并根據(jù)子樣本集C1信息結(jié)合有限元仿真響應(yīng)建立Kriging模型。

Step3：根據(jù)復(fù)相關(guān)系數(shù)R2判斷所建Kriging模型是否達(dá)到預(yù)期精度要求，符合則選用該模型，否則返回Step1重新選擇子樣本集C1重復(fù)以上步驟，直至滿足精度要求。

2 重力壩可靠指標(biāo)β計算模型

2.1可靠指標(biāo)β優(yōu)化模型

根據(jù)經(jīng)典的可靠度理論，定義可靠指標(biāo)β是在標(biāo)準(zhǔn)正態(tài)空間中，原點(diǎn)O距極限狀態(tài)曲面最短距離，驗(yàn)算點(diǎn)P*為極限狀態(tài)曲面一點(diǎn)，如圖1所示。對于任一觀測點(diǎn)x都唯一對應(yīng)一響應(yīng)值y=g(x)，若在保證y*=0(極限狀態(tài)曲面)前提下找到某一觀測點(diǎn)x*在標(biāo)準(zhǔn)正態(tài)空間中距原點(diǎn)最近，則觀測點(diǎn)x*即為驗(yàn)算點(diǎn)P*。則可靠指標(biāo)β的計算可轉(zhuǎn)化為一數(shù)值優(yōu)化問題，即：

(6)

式中：x*為標(biāo)準(zhǔn)正態(tài)空間設(shè)計驗(yàn)算點(diǎn)；y=g(x)=0為Kriging模型擬合得到極限狀態(tài)方程；ε為允許誤差，取0.000 1?？紤]到Kriging模型只是近似模型，因此人為設(shè)定一個極小的允許誤差ε=0.000 1。對于式(6)可在MATLAB中調(diào)用fmincon函數(shù)優(yōu)化求解，具體操作步驟為：①根據(jù)計算參數(shù)建立目標(biāo)函數(shù)my_fun；②對擬合Kriging模型全局化(global dmodel，dmodel為擬合Kriging模型)，并根據(jù)允許誤差大小構(gòu)建約束函數(shù)my_con；③調(diào)用fmincon函數(shù)對式(6)求解，調(diào)用形式為[x，fval]=fmincon(@my_fun，x0，A，b，Aeq，beq，lb，ub，@my_con)，其中，x0以所有計算參數(shù)均值作為初始值；A和Aeq為不等式約束和等式約束的系數(shù)矩陣；b，beq，lb和ub為線性不等式約束的上、下界，在此問題中A，b，Aeq，beq，lb，ub=[]。

圖1 可靠指標(biāo)幾何意義示意圖

2.2 重力壩可靠指標(biāo)β的計算方法

本文特根據(jù)第2節(jié)所述重力壩Kriging模型建立及校驗(yàn)方法，結(jié)合重力壩實(shí)測資料信息，對不具備實(shí)測資料的參數(shù)按衰減函數(shù)確定其服役后統(tǒng)計信息。在此基礎(chǔ)上，結(jié)合1.2及1.3節(jié)所述方法，建立最佳的Kriging模型擬合重力壩的極限狀態(tài)方程，根據(jù)式(6)構(gòu)建重力壩可靠指標(biāo)優(yōu)化模型，對重力壩的可靠性進(jìn)行分析計算。

重力壩可靠指標(biāo)的具體計算步驟如下：

Step1：根據(jù)重力壩基本資料確定參數(shù)統(tǒng)計信息，抽樣獲取3 000個樣本點(diǎn)組成總樣本點(diǎn)集U；

Step2：隨機(jī)選取40個樣本點(diǎn)形成子樣本集C1，利用有限元ABAQUS建立重力壩模型進(jìn)行結(jié)構(gòu)分析，獲得各樣本點(diǎn)處的結(jié)構(gòu)響應(yīng)(應(yīng)力)；

Step3：依據(jù)應(yīng)力計算成果，結(jié)合式(2)和(3)分別計算重力壩壩體單元強(qiáng)度響應(yīng)值及抗滑穩(wěn)定響應(yīng)值；

Step4：根據(jù)樣本點(diǎn)的響應(yīng)值，利用MATLAB對其進(jìn)行擬合得到近似極限狀態(tài)方程y=g(x)，并根據(jù)1.3節(jié)所述檢驗(yàn)方法判斷其是否滿足要求，符合精度要求則采用此y=g(x)作為極限狀態(tài)方程，否則返回Step2；

Step5：基于擬合所得近似極限狀態(tài)方程y=g(x)，根據(jù)2.1節(jié)所述方法計算重力壩單元強(qiáng)度破壞和抗滑穩(wěn)定的可靠指標(biāo)。

3 工程應(yīng)用

3.1 工程資料及模型建立

某混凝土重力壩總壩長308.50 m，從左岸至右岸共6個壩段，本文選取4#非溢流壩段作為研究對象。根據(jù)工程地質(zhì)勘測資料，該混凝土重力壩所處基巖主要為微風(fēng)化狀態(tài)，未勘測出深層的裂縫分布。因此，本文重點(diǎn)基于本文方法分析4#壩段的單元強(qiáng)度破壞和沿建基面滑動失穩(wěn)的可靠性。該壩段全長50.00 m，壩頂高程179.00 m，最大壩高96.00 m，最小壩高56.00 m，正常蓄水位173.00 m。根據(jù)該壩段實(shí)際情況，本文利用ABAQUS進(jìn)行建立有限元模型，其中，基巖尺寸為壩踵和壩趾分別向外延伸1.5倍最大壩高，壩基深度取100.00 m。4#壩段三維有限元模型如圖2所示。其中，該模型總節(jié)點(diǎn)數(shù)27 679個、總單元數(shù)23 534個。

圖2 4#壩段網(wǎng)格劃分圖

3.2 壩體和壩基參數(shù)的確定

重力壩大多參數(shù)是建壩期間所實(shí)測的，重力壩服役多年后其結(jié)構(gòu)參數(shù)必將發(fā)生演化改變。據(jù)此，本文考慮材料參數(shù)的老化衰減的影響，對現(xiàn)階段實(shí)測數(shù)據(jù)不足的參數(shù)，根據(jù)這些參數(shù)的設(shè)計標(biāo)準(zhǔn)值和經(jīng)驗(yàn)衰減函數(shù)計算其現(xiàn)在確定值。壩體密度和揚(yáng)壓力衰減函數(shù)分別為φ1=e-0.000 5 t和φα=e0.005 t(t為大壩服役年數(shù)，下同)；混凝土和壩基的彈性模量、抗壓強(qiáng)度、抗拉強(qiáng)度、摩擦系數(shù)和粘聚力的衰減函數(shù)均為φ2=e-0.005 t。同時根據(jù)重力壩物理與力學(xué)參數(shù)變異性對其可靠性計算結(jié)果影響大小，影響較小的參數(shù)視為定值考慮，而影響較大的參數(shù)視為不確定變量，其定值參數(shù)統(tǒng)計見表1，不確定參數(shù)界限見表2。材料力學(xué)模型采用D-P準(zhǔn)則。計算施加的荷載組合為：自重+上游水壓力+揚(yáng)壓力。

3.3 重力壩可靠指標(biāo)計算結(jié)果

為直觀分析重力壩的可靠狀態(tài)，特根據(jù)2.2節(jié)計算步驟計算該壩段單元可靠指標(biāo)β。壩體單元指標(biāo)分布云圖如圖3所示，壩體抗滑面和巖基抗滑面單元指標(biāo)分布云圖分別如圖4(a)、(b)所示，以及抗滑面單元可靠指標(biāo)統(tǒng)計如圖5所示。根據(jù)GB50199—2013《水利水電結(jié)構(gòu)可靠度設(shè)計統(tǒng)一標(biāo)準(zhǔn)》相關(guān)規(guī)定，對于Ⅰ級水利樞紐可靠指標(biāo)需滿足4.2以上。由圖3可知，壩體大部分單元可靠指標(biāo)處于5以上，僅在壩踵附近出現(xiàn)少量的可靠指標(biāo)較小的單元，不足以威脅大壩正常運(yùn)行，其中最小可靠指標(biāo)為2.79，位于壩踵處7 113號單元。再由圖4可知，抗滑面也僅在靠近壩踵區(qū)域單元可靠指標(biāo)相對較小，靠近上游單元可靠指標(biāo)較下游單元大。結(jié)合圖5可知，壩體和巖基抗滑面共6個單元可靠指標(biāo)小于4，占總滑面單元(734個)0.8%，該壩段沿建基面抗滑穩(wěn)定可靠度亦較安全。該重力壩實(shí)際運(yùn)行良好，暫未發(fā)現(xiàn)影響大壩正常運(yùn)行的安全隱患，這也為本文方法計算結(jié)果提供了一定的保障。

表1 定值參數(shù)統(tǒng)計

注：壩體密度初始設(shè)計值2 450，衰減系數(shù)0.991 8；揚(yáng)壓力初始設(shè)計值0.30，衰減系數(shù)1.088 4。

表2 不確定參數(shù)界限

注：該壩已服役17年，由衰減函數(shù)得衰減系數(shù)為0.918 5。

圖3 壩體單元可靠指標(biāo)β分布

圖4 巖基單元可靠指標(biāo)指標(biāo)β分布

4 結(jié) 論

(1)考慮到重力壩概率可靠性分析時其參數(shù)實(shí)時獲取難特點(diǎn)，本文基于參數(shù)的衰減函數(shù)獲取重力壩服役后的參數(shù)統(tǒng)計信息，可更加真實(shí)反映重力壩的可靠狀態(tài)。

(2)本文基于Kriging模型通過合理選取并確定參數(shù)概率分布信息，擬合并校驗(yàn)得到最佳極限狀態(tài)方程以代替重力壩極限狀態(tài)方程，有效地適應(yīng)了重力壩極限狀態(tài)方程高度非線性甚至非顯式表達(dá)的特點(diǎn)。

圖5 抗滑面單元可靠指標(biāo)指標(biāo)β統(tǒng)計

(3)本文將重力壩可靠指標(biāo)求解問題轉(zhuǎn)化為優(yōu)化問題，構(gòu)建了重力壩可靠指標(biāo)優(yōu)化模型，實(shí)現(xiàn)了重力壩所有單元強(qiáng)度破壞及沿壩基面抗滑穩(wěn)定可靠指標(biāo)高效計算，為深入分析重力壩可靠狀態(tài)提供了依據(jù)。