王博臣，侯玉亮，夏涼,*，史鐵林

1.華中科技大學(xué) 數(shù)字制造裝備與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074 2.鄭州大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，鄭州 450001

斷裂是材料的主要破壞形式之一，斷裂的模擬在科學(xué)研究和工程應(yīng)用中具有十分重要的意義，同時(shí)也極具挑戰(zhàn)性。在航空領(lǐng)域，飛行器復(fù)雜的服役環(huán)境以及設(shè)計(jì)、制造過(guò)程中形成的缺陷、微裂紋導(dǎo)致結(jié)構(gòu)強(qiáng)度降低。為檢驗(yàn)結(jié)構(gòu)強(qiáng)度是否滿(mǎn)足設(shè)計(jì)性能指標(biāo)，需要對(duì)在研飛行器進(jìn)行強(qiáng)度校驗(yàn)，包括地面強(qiáng)度試驗(yàn)和斷裂數(shù)值模擬等驗(yàn)證性措施；其中，對(duì)斷裂數(shù)值模擬的有效性提出了苛刻的要求。材料斷裂的相關(guān)研究經(jīng)歷了多個(gè)階段的發(fā)展，最早可以追溯到Griffith的相關(guān)工作，他通過(guò)彈性應(yīng)變能和裂紋表面能之間的平衡關(guān)系研究裂紋的擴(kuò)展。在Griffith相關(guān)工作的基礎(chǔ)上，Irwin建立了線(xiàn)彈性斷裂力學(xué)，借助連續(xù)介質(zhì)力學(xué)的方法研究材料的斷裂機(jī)制。近幾十年來(lái)，隨著有限元方法和計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值斷裂模擬方法被廣泛應(yīng)用于斷裂研究?；谟邢拊椒ǖ臄嗔涯Ｐ鸵话憧煞譃?種：分離式裂紋模型和彌散式裂紋模型。

單元?jiǎng)h除法是一種典型的分離式裂紋模型，該方法根據(jù)單元的應(yīng)力水平或彈性應(yīng)變能水平進(jìn)行判斷，將超過(guò)材料損傷閾值的單元進(jìn)行刪除；然而，該方法難以模擬裂紋的分叉和復(fù)雜裂紋形式。界面單元法在單元之間插入內(nèi)聚力界面單元，當(dāng)界面單元的內(nèi)聚力達(dá)到斷裂條件時(shí)發(fā)生失效，裂紋在該界面單元處相應(yīng)成核或擴(kuò)展。在界面單元法中，裂紋成核與擴(kuò)展只允許發(fā)生在單元的界面處，模擬結(jié)果存在較強(qiáng)的有限元網(wǎng)格依賴(lài)性，難以模擬復(fù)雜裂紋形式。擴(kuò)展有限元法采用擴(kuò)展形函數(shù)表示裂紋兩側(cè)位移場(chǎng)的間斷，裂紋可以在自由度被擴(kuò)展的單元內(nèi)部成核與擴(kuò)展，因此一定程度上避免了網(wǎng)格依賴(lài)性問(wèn)題。由于需要顯式地追蹤裂紋尖端，擴(kuò)展有限元法較難被應(yīng)用于多裂紋擴(kuò)展、三維問(wèn)題等復(fù)雜斷裂問(wèn)題的模擬。

