考慮接觸問(wèn)題的多邊形有限元重力壩抗震斷裂研究

徐強(qiáng) 董金凱 陳健云 李靜

摘要：

抗震斷裂問(wèn)題是混凝土重力壩的重要研究?jī)?nèi)容，數(shù)值模擬則是針對(duì)此類問(wèn)題的主要研究方法，現(xiàn)有斷裂算法存在計(jì)算效率較低，大多未考慮接觸等問(wèn)題。為此，構(gòu)造了高效的平均值坐標(biāo)插值多邊形單元和接觸薄層單元，并將二者耦合后用于結(jié)構(gòu)裂縫擴(kuò)展的數(shù)值模擬。該方法將具有奇異特性以及復(fù)雜幾何區(qū)域網(wǎng)格自適應(yīng)性強(qiáng)的多邊形有限元方法與平均值坐標(biāo)插值方法、局部多邊形網(wǎng)格策略結(jié)合起來(lái)模擬裂縫擴(kuò)展；通過(guò)等參子三角形方法進(jìn)行數(shù)值積分，求解剛度矩陣；基于最大周向應(yīng)力準(zhǔn)則，采用混合斷裂強(qiáng)度因子作為起裂判據(jù)；針對(duì)裂縫面間的接觸問(wèn)題，引入張剪狀態(tài)下的黏聚力模型以及壓剪狀態(tài)下的光滑接觸模型。基于上述方法，通過(guò)Matlab編程的方式開發(fā)出一套能模擬裂縫擴(kuò)展的程序，并通過(guò)對(duì)含邊裂縫的PMMA板試驗(yàn)和含預(yù)設(shè)裂縫的四點(diǎn)剪切梁試驗(yàn)的模擬，驗(yàn)證了程序的準(zhǔn)確性與有效性。最后，以含有預(yù)設(shè)裂縫的柯伊納混凝土重力壩為例，對(duì)其在是否考慮接觸問(wèn)題下裂紋擴(kuò)展的情況進(jìn)行了對(duì)比分析，結(jié)果表明該程序模擬裂縫真實(shí)準(zhǔn)確，具有計(jì)算速度快、結(jié)果文件小等獨(dú)特優(yōu)勢(shì)。

關(guān) 鍵 詞：

混凝土重力壩； 平均值坐標(biāo)插值多邊形有限元； 裂縫擴(kuò)展； 接觸算法

中圖法分類號(hào)： TV642.3；TV313

文獻(xiàn)標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2023.08.030

0 引 言

中國(guó)水能資源十分豐富，位居世界首位。水電具有循環(huán)再生、長(zhǎng)期成本低、綜合效益大等優(yōu)點(diǎn)。重力壩是中國(guó)高壩建設(shè)的主要壩型，其多位于西南、西北等強(qiáng)震頻發(fā)的地區(qū)，而地震則是造成其損傷斷裂破壞的主要潛在危險(xiǎn)。從大壩的斷裂過(guò)程來(lái)看，結(jié)構(gòu)損傷裂縫往往從壩體表面開始，在施工階段已經(jīng)形成。在地震作用下，表面裂縫很有可能繼續(xù)發(fā)展，深入壩體。因此，大壩安全評(píng)價(jià)的首要問(wèn)題是確定已有裂縫是否會(huì)擴(kuò)展，預(yù)測(cè)在各種可能條件下的擴(kuò)展路徑以及危險(xiǎn)程度，以便在它們進(jìn)一步擴(kuò)展前采取控制措施［1］。

