基于混凝土塑性損傷本構的高拱壩損傷開裂分析

張建偉，劉鵬飛，王 濤，焦延濤，李兆恒

(1.華北水利水電大學 水利學院，鄭州 450046；2.水資源高效利用與保障工程河南省協同創新中心，鄭州 450046；3.河南省水工結構安全工程技術研究中心，鄭州 450046；4.廣東省水利水電科學研究院，廣州 510635)

拱壩體型單薄，節約建筑材料，又具備承載力強、經濟性好等優點，在我國水能資源豐富的西南和西北地區選用率很高。歷史上已經發生的拱壩震害實例說明拱壩的抗震性能良好，但這些實例僅限于較低的拱壩，目前還沒有高拱壩受到強震的震害實例。混凝土作為拱壩的主要建筑材料可以充分發揮其抗壓強度，但其抗拉強度低，容易造成結構剛度退化進而導致結構失穩。特別是在地震頻發的西南、西北地區，壩體在地震循環荷載作用下可能會發生較大變形進而導致結構失效，因而拱壩的抗震安全問題受到了廣泛關注[1]。

基于線性彈性力學的思想結合拉應力規范控制標準在有限元分析中常用來模擬高拱壩應力-應變等物理參數指標，且拱壩壩體可視為素混凝土結構，在合理的荷載控制標準內，其力學特性可視為線彈性。若將拱壩在地震中所受應力考慮為逐漸增大的拉壓循環荷載，相比于壓應力，混凝土的拉應力會先導致結構內部出現裂縫，結構剛度顯著降低并出現應變軟化現象。若采取線性彈性力學的思想來考慮高拱壩混凝土的應變軟化特性顯然是不合理的，應客觀基于實際情況選擇合適的混凝土塑性本構模型作為考慮影響大壩安全的關鍵因素，以便確定地震動慣性力作用下大壩結構的薄弱部位以及極限地震動慣性力作用下的大壩結構的失穩問題[2-4]。

張楚漢等[5,6]對當前混凝土細觀力學的發展進行綜述，介紹了細觀層次的混凝土仿真預處理技術、仿真數值模型和方法等方面的科研學術趨勢，給出了混凝土材料力學特性上有待挖掘的研究內容。賈明曉[7]對國內外宏細觀本構模型的研究進展進行了總結,比較了混凝土宏觀和細觀研究各自的優缺點,通過對比宏細觀多尺度建模和分析方法,從不同角度模擬混凝土宏觀非線性特征和細觀損傷破壞的全過程,形象地闡述了混凝土損傷破壞的本質。Nguyen等[8-10]總結諸多非均勻材料的建模及發展現狀，考慮了水泥砂漿的力學模型,提出了混凝土的多尺度破壞模型,進而提高了對混凝土結構的多尺度分析效率。Rui Faria等[11]針對地震動慣性力作用下混凝土結構呈現出的拉-壓循環荷載的特點提出針對性較強的混凝土塑性損傷本構，并在拱壩結構進行試驗應用；林皋院士[12]也曾提出對于混凝土拱壩，特別是為了對三維拱壩-無限地基體系的抗震安全性做出更為科學的評價，進行混凝土拱壩地震損傷破壞發展全過程的數值模擬是十分必需的。

本次仿真試驗基于Lee和Fenve[13]針對準脆性材料提出的彈塑性損傷模型，根據實際資料構造地基-大壩-庫水三維有限元模型，且在動力分析步之前設置靜力分析步計算拱圈層層施加的壩體初始應力場，依據結構設計抗震等級和所處壩址處的地震烈度劃分區域選取天然頻譜特性的地震波，同時應用SeismoSignal軟件處理地震波，梯度設計地震動響應工況(PGA=0.2、0.4、0.6 g，PGA為地震動峰值加速度)，以便獲取更加貼合實際情況的拱壩損傷結果。

1 基本原理

1.1 塑性損傷原理