相場(chǎng)斷裂模型是一種彌散式裂紋模型，近年來(lái)得到了快速發(fā)展。與分離式裂紋模型相比，相場(chǎng)斷裂模型通過(guò)一個(gè)彌散分布的標(biāo)量場(chǎng)近似描述裂紋，通過(guò)求解偏微分方程模擬裂紋成核與擴(kuò)展，避免了不連續(xù)位移場(chǎng)中裂紋幾何描述和尖端追蹤的難題，極大地簡(jiǎn)化了算法實(shí)現(xiàn)的復(fù)雜度，適用于多種復(fù)雜斷裂問(wèn)題。相場(chǎng)斷裂模型起源于1921年提出的Griffith能量理論。該理論突破了傳統(tǒng)的強(qiáng)度理論，通過(guò)能量平衡研究裂紋的擴(kuò)展。1998年，F(xiàn)rancfort和Marigo在Griffith能量理論的基礎(chǔ)上提出了斷裂變分原理，將彈性體的總能量表述為關(guān)于彈性體位移場(chǎng)和裂紋分布的泛函，可通過(guò)最小化總能量來(lái)表征裂紋的成核與擴(kuò)展。斷裂變分原理是對(duì)Griffith能量理論的繼承和發(fā)展，實(shí)現(xiàn)了裂紋起裂、分叉等現(xiàn)象的模擬。在斷裂變分原理的基礎(chǔ)上，Bourdin等進(jìn)一步引入了一個(gè)標(biāo)量場(chǎng)(相場(chǎng))對(duì)裂紋進(jìn)行正則化描述，數(shù)值上實(shí)現(xiàn)了斷裂變分原理，稱(chēng)為相場(chǎng)斷裂模型。目前，相場(chǎng)斷裂模型已被廣泛應(yīng)用于動(dòng)態(tài)斷裂、塑性材料、各向異性材料、巖石材料、復(fù)合材料、結(jié)構(gòu)優(yōu)化等領(lǐng)域。

在相場(chǎng)斷裂模型中，裂紋彌散寬度由定義的正則化參數(shù)()控制，該參數(shù)亦可被視為材料相關(guān)的本構(gòu)參數(shù)。相場(chǎng)斷裂模型要求在裂紋附近進(jìn)行精細(xì)化的有限元網(wǎng)格剖分。Miehe等研究表明裂紋處有限元網(wǎng)格特征尺寸()和正則化參數(shù)需要滿(mǎn)足關(guān)系：≤2。通常情況下，正則化參數(shù)與計(jì)算區(qū)域整體尺寸相比是一個(gè)較小的值，將導(dǎo)致規(guī)模龐大的有限單元數(shù)量。對(duì)于裂紋分布已知或可預(yù)測(cè)的簡(jiǎn)單斷裂問(wèn)題，通?？稍陬A(yù)測(cè)的裂紋擴(kuò)展路徑區(qū)域進(jìn)行網(wǎng)格細(xì)化，而其他區(qū)域采用較粗的網(wǎng)格，控制有限元模型的計(jì)算規(guī)模。對(duì)于裂紋軌跡無(wú)法預(yù)測(cè)甚至可能出現(xiàn)多裂紋交匯的復(fù)雜斷裂問(wèn)題，例如復(fù)合材料、周期性點(diǎn)陣填充結(jié)構(gòu)等，難以對(duì)其裂紋擴(kuò)展軌跡進(jìn)行有效先驗(yàn)預(yù)測(cè)。針對(duì)該類(lèi)斷裂模擬問(wèn)題，需要對(duì)全區(qū)域進(jìn)行精細(xì)化網(wǎng)格剖分，這對(duì)計(jì)算能力提出了極高的挑戰(zhàn)。針對(duì)于此，Burke等發(fā)展了殘差驅(qū)動(dòng)的自適應(yīng)網(wǎng)格細(xì)化算法，并證明基于該算法的自適應(yīng)網(wǎng)格使得能量泛函收斂到局部極小值。Patil等基于多尺度有限元方法實(shí)現(xiàn)了自適應(yīng)網(wǎng)格剖分算法，根據(jù)斷裂相場(chǎng)數(shù)值及其增量判斷裂尖和裂紋位置，在細(xì)網(wǎng)格尺度對(duì)裂紋和裂尖進(jìn)行分析，在粗網(wǎng)格尺度對(duì)其余位置進(jìn)行分析。Ziaei-Rad和Shen實(shí)現(xiàn)了相場(chǎng)斷裂模型的GPU(Graphical Processing Unit)并行計(jì)算，顯著降低了求解時(shí)間。