由于重力壩具有體積大、受力情況復(fù)雜等特點(diǎn)，現(xiàn)階段主要利用數(shù)值模擬的方法對(duì)其進(jìn)行研究。近年來(lái)，很多學(xué)者針對(duì)混凝土重力壩斷裂問(wèn)題進(jìn)行了多方面研究［2-4］。王超等［5］和Siva等［6］探討了在地震作用下初始裂縫的位置以及深度對(duì)壩體的影響。在對(duì)結(jié)構(gòu)進(jìn)行損傷斷裂數(shù)值模擬時(shí)，有限元方法作為一種重要手段，已經(jīng)廣泛應(yīng)用于工程領(lǐng)域和學(xué)術(shù)研究中。二維有限元方法中使用的基本單元形狀為三角形和四邊形。三角形單元為常應(yīng)變單元，而四邊形單元對(duì)復(fù)雜幾何區(qū)域離散比較困難，使得傳統(tǒng)有限元對(duì)于一些特定的問(wèn)題求解有明顯的缺陷。將多邊形單元用于有限元分析中能克服傳統(tǒng)有限元的某些缺陷。Wachspress［7］首次在凸多邊形單元上構(gòu)造了有理多項(xiàng)式插值函數(shù)，用于在有限元分析中計(jì)算多邊形單元的形函數(shù)；Floater［8］基于多邊形單元的重心坐標(biāo)提出平均值坐標(biāo)法，適用于計(jì)算任意凸多邊形和凹多邊形的插值函數(shù)；Sukumar［9］基于最大熵方法推導(dǎo)了任意多邊形單元更為一般的插值函數(shù)；多邊形有限元各向異性的特點(diǎn)使它在解決裂縫擴(kuò)展［10］、斷裂［11］等方面更加方便有效。Tabarraei和Sukumar［12］將自然鄰點(diǎn)形函數(shù)應(yīng)用于擴(kuò)展有限元（X-FEM）中以模擬裂縫的擴(kuò)展；Manzini等［13］研究了各種凸多邊形有限元以及凹多邊形有限元的性能，使其能用于裂尖以及裂縫的擴(kuò)展的模擬；Natarajan等［14］、Mousavi等［15］、Talischi和Paulino［16］對(duì)凸多邊形有限元以及凹多邊形有限元數(shù)值積分規(guī)則進(jìn)行了許多研究，使之可以用于模擬裂縫兩側(cè)的單元積分；Ooi等［17-18］提出用基于多邊形單元的比例邊界有限元法模擬結(jié)構(gòu)裂縫擴(kuò)展；Spring［19］和Leon［20］等人將多邊形有限元方法用于黏彈性動(dòng)態(tài)斷裂問(wèn)題中。施明光等［21］采用多邊形比例邊界有限元，結(jié)合基于拓?fù)涞木植烤W(wǎng)格重剖分方法，首次模擬了復(fù)合材料裂縫擴(kuò)展。在裂縫擴(kuò)展過(guò)程中，由于受力狀態(tài)不同，如果不考慮裂縫面上的接觸條件，則不能反映裂縫的真實(shí)開裂狀態(tài)以及導(dǎo)致有限元分析時(shí)網(wǎng)格相互嵌入的問(wèn)題。針對(duì)張剪狀態(tài)，F(xiàn)ang等［22］和伍江飛［23］對(duì)黏聚裂縫模型進(jìn)行了研究，并應(yīng)用于混凝土斷裂的模擬。針對(duì)壓剪狀態(tài)，薛冰寒等［24］提出了求解彈性摩擦接觸問(wèn)題的IGA-B可微方程組；Zhou等［25］和喻葭臨等［26］基于擴(kuò)展有限元方法，構(gòu)造了與之相適應(yīng)的接觸算法；章鵬等［27］基于比例邊界有限元，提出了一種在裂尖單元的裂縫面上采用邊對(duì)邊約束方法處理接觸問(wèn)題的方法。

由于常規(guī)單元在裂縫尖端不能反映裂尖的奇異性，而基于平均值坐標(biāo)插值的多邊形有限元具有1/4點(diǎn)的奇異性，同時(shí)平均值坐標(biāo)插值的多邊形有限元具有突出的幾何拓?fù)湫阅埽芨鼮闇?zhǔn)確地描述裂縫擴(kuò)展路徑。因而，平均值坐標(biāo)插值多邊形有限元常被用于模擬裂尖擴(kuò)展。但是現(xiàn)階段平均值坐標(biāo)插值多邊形有限元模擬裂尖仍存在以下問(wèn)題：平均值坐標(biāo)插值多邊形函數(shù)復(fù)雜，應(yīng)變矩陣計(jì)算復(fù)雜度高，剛度矩陣合成效率低；平均值坐標(biāo)插值多邊形單元，在模擬裂縫擴(kuò)展時(shí)，大多沒(méi)有考慮裂縫閉合的接觸問(wèn)題；裂縫擴(kuò)展重剖分算法一般都比較復(fù)雜，計(jì)算效率較低。基于以上問(wèn)題，本文采用等參元方法解決剛度矩陣合成效率低的問(wèn)題；開發(fā)接觸薄層單元并耦合平均值坐標(biāo)插值多邊形單元解決裂縫閉合的接觸問(wèn)題；使用一次裂縫跨越一個(gè)單元解決重剖分的幾何復(fù)雜度的問(wèn)題；最后通過(guò)Matlab編程實(shí)現(xiàn)了上述算法過(guò)程。

