基于統(tǒng)一相場(chǎng)理論的多邊形有限元模擬研究

摘要：固體材料裂縫擴(kuò)展一直是工程中最普遍的破壞方式之一。統(tǒng)一相場(chǎng)理論在模擬裂縫擴(kuò)展方面具有一定優(yōu)勢(shì)，但傳統(tǒng)有限元方法在相場(chǎng)復(fù)雜區(qū)域離散時(shí)存在較大困難，為此將多邊形有限元方法與統(tǒng)一相場(chǎng)損傷模型相結(jié)合，基于Matlab軟件開(kāi)發(fā)出多邊形有限元統(tǒng)一相場(chǎng)損傷模型。該模型采用多邊形離散方式并通過(guò)子問(wèn)題交錯(cuò)迭代算法進(jìn)行求解，即通過(guò)固定相場(chǎng)求解位移場(chǎng)，然后基于位移場(chǎng)結(jié)果求解相場(chǎng)，反復(fù)循環(huán)計(jì)算直至兩者結(jié)果收斂。通過(guò)對(duì)比四邊形單元與多邊形單元離散下L形板的裂縫擴(kuò)展結(jié)果，發(fā)現(xiàn)多邊形相場(chǎng)計(jì)算結(jié)果與傳統(tǒng)有限元計(jì)算和實(shí)驗(yàn)結(jié)果基本一致，證明了此方法的可靠性。提出的統(tǒng)一相場(chǎng)多邊形有限元模擬方法有望在對(duì)工程結(jié)構(gòu)復(fù)雜區(qū)域進(jìn)行離散時(shí)發(fā)揮重要作用。

關(guān) 鍵 詞：裂縫擴(kuò)展； 固體材料破壞； 多邊形有限元； 統(tǒng)一相場(chǎng)損傷模型； 交錯(cuò)求解算法

中圖法分類號(hào)： TV642

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

DOI：10.16232/j.cnki.1001-4179.2025.02.027

0 引 言

在實(shí)際生產(chǎn)生活中，固體結(jié)構(gòu)的各種損傷可以近似看作裂縫，由裂縫所產(chǎn)生的斷裂是固體結(jié)構(gòu)最主要的失效方式之一。斷裂可分成脆性斷裂和準(zhǔn)脆性斷裂，前者指的是脆性固體材料發(fā)生的斷裂破壞，由于其屈服強(qiáng)度高于斷裂應(yīng)力，因此，脆性材料發(fā)生的斷裂又稱之為低應(yīng)力斷裂。準(zhǔn)脆性斷裂是指混凝土等準(zhǔn)脆性固體材料發(fā)生的斷裂破壞，混凝土等準(zhǔn)脆性材料是工程中普遍存在的，而混凝土又是土木工程中用量最大、應(yīng)用范圍最廣的材料，因此，對(duì)現(xiàn)有固體材料結(jié)構(gòu)破壞理論、數(shù)值模型的創(chuàng)新研究十分必要。

早年間，國(guó)內(nèi)外專家學(xué)者針對(duì)斷裂力學(xué)與損傷力學(xué)領(lǐng)域做了大量的研究。Griffith^［1］提出了有關(guān)陶瓷、玻璃等脆性固體材料的裂縫擴(kuò)展研究方法。Irwin^［2］隨后對(duì)應(yīng)力強(qiáng)度因子準(zhǔn)則進(jìn)行了相關(guān)研究，提出了能量準(zhǔn)則與應(yīng)力強(qiáng)度因子準(zhǔn)則的關(guān)系。在1959年，Barenblatt^［3］與Dugdale^［4］研究出了黏聚裂縫模型，此模型將斷裂力學(xué)研究方法由線性領(lǐng)域擴(kuò)展至非線性領(lǐng)域。后來(lái)，Willis^［5］將Griffith的能量準(zhǔn)則和黏聚裂縫模型兩者的關(guān)系進(jìn)行了關(guān)聯(lián)性討論；隨后，Rice^［6］提出了J積分的概念用來(lái)表征裂縫尖端區(qū)域的應(yīng)變場(chǎng)強(qiáng)度。另一方面，相關(guān)學(xué)者引入了固體材料裂縫帶理論^［7］，使之能夠處理含應(yīng)變軟化階段的本構(gòu)關(guān)系，后來(lái)研究人員又分別發(fā)展出固體損傷斷裂分析的連續(xù)方法^［8］和非連續(xù)方法^［9］。

