考慮壩體材料參數(shù)隨機性的拱壩易損性分析

徐 強,姚 文 斌,陳 健 云,張 天 然,曾 祥 興

(1.大連理工大學 海岸與近海工程國家重點實驗室,遼寧 大連 116024; 2.大連理工大學 工程抗震研究所,遼寧 大連 116024)

0 引 言

中國擁有豐富的水能資源,但是其空間分布十分不均勻,多分布于西南地區(qū),而西南地區(qū)多為河谷狹窄的山區(qū),因此,拱壩成為該地區(qū)大壩修建的主要壩型之一,庫容往往能達到數(shù)十億甚至上百億m3[1]。然而,中國西部地區(qū)位于喜馬拉雅火山地震帶上,地震較為頻繁,強地震發(fā)生時,壩體材料在較大的應(yīng)力作用下易失效。一旦拱壩失事,將給下游造成巨大的生命財產(chǎn)損失。針對大壩破壞失事方面的模擬,隨著計算機水平的不斷提高和數(shù)值分析方法的不斷完善,有限元分析方法逐漸成為研究壩體在不同工況下?lián)p傷破壞情況的重要手段。

由于運輸、施工和筑壩材料屬性等因素的影響,壩體不同部位的混凝土材料參數(shù)存在一定的空間隨機性,然而,傳統(tǒng)的有限元分析中,假定大范圍的筑壩材料參數(shù)為一個確定的值,無法反映出壩體實際的材料分布情況,造成分析結(jié)果與實測值有一定差距。如果能將這一空間隨機性在數(shù)值分析過程中充分模擬,分析結(jié)果會更能反映壩體的實際損傷分布,從而為拱壩的結(jié)構(gòu)設(shè)計和加固提供依據(jù)。

目前,對混凝土材料隨機性進行模擬的理論主要有隨機變量理論和隨機場理論。隨機變量理論假定空間內(nèi)任意點之間的材料參數(shù)相互獨立,沒有相關(guān)性,該理論的優(yōu)點在于簡單易于實現(xiàn),缺點是忽略了結(jié)構(gòu)不同部位之間的相關(guān)性。而隨機場理論則充分考慮了材料參數(shù)的隨機性和空間位置上的相關(guān)性,更能反映實際的結(jié)構(gòu)參數(shù)特性。Ardebili等[2]以Koyna重力壩為研究對象,將混凝土彈性模量、質(zhì)量密度和抗拉強度均假定為隨機場,并采用協(xié)方差矩陣分解和中點離散技術(shù)生成有限元模型展開研究,研究結(jié)果發(fā)現(xiàn)混凝土的非均質(zhì)性影響著結(jié)構(gòu)的抗震性能評價;Carl和 Padgett等[3]在對大型混凝土重力壩的可靠度分析中,同時考慮了地震動和壩體材料的隨機性;陳健云等[4]在對Hardfill壩的有限元分析中同時考慮材料密度、拉伸彈性模量和抗拉強度等參數(shù)的空間隨機性,并通過蒙特卡洛模擬方法進行分析,驗證了材料參數(shù)的隨機性對結(jié)構(gòu)損傷的影響。

在研究壩體在不同地震動強度下的非線性響應(yīng)時,需要對加速度時程反復(fù)進行調(diào)幅計算,計算量大且計算結(jié)果可能不收斂。國外學者Estekanchi提出了耐震時程法(ETA)[5],且詳細闡述了耐震時程加速度曲線的合成方法,并驗證了該方法在保證計算精度的同時可有效減少了計算量。白久林等[6]合成了基于中國抗震反應(yīng)譜的ETA時程曲線,與增量動力分析方法(IDA)對比后,驗證了ETA方法的可行性。

2003年,由中國學者李杰、陳建兵等[7]提出的概率密度演化方法(Probability Density-Evolution Method)近些年來逐漸發(fā)展完善,并得到很好的應(yīng)用。該方法通過推導出解耦的廣義概率密度演化方程,求解后得到響應(yīng)全過程的概率信息[8]。概率密度演化方法對于非線性強、自由度數(shù)量高的復(fù)雜結(jié)構(gòu)也有較好的適用性。陳健云等[9]結(jié)合廣義概率密度演化理論,提出了一種混凝土壩全周期抗震性能分析方法,并與蒙特卡洛模擬方法的計算結(jié)果進行對比,驗證了概率密度演化方法的可行性。