Lee和Fenves[13]提出的塑性損傷模型(Concrete Damage Plastic，CDP)能夠模擬由準脆性材料構成的各種結構類型，并在混凝土塑性損傷方面得到廣泛的認可。其原理是采用各向同性的彈性損傷理論，結合各向同性拉伸、壓縮塑性理論來表征混凝土的非彈性行為，用非關聯多重硬化塑性和各向同性彈性損傷理論來表征材料破壞過程發生的不可逆損傷行為。

彈塑性增量理論根據應力-應變相關準則將應變率可視為彈性和塑性兩部分相加之和：

(1)

(2)

式中：E0為初始彈性模量。

為了對應混凝土出現軟化現象后呈現出的非線性，應力可表示為：

(3)

式中：d為損傷因子變量dt和dc的函數；dc為混凝土在受壓損傷引起的剛度退化；dt為混凝土受拉損傷引起的剛度退化。

混凝土受單軸循環力的情況下，d可視為混凝土損傷變量，d=0表示混凝土表現為直線形態的線彈性階段，混凝土未發生破壞；d=1表示混凝土結構失效，有影響結構完整性的裂縫出現。故可引入以下假定：

1-d=(1-stdc)(1-scdt)

(4)

與應力反向相關的剛度復原應力狀態的函數分別用st，sc表示，其可以用方程表示為：

(5)

在這個本構模型中，混凝土的剛度削弱有如下假定：

Dc=(1-dc)E0

(6)

Dt=(1-dt)E0

(7)

圖1表示該本構模型中材料在單軸受拉和受壓狀態下的應力-應變曲線。

圖1 混凝土單軸受拉與受壓下結構應力-應變軟化示意圖Fig.1 Schematic diagram of stress-strain softening of concrete structures under uniaxial tension and compression

混凝土結構在單軸受拉或是受壓狀態下結構破壞時呈現的體積膨脹現象，選用塑性流動勢函數G的雙曲函數為：

(8)

(9)

(10)

I=diag(111)

(11)

式中：s為有效應力偏量；φ為在高強度圍壓情況下p-q平面的膨脹角；σt0為混凝土結構單軸應力狀況下達到最大情況的極值；ξ為這個雙曲函數逼近漸近線時的偏心度。

通過有關塑性流動法則控制的塑性勢函數,可用來表示有效應力空間上的塑性應變：

(12)

Lubliner提出的CDP模型遵從了不同拉、壓條件下屈服強度也互不相同的準則，經過Lee和Fenves進一步優化后，采用有效應力、內部狀態變量為自變量函數表達式為：

(13)

1.2 模型損傷因子d取值

在CDP數值模型本構模型中，模型以單軸狀態的損傷演化方程擴展到三維狀態，確定單軸的損傷因子d是關鍵環節。在規范[14]闡述的混凝土應力形變等物理指標的基礎上，關聯到損傷因子d，詳盡描述混凝土在過度荷載的情況下材料應變發生軟化的情況。其中通過混凝土非線性階段的相關物理指標如應力、形變，按照以下方式得到損傷因子的具體參數：

(14)

1.3 考慮壩體-水體動力荷載下的附加質量法

在水工建筑物抗震規范規定采用動力時程法時，除了要考慮壩前靜水壓力的影響因素外，根據牛頓第二定律可得壩前的水體也會在地震動慣性力的影響下獲得除靜水壓力外的動水壓力荷載，即可將水體看成質點并乘以相應的地震加速度獲得相應的動水壓力荷載，而該部分附加的動水壓力僅以水平向方式施加在壩體表面上，故對應工況下的水深處附加質量公式如下：

(15)

式中：h為壩前水位至計算點水位的距離；H0為對應計算水頭下庫水深度；ρw為水體的密度；ma(h)為不同計算點深度處動水所需附加的質量。

2 仿真結果分析

2.1 工程主體建模及地震波選取