1 平均值坐標(biāo)插值多邊形有限元模型

1.1 多邊形單元插值形函數(shù)的一般性質(zhì)

4 數(shù)值算例驗(yàn)證

為驗(yàn)證上述裂縫擴(kuò)展程序的合理性，本文選用文獻(xiàn)［33］中所述的數(shù)值模擬中的兩個(gè)算例進(jìn)行對(duì)比分析，論證程序的正確性。

第一個(gè)算例為PMMA邊緣裂縫板試驗(yàn)，其幾何尺寸、邊界條件、荷載條件如圖7所示。材料彈性模量為3.86 GPa，泊松比為0.31，厚度為1 m，斷裂韌度為1 MPa·m1/2。尺寸a=4 mm，c=5 mm，w=30 mm，h=45 mm。分別改變裂縫位置d=6，10 mm和14 mm，研究裂縫位置對(duì)裂縫路徑的影響。假定為平面應(yīng)變條件。本文所述程序模擬裂縫結(jié)果與文獻(xiàn)［33］模擬裂縫結(jié)果對(duì)比情況如圖8所示。

第二個(gè)算例為含預(yù)制裂縫的四點(diǎn)剪切梁試驗(yàn)，該試驗(yàn)已成為驗(yàn)證混合型裂縫擴(kuò)展模型的基準(zhǔn)。材料彈性模量為24.8 GPa，泊松比為0.18，厚度為0.152 m，斷裂韌度為1.65 MPa·m1/2。假定為平面應(yīng)力狀態(tài)。具體尺寸、荷載條件等情況如圖9所示。應(yīng)用本程序計(jì)算的裂縫路徑及擴(kuò)展情況與文獻(xiàn)［33］所述情況對(duì)比如圖10所示。

經(jīng)過(guò)兩個(gè)算例的對(duì)比分析可見，兩個(gè)算例的裂縫路徑以及開裂情況均與本文所用程序計(jì)算結(jié)果相吻合。因此，論證了該斷裂程序搜索斷裂路徑的正確性與合理性。

5 重力壩裂縫擴(kuò)展模擬

印度的柯伊納混凝土重力壩是壩體在地震作用下遭受到破壞的典型實(shí)例之一，可在一定程度上反映壩體的裂縫擴(kuò)展情況。國(guó)內(nèi)外很多學(xué)者對(duì)柯伊納混凝土重力壩在地震作用下的動(dòng)力響應(yīng)進(jìn)行了深入研究。因此，本文選用該混凝土重力壩，利用本文方法編制的程序，對(duì)其進(jìn)行裂縫擴(kuò)展模擬的研究。

柯伊納混凝土重力壩，壩高103 m，壩頂寬14.8 m，壩底寬70 m，壩高66.5 m處下游壩面出現(xiàn)折坡。為使計(jì)算結(jié)果具有較高的精度，采用在壩基與壩體邊緣處加密的有限元模型，如圖11所示。斷裂分析過(guò)程需要設(shè)置初始裂縫，所以斷裂過(guò)程就是裂縫擴(kuò)展過(guò)程，初始裂縫設(shè)置在下游折坡處下，裂縫長(zhǎng)度0.5 m。壩體混凝土材料的彈性模量為20 GPa，泊松比為0.167，密度為2 400 kg/m3，抗拉強(qiáng)度為2 MPa，斷裂能為200 N/m，混凝土斷裂韌度取為2.747 MPa·m1/2。

外荷載主要考慮重力、靜水壓力（靜水壓力水位為91.75 m），以及地震作用。地震采用柯伊納地震實(shí)測(cè)波，考慮峰值加速度約為0.5g的水平方向地震動(dòng)與峰值加速度約為0.3g的垂直方向地震動(dòng)同時(shí)作用（見圖12）。

5.1 基本步驟

基于平均值插值多邊形有限元方法，考慮裂縫面接觸問(wèn)題，給出了工程實(shí)例在實(shí)測(cè)地震動(dòng)作用下裂縫擴(kuò)展模擬的情況，程序具體實(shí)施過(guò)程如下：