相場(chǎng)損傷斷裂模型所參考的是Francfort等^［10］提出的線彈性固體能量變分原理，即固體材料產(chǎn)生的實(shí)際裂縫能夠使得整體系統(tǒng)的總能量達(dá)到最小值。2010年Miehe等^［11］根據(jù)熱力學(xué)方法推導(dǎo)出脆性固體材料破壞相場(chǎng)損傷模型。為了使相場(chǎng)模型適用于準(zhǔn)脆性材料，2017年Wu將斷裂力學(xué)與損傷力學(xué)結(jié)合，引入能量衰減函數(shù)與裂縫幾何函數(shù)，并根據(jù)熱力學(xué)的基本方法建立了適合脆性固體材料與混凝土等準(zhǔn)脆性固體材料的統(tǒng)一相場(chǎng)損傷理論^［12］，應(yīng)用于各類固體材料^［13-16］和結(jié)構(gòu)損傷分析^［17-21］。近年來(lái)，針對(duì)模擬結(jié)構(gòu)裂縫擴(kuò)展，比例邊界有限元法^［22］和無(wú)網(wǎng)格法^［23］也具有一定優(yōu)勢(shì)，比如前者對(duì)裂縫應(yīng)用多邊形網(wǎng)格重剖分方法，也能夠?qū)崿F(xiàn)裂縫擴(kuò)展的自動(dòng)模擬，后者則主要體現(xiàn)在建模和計(jì)算效率上的優(yōu)勢(shì)。但是，比例邊界有限元法需要在裂縫處進(jìn)行網(wǎng)格重構(gòu)，無(wú)網(wǎng)格法則需要對(duì)裂縫邊界進(jìn)行處理，而統(tǒng)一相場(chǎng)損傷模型不需要進(jìn)行網(wǎng)格重構(gòu)以及對(duì)裂縫邊界的處理，因此，統(tǒng)一相場(chǎng)損傷模型具有一定優(yōu)勢(shì)。

基于統(tǒng)一相場(chǎng)損傷模型的有限元方法對(duì)二維問(wèn)題進(jìn)行分析時(shí)，采用的是三角形與四邊形單元作為基本分析單元，但三角形在有限元分析中是常應(yīng)變單元，同時(shí)，作為有限元分析中常用的四邊形單元在計(jì)算復(fù)雜幾何區(qū)域時(shí)非常復(fù)雜。近年來(lái)，多邊形有限元離散方法受到學(xué)者的重視，Wachspress^［24］首先推導(dǎo)出多項(xiàng)式插值函數(shù)用于凸多邊形單元的形函數(shù)；隨后，F(xiàn)loater^［25］提出能夠計(jì)算任意多邊形插值函數(shù)的平均坐標(biāo)法。2004年Sukumar^［26］構(gòu)造了可以廣泛應(yīng)用多邊形單元的插值形函數(shù)。Biabanaki等^［27］提出了任意界面建模多邊形有限元的技術(shù)。隨后Khoei等^［28］又利用多邊形有限元解決了非共性的網(wǎng)格接觸碰撞問(wèn)題。2015年Khoei^［29］對(duì)滑動(dòng)接觸采用多種多邊形插值形函數(shù)進(jìn)行研究。2018年徐強(qiáng)等^［30］詳細(xì)介紹了多邊形單元的插值形函數(shù)基本的性質(zhì)，且對(duì)常用的三種插值形函數(shù)進(jìn)行具體對(duì)比分析。

由于多邊形有限元具有靈活離散的特性，這使得它在解決裂縫擴(kuò)展、損傷、斷裂等方面更加便捷精確，而統(tǒng)一相場(chǎng)損傷模型在復(fù)雜幾何區(qū)域離散時(shí)較為困難。因此，本文將多邊形有限元與統(tǒng)一相場(chǎng)損傷模型結(jié)合，并通過(guò)Matlab軟件開(kāi)發(fā)出多邊形有限元統(tǒng)一相場(chǎng)損傷程序，該程序基于多邊形有限元離散，采用子問(wèn)題交錯(cuò)迭代算法進(jìn)行求解。

1 統(tǒng)一相場(chǎng)損傷理論基礎(chǔ)

相場(chǎng)損傷理論從Griffith的能量斷裂理論和熱力學(xué)的基本理論出發(fā)，通過(guò)引入一個(gè)標(biāo)量自動(dòng)模擬裂縫擴(kuò)展，使得計(jì)算和模擬得到簡(jiǎn)化。Miehe等^［31］用該理論成功計(jì)算出了脆性固體材料的裂縫擴(kuò)展過(guò)程。

相場(chǎng)損傷理論雖然可以計(jì)算脆性固體材料，但是在混凝土等準(zhǔn)脆性固體材料上無(wú)法計(jì)算。因此，Wu^［12］整合了內(nèi)聚裂縫模型與相場(chǎng)損傷理論，推導(dǎo)出適用于混凝土等準(zhǔn)脆性材料的特征函數(shù)的表達(dá)形式，成功提出適用于混凝土等準(zhǔn)脆性固體材料的統(tǒng)一相場(chǎng)損傷理論。

1.1 特征函數(shù)的規(guī)范化

統(tǒng)一相場(chǎng)損傷理論基于傳統(tǒng)相場(chǎng)模型，進(jìn)行規(guī)范化后得出同時(shí)適用于脆性材料和準(zhǔn)脆性材料的特征函數(shù)一般表達(dá)形式，其中，裂縫幾何函數(shù)和能量衰減函數(shù)分別如下：

