陸 軍，武清璽

（河海大學 土木工程學院，南京 210098）

拱壩是一種重要的壩型，以其結構合理和體形優(yōu)美而著稱，其所有的優(yōu)良性能建立在拱壩穩(wěn)定可靠的基礎上。壩體一旦失事，將嚴重影響工農業(yè)生產(chǎn)，因此拱壩的可靠度研究和安全分析早已得到人們的重視。本文將考慮基巖材料的非線性性質，研究非線性方程的求解方法，并將材料的力學特性、拱壩受到的外來作用等作為隨機變量，以有限元計算的結果為基礎，利用二次響應面法計算拱壩結構重要部位單元的可靠度，通過整理和分析得出壩體可靠指標的變化規(guī)律。

1 結構可靠度計算的響應面法

目前，結構的極限狀態(tài)方程一般都基于抗力－荷載效應模型?，F(xiàn)有可靠度計算方法多是以極限狀態(tài)方程具有明確的解析表達式為基礎的。但是對于某些復雜結構系統(tǒng)，由于結構本身的復雜性，基本隨機變量的輸入與輸出量之間的函數(shù)關系常常是高度非線性的，有時甚至得不到明確的解析表達式，在計算這類復雜結構的可靠度時，采用一次二階矩法、JC法等就存在一定的困難。Box和Wilson于1951年首先提出響應面法，當時對于響應面法的研究僅限于如何用統(tǒng)計的方法得到一個顯式函數(shù)，用來逼近一個復雜的隱式函數(shù)。直到1985年，文獻［1］首先用一次響應面法分析研究了土坡穩(wěn)定的可靠度問題，而后文獻［2－4］研究才逐步拓展了這一領域。響應面法是統(tǒng)計學的綜合試驗技術，用于處理復雜系統(tǒng)的輸入（隨機變量）與輸出（系統(tǒng)響應）的轉換關系問題。該方法采用有限的試驗，通過回歸擬合解析表達式Z′＝g′（X）代替真實曲面Z＝g（x）。

在可靠度計算中，設極限狀態(tài)方程為：

該方程將n維隨機空間分為失效區(qū)和安全區(qū)兩部分。結構的失效概率為：

式中 fX（x1，x2，…，xn）為n維隨機向量的聯(lián)合概率密度函數(shù)。在多數(shù)實際問題中，不但極限狀態(tài)方程是非線性的，而且n維隨機向量的聯(lián)合概率密度函數(shù)也難以得到，所以精確求解式（2）非常困難。對于大型、復雜的工程結構，可靠度計算常采用數(shù)值方法，如隨機有限元法等。由于基巖的本構關系是非線性的，利用隨機有限元法求解要經(jīng)過多次迭代，計算繁瑣且結果的誤差也較大，難以滿足精度要求。近年來得到迅速發(fā)展的結構可靠度計算的響應面法，可以應用確定性有限元計算程序對結構進行可靠度分析，使這一方法展現(xiàn)出良好的應用前景。對于非線性極限狀態(tài)方程，采用在設計驗算點處的一次Taylor展開式求得結構的可靠指標的精度不理想。二次可靠度分析法的計算精度雖有所提高，但工作量成倍增加。結構可靠度計算的二次響應面法與一次二階矩法相比有2個優(yōu)點：①用一個待定的顯式響應面函數(shù)逼近實際的隱式極限狀態(tài)函數(shù)，使可靠度計算得到簡化；②計算結構可靠指標時以二次響應面取代一次近似式，使計算精度有所提高。

設二次響應面函數(shù)為：

式中 a，bi，ci為待定系數(shù)，需由2n＋1個樣本點確定。 這一方法的計算過程為：第1步，以均值μX點為中心，樣本點在區(qū)間（μX-fσX，μX+fσX）中選取，建議f的取值為1～3；用g（X）在2n＋1個樣本點處的函數(shù)值，可以確定響應面函數(shù)中的待定系數(shù)；響應面函數(shù)得到以后，即可求出原極限狀態(tài)面上設計驗算點的近似值XD。第2步，選取新的樣本中心點，新的樣本中心點XM選在由均值點mX到XD的直線上，并保證在新的中心點處滿足：

這樣選取新中心點的目的是為了使所選樣本點包含原極限狀態(tài)面更多的信息。第3步是以為中心點選取新的一組樣本點，重復上述過程，即可得到新的響應面函數(shù)，并近似求出相應的可靠指標。整個過程需求出4n＋3個g（X）值，這一方法的逼近過程如圖1所示［5-8］。

2 非線性有限元

2.1 材料的破壞準則與本構模型

基于剪壓破壞的理想彈塑性材料模擬基巖，屈服準則采用Mohr－Coulomb準則［9］：

壩體混凝土材料采用彈性模型，采用拉裂壓碎破壞準則：