針對(duì)復(fù)雜斷裂問(wèn)題相場(chǎng)模擬計(jì)算規(guī)模龐大的瓶頸，本文計(jì)劃基于子結(jié)構(gòu)法和損傷識(shí)別，發(fā)展適用于周期性結(jié)構(gòu)脆性斷裂的高效相場(chǎng)模擬方法。首先，根據(jù)周期性特征對(duì)結(jié)構(gòu)進(jìn)行分區(qū)，選取特征子域結(jié)構(gòu)進(jìn)行內(nèi)部節(jié)點(diǎn)的自由度凝聚，生成相應(yīng)的超單元，并由其組裝周期性結(jié)構(gòu)的整體剛度矩陣,降低計(jì)算規(guī)模。然后，離線(xiàn)(Off-line)對(duì)特征子域結(jié)構(gòu)施加拉伸和剪切測(cè)試載荷并進(jìn)行相場(chǎng)斷裂模擬，得到子域內(nèi)部損傷與子域彈性應(yīng)變能之間的定量關(guān)系，標(biāo)定合適的能量判定閾值。在線(xiàn)(On-line)斷裂模擬過(guò)程中將僅對(duì)閾值以上子域的內(nèi)部位移和斷裂相場(chǎng)進(jìn)行解耦分析，避免非裂紋擴(kuò)展子域的重復(fù)計(jì)算消耗。

1 相場(chǎng)斷裂模型

1.1 裂紋正則化

圖1(a)為一含裂紋彈性體的示意圖：∈、?∈-1分別為彈性體區(qū)域及其邊界，∈-1為裂紋；∈{1,2,3}表示空間維度，、、分別表示一維、二維、三維空間；邊界分為?和?，分別對(duì)應(yīng)Dirichlet和Neumann邊界。如圖1(b)所示，通過(guò)定義一個(gè)取值在區(qū)間[0,1]的標(biāo)量相場(chǎng)(,)近似描述裂紋的分布，即裂紋的正則化描述；其中，為該材料點(diǎn)的位置，為加載時(shí)間。當(dāng)相場(chǎng)取值為0時(shí)表示材料完好(無(wú)裂紋)；相場(chǎng)取值為1時(shí)表示材料完全破壞(裂紋)；當(dāng)相場(chǎng)取值介于0～1時(shí)，材料處于完好和完全破壞之間的過(guò)渡狀態(tài)，即彈性體中彌散分布著裂紋。

圖1 裂紋的正則化表示Fig.1 Regularized representation of crack

滿(mǎn)足以上斷裂相場(chǎng)條件的標(biāo)量場(chǎng)可通過(guò)求解以下邊值問(wèn)題得到：

(1)

(2)

邊值問(wèn)題式(1)的歐拉變分形式為

(3)

式中：(,) 滿(mǎn)足Dirichlet邊界條件

={|(,)=1,在上}

(4)

為裂紋表面函數(shù)，計(jì)算公式為

(5)

式中:為彈性體中彌散分布裂紋的密度函數(shù)，計(jì)算公式為

(6)

1.2 控制方程強(qiáng)形式

Francfort與Marigo提出了斷裂變分原理，該原理是相場(chǎng)斷裂模型的理論基礎(chǔ)。斷裂變分原理將產(chǎn)生裂紋彈性體的自由能分為彈性應(yīng)變能和裂紋表面能2部分，裂紋成核與擴(kuò)展的同時(shí)伴隨著彈性應(yīng)變能向裂紋表面能的轉(zhuǎn)換，即彈性應(yīng)變能和裂紋表面能分別驅(qū)動(dòng)和阻礙著裂紋的成核與擴(kuò)展：

(7)

式中：為位移；為應(yīng)變張量；()為彈性應(yīng)變能密度函數(shù)，其在彈性體上積分為彈性應(yīng)變能；為臨界能量釋放率,其在裂紋表面上積分為裂紋表面能。由于裂紋位移場(chǎng)的不連續(xù)性，且裂尖處存在奇異的應(yīng)力場(chǎng)，裂紋表面能難以計(jì)算。

通過(guò)引入相場(chǎng)標(biāo)量(,)對(duì)裂紋進(jìn)行正則化描述，即由全域彌散分布的裂紋近似描述真實(shí)裂紋，彈性體自由能可近似表述為

(8)

(9)

為防止材料在受壓狀態(tài)下出現(xiàn)裂紋，需要對(duì)應(yīng)變能進(jìn)行拉壓狀態(tài)分解

()=()+()

(10)