本文根據(jù)隨機場理論建立了壩體材料參數(shù)隨機的白鶴灘拱壩模型,根據(jù)標準設(shè)計反應(yīng)譜生成了50組ETA人工加速度時程曲線,并作用于模型底部,分析了相關(guān)長度和參數(shù)平均值對壩體結(jié)構(gòu)損傷分布的影響。選取壩頂位移、損傷體積比和橫縫開度為損傷評價指標,對50組分析結(jié)果進行整理并建立其概率密度演化模型,根據(jù)概率密度曲線得出拱壩的可靠度,以評價白鶴灘拱壩在地震動作用下的安全性。

1 理論方法

1.1 混凝土材料空間變異性模擬

采用隨機場對拱壩進行模擬時,首先要對隨機場進行離散。合理的離散有助于提高模擬的精度。本文以斷裂能和抗拉強度為隨機變量,基于協(xié)方差矩陣分解法和中心點法生成隨機場[2]。

假設(shè)隨機場單元和有限元模型網(wǎng)格劃分單元一致,共生成n組隨機場。采用指數(shù)型自相關(guān)函數(shù):

(1)

式中:ζ為兩單元之間的距離,lcorr為材料的相關(guān)長度。壩體不同部位材料的相關(guān)性取決于相關(guān)長度lcorr,相關(guān)長度的含義是若壩體上任意兩點之間的距離大于該長度,則這兩點的材料參數(shù)相互獨立,反之,如果結(jié)構(gòu)上兩點之間距離小于相關(guān)長度,則兩點的材料參數(shù)值具有相關(guān)性,即相關(guān)長度越大,壩體材料越趨向于均值。

由離散化的自相關(guān)函數(shù)構(gòu)造相關(guān)矩陣C:

C=(ρij)m×m

(2)

式中:m為模型單元數(shù)。

構(gòu)造矩陣Q,并求其正態(tài)分布函數(shù)值:

(3)

V=F-1[Φ(Q)]

(4)

需要說明的是,V(i,j)的具體含義為假設(shè)材料參數(shù)服從標準正態(tài)分布,其值為負無窮到第i組隨機場第j個單元的材料參數(shù)的分布函數(shù)值。

假設(shè)壩體材料參數(shù)服從韋伯分布,則i組隨機場各單元參數(shù)的分布函數(shù)值取V(i,j),求出其逆累積分布函數(shù),得出參數(shù)值。

1.2 耐震時程法(ETAM)

耐震時程法作為一種動力推覆過程,其優(yōu)點在于保證預(yù)測結(jié)構(gòu)響應(yīng)精度的同時,可以很大程度減少計算量[10]。其思路是生成隨時間逐漸遞增的ETA加速度時程曲線,且不同周期對應(yīng)的反應(yīng)譜呈既定的倍數(shù)擬合關(guān)系。

ETA加速度時程的生成思路是先生成平穩(wěn)加速度時程,經(jīng)反復(fù)迭代優(yōu)化后,使之各時段加速度對應(yīng)的反應(yīng)譜與目標反應(yīng)譜有更好的擬合。在地震動持續(xù)時間內(nèi),不同時間段對應(yīng)的目標反應(yīng)譜與目標反應(yīng)譜為倍數(shù)關(guān)系:

(5)

式中:TTarget為目標時間;SaC(T)為標準設(shè)計反應(yīng)譜;t為擬合反應(yīng)譜的時間;SaT(T,t)是0～t時刻的目標加速度反應(yīng)譜。

根據(jù)加速度反應(yīng)譜與位移譜的數(shù)學關(guān)系可知,目標位移反應(yīng)譜為

(6)

式中:SuT(T,t)為時刻t的目標位移反應(yīng)譜。

想要讓目標加速度時程在每一點都能同時較好地擬合加速度反應(yīng)譜和位移反應(yīng)譜幾乎是不可能的,但可以采用式(7)對加速度時程進行無約束優(yōu)化,從而減小同目標反應(yīng)譜之間的誤差:

α[Su(T,t)-SuT(T,t)]2}dtdT

(7)

式中:ag是初始生成的地震動時程曲線,Tmax是地震動時間最大值;α是位移譜的權(quán)重系數(shù),本文只考慮加速度反應(yīng)譜的影響,因此α值取0;Sa(T,t)和Su(T,t)分別表示周期T下0～t時刻的加速度反應(yīng)譜和位移反應(yīng)譜。

圖1為ETA時程和IDA時程的對比圖。其中,ETA時程對0～5 s、0～10 s、0～15 s和0～20 s四個時段對應(yīng)的反應(yīng)譜進行優(yōu)化擬合,并計算不同時段時程對應(yīng)的反應(yīng)譜如圖2所示。得到的ETA時程曲線具有兩方面的特點:① 隨著時間的推移,其峰值加速度不斷增大;② 不同時段對應(yīng)的加速度反應(yīng)譜與標準譜成倍數(shù)關(guān)系。

注:左為ETA法,右為IDA時程法。圖1 ETA法與IDA法時程對比Fig.1 Comparison of ETA Method and IDA Method

圖2 ETA時程加速度反應(yīng)譜Fig.2 ETA time history acceleration response spectrum

1.3 概率密度演化法(PDEM)

通常,壩體在地震激勵下的動力方程可以表示為

(8)

Z(t)=H(θ,t)

(9)

求導后,廣義速度向量可表示為

(10)

由于沒有其他隨機因素的加入和消失,因此根據(jù)概率守恒原理經(jīng)推導可得[12]:

(11)

當僅考慮系統(tǒng)中某一個分量時,式(11)的概率密度演化方程簡化為

(12)

式中:pzθ(z,θ,t)為Z(t)與θ的聯(lián)合概率密度函數(shù)。

根據(jù)初始條件和邊界條件對偏微分方程進行求解,得到感興趣變量的概率密度函數(shù):

(13)

對于式(13)這種高維積分問題,很難得出解析解,因此,本文分別采用MOL算法和具有TVD性質(zhì)的有限差分格式進行求解[13]。

1.4 直接概率積分法(DPIM)

若以θ和Y分別表示輸入隨機向量和輸出隨機向量,根據(jù)概率守恒原理可知:

(14)

其中,Pθ(θ)和PY(Y)分別表示θ和Y的概率密度函數(shù),Ωθ和ΩY分別表示輸入隨機向量和輸出隨機向量的樣本空間。通過推導,動力系統(tǒng)中任意時刻t的輸出隨機響應(yīng)的概率密度函數(shù)為

(15)

方程(15)稱為概率密度積分方程,由于該函數(shù)的求解過程中涉及到狄拉克函數(shù)以及高維度積分,導致該方程地求解極為困難。利用直接概率積分法可以高效地求解概率密度積分方程,得到感興趣參數(shù)的概率密度函數(shù)。直接概率積分法包含兩方面的重要技術(shù):① 基于GF-偏差法選取代表點集并對概率空間進行剖分[14];② 將原方程中不連續(xù)的高斯函數(shù)替換為連續(xù)高斯函數(shù),進行光滑化處理。運用下式計算輸出隨機響應(yīng)概率密度函數(shù)[15]:

(16)

進一步推導后可得:

(17)

式中:θq表示第q個代表點;N表示所選代表點數(shù)量;Pq為該代表點的賦得概率;σ為光滑化參數(shù)。

2 隨機參數(shù)的選取

考慮到拱壩壩體材料和地震動都具有隨機性,因此以相關(guān)長度表征壩體材料的隨機性,結(jié)構(gòu)上距離超過相關(guān)長度的兩點認為其隨機材料屬性互相獨立,本文認為相關(guān)長度服從壩高10%～60%的均勻分布[1];以標準地震動反應(yīng)譜的特征周期表征地震動的隨機性,認為其服從0.20～0.35 s的均勻分布。