拉西瓦水電站所選壩型為混凝土雙曲拱壩，最大壩高250 m，壩頂高程2 460 m，拱頂寬10 m，建基面高程2 210 m。建立拱壩有限元模型，圖2為模型的壩體俯視圖，圖3為拉西瓦拱壩整體有限元模型及拱冠梁切片。

圖2 雙曲拱壩壩體有限元模型Fig.2 Simulation model of double curved arch dam

拱壩有限元建模基于大型軟件ABAQUS，壩體模型以拱圈中心線將壩體分為左半拱與右半拱，對于不同高程處的水平拱圈采取輸入關鍵點控制模型精度，同時利用曲線擬合關鍵點，再閉合控制高程處的曲線后向上掃略形成層層拱圈，有限元模型壩體最大壩高處為250 m。模型壩基以最大壩高處為計算依據，豎直向下延伸1.5倍壩高，即375m；為保證地震作用下庫水因素的影響向上游延伸取3倍的壩高，即750 m；同時對于影響因素較小的左岸、右岸以及壩后同取1.5倍的壩高，拱壩壩體-地基有限元模型的整體計算范圍如圖3所示。

圖3 據設計資料建立的整體有限元模型及壩體中部切片Fig.3 Integral finite element model and mid-section of dam based on design data

拱壩仿真模型的壩體混凝土依據設計資料取強度為C35相對應下的材料本構模型(參數見表1)，壩體混凝土的彈性模量取值為31 GPa，泊松比的取值為0.167，壩體混凝土采取均質密度2 400 kg/m3；計算模型的地基選取無質量地基模型，地基的彈性模量依據設計資料中取值29.25 GPa，泊松比取值0.25。仿真試驗進行地震動力分析時，壩體材料與地基的彈性模量均已依據規范取值為靜態彈性模量的1.3倍。此外，地基約束條件為地基四周采取法向約束，地基底部采取三向全約束。

查詢中國地震動區規劃圖，得到壩址處設計的基本烈度為Ⅶ度。水工建筑物抗震規范中依據壩體結構重要性(拉西瓦拱壩為甲類設防工程)在壩址所處地震基本烈度的基礎上，對于壩體結構的抗震等級有提高一級的規定，由此得到拱壩結構的抗震烈度為Ⅷ度。同時拱壩設計資料中在遭遇罕遇地震動下的100年基準期超越概率為2%時的基巖水平峰值加速度為0.23 g，根據抗震設計規范計算得到結構的反應譜的相關曲線參數(如圖4中的黑色圖線)，并將對應的參數曲線提交至甄選地震波的網站，選擇出如圖5所示自然頻譜特性下的三向地震波，圖4為所選地震波三方向反應譜與對應拱壩的規范反應譜的對比圖。

表1 C35混凝土對應下的材料本構參數Tab.1 Material constitutive parameters of C35 concrete corresponding

圖4 所選地震波三方向反應譜與對應拱壩的規范反應譜Fig.4 Three-directional response spectrum of seismic waves and normalized response spectrum of dam

圖5 所選地震波的加速度時程Fig.5 Acceleration time history of selected seismic waves

該條記錄到地震波持續時間為101.11 s，并按相關規范要求選取前10 s地震波數據進行拉西瓦拱壩動力時程計算。

2.2 壩體有限元模型模態校驗分析

在進行有限元動力仿真試驗之前，為保證所建模型及參數的正確可靠性，即能夠真實反應結構固有的動力特性，主要通過提取模型的模態陣型和自振頻率進行驗證。本次試驗主要通過計算壩前無水工況與以附加質量形式考慮正常蓄水位工況下結構的模態信息進行對比論證，以保證所建模型的精度。限于篇幅僅列出與無水工況下各階模態頻率對比情況，見表2。

