基于預(yù)插黏性界面單元的Koyna重力壩強(qiáng)震破壞過程分析

徐海濱，杜修力，楊貞軍

（1．北京工業(yè)大學(xué) 城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100124；2．浙江大學(xué) 建筑工程學(xué)院，杭州 310058）

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

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

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

本文基于黏結(jié)裂縫模型的模擬方法，在通用有限元軟件ABAQUS平臺(tái)上，將文獻(xiàn)［11］中所開發(fā)的二維算法和程序，擴(kuò)展運(yùn)用到混凝土重力壩在強(qiáng)震作用下的斷裂破壞研究中，以Koyna重力壩為研究對(duì)象，對(duì)混凝土重力壩的斷裂破壞過程和宏觀力學(xué)性能進(jìn)行研究，探索其破壞機(jī)制。

1 黏性裂縫模型

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

式中：t0為黏結(jié)強(qiáng)度，δsep為開口位移，G為斷裂能

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

1．1 彈性階段

單元頂面和底面受到的應(yīng)力與其相對(duì)位移滿足線性關(guān)系：

式中：t為兩個(gè)方向上的應(yīng)力分量，tn裂縫面法向應(yīng)力，ts為裂縫面內(nèi)的剪應(yīng)力。δ為對(duì)應(yīng)的位移分量。K為單元的剛度矩陣，在通常的數(shù)值模擬中，一般取對(duì)角線外數(shù)值為零。

圖1 Cohesive單元的應(yīng)力－相對(duì)位移關(guān)系［11］Fig．1 The stress-relative displacement relationship of cohesive element

1．2 軟化階段

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

〈〉是Macaulay括號(hào)：

如圖1中的線性軟化準(zhǔn)則，損傷指標(biāo)如下：

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

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

應(yīng)力為：

本文采用名義應(yīng)力平方準(zhǔn)則為起裂準(zhǔn)則，其表達(dá)式如下：

其中：tn0、ts0分別為各個(gè)方向單獨(dú)作用的起裂應(yīng)力。

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

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

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

3 算例驗(yàn)證

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

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

圖4 試驗(yàn)測(cè)得的沖擊荷載曲線Fig．4 Load history used for the impact test

由于該三點(diǎn)加載梁結(jié)構(gòu)和邊界條件關(guān)于中軸對(duì)稱，荷載位于中軸線上，因此，裂縫將沿著試塊的中軸線展開?；谝陨峡紤]，本次模擬中，只在裂縫開展路徑上嵌入裂縫單元。圖6為荷載加載完畢時(shí)，有限元網(wǎng)格的變形圖。其中紅色單元為已經(jīng)開裂的裂縫單元，損傷指標(biāo)（SDEG）均大于0．99。

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

圖6 三點(diǎn)加載梁受沖擊荷載后斷裂形態(tài)Fig．6 Deformed mesh for the impact test

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

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

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

4 Koyna重力壩強(qiáng)震破壞分析

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

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

圖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

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

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

圖13 Koyna重力壩的斷裂擴(kuò)展過程（細(xì)網(wǎng)格）Fig．13 The fracture extension process of Koyna gravity dam（fine mesh）

圖14 壩頂水平位移響應(yīng)Fig．14 The horizontal displacement response of dam crest

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

圖15 壩頂豎向位移響應(yīng)Fig．15 The vertical displacement response of dam crest

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

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

5 結(jié) 論

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

［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］方修君，金峰，王進(jìn)廷．基于擴(kuò)展有限元法的Koyna重力壩地震開裂過程模擬［J］．清華大學(xué)學(xué)報(bào)（自然科學(xué)版），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］杜效鵠，段云嶺，王光綸．重力壩斷裂數(shù)值分析研究［J］．水利學(xué)報(bào)，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，F(xiàn)aria 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，F(xiàn)inite 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，F(xiàn)enves G L．A plastic-damage concrete model for earthquake analysis of dams［J］．Earthquake Engineering＆Structural Dynamics．1998，27（9）：937－956．

［19］張社榮，王高輝，孫博，等．基于擴(kuò)展有限元法的重力壩強(qiáng)震潛在失效模式分析［J］．水利學(xué)報(bào)，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．