（1） 首先利用ANSYS建立壩體有限元模型，施加荷載，約束節(jié)點(diǎn)。求解后，讀出單元屬性信息，并導(dǎo)入Matlab中，通過(guò)Matlab編程對(duì)節(jié)點(diǎn)施加靜力和動(dòng)力荷載。

（2） 根據(jù)輸入實(shí)測(cè)地震動(dòng)進(jìn)行時(shí)程計(jì)算，計(jì)算得到初始裂縫擴(kuò)展角，轉(zhuǎn)入動(dòng)力裂縫求解程序，對(duì)裂縫單元重剖分，剖分為一個(gè)五邊形單元、一個(gè)三邊形單元加一個(gè)五邊形單元，或者剖分為一個(gè)五邊形單元與兩個(gè)四邊形單元。

（3） 動(dòng)力裂縫求解過(guò)程中，考慮裂縫面間的接觸問(wèn)題，因此，利用接觸算法求解附加剛度矩陣，與整體剛度矩陣疊加后進(jìn)行迭代，計(jì)算真實(shí)開度。此算法既考慮了張開時(shí)的黏聚力影響，又防止出現(xiàn)閉合時(shí)單元相互嵌入的問(wèn)題。

（4） 通過(guò)時(shí)程計(jì)算，利用動(dòng)力裂縫求解程序，判斷裂縫擴(kuò)展方向，確定裂縫擴(kuò)展路徑，完成對(duì)裂縫擴(kuò)展情況的模擬。

5.2 裂縫擴(kuò)展模擬情況分析

由于開裂，將裂縫附近單元重新剖分為獨(dú)立單元，在地震加載過(guò)程中，單元會(huì)存在相互嵌入的問(wèn)題。因此本文采用前述接觸算法應(yīng)用于柯伊納混凝土重力壩裂縫擴(kuò)展模擬分析，成功解決了模擬過(guò)程中相互嵌入的問(wèn)題（見圖13）。

對(duì)是否考慮接觸問(wèn)題的裂縫擴(kuò)展情況進(jìn)行模擬分析，開裂過(guò)程取初始裂縫處為開度時(shí)程監(jiān)測(cè)點(diǎn)，如圖14所示，其縫面張開閉合時(shí)程曲線如圖15所示。具體開裂過(guò)程如圖16所示，最終裂縫路徑如圖17所示。可以看出，本文提出的接觸算法，消除了地震過(guò)程中縫面的負(fù)開度，既解決了縫面間的嵌入問(wèn)題，同時(shí)又提高了裂縫擴(kuò)展模擬的準(zhǔn)確性。

通過(guò)對(duì)應(yīng)力強(qiáng)度因子時(shí)程曲線（見圖18）進(jìn)行分析，可以發(fā)現(xiàn)，裂縫均于0.22 s時(shí)開裂，裂縫擴(kuò)展速度很快，均在0.5 s內(nèi)開裂完成。開裂過(guò)程為Ⅰ型和Ⅱ型耦合作用。考慮接觸情況下，開裂速度略快，等效應(yīng)力強(qiáng)度因子峰值出現(xiàn)略晚且小于不考慮接觸情況。

由于豎直方向位移響應(yīng)非常微小，故僅對(duì)水平方向位移響應(yīng)進(jìn)行分析，如圖19所示。可以發(fā)現(xiàn)，考慮接觸的情況下，大大減小了x軸正向位移，也反映出接觸算法的作用。

6 結(jié) 語(yǔ)

本文基于多邊形有限元理論，考慮裂縫接觸面模型，利用Matlab開發(fā)了裂縫擴(kuò)展模擬的程序。并通過(guò)對(duì)PPMA邊裂縫板試驗(yàn)以及含預(yù)設(shè)裂縫的四點(diǎn)剪切梁試驗(yàn)的裂縫擴(kuò)展情況進(jìn)行模擬，驗(yàn)證了該程序的正確性。隨后，利用該程序?qū)κ欠窨紤]接觸問(wèn)題情況下的柯伊納混凝土重力壩裂縫擴(kuò)展問(wèn)題進(jìn)行對(duì)比分析，研究考慮接觸問(wèn)題后如何解決單元相互嵌入的問(wèn)題。