采用matlab編程在二維空間中選取代表點,并基于GF-偏差優(yōu)化方法對選取點集進行優(yōu)化。點集的GF值表征點集與理論分布的偏差程度。通過優(yōu)化,點集的GF值從0.32下到0.28,優(yōu)化效果良好。采用維諾圖計算各代表點的賦得概率,點集分布如圖3所示。

圖3 代表點選取結(jié)果Fig.3 Results of representative point selection

3 拱壩有限元計算

3.1 有限元計算模型

以白鶴灘高拱壩為研究對象,建立其有限元模型如圖4所示。其最大壩高為289 m,壩頂和壩底最大厚度分別為13 m和72 m,白鶴灘壩體-地基有限元模型由拱壩和具有輻射阻尼的巖石地基組成;壩體共設(shè)置了30條橫縫。模型由25 776個六邊形實體單元組成,其中包括4 656個壩體單元和21 120個地基單元。

圖4 白鶴灘拱壩有限元模型Fig.4 Finite element model of Baihetan arch dam

模型考慮斷裂能和抗拉強度服從威布爾分布,其他參數(shù)設(shè)置為定值?？紤]自重影響,壩體混凝土密度為2 400 kg/m3,彈性模量為36 GPa,泊松比0.167;基巖密度為2 800 kg/m3,彈性模量26 GPa,泊松比0.24;靜水壓力方面,上游水位825 m,下游水位604 m;動水壓力按照Westergaard附加質(zhì)量進行施加,來模擬動水壓力。

3.1.1混凝土本構(gòu)模型

本文采用塑性損傷本構(gòu)模型,其應(yīng)力-應(yīng)變關(guān)系如圖5所示,在不同階段,其關(guān)系式為

圖5 塑性損傷本構(gòu)的應(yīng)力應(yīng)變關(guān)系曲線Fig.5 Stress-strain curve of plastic damage constitutive

(18)

3.1.2拱壩壩體材料隨機性的實現(xiàn)

根據(jù)不同參數(shù)對壩體抗震性能的影響程度,選取拱壩混凝土材料的抗拉強度和斷裂能為隨機變量,其他參數(shù)為定值。假定抗拉強度服從1.5～3.5 MPa的均勻分布,斷裂能服從100～300 N/m的均勻分布,根據(jù)選取的50組代表點對應(yīng)的不同相關(guān)長度,采用協(xié)方差矩陣分解法和中心點法生成50組隨機場,按相關(guān)長度從小到大排列。圖6(a)展示了由上述理論生成的白鶴灘拱壩隨機場模型,從上至下相關(guān)長度在變化范圍內(nèi)逐漸增大。與圖6(b)[2]所示的國外學者基于同一理論生成的Koyna重力壩隨機場模型進行對比可以發(fā)現(xiàn),其隨機性分布較為一致。隨著相關(guān)長度的增大,壩體材料分布逐漸趨向于均勻化,在較大范圍內(nèi)材料參數(shù)趨于一致,說明應(yīng)用隨機場理論能較好地對壩體材料的隨機性進行模擬。

圖6 不同相關(guān)長度對應(yīng)的隨機場Fig.6 Random fields corresponding to different correlation lengths

3.1.3橫縫模型

白鶴灘拱壩共設(shè)置30條橫縫,橫縫分布如圖7所示。在本研究中,采用接觸邊界模型模擬混凝土拱壩的橫縫,橫縫的接觸關(guān)系[16]如圖8所示,法向接觸壓力的計算公式為

圖7 壩體橫縫分布Fig.7 Distribution of dam transverse joints

(19)

式中:初始接觸距離c取100 mm,初始壓力p0取0.3 MPa。

每個橫縫接觸面包含主面和從面,主從面的接觸關(guān)系如圖8所示。

圖8 主從面接觸模型Fig.8 Dam transverse joint model

3.1.4黏彈性人工邊界

黏彈性人工邊界可以較為合理地模擬地基與結(jié)構(gòu)的相互作用關(guān)系,具有多方面優(yōu)點,不僅可以吸收邊界輻射能量,同時可以反映地基的彈性恢復(fù)能力[17],其如圖9所示。

圖9 黏彈性人工邊界模型Fig.9 Viscoelastic artificial boundary model