式中：()、()分別為拉伸和壓縮部分的應(yīng)變能。應(yīng)變能的分解方法有多種，本文基于應(yīng)變能譜分解方法對(duì)其進(jìn)行拉壓分解

(11)

=+

(12)

(13)

(,)=[()+]()+()

(14)

式中：為引入的一個(gè)極小值，用于保證材料完全破壞時(shí)單元?jiǎng)偠染仃嚨姆瞧娈愋裕?)為彈性應(yīng)變能退化函數(shù)，描述斷裂相場(chǎng)對(duì)彈性應(yīng)變能拉伸部分的影響，需滿(mǎn)足：

(15)

通常選取

(16)

拉壓分解后的彈性體自由能表述為

(17)

對(duì)式(17)進(jìn)行位移場(chǎng)和相場(chǎng)的變分可得

(18)

相應(yīng)邊值問(wèn)題為

(19)

=-Δ

(20)

(21)

由于斷裂的不可逆性，為防止卸載后裂紋發(fā)生愈合，引入歷史狀態(tài)變量H

(22)

替代邊值問(wèn)題式(19)中的。

1.3 控制方程弱形式

忽略體力的作用，位移場(chǎng)的弱形式為

(23)

根據(jù)彈性應(yīng)變能密度的定義，應(yīng)力計(jì)算公式為

(24)

式中：為二階單位張量。

將相場(chǎng)控制方程乘上試函數(shù)在全域上積分，得到其等效積分形式

(25)

根據(jù)分部積分法則、散度定理，并引入邊界條件，可得相場(chǎng)控制方程的弱形式

(26)

1.4 有限元離散

對(duì)位移場(chǎng)和相場(chǎng)進(jìn)行有限元離散，可得位移場(chǎng)平衡方程和相場(chǎng)演化方程

(27)

式中：、分別為位移場(chǎng)和相場(chǎng)的整體剛度矩陣；、分別為位移向量和相場(chǎng)向量；、分別為位移場(chǎng)和相場(chǎng)的外載荷向量。

(28)

(29)

2 子結(jié)構(gòu)建模及損傷識(shí)別

對(duì)于周期性結(jié)構(gòu)的斷裂模擬，可利用結(jié)構(gòu)的周期性特征，將結(jié)構(gòu)劃分為若干幾何構(gòu)型一致的子域，通過(guò)子結(jié)構(gòu)法對(duì)子域的內(nèi)部計(jì)算節(jié)點(diǎn)進(jìn)行靜態(tài)凝聚降低計(jì)算規(guī)模；其后，結(jié)合線(xiàn)下測(cè)試標(biāo)定的應(yīng)變能損傷閾值，判定是否需進(jìn)一步求解子域的內(nèi)部響應(yīng)，避免非必要的計(jì)算消耗。

2.1 子結(jié)構(gòu)建模

根據(jù)結(jié)構(gòu)的周期性特征，可將其劃分成若干幾何構(gòu)型一致的子域。數(shù)值上，具有相同幾何特征的子域結(jié)構(gòu)的剛度矩陣一致，因此以下某一代表性子域結(jié)構(gòu)的位移或相場(chǎng)的有限元列式為

=

(30)

根據(jù)節(jié)點(diǎn)是否位于子域的邊界上，式(30)可分區(qū)域表示為

(31)

式中：下標(biāo)i、b分別表示該自由度對(duì)應(yīng)節(jié)點(diǎn)位于子域內(nèi)部和邊界。由式(32)可得

(32)

將式(32)代入式(31)可得

(33)

定義以下矩陣和向量：

(34)

可得邊界節(jié)點(diǎn)響應(yīng)的有限元列式：

=

(35)

以上過(guò)程稱(chēng)為子結(jié)構(gòu)法靜態(tài)凝聚。子結(jié)構(gòu)建模的過(guò)程如圖2所示，具體步驟為：① 對(duì)特征子域進(jìn)行靜態(tài)凝聚；② 由凝聚的剛度矩陣組裝整體剛度矩陣；③ 求解整體結(jié)構(gòu)有限元問(wèn)題；④ 根據(jù)節(jié)點(diǎn)響應(yīng)解耦求解內(nèi)部響應(yīng)。