如表2所示為正常蓄水位工況與無水工況下拱壩的前五階頻率，計算得到的拉西瓦拱壩正常蓄水位工況的模態基頻為1.377 Hz，與姚栓喜等[16]計算得到的拉西瓦拱壩在正常蓄水位下的基頻1.49 Hz較為相近，并對應其提出的拉西瓦拱壩第一階振型為反對稱振型的論斷。同時正常蓄水位工況對比無水工況下拉西瓦拱壩的各階模態一致偏低，驗證了用附加質量法考慮正常蓄水位的工況會降低了壩體模態頻率的結論，說明了建模和參數選擇的合理性，以及考慮正常蓄水位工況下附加質量方法使用的正確性。

表2 不同工況下拱壩的各階模態頻率 Hz

2.3 壩體仿真試驗設計及試驗結果分析

本次試驗設計思路旨在壩體施工期伊始到壩體整體完建，然后施加壩前靜水壓力到遭遇設計不同梯度地震峰值加速度下的壩體損傷開裂差異性質的全過程分析，實施的操作主要步驟如下：①對劃分好網格的壩體沿高程上均勻進行分層，對劃分好的層層拱圈設置單元集，靜力分析步初始時刻先殺死壩體整體單元集，考慮層層疊加拱圈時的施工應力場時激活對應層層拱圈單元集，直至頂拱拱圈激活完成后施加壩體自重；②然后加載靜水壓力荷載，計算正常蓄水位工況下的拱壩應力分布場；③最后進行地震動時程分析步計算獲取壩體損傷開裂結果。

本文計算輸入的地震動記錄為San Fernando地震的記錄，因水工建筑物抗震規范中放大系數為2.5倍的要求，梯度設置峰值加速度為0.2、0.4、0.6 g時3個地震峰值工況下的動力時程分析，計算時長為10s，步長0.02 s。

(1)壩頂中點位移結果分析。距離拱端最遠的拱壩中部是壩體結構高程差最大的壩段，除去拱向作用若能將中部壩段看作是高程差最大的懸臂梁，則壩體的位移動力特性會集中體現在該區域內且位移動力響應最大值會出現在梁端，故選取拱壩的中點為控制點，以建模的3個方向為正方向來表述拱壩對應方向上的位移特征，圖6為設計梯度峰值工況下的壩頂中點位移時程。

圖6中，0～13 s為靜力分析步，后面13 s末到23 s為地震動力分析步，前12 s靜力步為有限元生死單元技術進行壩體分層加載水平拱圈，在第12步施加完最后的頂拱單元集后激活壩體自重，在0～12 s的靜力分析步時段內，對比3個方向的位移，只有順河向和垂直向的壩體位移在第12步末發生了輕微的變化，特別是在垂直方向上，因壩體自重因素的影響在垂直負方向上有了較大的位移突變，在第13靜力步施加正常蓄水位下的靜水壓力荷載，除橫河向受兩岸巖體約束外，豎直方向和順河向均產生較大突變，尤其是在豎直方向上結構施加靜水壓力后由負向變化為正向。在進入地震動力分析步之后，在地震動慣性力作用下前期各工況下壩頂位移的差異性并不明顯，但隨著后期地震峰值的不斷增大，特別是在0.6 g的工況下，參照后面的損傷開裂分析結果可知因壩體出現較大的塑性變形，使得壩體結構震動的相對平衡位置也發生了改變。

圖6 設計梯度峰值工況下的壩頂中點位移時程Fig.6 Midpoint displacement history of dam top under design gradient peak conditions

(2)壩體應力結果分析。為更清楚表達結構地震過程中因應力分布狀態導致的損傷分布，分別提取了壩體結構的最大主應力和最小主應力的包絡圖。

工況1(峰值加速度為0.2 g)計算得到的壩體應力包絡圖(圖7)。

圖7 拱壩壩體主應力包絡圖(PGA=0.2 g)Fig.7 Envelope diagram of principal stress of dam (PGA=0.2 g)

除去壩踵部位的應力分布情況，在0.2 g峰值加速度下地震計算中，壩體上下游面第一、三主應力的大小分別為3.692和7.992 MPa，壩體上出現拉應力最大的位置為上游壩踵部位，對應后面的損傷分析中的0.2g峰值加速度下僅壩踵出現損傷的工況。

