基于預插黏性界面單元的Koyna重力壩強震破壞過程分析

徐海濱，杜修力，楊貞軍

（1．北京工業大學 城市與工程安全減災教育部重點實驗室，北京 100124；2．浙江大學 建筑工程學院，杭州 310058）

強震作用下的壩體開裂過程模擬對于大壩抗震安全分析具有重要意義。目前有兩種模擬混凝土壩體損傷開裂的數值分析方法得到了重視和快速的發展，一種是從斷裂力學的角度出發，研究大壩的斷裂響應。從宏觀角度來看混凝土是一種準脆性材料，易發生脆性斷裂，適于采用斷裂力學方法。如Skrikerud等［1］和Ayari等［2］用離散裂縫模型分析大壩的動力開裂和閉合。Bhattacharjee等［3］和 Calayir等［4］用改進的彌散裂縫模型—旋轉裂縫模型在考慮非正交裂紋的影響下，對Koyna重力壩進行了動力開裂分析。方修君等［5］對質量矩陣進行了修正，將擴展有限元法應用于動態斷裂問題的求解，對Koyna重力壩進行了地震響應分析。杜效鵠等［6］在分析中引入黏聚裂紋模型，并對重力壩模型的復合開裂過程進行了模擬。另一種為基于連續介質損傷力學的裂縫計算模型。Cevera等［7］建立了各向同性損傷模型，研究了Koyna重力壩在地震作用下的損傷破壞情況。Yazdchi等［8］用有限元與邊界元相互結合的方法分析了Koyna壩的縮尺模型在地震作用下壩體的損傷發展。

但是離散裂縫模型和彌散裂縫模型都不同程度地存在著網格敏感性。彌散裂縫模型由于缺乏不連續的位移模式，在網格邊界不與主應力方向平行時非常容易產生應力鎖死現象。擴展有限元法（XFEM）是美國西北大學Belytschko等［9］首先提出來的，XFEM引入非連續的階躍函數來表征裂縫兩側的非連續位移場，不需要網格重分，但是該方法裂紋貫穿單元需要改變裂紋單元的積分格式，且程序實現很復雜、難度較大。損傷力學的裂縫模型只是將斷裂過程視為一個損傷積累的過程，不能直接描述裂縫的開裂和擴展的過程。

近年來，黏性裂縫模型（Cohesive Crack Model，CCM）受到學者們的持續關注并發展迅速［10］，這種較成熟的模擬裂縫擴展的模型提出了在裂縫尖端存在一個過渡區（Fracture process zone），假定在這個區域中存在內聚力（法向、切向或者混合），以此來分擔集中于尖端的應力，從而避免在解有限元方程過程中的奇異問題，而且可以高效簡便地在各種數值計算方法中得以實現，如有限元和離散元等。Yang、Su等［11－13］開發了在實體有限元網格中靈活插入黏性界面單元的算法和程序，對混凝土軸心受拉試件進行了二維和三維復雜隨機斷裂的蒙特卡羅模擬，隨后又對四個典型的混凝土靜態和動態三維斷裂例子進行模擬，驗證了此方法模擬斷裂的可行性。

本文基于黏結裂縫模型的模擬方法，在通用有限元軟件ABAQUS平臺上，將文獻［11］中所開發的二維算法和程序，擴展運用到混凝土重力壩在強震作用下的斷裂破壞研究中，以Koyna重力壩為研究對象，對混凝土重力壩的斷裂破壞過程和宏觀力學性能進行研究，探索其破壞機制。

1 黏性裂縫模型

20世紀50年代末，Barenblatt［14］和 Dugdale［15］首先提出適用于金屬延性材料的黏結裂縫模型，然后Hiller-borg等［16］提出適用于混凝土等材料的虛擬裂縫模型。黏結裂縫模型和虛擬裂縫模型開啟了對斷裂過程區中的能量進行數值模擬。黏結裂縫模型在斷裂過程區內考慮：骨料咬合作用、裂縫面間的摩擦作用和材料的黏結作用。這些作用在數學模型中表示為：一個垂直裂縫面的拉應力tn和裂縫平面內的剪應力ts，而且考慮應變軟化作用，也就是說，隨著裂縫面相對法向位移δn和切向位移δs的增大，相應的應力逐漸減小。這種應力隨著裂縫面相對位移的增大而減小的規律，可用應力－相對位移曲線表示，如圖1所示：在彈性階段，裂縫并未出現，應力隨著相對位移的增大而線性增大；在軟化階段，當應力達到起裂水平，材料開始軟化，軟化曲線與坐標軸包圍的面積稱為材料的Ⅰ型斷裂能Gf和Ⅱ型斷裂GfII。需要說明的是，由于混凝土的抗壓強度遠比抗拉強度大，因此，在裂縫的法向應力－相對位移曲線（圖1）中，當應力為負，即為壓應力時，不存在應變軟化。斷裂能與黏結強度的關系如式（1）所示