相場(chǎng)斷裂模擬為非線(xiàn)性問(wèn)題，需多次迭代施加載荷，重復(fù)以上計(jì)算步驟。需要注意的是，若子結(jié)構(gòu)內(nèi)部出現(xiàn)裂紋，需對(duì)該子結(jié)構(gòu)的位移和相場(chǎng)凝聚剛度矩陣進(jìn)行更新，并重新組裝結(jié)構(gòu)的整體剛度矩陣。

圖2 周期性結(jié)構(gòu)的子結(jié)構(gòu)建模Fig.2 Substructuring of periodic structure

2.2 損傷識(shí)別

在利用子結(jié)構(gòu)法對(duì)周期性結(jié)構(gòu)進(jìn)行斷裂模擬時(shí)，若子結(jié)構(gòu)內(nèi)部無(wú)裂紋產(chǎn)生，則無(wú)需求解子結(jié)構(gòu)內(nèi)部的位移場(chǎng)和相場(chǎng)，此類(lèi)子結(jié)構(gòu)僅需執(zhí)行一次靜態(tài)凝聚；而內(nèi)部產(chǎn)生裂紋的子結(jié)構(gòu)，需進(jìn)一步求解子結(jié)構(gòu)內(nèi)部的位移場(chǎng)和裂紋相場(chǎng)。由于裂紋的成核與擴(kuò)展會(huì)引起材料性能的退化，故在每一個(gè)加載步中均需對(duì)含裂紋子結(jié)構(gòu)重新進(jìn)行靜態(tài)凝聚。為判別子結(jié)構(gòu)內(nèi)部是否有潛在裂紋形成的可能性，需進(jìn)行子結(jié)構(gòu)的損傷識(shí)別。由于每次加載時(shí)位移場(chǎng)平衡方程的邊界條件發(fā)生改變，而裂紋相場(chǎng)演化方程的邊界條件保持不變，所以選擇利用子結(jié)構(gòu)邊界位移響應(yīng)進(jìn)行損傷識(shí)別。同時(shí)，斷裂的產(chǎn)生伴隨著彈性應(yīng)變能和裂紋表面能之間的轉(zhuǎn)換，故將通過(guò)子結(jié)構(gòu)邊界位移響應(yīng)計(jì)算得到的彈性應(yīng)變能作為損傷識(shí)別的指標(biāo),計(jì)算公式為

(36)

以二維情況為例，為研究某一代表性子域結(jié)構(gòu)內(nèi)部損傷與邊界位移響應(yīng)關(guān)系，對(duì)其施加如圖3所示6種典型離線(xiàn)測(cè)試工況，其中為位移載荷。若子域結(jié)構(gòu)具有幾何對(duì)稱(chēng)性，可相應(yīng)減少測(cè)試工況的類(lèi)型。根據(jù)離線(xiàn)測(cè)試結(jié)果標(biāo)定該類(lèi)型子域結(jié)構(gòu)內(nèi)部發(fā)生潛在裂紋損傷的能量識(shí)別閾值。

圖3 測(cè)試工況：橫向/縱向的拉伸和壓縮,面內(nèi)剪切Fig.3 Testing loads: Horizontal/vertical tension and compression, and in-plane shear

2.3 算法流程

基于子結(jié)構(gòu)建模和損傷識(shí)別，本文發(fā)展了一種適用于周期性結(jié)構(gòu)脆性斷裂模擬的高效相場(chǎng)方法，其計(jì)算流程歸納如下：

將待分析周期性結(jié)構(gòu)進(jìn)行有限元網(wǎng)格剖分，保證相同特征子域分區(qū)的網(wǎng)格形式一致。

對(duì)特征子域結(jié)構(gòu)的位移和相場(chǎng)的剛度矩陣進(jìn)行靜態(tài)凝聚，組裝結(jié)構(gòu)整體剛度矩陣。

對(duì)特征子域結(jié)構(gòu)進(jìn)行離線(xiàn)典型工況測(cè)試，標(biāo)定基于應(yīng)變能的損傷識(shí)別閾值。