工況2(峰值加速度為0.4 g)計算得到的壩體應力包絡圖(圖8)。

圖8 拱壩壩體主應力包絡圖(PGA=0.4 g)Fig.8 Envelope diagram of principal stress of dam (PGA=0.4 g)

在0.4 g的峰值加速度下，壩體的第一、三主應力極值分別達到了3.935和10.16 MPa，第一主應力包絡圖中拉應力較大的部位為拱冠梁中部頂端以及2/3壩高的位置處，參照比色卡對應的取值，可得到，壩體對應的拉應力的值得范圍為3.1～3.5 MPa，通對比后面壩體出現的損傷部位與第一主應力拉應力較大的部位相對應。

工況3(峰值加速度為0.6g)計算得到的壩體應力包絡圖(圖9)。

圖9 拱壩壩體主應力包絡圖(PGA=0.6 g)Fig.9 Envelope diagram of principal stress of dam (PGA=0.6 g)

在0.6 g峰值加速度下，壩體的第一、三應力的極值分別達到了6.703和14.44 MPa，第一主應力包絡圖中在壩高2/3處的中間部位出現了拉應力較大的區域并有沿四周均勻分布的趨勢，對照文章后面0.6 g峰值加速度下出現的損傷區域，可判定壩體出現局部貫穿裂縫，并向四周擴展的發展趨勢。

從圖8～圖10并對照后面的損傷結果可得，第一主應力在很小的變化范圍內便決定了結構的受損程度。對照混凝土特性為抗壓不抗拉，特別是對于由素混凝土構成的拱壩而言，拱壩結構在地震循環荷載作用下產生的拉應力和壓應力均會導致結構受損，但一般是由于張拉出現損傷，先導致結構剛度降低，一定程度上降低了結構的抗壓強度，循環往復，使得結構損傷區域不斷擴展。

(3)損傷開裂結果提取。

當輸入0.2 g峰值加速度的地震動時，壩體損傷如圖10所示。巖基和壩體的膠結面屬于兩種不同材料的結合面，可視為結構的薄弱破壞面，在地震強度未超越壩體混凝土材料強度時，損傷區域也僅限于壩踵部位，壩體其余部位并未出現損傷。

圖10 0.2 g峰值工況下的壩面受損區域Fig.10 Damaged area of dam surface under 0.2 g seismic condition

當在動力分析步設計峰值達到0.4 g的地震工況時，除去壩踵部位的受損，在拱冠梁處出現了局部的受損區域，對照前面針對特征部位提取的壩體中部拱冠梁梁端的位移可知，在地震峰值為0.4 g的工況下，壩體出現損傷的區域并不顯著，壩頂位移振動在地震動慣性力作用結束后基本回歸至平衡位置。0.4 g工況下的損傷并未導致壩體結構失穩，壩體未出現宏觀的貫穿性裂縫，大壩仍能正常運行，如圖11所示。

圖11 0.4 g峰值工況下的壩面受損區域Fig.11 Damaged area of dam surface under 0.4 g seismic condition

當輸入地震峰值為0.6 g時，壩基膠結面處損傷范圍不斷加大，同時壩體上也出現了較為顯著損傷區域，特別是拱壩中部，上游面的損傷與下游面相對應，可視為壩體結構出現宏觀貫穿性質的裂縫，壩體結構受損嚴重，可視為在該地震峰值下，壩體結構部分失效，失效部分退出工作機制(見圖12)。

圖12 0.6 g峰值工況的壩體損傷圖Fig.12 Damaged area of dam surface under 0.6 g seismic condition

上述系列試驗觀測到的損傷區域直觀上只能觀察到大壩結構的表面，為了了解結構內部的損傷情況及損傷發展過程，特提取了拱壩中軸線處的拱冠梁損傷切片圖來表述損傷量的一般發展規律。在圖13(a)中，0.2 g峰值加速度的地震工況下，損傷區域僅局限于壩踵部位，壩體上部結構未出現明顯損傷；圖13(b)中，0.4 g的工況下壩體的拱冠處出現了少量的塑性損傷區域；圖13(c)中，0.6 g峰值加速度的工況下，在拱壩壩高的2/3處下游面出現了部分的塑性損傷區域，并有向上游面擴展的趨勢。