式中：t0為黏結強度，δsep為開口位移，G為斷裂能

這就決定了每個裂縫單元的極限位移為

1．1 彈性階段

單元頂面和底面受到的應力與其相對位移滿足線性關系：

式中：t為兩個方向上的應力分量，tn裂縫面法向應力，ts為裂縫面內的剪應力。δ為對應的位移分量。K為單元的剛度矩陣，在通常的數值模擬中，一般取對角線外數值為零。

圖1 Cohesive單元的應力－相對位移關系［11］Fig．1 The stress-relative displacement relationship of cohesive element

1．2 軟化階段

進入軟化階段后，裂縫單元出現損傷，由于材料的損傷單元的剛度將減小。損傷變量D介于0～1之間，它是一個計算變量，后處理中用SDEG表征，其又是有效相對位移δm的函數：

〈〉是Macaulay括號：

如圖1中的線性軟化準則，損傷指標如下：

其中：δm，max是加載歷史中的最大有效相對位移。δmo和δmf分別是裂縫起裂和完全破壞時的有效相對位移。

用初始剛度knn和kss表示退化后的剛度kn和ks：

應力為：

本文采用名義應力平方準則為起裂準則，其表達式如下：

其中：tn0、ts0分別為各個方向單獨作用的起裂應力。

2 嵌入黏性界面單元 CIEs（Cohesive Inter-face Elements）

本文的模擬方法是將裂縫單元嵌入各實體單元邊界上［11］，對初始有限元網格進行處理，圖2為嵌入裂縫單元前后的有限元網格以及裂縫單元的示意圖，需要說明的是，插入的裂縫單元在幾何上厚度為零，為了示意裂縫位置，在圖中顯示為帶一定厚度的單元。

圖2 嵌入cohesive單元前后網格示意圖Fig．2 Inserting cohesive elements in the initial mesh

3 算例驗證

Du等［17］開展了承受沖擊荷載的三點混凝土加載梁。試件的幾何尺寸和邊界條件見圖3，三點加載梁在梁頂面中點位置受到沖擊荷載作用，在底面中點位置有一個12．7 mm的預制缺口。試驗中測到的沖擊荷載（見圖4）作為本次模擬的外部荷載，施加在作用點上。試件的混凝土彈模E＝34 480 MPa，泊松比 ν＝0．2，密度 ρ＝2 500 kg／m3。根據 Du等［17］的建議，本文采用如圖5所示的極數軟化曲線表示為混凝土的軟化屬性，其斷裂能為152 N／m。

圖3 三點加載梁的幾何尺寸和加載條件Fig．3 Impact test of a concrete beam

圖4 試驗測得的沖擊荷載曲線Fig．4 Load history used for the impact test

由于該三點加載梁結構和邊界條件關于中軸對稱，荷載位于中軸線上，因此，裂縫將沿著試塊的中軸線展開?；谝陨峡紤]，本次模擬中，只在裂縫開展路徑上嵌入裂縫單元。圖6為荷載加載完畢時，有限元網格的變形圖。其中紅色單元為已經開裂的裂縫單元，損傷指標（SDEG）均大于0．99。

圖5 三點加載梁采用的軟化曲線Fig．5 Exponential softening law for the impact test

圖6 三點加載梁受沖擊荷載后斷裂形態Fig．6 Deformed mesh for the impact test

圖7 三點加載梁的荷載－位移曲線比較Fig．7 Load-displacement curves for the impact test

圖7 為本次模擬所得荷載－位移曲線和試驗測得的曲線。從圖中可以看出，兩者在初始階段吻合的較好，在末端稍有差別，但偏差不大。另外，提取了裂縫尖端位置隨時間的變化曲線，其中裂縫尖端是指損傷指數（SDEG）介于0．99上下的兩裂縫單元共用結點處。圖8為模擬結果與試驗測得的裂縫尖端的發展歷程，試驗中只得到三個時間點的裂縫尖端位置，從比較可見，模擬結果與試驗數據以及蘇項庭等三維模擬結果［13］均較吻合。

圖8 裂縫尖端的位置變化比較Fig．8 Crack-tip extension history for the impact test

4 Koyna重力壩強震破壞分析