該程序采用平均值坐標(biāo)插值多邊形有限元在關(guān)心的復(fù)雜幾何區(qū)域進(jìn)行網(wǎng)格劃分，具有網(wǎng)格自適應(yīng)性強(qiáng)的優(yōu)點(diǎn)；由于平均值坐標(biāo)插值多邊形有限元節(jié)點(diǎn)應(yīng)力值為奇異解的特性，將其應(yīng)用于壩體斷裂分析具有較高的精度；采用等參元方法使平均值坐標(biāo)插值多邊形有限元的計(jì)算效率大幅提升；采用裂縫面不同狀態(tài)下的接觸算法，很好地模擬了接觸界面上的粘結(jié)與分離，更真實(shí)地模擬了裂縫開裂情況。同時(shí)該程序還具有計(jì)算速度快、結(jié)果文件小、適用性強(qiáng)等優(yōu)點(diǎn)。

參考文獻(xiàn)：

［1］ 杜效鵠，段云嶺，王光綸.重力壩斷裂數(shù)值分析研究［J］.水利學(xué)報(bào)，2005（9）：1035-1042.

［2］ 卿龍邦，喻渴來(lái)，徐東強(qiáng).基于擴(kuò)展有限元法的混凝土重力壩宏細(xì)觀斷裂數(shù)值分析［J］.水力發(fā)電學(xué)報(bào)，2017，36（6）：94-102.

［3］ 暴艷利，鐘紅，林皋.基于多邊形比例邊界有限元的重力壩地震斷裂模擬［J］.水電能源科學(xué)，2015，33（4）：72-75，42.

［4］ 鐘紅，牟昊，張文宣.基于有限斷裂法和比例邊界有限元法的裂紋擴(kuò)展模擬［J］.計(jì)算力學(xué)學(xué)報(bào)，2017，34（2）：168-174.

［5］ 王超，張社榮，王高輝.初始裂縫對(duì)重力壩地震響應(yīng)特性的影響［J］.天津大學(xué)學(xué)報(bào)（自然科學(xué)與工程技術(shù)版），2016，49（4）：392-399.

［6］ SIVA P，MAHESH M，RAJ K D，et al.Critical crack lengths of concrete gravity dam by using fracture mechanics［J］.Materials Today：Proceedings，2021，38（5）：3149-3159.

［7］ WACHSPRESS E L.A Rational finite element basis［J］.Journal of Lubrication Technology，1975，98（4）：635.

［8］ FLOATER M S.Mean value coordinates［J］.Computer Aided Geometric Design，2003，20（1）：19-27.

［9］ SUKUMAR N.Construction of polygonal interpolants：a maximum entropy approach［J］.International Journal for Numerical Methods in Engineering，2004，61（12）：2159-2181.

［10］ 張慧華，祝晶晶.復(fù)雜裂縫問(wèn)題的多邊形數(shù)值流形方法求解［J］.固體力學(xué)學(xué)報(bào)，2013，34（1）：38-46.

［11］ BISHOP J E.Simulating the pervasive fracture of materials and structures using randomly close packed Voronoi tessellations ［J］.Computational Mechanics，2009，44（4）：455-471.

［12］ SUKUMAR N，TABARRAEI A.Conforming polygonal finite elements［J］.International Journal for Numerical Methodsin Engineering，2004，61（12）：2045-2066.

［13］ GIANMARCO M，ALESSANDRO R，SUKUMAR N.New perspectives on polygonal and polyhedral finite element methods［J］.Mathematical Models and Methods in Applied Sciences，2014，24（8）：1665-1699.

［14］ NATARAJAN S，BORDAS S，MAHAPATRA D R.Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping［J］.International Journal for Numerical Methods in Engineering，2009，80（80）：103-134.

［15］ MOUSAVI S E，XIAO H，SUKUMAR N.Generalized Gaussian quadrature rules on arbitrary polygons［J］.International Journal for Numerical Methods in Engineering，2010，82（1）：99-113.

［16］ CAMERON T，GLAUCIO H P.Addressing integration error for polygonal finite elements through polynomial projections：A patch test connection［J］.Mathematical Models and Methods in Applied Sciences，2013，24（8）：1701-1727.

［17］ OOI E T，SONG C，TIN-LOI F，et al.Automatic modelling of cohesive crack propagation in concrete using polygon scaled boundary finite elements［J］.Engineering Fracture Mechanics，2012，93：13-33.

［18］ OOI E T，SONGC，TIN-LOI F，et al.Polygon scaled boundary finite elements for crack propagation modelling［J］.International Journal for Numerical Methods in Engineering，2012，91（3）：319-342.