人工邊界上,節(jié)點上的參數(shù)由式(20)計算給出:

(20)

式中:Kn為彈簧-阻尼元件法向剛度系數(shù);Cn為彈簧-阻尼元件法向阻尼系數(shù);Ks為彈簧-阻尼元件切向剛度系數(shù);Cs彈簧-阻尼元件切向阻尼系數(shù);Cp,Cs分別為P波和S波波速;r為結(jié)構(gòu)幾何中心到人工邊界點的等效距離;ρ為介質(zhì)密度;G為剪切模量;λ為第一拉梅常數(shù);a,b為參數(shù)。

3.2 ETA地震動時程

根據(jù)選取的特征周期生成標準設(shè)計反應(yīng)譜作為目標譜,通過與0～5 s、0～10 s、0～15 s和0～20 s等四個不同時間段內(nèi)目標反應(yīng)譜進行擬合,生成3條ETA人工地震動時程,作為3個方向的輸入地震動。圖10給出了Tg=0.322 s時,3個方向的ETA加速度時程,隨著時間的增加,其加速度幅值不斷增大。圖2顯示加速度時程對應(yīng)的加速度反應(yīng)譜在標準反應(yīng)譜周圍振蕩,擬合效果良好。

圖10 ETA時程Fig.10 ETA time history

值得注意的是,圖10所示的ETA時程在不同時間點對應(yīng)不用的峰值加速度,進行一次時程分析便可以得到不同峰值加速度下的動力響應(yīng),相較IDA分析方法,無需進行大量調(diào)幅計算,顯著降低了計算量,適用于對大型結(jié)構(gòu)的分析。

3.3 動力響應(yīng)分析

3.3.1響應(yīng)分析

通過ABAQUS軟件計算后,發(fā)現(xiàn)壩體下游面損傷主要集中的中上部,部分工況下會形成貫穿性損傷,上游面損傷主要集中在拱端,損傷符合一般規(guī)律?？紤]壩體材料的隨機性后,損傷的分布更加離散,且壩體的損傷情況主要會受到壩體材料參數(shù)均值和相關(guān)長度的影響。圖11(a)展示了3個具有相近均值,但相關(guān)長度不同的3個拱壩有限元模型的材料參數(shù)分布情況,圖11(b)展示了峰值加速度為0.6g時,以上3種工況的壩體損傷情況。通過對比,在壩體材料均值接近的情況下,相關(guān)長度越小,壩體的損傷越嚴重,隨著相關(guān)長度的增大,壩體的抗震能力逐漸提高。同樣,圖12對比了相關(guān)長度增大,但均值逐漸減小的3組工況,經(jīng)對比,雖然相關(guān)長度不斷增大,但隨著材料參數(shù)的下降,壩體的破壞逐漸加重,說明壩體材料參數(shù)的均值也是決定壩體抗震能力的關(guān)鍵因素。

圖11 參數(shù)均值相近,相關(guān)長度不同工況的計算結(jié)果Fig.11 Results under similar mean parameters and different correlation lengths

圖12 相關(guān)長度增大,參數(shù)均值減小工況的計算結(jié)果Fig.12 Results under decreasing mean parameters and increasing correlation lengths

3.3.2概率密度分析

采用橫縫開度、損傷體積比和相對位移作為損傷評價指標。計算結(jié)果顯示,隨著地震動峰值加速度(PGA)值增大,各指標呈現(xiàn)明顯的增大趨勢,不同工況的分析結(jié)果離散性較高,考慮為壩體材料的隨機性導致。采用PDEM法和DPIM法分別計算了各指標概率演化過程。圖13展示了隨PGA值的增大,不同指標的變化趨勢。各指標均隨時間的推移逐漸增大,且分布范圍逐漸變大。圖14展示了橫縫開度的概率分布情況。從等高線圖中可以發(fā)現(xiàn),各時刻的概率分布近似服從正態(tài)分布,均值約等于該時刻的平均橫縫開度,且PGA越大,方差越大,概率分布曲線由“高瘦型”向“矮胖型”轉(zhuǎn)變,說明PGA越大,離散型越高。

圖13 各評價指標隨PGA變化曲線Fig.13 Curve of evaluation index with the change of PGA