Koyna重力壩作為少數幾個在強震中破壞且有較完整記錄的重力壩之一，為混凝土重力壩動力分析常用典型大壩。Koyna大壩位于印度的Koyna河上，壩高103 m，壩頂寬度 14．8 m，壩底寬70 m，壩段厚16 m，在壩高66．5 m處，下游壩面的坡面發生突然改變。1967年該壩壩址區域遭受一次6．5級強烈地震作用，在91．75 m高的庫前水位和0．474 g水平和0．312 g的豎向加速度峰值地震作用下，壩體發生開裂，并在裂縫處發生滲漏。壩體混凝土材料彈性模量 E＝31 027 MPa，μ＝0．15，ρ＝2 463 kg／m3，動態拉伸屈服強度 σt＝2．9 MPa，抗壓強度 σc＝24．1 MPa，斷裂能 250 N／m，Rayleigh阻尼因子根據線彈性分析得到的前兩階頻率計算。

計算模型如圖9所示，取Koyna最大壩高壩段為研究對象（壩高103 m），地面實測地震波如圖10所示。本文對重力壩模型劃分了粗細兩套網格，有限元模型采用三角形單元，實體單元邊界上均嵌入黏性裂縫單元，如圖11所示。粗網格含有6 336個結點和1 858個裂縫單元，細網格含有13 707個結點和3 577個裂縫單元。

圖9 Koyna重力壩模型尺寸Fig．9 Koyna gravity dam model

圖10 Koyna地震地面加速度記錄Fig．10 Seismic acceleration records of Koyna

圖11 有限元模型及插入的cohesive單元模型Fig．11 Finite element model and inserting cohesive elements

基于黏結裂縫模型理論，在壩體實體單元邊界處插入cohesive單元，采用名義應力平方準則，得到Koyna實測地震波作用下的大壩最終破壞模式。圖12和圖13是Koyna重力壩在地面強震作用下斷裂和擴展的過程，裂縫的動態開展過程主要是靠SDEG（Scalar Dam-age Variable）來表征單元的破壞程度，SDEG是一個無量綱量，大小從0到1，開裂的裂縫單元均為損傷指標（SDEG）大于0．99。從圖12中可以看出裂紋從坡度突變處最先萌生并以一定角度斜向壩體內部擴展，主裂縫逐漸發展為兩條，一條主裂縫沿著水平向擴展，最終擴展至上游面形成貫穿性裂縫；另一條主裂縫仍舊沿著一定角度斜向上游面擴展形成貫穿性裂縫；同時兩條主裂縫周邊衍生出許多微裂紋。圖13得到的前期裂縫發展軌跡與圖12非常相似，只是斜向的裂縫沒有形成貫穿，但是發展趨勢是一樣的，這個主要是裂縫局部化導致的，粗細網格的劃分導致插入裂縫單元的密度不同，擴展路徑稍有不同，但趨勢是相似的，可見此模型具有一定的網格依賴性，但是網格依賴性不大。

圖12 Koyna重力壩的斷裂擴展過程（粗網格）Fig．12 The fracture extension process of Koyna gravity dam（coarse mesh）

圖13 Koyna重力壩的斷裂擴展過程（細網格）Fig．13 The fracture extension process of Koyna gravity dam（fine mesh）

圖14 壩頂水平位移響應Fig．14 The horizontal displacement response of dam crest

圖14 和圖15為兩套網格在地震作用下壩頂A點的水平和豎向位移響應，由于大壩在2．5 s左右時完全斷裂，所以位移時程只有2．5 s。水平向位移響應前段吻合的很好，后面稍有差別，但是趨勢一樣。豎向位移響應整體吻合的很好。由此也可以看出本方法的網格依賴性不大。

圖15 壩頂豎向位移響應Fig．15 The vertical displacement response of dam crest

圖16 為Koyna重力壩的幾種斷裂破壞形態，圖（a）為Lee等［18］提出的混凝土塑性損傷模型，并依此對Koyna大壩進行了地震反應分析，得到壩體損傷結果；圖（b）為張社榮等［19］基于擴展有限元法（XFEM）對Koyna重力壩地震破壞過程進行分析，得到的大壩開裂破壞分布；圖（c）為振動臺試驗的結果［20］；圖（d）、（e）為本文方法得到的模擬結果。通過比較可以看出，本文方法得到的最終斷裂破壞形態，其中斜向的主裂縫與圖（a）、（b）、（c）基本一致，水平向貫穿主裂縫與圖（a）描述的基本一致，說明采用cohesive單元可以較好地模擬大壩在強震作用下的失效模式。

圖16 Koyna重力壩斷裂破壞形態Fig．16 Fracture failure patterns of Koyna gravity dam

5 結 論