α（d）=ξ1d+ξ2d2+…+ξmdm（1）

ω（d）=11+φ（d）=（1-d）p（1-d）p+Q（d）（2）

其中，

φ（d）=Q（d）（1-d）p

式中：Q（d）為連續(xù)函數(shù)，d為引入的裂縫相場(chǎng)，為了滿足α（d）∈［0，1］且單調(diào)遞增，則參數(shù)ξ∈［0，2］。

通過(guò)與傳統(tǒng)脆性相場(chǎng)模型^［32］的特征函數(shù)表達(dá)形式對(duì)比發(fā)現(xiàn)，脆性相場(chǎng)模型只是統(tǒng)一相場(chǎng)損傷理論的一種特殊表達(dá)形式。在引入新的特征函數(shù)后，裂紋密度函數(shù)表達(dá)形式如下：

γ（d，SymbolQC@d）=1cα1bα（d）+bSymbolQC@d2（3）

其中，

cα=4∫10α（β）dβ

式中：SymbolQC@d表示梯度；b表示裂縫尺度bgt;0；cα代表縮放系數(shù)。

1.2 熱力學(xué)下的本構(gòu)關(guān)系

根據(jù)熱力學(xué)第一、第二定律要求，總系統(tǒng)能量耗散需滿足如下不等式：

=∫Ω（σ∶ε·-ψ·）dV≥0（4）

式中：ψ為固體自由能；σ應(yīng)力張量與ε應(yīng)變張量共軛，且考慮應(yīng)變率ε·的任意性； ∶為雙點(diǎn)積運(yùn)算符號(hào)。

材料發(fā)生開(kāi)裂后，Helmholtz自由能衰減，表達(dá)形式如下：

ψ（ε，d）=ω（d）ψ0（ε）ψ0（ε）=12ε∶E0∶ε（5）

式中：E0為彈性剛度張量。

將Helmholtz自由能函數(shù)代入式（4）中得到固體的應(yīng)力應(yīng)變關(guān)系如下：

σ=ψεY=-ψd=-ω′（d）Y（6）

式中：Y=ψ/ω，ω′（d）=ω/d分別為有效能量釋放率和能量衰減函數(shù)的導(dǎo)數(shù)。

將式（6）代入規(guī)范化后的特征函數(shù)中，為了方便表示，對(duì)得到的本構(gòu)關(guān)系進(jìn)行形式上的定義，給出統(tǒng)一相場(chǎng)模型的本構(gòu)關(guān)系如下：

σ=ω（d）σ，其中σ=ψ0ε=E0∶εY=-ω′（d）Y，其中Y=ψω=12ε∶E0∶ε（7）

1.3 應(yīng)變局部化及內(nèi)聚裂縫模型

在非彈性材料中，正則化后的應(yīng)變局部化裂縫帶寬如圖1所示，圖中D表示裂縫半帶寬。

對(duì)于式（7）給出的應(yīng)力-應(yīng)變本構(gòu)關(guān)系，非彈性應(yīng)變能表示如下：

εin=φ（d）C0∶σ（8）

式中：εin為非彈性應(yīng)變，φ（d）為單調(diào)遞增的函數(shù)，表達(dá)式見(jiàn)（2）式；C0為柔度張量。

固體材料在局部化帶內(nèi)產(chǎn)生的裂縫需要達(dá)到麥克斯韋變形協(xié)調(diào)條件，所以局部化帶內(nèi)所產(chǎn)生的裂縫跳躍位移如下：

u*=∫D-Dφ（d）C0∶σdx=2σE∫D0φ（d）dx（9）

裂縫跳躍位移與局部化帶內(nèi)應(yīng)力是一種內(nèi)聚力關(guān)系。假設(shè)足夠長(zhǎng)的桿，完全破壞時(shí)令相場(chǎng)d=d*，其左右兩端作用大小相等方向相反單調(diào)遞增的受拉位移，不考慮體積力，在開(kāi)始階段，應(yīng)變均勻分布，定義如下應(yīng)力峰值、參數(shù)和抗拉強(qiáng)度：

σc=βc2E0Gfcαb，其中βc=α′（0）φ′（0）ft=2E0Gfcαb·α′（0）φ′（0）（10）

式中：Gf為材料的斷裂能；b為裂縫尺度；α′，φ′為α函數(shù)和φ函數(shù)的導(dǎo)數(shù)。

為了更加方便說(shuō)明裂縫擴(kuò)展和準(zhǔn)脆性破壞的問(wèn)題，對(duì)于特征函數(shù)，簡(jiǎn)化形式如下：

α（d）=ξd+（1-ξ）d2ω（d）=（1-d）n（1-d）n+a1d（1+a2d+a3d）（11）

式中：ξ為裂縫幾何函數(shù)α的關(guān)鍵系數(shù)。

在完全破壞時(shí)，當(dāng)d=d*時(shí)，代入式（11），桿內(nèi)的應(yīng)力及裂縫帶兩側(cè)的非彈性位移值可以表示為d*的函數(shù)：