施加載荷增量進(jìn)行在線(xiàn)斷裂模擬，根據(jù)子結(jié)構(gòu)邊界響應(yīng)計(jì)算的彈性應(yīng)變能進(jìn)行損傷識(shí)別判定：若低于閾值，無(wú)需求解內(nèi)部響應(yīng)；若高于閾值，求解內(nèi)部響應(yīng)，并重新進(jìn)行靜態(tài)凝聚。

重新組裝位移和相場(chǎng)的整體剛度矩陣，重復(fù)步驟4，直至分析完成。

3 數(shù)值算例

考慮由方孔單胞周期性排列構(gòu)成結(jié)構(gòu)的斷裂模擬。首先，對(duì)方孔單胞進(jìn)行離線(xiàn)測(cè)試，標(biāo)定損傷識(shí)別的能量閾值；隨后，基于子結(jié)構(gòu)建模和損傷識(shí)別對(duì)兩類(lèi)周期性結(jié)構(gòu)進(jìn)行在線(xiàn)的相場(chǎng)斷裂模擬，并與常規(guī)相場(chǎng)法的模擬結(jié)果和計(jì)算規(guī)模進(jìn)行比較。

3.1 方孔單胞離線(xiàn)測(cè)試

如圖4所示，考慮2種不同尺寸的方孔單胞，分別定義為類(lèi)型A和B。A和B兩類(lèi)單胞均采用邊長(zhǎng)=0.005 mm的平面應(yīng)力四節(jié)點(diǎn)正方形單元進(jìn)行網(wǎng)格剖分。參照文獻(xiàn)[31]，脆性材料的參數(shù)設(shè)置為：=121.5 kN/mm,=80.77 kN/mm，=2.7×10kN/mm；相場(chǎng)模型中的正則化參數(shù)設(shè)置為=0.015 mm。

由于方孔單胞的幾何對(duì)稱(chēng)性，可相應(yīng)簡(jiǎn)化典型工況的數(shù)目；此外，由于相場(chǎng)模型對(duì)彈性應(yīng)變能進(jìn)行了拉壓狀態(tài)分解，壓縮狀態(tài)的材料難以退化形成裂紋。綜合以上因素，僅需考慮縱向拉伸和面內(nèi)剪切2種工況進(jìn)行損傷能量閾值的標(biāo)定。單胞A在加載過(guò)程中彈性應(yīng)變能-位移曲線(xiàn)和裂紋軌跡如圖5所示。單胞彈性應(yīng)變能在加載過(guò)程中先逐漸增大，隨后在裂紋成核引起材料退化時(shí)發(fā)生驟降。單胞B在典型工況加載過(guò)程中曲線(xiàn)變化趨勢(shì)和單胞A一致，僅應(yīng)變能的幅值成比例相應(yīng)增大。

圖4 2種不同尺寸單胞的幾何構(gòu)型和網(wǎng)格剖分Fig.4 Geometrical configurations and FE meshes of two cells of different sizes

選擇拉伸和剪切加載過(guò)程中方孔單胞相場(chǎng)最大值分別為0.05、0.10、0.20時(shí)的彈性應(yīng)變能作為損傷判定的閾值，并進(jìn)行編號(hào)，如表1所示。對(duì)

圖5 單胞A測(cè)試工況下彈性應(yīng)變能-位移曲線(xiàn)和裂紋軌跡Fig.5 Elastic energy-displacement curves and crack paths of cell A under testing loads

表1 單胞的損傷識(shí)別能量閾值Table 1 Energy thresholds for cell damage identification

比2種工況下標(biāo)定的能量閾值，剪切加載的閾值均小于拉伸加載的閾值。后續(xù)算例驗(yàn)算中，選取剪切工況下標(biāo)定的能量閾值作為單胞損傷的判定閾值。

3.2 Ⅱ 型斷裂模擬