本文基于黏性裂縫模型理論，在壩體實體單元邊界處插入黏性界面單元，對Koyna重力壩在強震作用下的破壞過程進行數值模擬，得到的地震破壞模式與實際震害、損傷破壞、擴展有限元破壞和模型試驗基本一致，說明采用黏性界面單元可較好地描述裂縫的起裂和擴展，能直觀地模擬大壩在強震作用下動力破壞過程。本文的優勢在于不需預設開裂路徑，避免了網格重劃分，但是預插的黏性界面單元，增加了自由度的數目，計算效率是個值得關注的問題，有待于和其他方法進行比較。

［1］Skrikerud P E，Bachmann H．Discrete crack modeling for dynamically loaded，unreinforced concrete structures［J］．Earthqu．Engng．Struct．Dynam，1986，14：297－315．

［2］Ayari L M，Saouma V E．A fracture mechanics based seismic analysis of concrete gravity dams using discrete cracks［J］．Engng．Fract．Mech，1990，35：587－598．

［3］Bhattacharjee S S，Leger P．Seismic cracking and energy dissipation in concrete gravity dams［J］．Earthqu．Engng．Struct．Mech，1993，22：991－1007．

［4］Calayir Y，KaratonM．Seismic fracture analysis of concrete gravity dams including dam-reservoir interaction［J］．Comput．Struct，2005，83：1595－1606．

［5］方修君，金峰，王進廷．基于擴展有限元法的Koyna重力壩地震開裂過程模擬［J］．清華大學學報（自然科學版），2008，48（12）：2065－2069．FANG Xiu-jun，JIN Feng，WANG Jin-ting．Seismic fracture simulation of the Koyna gravity dam using an extended finite element method［J］．Joumal of Tsinghua University（Science and Technonlogy）．2008，48（12）：2065－2069．

［6］杜效鵠，段云嶺，王光綸．重力壩斷裂數值分析研究［J］．水利學報，2005，36（9）：1035－1041．DU Xiao-hu， DUAN Yun-ling， WANG Guang-lun．Numerical analysis of fracture in gravity dam［J］．Journal of Hydraulic Engineering，2005，36（9）：1035－1042．

［7］Cevera M，Oliver J，Faria R．Seismic evaluation of concrete dams via continuum damage models［J］． Earthquake Engineering and Structural Dynamics．1995，24（9）：1225－1245．

［8］Yazdehi M，Khalili N，Valliappan S．Non-linear seismic behaviour of concrete dams using coupled finite element-boundary element technique［J］．International Journal for Numerical Methods in Engineering，1999，44：101－130．

［9］Belytschko T，Black T．Elastic crack growth in finite element with minimal remeshing［J］． International Journal for Numerical Methods in Engineering，1999，45（5）：610－620．

［10］Shet C，Chandra N．Analysis of energy balance when using cohesive zone models to simulate fracture processes［J］．Journal of Engineering Materials and Technology Transaction of the ASME．2002，124，440－450．

［11］Yang Z J，Su X T，Chen J F，et al，Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials［J］．International Journal of Solids and Structures，2009，46（17）：3222－3234．

［12］Su X T，Yang Z J，Liu G H，Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials：A 3D study［J］．International Journal of Solids and Structures，2010，47（17）：2336－235．

［13］Su X T，Yang Z J，Liu G H，Finite element modelling of complex 3D static and dynamic crack propagation by embedding cohesive elements in abaqus［J］．Acta Mechanica Solida Sinica，2010，23（3）：271－282．

［14］ Barenblatt G I．The formation of equilibrium cracks during brittles fracture： general ideas and hypothesis， axially symmetric cracks［J］．Applied Mathematics and Mechanics，1959，23：622－636．

［15］Dugdale D S．Yielding of steel containing slits［J］．Journal of the Mechanics and Physics of Solids，1960，8（2）：100－104．

［16］Hillerborg A，Modeer M，Petersson P E．Analysis of crack formulation and crack growth in concrete by means of fracture mechanics and finite elements［J］．Cement and Concrete Research，1976，6：773－782．

［17］Du J，Yon J H，Hawkins N M，et al．Fracture process zone for concrete for dynamic loading［J］．ACI Materials Journal，1992，89（3）：252－258．

［18］Lee J，Fenves G L．A plastic-damage concrete model for earthquake analysis of dams［J］．Earthquake Engineering＆Structural Dynamics．1998，27（9）：937－956．

［19］張社榮，王高輝，孫博，等．基于擴展有限元法的重力壩強震潛在失效模式分析［J］．水利學報，2012，43（12）：1431－1439．ZHANG She-rong，WANG Gao-hui，SUN Bo，et al．Seismic potential failure mode analysis of concrete gravity dam based on extended finite element method［J］．Journal of Hydraulic Engineering，2012，43（12）：1431－1439．

［20］National Research Council． Earthquake engineering for concrete dams：design，performance，and resaerch needs［M］．National Academics Press，USA，1990：99－100．