σ（d*）=ft［ξ+（1-ξ）d*］（1-d*）nξP（d）u*（d*）=4Gfξc0ft×∫^d*0［P（d）（1-d*）n·1ξ+（1-ξ）d*d-P（d）（1-d）n］·d·P（d）（1-d）nd（d）（12）

式中：P（d）=1+a2d+a3d2，n≥2，并同時(shí)定義正數(shù)a1，a2，a3如下：

a1=2ξcα·lchb（13）

a2=1ξ-4πξ2cα·Gfft2·k0²³+1-（n+1）（14）

a3=0，ngt;2

1ξ（cαu*ft2πGf）2-（1+a2），n=2

（15）

式中：lch為特征長(zhǎng)度，k0為初始斜率。

由式（12）的形式可知，其結(jié)果與裂縫尺度b無(wú)關(guān)，僅與材料彈性模量等相關(guān)參數(shù)有關(guān)。而確定了固體材料的基礎(chǔ)屬性，a1也就根據(jù)材料屬性確定，至于其他參數(shù)，則與局部損傷帶內(nèi)的軟化法則相關(guān)。

對(duì)于常見(jiàn)模型軟化有線性軟化、指數(shù)軟化、雙曲軟化等法則，而針對(duì)混凝土等準(zhǔn)脆性固體，則采用Cornelissen^［33］軟化法則更為準(zhǔn)確。

1.4 統(tǒng)一相場(chǎng)理論控制方程

由式（4）和（6）可以將能量耗散不等式簡(jiǎn)化如下：

=∫BYd·dV≥0（16）

由于裂縫相場(chǎng)擴(kuò)散的不可逆，即：

d·≥0（17）

所以裂縫相場(chǎng)能量耗散準(zhǔn)則滿足如下關(guān)系：

δ=∫BYδddV≤GfδAS=∫BGfδγdV（18）

式中：B為裂縫帶，AS為裂縫面積，γ為裂縫面積密度函數(shù)，δ為變分符號(hào)。

利用高斯-斯托克斯散度定理可給出如下裂縫相場(chǎng)擴(kuò)散能量準(zhǔn)則：

g（Y，d）=Y-Gfδd≤0（19）

裂縫面密度函數(shù)的變分導(dǎo)數(shù)δγ表示為

δγ=1cα1bα′（d）-2bΔd（20）

邊界條件如下：

q·nB≥0q=2bcαGfSymbolQC@d（21）

式中：q為通量，SymbolQC@d為相場(chǎng)梯度，兩者滿足線性本構(gòu)關(guān)系，nB為單位法向量。

根據(jù)式（19）和邊界條件（21）可以得到相場(chǎng)通量q與相場(chǎng)源Q（ε，d）的平衡關(guān)系如下：

Q（ε，d）=-ω′（d）Y-1cαbGfα′（d）（22）

位移場(chǎng)u也滿足經(jīng)典力平衡方程：

SymbolQC@σ+b*=0σ·n=t*（23）

式中：b*為體積均布力，t*為單位面積的面力。

可以看出，上述偏微分方程構(gòu)成了統(tǒng)一相場(chǎng)損傷理論的控制方程，其中位移場(chǎng)-相場(chǎng)問(wèn)題與系統(tǒng)總能量最小原理等效。

（u，d）=arg{min（（u，d））}（24）

系統(tǒng)總能量（u，d）表示如下：

（u，d）=∫Ωψ（ε（u，d））dV+∫BGfγ（d，SymbolQC@d）dA-（∫Ωb*·udV+∫Ωtt*·udA）（25）

其中系統(tǒng)總能量的三部分分別表示應(yīng)變能、裂縫表面能和總勢(shì)能，即在所有可能的裂縫路徑中，實(shí)際裂縫發(fā)展方向使系統(tǒng)總體能量最小，此時(shí)不需要已知或預(yù)設(shè)裂縫路徑。

2 平均值插值多邊形有限元模型

2.1 平均值插值形函數(shù)

多邊形插值形函數(shù)的一般性質(zhì)是多邊形有限元方法的基礎(chǔ)理論^［34］，平均值坐標(biāo)插值函數(shù)是基于平均值定理提出的平均值坐標(biāo)，采用凹多邊形和邊界節(jié)點(diǎn)的多邊形單元，并給出如下定義：

Nj（x）=ωj（x）ni=1ωi（x）（26）

ωj（x）=tan（αj-1/2）+tan（αj/2）‖x-xj‖（27）

式中：ωj（x）表示權(quán)函數(shù)，‖x-xj‖為點(diǎn)q與qj之間的距離，角度為αj，如圖2所示。

通過(guò)式（26）～（27）所定義的坐標(biāo)以及多邊形頂點(diǎn)可以定義一種多邊形單元坐標(biāo)格式。根據(jù)相關(guān)理論，平均值插值在多邊形有限元的邊界上是線性的，當(dāng)q處于多邊形邊界qjqj+1上時(shí)，0lt;αilt;π，i≠j，αj=π，所以tan（αj/2）=∞，形函數(shù)如下所示：