式中 c，φ分別為材料的抗剪斷粘聚力和內摩擦角；σ1，σ3分別為單元的最大、最小主應力；ft，fc分別為材料的抗拉、抗壓強度的設計值。

2.2 非線性有限元求解方法

考慮巖體為彈塑性材料，材料的性質與應力和變形的歷史有關，本構方程用增量形式表達。按載荷作用的實際情況劃分為若干個子步，在小的載荷增量下逐步計算。在某一載荷增量{Δδ}n的作用下，有限元的平衡條件為：

式中 {ΔR}n為載荷增量，求解式（8）得位移增量{Δδ}n。 在載荷增量{ΔR}n施加之前，已作用有累積載荷{R}n-1，結構中的應力、應變和位移均已知道，分別為{σ}n-1，{ε}n-1和{δ}n-1。 那么在新的累積載荷{R}n={R}n-1+{ΔR}n作用下，其總的應力、應變和位移分別見（9）式［10］：

3 工程算例分析

3.1 工程概況

福建省永泰縣青龍溪大洋水庫大壩壩型設計為實體混凝土三心變厚雙曲拱壩，壩頂高程為236.0m，最大壩高87.2m，大壩為Ⅲ等3級建筑物。水庫的正常蓄水位230m，校核洪水位236.4m。工程區(qū)位于區(qū)域地質構造相對穩(wěn)定地段，未見區(qū)域性的斷裂構造通過壩址，河槽段無斷層構造通過，局部節(jié)理裂隙發(fā)育地帶采用加密固結灌漿孔工程處理。

3.2 計算模型

本次計算采用三維模型進行計算，坐標系以壩頂拱冠梁上游點為坐標原點，左岸為X軸正方向，下游為Y軸正方向，豎直向下為Z軸正方向。模擬范圍為上游0.8倍壩高，下游1.5倍壩高，左右兩岸1.5倍壩高，壩基深度為1倍壩高，壩頂高程以上模擬20m山體。

網(wǎng)格的劃分依據(jù)福建省水利水電科學研究所提供的拱壩勘察設計資料進行，網(wǎng)格采用八節(jié)點六面體單元，節(jié)點總數(shù)152208，單元總數(shù)為141064，其中壩體單元總數(shù)為5244。離散后的整體模型如圖2所示。計算中對壩基巖體材料應力應變關系采用基于增量法的彈塑性模型，采用對巖體強度參數(shù)進行折減的方法分別模擬右岸和壩體河槽段巖體的兩條軟弱帶如圖2，壩體混凝土材料采用彈性模型。

3.3 計算參數(shù)及計算結果

由地質勘查資料，壩基網(wǎng)格根據(jù)材料物理參數(shù)的不同分為河槽段的巖石1和兩岸的巖石2單元組；右岸和河槽段的軟弱帶分別定義為軟弱1和軟弱2單元組；壩體單元為混凝土單元組。

計算擬定的荷載組合為：自重＋校核洪水位（236.4m）＋設計溫升（升溫溫度場與壩體穩(wěn)定溫度場的差值）。為了模擬大壩施工和逐級蓄水的過程，采用了分級加載的方式。其中，壩體重力荷載分6級加載；水荷載分9級加載；溫度荷載分6級加載。不考慮地應力，僅計入巖體自重所形成的初始應力場，采用定值考慮的材料計算參數(shù)如表1所示。

現(xiàn)對各隨機變量取均值時的計算結果加以分析。（應力以拉為正，壓為負）壩體上、下游面單元的安全性是最為關心的，所以本文首先分析了這些區(qū)域的應力計算結果，如圖3所示。由圖3（a）、（b）可以看出，壩體最小主應力分布規(guī)律是：上游面單元由拱冠梁處向兩側逐漸變小，由頂部向下部先逐漸變大再逐漸變小，在距壩頂約1/3壩高處，壓應力達到最大值；下游面單元的最小主應力由拱冠梁處向兩側逐漸變大，并在臨近壩肩處達到最大值。由于壩體形狀基本上左右對稱，所以得到的應力也大致為對稱分布。壩體拉應力分布區(qū)域較小，主要分布在壩踵附近，拉應力最大值為1.5MPa。

表1 材料計算參數(shù)表

圖3 擬定荷載組合下拱壩應力云圖

3.4 拱壩可靠度計算及分析

設計人員總是力求使混凝土設計參數(shù)取值合理、可行，盡可能使混凝土具有優(yōu)良的物理、力學性能。拱壩建成運行后上游面的水壓力、混凝土水化熱、上游面庫水和下游面氣溫的溫差都是壩體的主要荷載。其中水荷載作為主要荷載是無法避免的，溫度荷載可以通過精心組織施工減小其影響，所以本文將壩體上游水位、基巖、混凝土的某些物理力學參數(shù)作為隨機變量考慮，對可靠度計算結果影響甚微和缺乏統(tǒng)計資料的參數(shù)作為定值考慮。計算中采用的隨機變量的特征值如表2所示。