考慮如圖6所示的II型斷裂模擬算例，該結(jié)構(gòu)由900個(gè)方孔單胞A周期性排列構(gòu)成。該結(jié)構(gòu)邊長(zhǎng)6 mm，內(nèi)部預(yù)置3 mm長(zhǎng)度裂紋，下端固定，上端施加方向位移荷載，增量為5×10mm，共加載350步。全域離散的有限元模型含節(jié)點(diǎn)數(shù)111.8萬(wàn)；采用子結(jié)構(gòu)法對(duì)內(nèi)部節(jié)點(diǎn)進(jìn)行凝聚，計(jì)算節(jié)點(diǎn)數(shù)降低至7.4萬(wàn)。加載過(guò)程中，彈性應(yīng)變能大于損傷識(shí)別能量閾值的子結(jié)構(gòu)數(shù)目逐步增多,計(jì)算規(guī)模隨之相應(yīng)增大。分別采用常規(guī)相場(chǎng)法和本文方法進(jìn)行該結(jié)構(gòu)的斷裂模擬，比對(duì)不同損傷閾值對(duì)模擬精度的影響。圖7為采用2種相場(chǎng)法預(yù)測(cè)的載荷-位移曲線(xiàn)以及采用A2損傷閾值時(shí)的裂紋軌跡擴(kuò)展情況(圖6虛框區(qū)域)。圖8為采用2種相場(chǎng)法預(yù)測(cè)的最終裂紋軌跡(圖6 虛框區(qū)域)。圖7、圖8中紅框區(qū)域?yàn)榻?jīng)識(shí)別需求解內(nèi)部物理響應(yīng)和裂紋相場(chǎng)的子結(jié)構(gòu)。對(duì)比可見(jiàn)，采用A3損傷閾值時(shí)，本文方法與常規(guī)相場(chǎng)法的模擬結(jié)果偏差較大；而取值相對(duì)保守的A1、A2損傷閾值時(shí)，與常規(guī)相場(chǎng)法模擬結(jié)果的一致性較好。圖9所示為3種不同損傷閾值取

圖6 周期性結(jié)構(gòu)的Ⅱ型斷裂問(wèn)題Fig.6 Mode Ⅱ fracture of periodic structure

值情況下，加載過(guò)程中經(jīng)識(shí)別需進(jìn)一步求解內(nèi)部響應(yīng)子結(jié)構(gòu)數(shù)量的占比情況。即使是最為保守的A1損傷閾值情況下，最終經(jīng)識(shí)別需進(jìn)一步求解的子結(jié)構(gòu)數(shù)量的占比僅為16.7%(150/900)，計(jì)算規(guī)模顯著降低。

圖7 周期性結(jié)構(gòu)Ⅱ型斷裂問(wèn)題的載荷-位移曲線(xiàn)及裂紋擴(kuò)展情況Fig.7 Load-displacement curves and crack propagation of mode Ⅱ fracturing of periodic structure

圖8 周期性結(jié)構(gòu)Ⅱ型斷裂問(wèn)題的最終裂紋軌跡Fig.8 Ultimate crack paths of mode Ⅱ fracturing of periodic structure

圖9 Ⅱ型斷裂問(wèn)題過(guò)程中含裂紋子結(jié)構(gòu)占比Fig.9 Proportion of substructures with cracks during mode Ⅱ fracturing

3.3 L形梁斷裂

進(jìn)一步驗(yàn)證本文所提方法的適用性,考慮如圖10所示的由675個(gè)方孔單胞A周期性排列構(gòu)成的L形梁的斷裂模擬。將L形梁的下端固定，在圖示箭頭處豎直向上施加位移荷載，增量為1×10mm，共加載200步。根據(jù)前一算例分析結(jié)果，綜合考慮模擬精度和計(jì)算規(guī)模，選擇能量閾值A(chǔ)2進(jìn)行子結(jié)構(gòu)的損傷識(shí)別。分別若采用常規(guī)相場(chǎng)法和本文方法進(jìn)行該結(jié)構(gòu)的斷裂模擬。圖11為采用2種相場(chǎng)法預(yù)測(cè)的載荷-位移曲線(xiàn)。圖12為采用2種相場(chǎng)法預(yù)測(cè)的最終裂紋軌跡(圖10虛框區(qū)域)。2種方法模擬結(jié)果的一致性較好。圖12(b)中紅框區(qū)域?yàn)榻?jīng)識(shí)別最終需求解內(nèi)部響應(yīng)的子結(jié)構(gòu)，僅占子結(jié)構(gòu)總數(shù)的4.6%(31/675)，計(jì)算規(guī)模顯著低于常規(guī)相場(chǎng)法。