圖13 拱壩中軸線處拱冠梁損傷切片圖Fig.13 Slice diagram of arch crown beam damage at the central axis of arch dam

參考文獻[17]中的拉西瓦拱壩物理模型動力試驗，其裂縫分布和開裂狀態如圖14所示。將本文計算結果與該動力模型試驗的破壞形態進行對比可知，拱壩壩頂中部是拱壩抗震的薄弱部位，拱壩的上下游貫穿裂縫最初在該部位產生，隨后壩體發生應力重分布，在壩體上出現了大致與拱座平行的裂縫，最后裂縫趨于連通，壩頂中部被裂縫包圍的混凝土塊有脫離大壩主體的趨勢。但本文壩踵部位開裂損傷的情況在該動力模型上并沒有出現，原因是物理模型試驗著重考慮的是壩體結構的損傷開裂，同時物理模型試驗未考慮壩體與基巖接觸面屬于兩種材料的膠結面，結構上屬于薄弱易損區域，使得物理模型試驗與有限元試驗計算結果存在一定差距。

圖14 拉西瓦拱壩動力模型試驗破壞形態Fig.14 Failure patterns of arch dam by dynamic model test

3 結 論

基于有限元軟件建立耦聯體系壩體地基有限元模型，考慮壩前水體的動水壓力，并進行自振特性分析，檢驗模型建立的正確性以及參數選擇的合理性；同時針對拱壩結構素混凝土的結構特性確定相應塑性本構模型，并依據結構設計抗震等級和所屬抗震烈度區選擇相應自然地震波，然后設計合理的梯度地震工況，分析大壩結構損傷開裂及失穩問題：

(1)計算得到的拉西瓦拱壩正常蓄水位工況的主要模態頻率為1.377和1.587 5 Hz，相較于參考文獻[16]的工作基頻差異性不大，同時正常蓄水位工況對比與無水工況下拉西瓦拱壩的各階模態一致偏低，驗證了用附加質量法考慮正常蓄水位的工況會降低了壩體模態頻率的結論，說明了建模和參數選擇的合理性，以及考慮正常蓄水位工況下附加質量方法使用的正確性，為后續不同梯度峰值工況下的拱壩地震動仿真試驗提供先決條件。

(2)通過提取的應力結果比較分析可知，在0.2、0.4、0.6 g峰值加速度下，第一主應力極值分別為3.692、3.935、6.703 MPa，變化的幅值范圍較小，相較于第三主應力極值為7.992、10.16、14.44 MPa，第一主應力的變化范圍較小，對應了素混凝土結構抗壓不抗拉的特性，且壩體損傷區域對應了拉應力較大的區域，驗證了試驗地震動力分析的正確性。

(3)本次拉西瓦仿真試驗中在設計的不同峰值工況下，呈現出損傷開裂各有異同。首先在峰值加速度為0.2 g的工況下，大壩與基巖的膠結面最先出現損傷，作為兩種材料的膠結面屬于結構的薄弱面，故損傷區域最先出現在該區域；在0.4 g地震工況下，壩體開始出現損傷區域，但損傷區域面積較小；當地震峰值加速度達到0.6 g時，首先壩體的下游出現了大面積的塑性開裂區域，對應壩體的上游也出現了損傷，參照對應的拱冠梁切片可發現壩體出現了貫穿性質的裂縫，理論上結構已經整體破壞。

(4)將拱壩損傷開裂仿真結果與拉西瓦拱壩動力模型試驗進行對比可知，拱壩壩體損傷仿真模擬結果基本上與動力模型試驗結果相似，均在壩體的拱頂中部發生較大損傷，因此可判定該區域為大壩壩體的薄弱區域。

□