表2 計算中所采用的隨機變量統(tǒng)計特征值

壩體混凝土的抗裂或抗壓可靠度極限狀態(tài)方程可表示為：

式中 R為壩體混凝土的抗裂強度或抗壓強度，分別對應通過非線性有限元計算出的混凝土單元的第一和第三主應力；S為荷載效應，是關于隨機變量E1，E2，E3，tg準，c，Hu的隱式函數(shù)。右岸巖體的薄弱帶基于Mohr－Coulomb屈服準則的可靠度極限狀態(tài)方程可表示為：

式中 c，準分別為右岸巖體薄弱帶的粘聚力和摩擦角；σ1，σ3分別是通過非線性有限元計算出的該薄弱帶單元的第一和第三主應力，是關于隨機變量E1，E2，E3，tg準，c，Hu的隱式函數(shù)。

3.5 計算結果

利用有限元計算的應力結果，再結合響應面法可以得到壩體單元的可靠指標，得到足夠多單元的可靠指標后再繪制出可靠指標的等值線圖，比較直觀的反映出壩體結構的重要部位的可靠度情況。以48號單元為例給出某一單元的可靠指標計算過程。

（2）響應面函數(shù)確定后再利用梯度優(yōu)化法求出極限狀態(tài)面上設計驗算點的近似值XD為（14.9943，20.0042，19.9393，0.4659，0.7001，74.5734，1.7972）。 按照上文介紹的方法選取新的中心點XM（14.9929，20.0053，19.9243，0.4658，0.7002，76.9225，1.7966）。

（3）以新的中心點XM為中心開始選取新的一組樣本點，重復第一步的工作，即可求得極限狀態(tài)面上設計驗算點的值和相應的可靠指標 β=4.5。

本文取壩體主要區(qū)域的320個單元計算其可靠指標，繪制出抗壓可靠指標的等值線圖如圖4。該圖清晰地反映了拱壩結構的可靠指標分布情況。拱壩上游面拱冠梁自由端附近的單元、下游面壩肩附近的壩體單元可靠指標均較小。從圖4a、圖4b圖看出壩體單元抗壓可靠指標大小的變化趨勢與其主應力變化趨勢與圖3相同。由于拱壩結構內部以壓應力為主，拉應力主要分布在建基面附近。拉應力的最大值出現(xiàn)在壩踵處，此處單元的抗拉可靠指標的最小值為1.1，需采取措施增大材料的抗拉能力。軟弱1和軟弱2單元組中單元的可靠指標均在5.0以上，說明在給定的荷載下是比較安全的。

4 結語

（1）混凝土拱壩及壩基材料均具有多樣的物理、力學特性，其離散性強、變異性大，本構關系是非線性的。在對其作可靠度分析時，直接得到系統(tǒng)極限狀態(tài)方程的解析表達式往往是非常困難的，并且聯(lián)合概率密度函數(shù)也難以得到。響應面法結合有限元分析方法，為解決類似工程的可靠性分析提供了可供選用的一條新途徑。

（2）可靠度計算時采用響應面法解決拱壩類復雜結構問題穩(wěn)定性較好，計算可靠指標的梯度優(yōu)化法收斂速度快，一般僅需3～5次迭代計算。

（3）采用商用有限元計算軟件結合自編的FORTRAN程序實現(xiàn)結構可靠指標的計算方法在PC上解決本論文的計算具有實用價值。

［1］Wong F S.Slope reliability and response surface method［J］.Journal of Geotechnical Engineering，1985，b111（1）：32-53.

［2］Liu Y W，Moses F.A sequential response surface method and its application in the reliability analysis of aircraft structural systems［J］.Structural Safety，1994，16：39-46.

［3］武清璽，卓家壽.結構可靠度分析的變f序列響應面法及其應用［J］.河海大學學報，2001，29（2）：75-78.

［4］武清璽，俞曉正.混凝土面板堆石壩可靠度計算方法研究［J］.巖土工程學報，2004，26（4）：468-472.

［5］武清璽.結構可靠性分析及隨機有限元法［M］.北京：機械工業(yè)出版社，2005.

［6］結構可靠度理論［M］.北京：中國建筑工業(yè)出版社，2005.

［7］吳世偉.結構可靠度分析［M］.北京：人民交通出版社，1990.

［8］熊鐵華，常曉林.響應面法在結構體系可靠度分析中的應用［J］.工程力學，2006，23（4）：58-61.

［9］李詠偕，施澤華.塑性力學［M］.北京：中國水利水電出版社，1987.

［10］朱伯芳，有限單元法原理與應用［M］.北京：中國水利水電出版社，1998.