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

徐強 董金凱 陳健云 李靜

摘要：

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

關(guān) 鍵 詞：

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

中圖法分類號： TV642.3；TV313

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

DOI：10.16232/j.cnki.1001-4179.2023.08.030

0 引 言

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

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

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

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

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

4 數(shù)值算例驗證

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

第一個算例為PMMA邊緣裂縫板試驗，其幾何尺寸、邊界條件、荷載條件如圖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，研究裂縫位置對裂縫路徑的影響。假定為平面應(yīng)變條件。本文所述程序模擬裂縫結(jié)果與文獻(xiàn)［33］模擬裂縫結(jié)果對比情況如圖8所示。

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

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

5 重力壩裂縫擴展模擬

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

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

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

5.1 基本步驟

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

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

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

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

（4） 通過時程計算，利用動力裂縫求解程序，判斷裂縫擴展方向，確定裂縫擴展路徑，完成對裂縫擴展情況的模擬。

5.2 裂縫擴展模擬情況分析

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

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

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

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

6 結(jié) 語

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

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

參考文獻(xiàn)：

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

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

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

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

［5］ 王超，張社榮，王高輝.初始裂縫對重力壩地震響應(yīng)特性的影響［J］.天津大學(xué)學(xué)報（自然科學(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ù)雜裂縫問題的多邊形數(shù)值流形方法求解［J］.固體力學(xué)學(xué)報，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ù)合材料裂縫擴展模擬［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］ 伍江飛.基于黏聚裂縫模型的混凝土斷裂過程擴展有限元法模擬［J］.水利水電技術(shù)，2019，50（4）：205-211.

［24］ 薛冰寒，林皋，胡志強等.求解摩擦接觸問題的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］ 喻葭臨，于玉貞，張丙印，等.基于擴展有限元方法的界面接觸算法［J］.工程力學(xué)，2011，28（4）：13-17.

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

［28］ 李術(shù)才，王兆清，李樹忱.基于無理函數(shù)插值的多邊形有限元方法［J］.山東大學(xué)學(xué)報（工學(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)廷.用擴展有限元方法模擬混凝土的復(fù)合型開裂過程［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