［19］ SPRING D W，LEON S E，PAULINO G H.Unstructured polygonal meshes with adaptive refinement for the numerical simulation of dynamic cohesive fracture［J］.International Journal of Fracture，2014，189（1）：33-57.

［20］ LEON S E，SPRING D W，PAULINO G H.Reduction in mesh bias for dynamic fracture using adaptive splitting of polygonal finite elements［J］.International Journal for Numerical Methods in Engineering，2015，100（8）：555-576.

［21］ 施明光，徐艷杰，鐘紅，等.基于多邊形比例邊界有限元的復(fù)合材料裂縫擴(kuò)展模擬［J］.工程力學(xué)，2014（7）：1-7.

［22］ FANG W H，WU J F，YU T T，et al.Simulation of cohesive crack growth by a variable-node XFEM［J］.Frontiers of Structural and Civil Engineering，2020，14（1）：215-228.

［23］ 伍江飛.基于黏聚裂縫模型的混凝土斷裂過(guò)程擴(kuò)展有限元法模擬［J］.水利水電技術(shù)，2019，50（4）：205-211.

［24］ 薛冰寒，林皋，胡志強(qiáng)等.求解摩擦接觸問(wèn)題的IGA-B可微方程組方法［J］.工程力學(xué)，2016，33（10）：35-43.

［25］ ZHOU X P，CHEN J W，BERTO F.XFEM based node scheme for the frictional contact crack problem［J］.Computers and Structures，2020，231：106221.

［26］ 喻葭臨，于玉貞，張丙印，等.基于擴(kuò)展有限元方法的界面接觸算法［J］.工程力學(xué)，2011，28（4）：13-17.

［27］ 章鵬，杜成斌，江守燕.比例邊界有限元法求解裂縫面接觸問(wèn)題［J］.力學(xué)學(xué)報(bào)，2017，49（6）：1335-1347.

［28］ 李術(shù)才，王兆清，李樹忱.基于無(wú)理函數(shù)插值的多邊形有限元方法［J］.山東大學(xué)學(xué)報(bào)（工學(xué)版），2008，38（2）：66-70.

［29］ HORMANN K.Barycentric coordinates for arbitrary polygons in the plane［R］.Clausthal-Zellerfeld：Institute of Computer Science，Clausthal University of Technology，2005.

［30］ KHOEI A R，YASBOLAGHI R，BIABANAKI S.A polygonal finite element method for modeling crack propagation with minimum remeshing［J］.International Journal of Fracture，2015，194（2）：123-148.

［31］ ERDOGA F，SIH G C.On the crack extension in plates under plane loading and transverse Shear［J］.Journal of Basic Engineering，1963，85：519-525.

［32］ 方修君，金峰，王進(jìn)廷.用擴(kuò)展有限元方法模擬混凝土的復(fù)合型開裂過(guò)程［J］.工程力學(xué)，2007（增1）：46-52.

［33］ YANG Z J.Fully automatic modelling of mixed-mode crack propagation using scaled boundary finite element method［J］.Engineering Fracture Mechanics，2006，73（12）：1711-1731.

（編輯：黃文晉）

Abstract：

Seismic fracture research is an important research content of concrete gravity dams.The numerical modeling is the main research method.Current fracture algorithms are inefficient and do not account the contact problems.In this paper，high efficiency polygonal finite element of mean value coordinate interpolation and contact thin layer element are constructed.The polygonal finite element method which is singularity and adapts well to complex geometries，is combined with the mean value coordinate interpolation method and the local polygonal mesh strategy to simulate the crack propagation.The stiffness matrix is solved by numerical integration using isoparametric sub-triangle method.Based on the maximum circumferential stress criterion，the hybrid fracture strength factor is employed as the crack initiation criteria.To address fracture surface contact，the cohesive force model in tension-shear and the smooth contact model in compression-shear are used.Based on the above methods，we use the Matlab to compile a program for simulating crack propagation.The program is validated by simulating PMMA plate tests with edge fractures and four-point shear beam tests with preset cracks.We present a comparative analysis of crack expansion with and without considering the contact problem for the Koyna concrete gravity dam in India，which has preset cracks.The conclusion is that the program has distinctive benefits such as realistic and accurate crack simulation，rapid calculating speed and small result files.

Key words：

concrete gravity dam；polygon finite element of mean value coordinate interpolation；crack propagation；contact algorithm