Nj（x）=limαj→πtan（αj-1/2）+tan（αj/2）‖xj-x‖nitan（αj-1/2）+tan（αj/2）‖xi-x‖=‖xj+1-x‖‖xj-x‖+‖xj+1-x‖（28）

同理，

Nj+1（x）=‖xj-x‖‖xj-x‖+‖xj+1-x‖（29）

Nj（x）=0，i≠j，j+1（30）

式（28）～（30）表明平均值插值在多邊形的問(wèn)題上為邊界上兩個(gè)頂點(diǎn)值的線性組合。

插值函數(shù)僅僅和邊界上兩個(gè)節(jié)點(diǎn)相關(guān)，可以保證插值函數(shù)的連續(xù)性，并且有利于其邊界條件計(jì)算，基于重心坐標(biāo)法對(duì)形函數(shù)可做如下推導(dǎo)：

tan（αj/2）=1-cosαjsinαj=rj（q）rj+1（q）-Dj-1（q）2Aj（q）（31）

tan（αj-1/2）=rj-1（q）rj（q）-Dj-1（q）2Aj-1（q）（32）

ωj（x）=rj-1（q）-Dj-1（q）/rj（q）2Aj-1（q）+

rj+1（q）-Dj（q）/rj（q）2Aj（q）（33）

式中：Dj=〈（qj-q），（qj+1-q）〉，rj（q）=‖qj-q‖，Aj（q）是多邊形單元節(jié)點(diǎn)q和其相鄰的頂點(diǎn)即qj，qj+1組成的三角形面積。

2.2 多邊形有限元積分

多邊形有限元需要使用下述的數(shù)值積分方法來(lái)對(duì)剛度矩陣進(jìn)行計(jì)算。本文采取Hammer積分和等參子三角法。

（1） 坐標(biāo)變換雅可比矩陣為

J=xξ yξxη yη（34）

式中：ξ，η為參數(shù)坐標(biāo)系。

xξ=Niξxi，yξ=Niξyiyη=Niηyi，xη=Niηxi（35）

（2） 標(biāo)準(zhǔn)單元和實(shí)際單元應(yīng)變矩陣B的對(duì)應(yīng)關(guān)系為

NixNiy=［J］^-1NiξNiη（36）

（3） 計(jì)算單元?jiǎng)偠染仃?/p>

［K］e=∫1-1∫1-1［B］T［D］［B］tJdξdη（37）

在三角形單元，Hammer積分法面積坐標(biāo)為λi=Δi/Δ（i=1，2，3），它表示為中心點(diǎn)到對(duì)邊的兩點(diǎn)所構(gòu)成的三角形與整個(gè)三角形面積比值，公式如下：

∫10∫^1-λi0f（λ1，λ2，λ3）dλ2dλ1=nk=1f（λk1，λk2，λk3）ωk（38）

式中：［D］為彈性矩陣；n為三角形單元積分點(diǎn)個(gè)數(shù)；λk1，λk2，λk3指面積坐標(biāo)的Hammer求積點(diǎn)；ωk為權(quán)值。本文是1個(gè)積分點(diǎn)，所用權(quán)值為1，對(duì)應(yīng)的坐標(biāo)是（1/3，1/3，1/3）。

3 子問(wèn)題交錯(cuò)迭代算法

對(duì)于非線性方程組求解，較為常用的求解方法有牛頓-拉普森迭代法，對(duì)于非線性所解出的非線性方程組，位移寫(xiě)成平衡方程的形式：

K（u）u=F（39）

式中：K（u）為切向剛度陣，u為位移向量，F(xiàn)為所施加的荷載向量。

平衡迭代方程如下：

ui+1=ui［Ki］^-1（F-Fi）（40）

將式（39）和（40）進(jìn)行牛頓-拉普森迭代線性化，得到迭代平衡方程如下：

Kuu KudKdu KddΔuΔd=rurd（41）

其中：Kuu為位移互相耦合的子剛度陣，Kud為位移和相場(chǎng)耦合的子剛度陣，Kdd為相場(chǎng)互相耦合的子剛度陣；ru和rd分別為位移場(chǎng)的荷載分量和相場(chǎng)的荷載分量。

采用牛頓整體迭代法計(jì)算時(shí)，收斂性很差，基于此，對(duì)上式忽略對(duì)稱項(xiàng)，得到如下方程：

Kuu00KddΔuΔd=rurd（42）

子問(wèn)題交錯(cuò)迭代算法的思路是：首先初始化相場(chǎng)和位移場(chǎng)，固定相場(chǎng)變量，相場(chǎng)損傷變量d對(duì)于位移場(chǎng)來(lái)說(shuō)相當(dāng)于已知，此時(shí)式（41）變成式（42），接下來(lái)代入（43）式來(lái)解決位移場(chǎng)子問(wèn)題

Kuuδu=ru（43）