圖10 單胞A周期性排列構(gòu)成的L形梁Fig.10 L-beam periodically composed with cell A

圖11 單胞A組成的L形梁斷裂模擬載荷-位移曲線(xiàn)Fig.11 Load-displacement curves of fracture modeling for L-beam composed of cell A

圖12 單胞A組成的L形梁斷裂模擬最終裂紋軌跡Fig.12 Ultimate crack paths of fracture modeling for L-beam composed of cell A

接下來(lái)，進(jìn)一步驗(yàn)證方法的適用性，考慮由300個(gè)方孔單胞B周期性排列構(gòu)成的L形梁的斷裂模擬，如圖13所示。類(lèi)似地，選擇表1中對(duì)應(yīng)的能量閾值B2進(jìn)行子結(jié)構(gòu)的損傷識(shí)別。圖14為采用2種相場(chǎng)法預(yù)測(cè)的載荷-位移曲線(xiàn)。圖15為采用2種相場(chǎng)法預(yù)測(cè)的最終裂紋軌跡(圖13虛框區(qū)域)。與前一算例類(lèi)似，2種方法模擬結(jié)果的一致性較好，且計(jì)算規(guī)模得到顯著降低；圖15(b)中紅框區(qū)域識(shí)別出的子結(jié)構(gòu)數(shù)目占比為6.3%(19/300)。對(duì)比圖11、圖14的載荷-位移曲線(xiàn)及圖12、圖15的最終裂紋擴(kuò)展形式，周期性結(jié)構(gòu)的承載極限和裂紋擴(kuò)展的軌跡取決于結(jié)構(gòu)整體的幾

圖13 單胞B周期性排列構(gòu)成的L形梁Fig.13 L-beam periodically composed of cell B

圖14 單胞B組成的L形梁斷裂模擬載荷-位移曲線(xiàn)Fig.14 Load-displacement curves of fracture modeling for L-beam composed of cell B

圖15 單胞B組成的L形梁斷裂模擬最終裂紋軌跡Fig.15 Ultimate crack paths of fracture modeling for L-beam composed of cell B

何構(gòu)型和邊界條件，受單胞尺寸的影響不顯著。單胞尺寸較大時(shí)，裂紋間斷式成核擴(kuò)展，L形梁的載荷-位移曲線(xiàn)隨著裂紋的擴(kuò)展波浪式下降(圖14)，一定程度上提升了結(jié)構(gòu)的抗裂韌性(曲線(xiàn)下圍面積)。

4 結(jié) 論

基于子結(jié)構(gòu)建模和損傷識(shí)別，發(fā)展了一種適用于周期性結(jié)構(gòu)脆性斷裂的高效相場(chǎng)模擬方法。相較于傳統(tǒng)相場(chǎng)法，本文所提方法在維持?jǐn)嗔涯M精度的同時(shí)，能夠有效地避免非裂紋擴(kuò)展區(qū)域的重復(fù)計(jì)算消耗、顯著降低了有限元計(jì)算分析規(guī)模。后續(xù)研究工作中，計(jì)劃基于人工神經(jīng)網(wǎng)絡(luò)對(duì)單胞的承載類(lèi)型進(jìn)行區(qū)分，根據(jù)承載類(lèi)型標(biāo)定類(lèi)型相關(guān)的損傷閾值，以實(shí)現(xiàn)更復(fù)雜周期性結(jié)構(gòu)的斷裂模擬。針對(duì)三維周期性結(jié)構(gòu)斷裂模擬求解規(guī)模大，靜態(tài)凝聚耗時(shí)長(zhǎng)的問(wèn)題，后續(xù)計(jì)劃結(jié)合其他平臺(tái)或其他高效算法(如基于子結(jié)構(gòu)法的并行計(jì)算)以降低耗時(shí)，將本方法擴(kuò)展到三維斷裂問(wèn)題。此外，基于發(fā)展的高效斷裂分析方法，有望進(jìn)一步結(jié)合優(yōu)化設(shè)計(jì)，通過(guò)優(yōu)化單胞的幾何或拓?fù)錁?gòu)型，提升周期性結(jié)構(gòu)的抗裂韌性。