圖14 橫縫開度概率密度演化過程Fig.14 Evolution process of probability density of transverse joint openness

在PGA較小時,各指標概率分布范圍較為集中,隨著PGA的增大,分布范圍趨于離散。兩種概率密度方法計算結(jié)果高度一致,但DPIM法的平均計算速度是PDEM法的20倍左右,優(yōu)勢明顯。

為對比直接概率積分法和概率密度演化方法的計算結(jié)果是否吻合,取10 s、15 s和20 s為特定時刻,對各特定時刻的概率密度分布情況以及累計概率分布情況進行對比。結(jié)果顯示各指標使用兩種方法計算的結(jié)果差距較小,DPIM法計算結(jié)果在局部略有振蕩,但幅度較小。圖15展示了以壩體橫縫開度為指標的計算對比結(jié)果,各時刻均較為接近。

圖15 概率密度分布及累計概率密度在特定時刻計算結(jié)果對比Fig.15 Comparison of probability density and accumulated probability density at specific time

3.3.3失效概率分析

失效概率是表征結(jié)構(gòu)破壞程度的有效指標,當結(jié)構(gòu)參數(shù)達到閾值后,認為結(jié)構(gòu)發(fā)生破壞,反之,認為結(jié)構(gòu)是安全的。首先采用概率密度方法計算結(jié)構(gòu)的失效概率[18]。根據(jù)3.3.2節(jié)中所求概率密度曲線,通過式(21)積分后求得失效概率:

(21)

式中:XL為參數(shù)的閾值,閾值的選取可以參照文獻[19],計算結(jié)果如圖16所示。通過對比,直接概率積分法與概率密度演化法的計算結(jié)果較為吻合,再次驗證了該方法的正確性。各項指標計算的失效概率快速增大的區(qū)域均集中在0.2g～0.4g之間,且互相之間較為吻合,具有一定等效關(guān)系。

圖16 評價指標失效概率曲線Fig.16 Failure probability curve of evaluation index

通過對比兩種方法計算的各損傷評價指標失效概率曲線發(fā)現(xiàn),以損傷體積比為評價指標時,兩種方法的計算結(jié)果在初期有所差異;隨著PGA的增大,基于相對位移指標計算的失效概率曲線增速較其他兩種指標計算的失效概率曲線更加平緩,說明不同損傷評價指標具有一定離散性;使用橫縫開度作為評價指標計算的結(jié)果較為合理。

4 結(jié) 論

本文引入隨機場理論,模擬拱壩混凝土材料參數(shù)的隨機性。為避免偶然性,計算多組樣本,并基于維諾圖賦以相應(yīng)概率。基于耐震時程法(ETA)理論,根據(jù)不同特征周期對應(yīng)的標準譜生成ETA人工地震動時程。通過對拱壩有限元模型的分析研究了抗拉強度和斷裂能隨機分布對白鶴灘拱壩損傷分布的影響。將各工況結(jié)果乘以相應(yīng)賦得概率并求和后,作為最終結(jié)果,并運用概率密度方法建立了各指標隨時間變化的時變概率信息,得到其概率密度演化曲線,曲線的走向與參數(shù)均值保持一致。設(shè)定損傷閾值,基于概率密度曲線對超過閾值的部分進行積分得到失效概率,建立了峰值加速度與失效概率的關(guān)系。

研究結(jié)果表明:考慮混凝土抗拉強度和斷裂能的隨機性后,白鶴灘拱壩動力響應(yīng)具有明顯的離散性,不同工況的分析結(jié)果隨著PGA的增大其差異逐漸增大;分析可知,當壩體材料隨機參數(shù)均值相近時,隨著相關(guān)長度的增大,壩體的損傷情況呈減小趨勢,即壩體材料參數(shù)越趨向于均值,壩體抗震性能越好;而當壩體材料隨機參數(shù)均值不同時,即使相關(guān)長度較大,參數(shù)均值較小的破壞依然嚴重。經(jīng)對比,運用橫縫開度和損傷體積比作為損傷評價指標計算結(jié)構(gòu)的失效概率時,其計算結(jié)果更為接近。