然后，再基于上一步所解出的位移場(chǎng)u的結(jié)果，通過(guò)式（42）求解相場(chǎng)子問(wèn)題。依此類推，反復(fù)循環(huán)計(jì)算直至兩者收斂，再進(jìn)入下一增量步的計(jì)算。

與整體牛頓迭代算法相比，本文算法具有良好的數(shù)值穩(wěn)定性，但在某些復(fù)雜模型上，增量步需要幾千次迭代，收斂性較慢。

4 數(shù)值算例及多邊形相場(chǎng)驗(yàn)證

4.1 脆性破壞算例

一個(gè)左側(cè)帶寬度為0.5 mm缺口的正方形板如圖3所示，在板的頂部施加豎直向上的位移荷載。

計(jì)算所用參數(shù)如下：彈性模量E=210 GPa，泊松比μ=0.3，斷裂能Gf=2.7×10^-3 kN/mm，網(wǎng)格h=0.003 9 mm的均勻網(wǎng)格，裂縫加載步Δu=1×10^-5 mm，共計(jì)1 000次迭代。分別計(jì)算了正則化參數(shù)b=0.015 mm和b=0.030 mm兩種算例，圖4給出了部分加載步下的結(jié)果，其中模型頂部力-位移曲線如圖5所示。

對(duì)比文獻(xiàn)［35］中的結(jié)果可知，本文的計(jì)算結(jié)果與文獻(xiàn)中的結(jié)果在裂縫擴(kuò)展趨勢(shì)與速率方面基本一致，證實(shí)了本文計(jì)算的準(zhǔn)確性與合理性。

4.2 準(zhǔn)脆性破壞算例

4.2.1 四邊形準(zhǔn)脆性破壞算例

如圖6所示為混凝土準(zhǔn)脆性斷裂相場(chǎng)模型的示意圖。除了脆性模型中b和Gf參數(shù)外，在混凝土相場(chǎng)模型中需要增加混凝土的抗拉強(qiáng)度f(wàn)t。

L形板實(shí)驗(yàn)最早由Winkler^［36］所做。為了與其結(jié)果保持相對(duì)一致，邊界條件如下：模型底部采用固定端，并在右側(cè)距離邊界30 mm處施加向上的位移荷載。

計(jì)算所用參數(shù)如下：彈性模量取E=2.0×10⁴ MPa，泊松比μ=0.18，混凝土材料抗拉強(qiáng)度f(wàn)t=2.5 MPa，斷裂能取Gf=130 J/m²，采用h=0.25 mm的均勻網(wǎng)格，分別計(jì)算相場(chǎng)參數(shù)b=2.5，5，10，20 mm的4種算例，計(jì)算基于ABAQUS/Standard的umat接口。

通過(guò)對(duì)比實(shí)驗(yàn)結(jié)果可知（圖7～9），本文的計(jì)算模擬結(jié)果與其裂縫擴(kuò)展趨勢(shì)等方面的規(guī)律是一致的，證實(shí)了混凝土準(zhǔn)脆性材料模型的計(jì)算結(jié)果的合理性。

4.2.2 多邊形準(zhǔn)脆性破壞算例

多邊形有限元相場(chǎng)分析采用4.2.1節(jié)所示模型，計(jì)算參數(shù)同上。為了對(duì)比分析，多邊形和四邊形都采用均勻網(wǎng)格，計(jì)算采用Matlab軟件開(kāi)發(fā)出的多邊形有限元統(tǒng)一相場(chǎng)損傷程序，該程序基于多邊形有限元離散，使用子問(wèn)題交錯(cuò)迭代算法進(jìn)行求解。

通過(guò)對(duì)比相同時(shí)間兩種離散方式的結(jié)果可以發(fā)現(xiàn)，將模型剖分成多邊形網(wǎng)格單元，進(jìn)行相場(chǎng)計(jì)算的結(jié)果與四邊形網(wǎng)格單元裂縫擴(kuò)展趨勢(shì)基本一致，證實(shí)了此方法的可靠性（圖10）。

但是兩者的計(jì)算結(jié)果存在一定差距：首先由于兩者計(jì)算所采用的軟件不同，一個(gè)是基于matlab，一個(gè)是基于Abaqus的二次開(kāi)發(fā)；其次，多邊形網(wǎng)格與四邊形網(wǎng)格大小不同，并且兩者難以完全統(tǒng)一，使得結(jié)果有細(xì)微差別。

上述計(jì)算結(jié)果表明多邊形統(tǒng)一相場(chǎng)損傷模型有望在對(duì)實(shí)際工程模型的復(fù)雜幾何區(qū)域進(jìn)行離散時(shí)發(fā)揮重要作用。

5 結(jié)論與展望

統(tǒng)一相場(chǎng)損傷模型雖然在模擬裂縫萌生、擴(kuò)展等方面有較大優(yōu)勢(shì)，但在有限元離散方面有一定的局限性。相場(chǎng)程序?qū)τ谟邢拊W(wǎng)格離散較為苛刻，即在有限元離散時(shí)一旦存在四邊形單元和三角形單元交替出現(xiàn)的情況，程序就無(wú)法計(jì)算。尤其是對(duì)設(shè)置了初始裂縫的模型，對(duì)裂縫周圍（為了結(jié)果的精確需進(jìn)行網(wǎng)格加密處理）進(jìn)行離散時(shí)經(jīng)常會(huì)出現(xiàn)兩種單元交替出現(xiàn)的情況，這給計(jì)算造成了很大的麻煩，而多邊形有限元具有靈活的離散特性，可以避免上述問(wèn)題。

本文將多邊形有限元與統(tǒng)一相場(chǎng)損傷損傷模型結(jié)合，并通過(guò)Matlab軟件開(kāi)發(fā)出多邊形有限元統(tǒng)一相場(chǎng)損傷程序。該程序基于多邊形有限元離散，采用子問(wèn)題交錯(cuò)迭代算法進(jìn)行求解，提出的計(jì)算方法有望在相場(chǎng)復(fù)雜區(qū)域離散時(shí)發(fā)揮重要作用。

但是以下工作仍然需要進(jìn)一步的優(yōu)化：

① 目前，多邊形網(wǎng)格程序僅能用于二維模型，針對(duì)復(fù)雜三維模型，多面體單元剖分相場(chǎng)損傷程序仍需進(jìn)一步開(kāi)發(fā)。

② 為保證相場(chǎng)計(jì)算的精度，有限元網(wǎng)格需足夠精細(xì)，導(dǎo)致計(jì)算量較大，所以需要進(jìn)一步發(fā)展更為高效的數(shù)值求解算法。

參考文獻(xiàn)：

［1］ GRIFFITH A A.The phenomena of rupture and flow in solids［J］.Philosophical Transactions of the Royal Society of London，1921，221（582）：163-198.

［2］ IRWIN G R.Analysis of stresses and strains near the end of a crack traversing a plate［J］.Journal of Applied Mechanics，1957，24：361-364.

［3］ BARENBLATT G I.The formation of equilibrium cracks during brittle fracture.General ideas and hypotheses.Axially-symmetric cracks［J］.Journal of Applied Mathematics and Mechanics，1959，23：622-636.

［4］ DUGDALE D.Yielding of steel sheets containing slits［J］.Journal of Mechanics and Physics of Solids，1960，8：100-109.

［5］ WILLIS J R.A comparison of the fracture criteria of Griffith and Barenblatt［J］.Journal of the Mechanics and Physics of Solids，1967，15（3）：151-162.

［6］ RICE J R.A Path independent integral and approximate analysis of strain concentration by notches and cracks［J］.Journal of Applied Mechanics，1968，35（2）：379-386.

［7］ BAZANT Z P，OH B H.Crack band theory for fracture of concrete［J］.Materiaux et Construction，1983，16：155-177.

［8］ RASHID Y R.Ultimate strength analysis of prestressed concrete pressure vessels［J］.Nuclear Engineering and Design，1968，7：334-344.

［9］ NGO D，SCORDELIS A.Finite element analysis of reinforced concrete beams［J］.Journal of the American Concrete Institute，1967，64（14）：152-163.

［10］FRANCFORT G，MARIGO J.Revisiting brittle fracture as an energy minimization problem［J］.Journal of the Mechanics and Physics of Solids，1998，46（8）：1319-1342.

［11］MIEHE C，WELSCHINGER F，HOFACKER M.Thermodynamcally consistent phase-field models of fracture：variational principles and multifield FE implementations［J］.International Journal of Numerical Methods in Engineering，2010，83：1273-1311.

［12］WU J Y.A unified phase-field theory for the mechanics of damage and quasi-brittle failure in solids［J］.Journal of the Mechanics and Physics of Solids，2017，103：72-99.

［13］FENG D C，WU J Y.Phase-field regularized cohesive zone model（CZM） and size effect of concrete［J］.Engineering Fracture Mechanics，2018，197：66-79.

［14］WANG Q，F(xiàn)ENG Y T，ZHOU W，et al.A phase-field model for mixed-mode fracture based on a unified tensile fracture criterion［J］.Computer Methods in Applied Mechanics and Engineering，2020，370：113270.

［15］MARBOEUF A，BENNANI L，BUDINGER M，et al.Electromechanical resonant ice protection systems：Numerical investigation through a phase-field mixed adhesive/brittle fracture model［J］.Engineering Fracture Mechanics，2020，230：106926.

［16］WU J Y，HUANG Y L，NGUYEN V P，et al.Crack nucleation and propagation in the phase-field cohesive zone model with application to hertzian indentation fracture［J］.International Journal of Solids and Structures，2022，241：111462.

［17］WU J Y，NGUYEN V P.A length scale insensitive phase-field damage model for brittle fracture［J］.Journal of the Mechanics and Physics of Solids，2018，119：20-42.

［18］ZHANG P，HU X，WANG X，et al.An iteration scheme for phase field model for cohesive fracture and its implementation in Abaqus［J］.Engineering Fracture Mechanics，2018，204：268-287.

［19］ZHANG P，HU X，YANG S，et al.Modelling progressive failure in multiphase materials using a phase field method［J］.Engineering Fracture Mechanics，2019，209：105-124.

［20］WANG Q，ZHOU W，F(xiàn)ENG Y T.The phase-field model with an auto-calibrated degradation function based on general softening laws for cohesive fracture［J］.Mechanics of Materials，2020，142：103282.

［21］LOEW P J，POH L H，PETERS B，et al.Accelerating fatigue simulations of a phase-field damage model for rubber［J］.Computer Methods in Applied Mechanics and Engineering，2020，370：113247.

［22］蔣新新，鐘紅，李云途，等.重力壩地震斷裂的多邊形比例邊界有限元模型研究［J］.水利學(xué)報(bào)，2024，55（1）：115-125.

［23］聶新宇，汪小芳，劉勇，等.基于無(wú)網(wǎng)格法的風(fēng)力發(fā)電機(jī)主軸強(qiáng)度仿真研究［J］.機(jī)電工程，2023，40（11）：1760-1767.

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

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

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

［27］BIABANAKI S O R，KHOEI A R.A polygonal finite element method for modeling arbitrary interfaces in large deformation problems［J］.Computational Mechanics，2012，50（1）：19-33.

［28］BIABANAKI S O R，KHOEI A R，WRIGGERS P.Polygonal finite element methods for contact impact problems on non-conformal meshes［J］.Computer Methods in Applied Mechanics and Engineering，2014，269（2）：198-221.

［29］KHOEI A R，YASBOLAGHI R，BIABANAKI S O R.A polygonal-FEM technique in modeling large sliding contact on nonconformal meshes：a study on polygonal shape functions［J］.Engineering Computations，2015，32（5）：1391-1431.

［30］徐強(qiáng)，張桂彬，陳健云，等.有限元中多邊形單元的研究及應(yīng)用［J］.人民長(zhǎng)江，2018，49（12）：77-86.

［31］MIEHE C，HOFACKER M，WELSCHINGER F.A phase field model for rate-independent crack propagation：robust algorithmic implementation based on operator splits ［J］.Computer Methods in Applied Mechanics and Engineering，2010，199：2765-2778.

［32］吳建營(yíng).固體結(jié)構(gòu)損傷破壞統(tǒng)一相場(chǎng)理論、算法和應(yīng)用［J］.力學(xué)學(xué)報(bào)，2021，53（2）：301-329.

［33］CORNELISSEN H，HORDIJK D，REINHARDT H.Experimental determination of crack softening characteristics of normalweight and lightweight concrete［J］.Heron，1986，31（2）：45-56.

［34］徐強(qiáng)，董金凱，陳健云，等.考慮接觸問(wèn)題的多邊形有限元重力壩抗震斷裂研究［J］.人民長(zhǎng)江，2023，54（8）：211-220.

［35］劉國(guó)威，李慶斌，左正.相場(chǎng)斷裂模型分步算法在ABAQUS中的實(shí)現(xiàn)［J］.巖石力學(xué)與工程學(xué)報(bào)，2016，35（5）：1019-1030.

［36］WINKLER B J.Traglastuntersuchungen von unbewehrten und bewehrten betonstrukturen auf der grundlage eines objektiven werkstoffgesetzes für beton［D］.Innsbruck：Universitt Innsbruck，2001.

（編輯：鄭 毅）

Polygon finite element simulation based on unified phase field theory

XU Qiang^1，2，WANG Shaokang²，CHEN Jianyun^1，2，WANG Mingming³，LIU Jing⁴

（1.State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，China； 2.Institute of Engineering Earthquake，Dalian University of Technology，Dalian 116024，China； 3.School of Electric Power Engineering，Kunming University of Science and Technology，Kunming 650500，China； 4.Yantai Laolan Reservoir Resettlement Service Center，Yantai 264003，China）

Abstract： Crack propagation in solid materials has always been one of the most common failures in engineering.The unified phase field theory has certain advantages in simulating crack propagation，but the traditional finite element method has great difficulties in discretizing the complex region of the phase field.Therefore，we combine the polygonal finite element method with the unified phase field damage model for the first time，and develop a polygonal finite element unified phase field damage model based on Matlab software.This model adopts the polygon discrete method and solves it by the sub-problem interleaved iterative algorithm，that is，the displacement field is solved by the fixed phase field，and then the phase field is solved based on the displacement field results，and the cycle calculation is repeated until the two results converge.By comparing the crack propagation results of L-shaped plates under the discrete quadrilateral elements and polygon elements，it is found that the calculated results of the polygon phase field are basically consistent with the traditional finite element calculations and experimental results，which proves the reliability of the method.The polygon finite element method of unified phase field theory is expected to play an important role in the discretization of complex areas of engineering structures.

Key words： crack propagation； damage of solid material； polygonal finite element； unified phase field damage model； interleaved iterative algorithm