擴(kuò)展有限元法在斷裂力學(xué)相場(chǎng)模型中的應(yīng)用

張子瑜，郝林

上海飛機(jī)設(shè)計(jì)研究院 載荷部，上海 201210

由于裂紋可能威脅到民用飛機(jī)主傳力結(jié)構(gòu)的完整性，進(jìn)而威脅飛行安全，對(duì)裂紋擴(kuò)展的模擬是疲勞及損傷容限分析的重要內(nèi)容。擴(kuò)展有限元方法(X-FEM)采用單位分解的思路，將單元形函數(shù)空間加以豐富，從而在不要求有限元網(wǎng)格與內(nèi)部邊界吻合的前提下模擬大量的物理問(wèn)題，其在斷裂力學(xué)領(lǐng)域的應(yīng)用日趨廣泛。Hansbo A和Hansbo P提出的幻影節(jié)點(diǎn)法(Phantom-node Method)能夠在技術(shù)層面等效地執(zhí)行X-FEM，此方法將被裂紋切割的單元復(fù)制后僅在被復(fù)制單元的真實(shí)域內(nèi)對(duì)相關(guān)量(如剛度矩陣等)做局部數(shù)值積分。這種方法將傳統(tǒng)X-FEM中添加額外自由度的需求替換為對(duì)網(wǎng)格的更改。

相場(chǎng)法在模擬裂紋擴(kuò)展問(wèn)題時(shí)則采取另一種思路，將表達(dá)破損區(qū)域的相場(chǎng)方程與力學(xué)方程耦合，從而使得裂紋可以自然擴(kuò)展，使得計(jì)算裂紋的前進(jìn)方向和距離不再必要。然而，由于相場(chǎng)法無(wú)法構(gòu)建出清晰的裂紋斷面，故該方法存在一些局限。例如，它很難模擬裂紋斷面之間的接觸及摩擦效應(yīng)，也會(huì)使裂紋區(qū)域單元的局部剛度趨于0而導(dǎo)致求解困難。

在早期介紹水平集(Level-set)方法的文獻(xiàn)[9] 中，X-FEM被用來(lái)模擬完全破壞的區(qū)域(即空洞)。相場(chǎng)法方面的研究者受到啟發(fā)，使用幻影節(jié)點(diǎn)法將被相場(chǎng)等高線(xiàn)切割的單元做局部數(shù)值積分，并刪去位于等高線(xiàn)內(nèi)部的單元，從而模擬破損區(qū)域。本文同樣將X-FEM與相場(chǎng)法結(jié)合，但提出一種構(gòu)建清晰裂紋斷面的替代方法：通過(guò)一定的算法，找到能夠近似地連接相場(chǎng)值為1的所有點(diǎn)的單根軸線(xiàn)(即相場(chǎng)中線(xiàn))，并通過(guò)X-FEM來(lái)實(shí)現(xiàn)相場(chǎng)中線(xiàn)對(duì)網(wǎng)格的切割。

另外，在以往構(gòu)建相場(chǎng)中線(xiàn)的文獻(xiàn)[11-14]中，研究者常連接位于相場(chǎng)等高線(xiàn)內(nèi)部的最大球的球心來(lái)生成中線(xiàn)，但這些方法存在局限性，如過(guò)分耗時(shí)，或僅能解決靜態(tài)裂紋問(wèn)題等。鑒于此，本文利用相場(chǎng)等高線(xiàn)與輔助網(wǎng)格的交點(diǎn)，構(gòu)建一種新的、效率更高的中線(xiàn)生成方法。

1 擴(kuò)展有限元基本算法

幻影節(jié)點(diǎn)法不直接添加額外的形函數(shù)和自由度，而是在裂紋處創(chuàng)建重疊的單元。圖1說(shuō)明了幻影節(jié)點(diǎn)X-FEM的基本思路。當(dāng)1個(gè)單元被裂紋完全切開(kāi)后，創(chuàng)建2個(gè)子單元(或稱(chēng)為副本單元)，2個(gè)子單元將從原來(lái)的母單元繼承不同的物理域，稱(chēng)此物理域?yàn)椤八槠?即圖1右側(cè)單元的陰影部分)。注意，圖1中圓圈所表示的節(jié)點(diǎn)不屬于單元碎片，因此它們被稱(chēng)為“幻影節(jié)點(diǎn)”?；糜肮?jié)點(diǎn)處的數(shù)值解不重要，但可用于對(duì)物理域內(nèi)某點(diǎn)的數(shù)值解進(jìn)行插值。一般說(shuō)來(lái)，一個(gè)子單元擁有的幻影節(jié)點(diǎn)將與其“兄弟”單元完全脫開(kāi)。在2個(gè)子單元被創(chuàng)建后，所有新的自由度會(huì)被適當(dāng)分配，且母單元將被刪除。

圖1 幻影節(jié)點(diǎn)X-FEM中裂紋切割單元Fig.1 Elements cut with phantom-node X-FEM

原始的X-FEM法會(huì)采用分段積分方法對(duì)被裂紋切割的單元求積分。在幻影節(jié)點(diǎn)X-FEM中，只需在每一個(gè)子單元的物理域內(nèi)進(jìn)行積分。因?yàn)榛糜肮?jié)點(diǎn)法更易于被系統(tǒng)化地執(zhí)行，所以更適合商用代碼使用。在有裂紋分叉的情形出現(xiàn)時(shí)，此方法更顯其優(yōu)勢(shì)。注意，文獻(xiàn)[5]的方法只適用于整個(gè)單元被裂紋完全切割的情況。Rabczuk等則為幻影節(jié)點(diǎn)法開(kāi)發(fā)出一種裂紋尖端單元，使得在進(jìn)行靜態(tài)裂紋仿真時(shí)，裂紋尖端可以位于單元內(nèi)部。本文使用的幻影節(jié)點(diǎn)X-FEM切割單元的算法被稱(chēng)為“單元碎裂算法”，其在每一個(gè)時(shí)間步中的基本流程如下：① 更新相鄰單元的信息；② 標(biāo)記裂紋尖端處的單元；③ 標(biāo)記新的切割點(diǎn)；④ 創(chuàng)建碎片；⑤ 從碎片創(chuàng)建子單元；⑥ 連接相鄰的碎片和子單元；⑦ 清理母單元。

2 相場(chǎng)法在斷裂力學(xué)中的應(yīng)用

本文采用Miehe等提出的相場(chǎng)斷裂力學(xué)模型。通常情況下，可以通過(guò)將總勢(shì)能最小化，或使虛功率等于零，來(lái)獲得控制方程的強(qiáng)形式:

+-=0

(1)

式中：為總應(yīng)變能；為耗散能；為外力功?？倓?shì)能為位移和相場(chǎng)的函數(shù):

(2)

若不考慮防止裂紋愈合的方法，脆性材料的耗散能一般可被寫(xiě)為

(3)

(4)

外力功則遵循傳統(tǒng)的形式

(5)

式中：為體積力密度；為單位面積的邊界力;為位移場(chǎng)；Г為Neumann邊界。

將式(2)～式(5)代入式(1)，并使用分部積分，可得如下強(qiáng)形式的控制方程：

(6)

且滿(mǎn)足邊界條件

(7)

為了保證相場(chǎng)的不可逆條件，Miehe等對(duì)式(6)的第2個(gè)方程進(jìn)行了改進(jìn)，引入了1個(gè)的對(duì)偶變量及1個(gè)額外的約束方程。那么，相場(chǎng)斷裂問(wèn)題的控制方程可寫(xiě)為

(8)

(9)

(10)

式(10)可使用有限元方法進(jìn)行離散求解。

3 生成相場(chǎng)中線(xiàn)的方法

構(gòu)建相場(chǎng)中線(xiàn)算法的核心思想為：創(chuàng)建一個(gè)比有限元網(wǎng)格粗的輔助網(wǎng)格，使其能夠覆蓋求解域。之后，求出相場(chǎng)等高線(xiàn)與輔助網(wǎng)格的交點(diǎn)，并將基于這些交點(diǎn)生成能夠表征相場(chǎng)中線(xiàn)的線(xiàn)段。其算法可分為2步：

1) 對(duì)于每一個(gè)輔助網(wǎng)格單元，求出中線(xiàn)交點(diǎn)及其方向向量。其中，中線(xiàn)交點(diǎn)即為相場(chǎng)中線(xiàn)與輔助網(wǎng)格單元邊的交點(diǎn)，而其對(duì)應(yīng)的方向向量代表著相場(chǎng)中線(xiàn)的方向。在執(zhí)行本步的過(guò)程中需遍歷所有輔助網(wǎng)格單元。

2) 在每一個(gè)輔助網(wǎng)格單元中，適當(dāng)?shù)貙⒅芯€(xiàn)交點(diǎn)相連來(lái)創(chuàng)建單元內(nèi)的中線(xiàn)線(xiàn)段。那么，完整的相場(chǎng)中線(xiàn)將由所有單元的中線(xiàn)線(xiàn)段組成。

注意，相場(chǎng)等高線(xiàn)由常規(guī)的行進(jìn)正方形(Marching Squares)算法生成。此外，使用排序算法將組成等高線(xiàn)的小線(xiàn)段沿著逆時(shí)針?lè)较蚺帕写鎯?chǔ)。圖2為相場(chǎng)云圖及=0.9所對(duì)應(yīng)的等高線(xiàn)示意圖。值得一提的是，即使在行進(jìn)正方形算法中將相場(chǎng)值設(shè)置為1，通常情況下也無(wú)法得到1條單一的等高線(xiàn)。這是因?yàn)?，有限元方法的離散性使其無(wú)法捕捉到相場(chǎng)值嚴(yán)格等于1的尖點(diǎn)。

圖2 相場(chǎng)云圖及其對(duì)應(yīng)的等高線(xiàn)Fig.2 A phase-field and its isocontour

3.1 計(jì)算輔助網(wǎng)格單元的中線(xiàn)交點(diǎn)

在計(jì)算相場(chǎng)中線(xiàn)與輔助網(wǎng)格單元的交點(diǎn)時(shí)，通常有3種情況需要處理，現(xiàn)論述如下。

輔助網(wǎng)格單元的邊和相場(chǎng)等高線(xiàn)只有一個(gè)交點(diǎn)。

以圖3左下角的單元為例來(lái)說(shuō)明算法。當(dāng)遍歷此單元的所有邊后，可以發(fā)現(xiàn)右側(cè)的邊與相場(chǎng)等高線(xiàn)只有1個(gè)交點(diǎn)，且此單元右上角的節(jié)點(diǎn)距離這個(gè)交點(diǎn)最近。那么將與這個(gè)節(jié)點(diǎn)相連的所有單元邊(這些邊形成1個(gè)“十”字)與相場(chǎng)等高線(xiàn)相交，則可得到4個(gè)交點(diǎn)。最后，將這4個(gè)點(diǎn)的坐標(biāo)取平均后即得相場(chǎng)中線(xiàn)與單元邊的交點(diǎn)(即圖3中的點(diǎn)5)。這一過(guò)程被稱(chēng)為“節(jié)點(diǎn)塊搜索”。

圖3 情況1示意圖Fig.3 Illustration of Case 1

相場(chǎng)中線(xiàn)與輔助網(wǎng)格交點(diǎn)的方向向量可基于相場(chǎng)等高線(xiàn)與輔助網(wǎng)格交點(diǎn)的切向量來(lái)計(jì)算。如前所述，由于組成相場(chǎng)等高線(xiàn)的小線(xiàn)段以逆時(shí)針排列存儲(chǔ)，故可將與輔助網(wǎng)格單元邊相交的相場(chǎng)等高線(xiàn)小線(xiàn)段的方向作為等高線(xiàn)交點(diǎn)的切向量。取第1個(gè)等高線(xiàn)交點(diǎn)和最后一個(gè)等高線(xiàn)交點(diǎn)的切向量，將二者平均后即可得到相場(chǎng)中線(xiàn)與輔助網(wǎng)格交點(diǎn)的方向向量。在圖3中，相場(chǎng)中線(xiàn)交點(diǎn)的方向向量即為0.5(-)，其中，為等高線(xiàn)交點(diǎn)1的切向量，為等高線(xiàn)交點(diǎn)4的切向量。注意，相場(chǎng)中線(xiàn)交點(diǎn)方向向量的方向與等高線(xiàn)交點(diǎn)1切向量的方向一致，因?yàn)楹笳叩姆较蛲ǔ４砹肆鸭y擴(kuò)展的方向。

輔助網(wǎng)格單元的邊和相場(chǎng)等高線(xiàn)有多個(gè)交點(diǎn)。

以圖4為例，考察左上角的單元，其左側(cè)的邊與相場(chǎng)等高線(xiàn)有2個(gè)交點(diǎn)，故只需將2個(gè)交點(diǎn)的坐標(biāo)取平均即得到相場(chǎng)中線(xiàn)交點(diǎn)。對(duì)于右側(cè)的邊，可以找到3個(gè)等高線(xiàn)交點(diǎn)，標(biāo)記為編號(hào)1～3，如圖4所示。注意，編號(hào)1～3事實(shí)上也代表了相場(chǎng)等高線(xiàn)小線(xiàn)段在數(shù)組中的存儲(chǔ)順序。將這3個(gè)交點(diǎn)按照距離單元右側(cè)邊第1個(gè)節(jié)點(diǎn)由近及遠(yuǎn)的順序排序后，可得數(shù)組{交點(diǎn)2, 交點(diǎn)3, 交點(diǎn)1}。之后，遍歷交點(diǎn)對(duì)，依次計(jì)算交點(diǎn)2、3的中點(diǎn)CutNode23和交點(diǎn)3、1的中點(diǎn)CutNode13。鑒于CutNode23位于相場(chǎng)等高線(xiàn)圍成的區(qū)域之外，故只有CutNode13可作為相場(chǎng)中線(xiàn)交點(diǎn)。此外，等高線(xiàn)交點(diǎn)2未與其他交點(diǎn)形成交點(diǎn)對(duì)，故它被稱(chēng)作孤立交點(diǎn)，可調(diào)用節(jié)點(diǎn)塊搜索來(lái)計(jì)算其對(duì)應(yīng)的相場(chǎng)中線(xiàn)交點(diǎn)。與情況1類(lèi)似，CutNode13的方向向量可由0.5(-)計(jì)算。

圖4 情況2示意圖Fig.4 Illustration of Case 2

還可以判斷相場(chǎng)中線(xiàn)方向向量與輔助網(wǎng)格單元的關(guān)系，即，它是在“進(jìn)入”“離開(kāi)”或者“掠過(guò)”當(dāng)前的輔助網(wǎng)格單元。具體的判斷準(zhǔn)則為：若相場(chǎng)中線(xiàn)交點(diǎn)方向向量與輔助網(wǎng)格單元邊法向量的點(diǎn)乘為正，則方向向量“離開(kāi)”；若點(diǎn)乘為負(fù)，則方向向量“進(jìn)入”；若方向向量與2個(gè)單元邊都有點(diǎn)乘，且二者符號(hào)相反(如圖4左下角單元中右上節(jié)點(diǎn)附近的相場(chǎng)中線(xiàn)交點(diǎn))，則方向向量“掠過(guò)”此單元。

在圖4中，雖然輔助網(wǎng)格單元邊與相場(chǎng)等高線(xiàn)只有3個(gè)交點(diǎn)，但情況2的算法可推廣到更多交點(diǎn)的情況。

找到輔助網(wǎng)格單元內(nèi)部的相場(chǎng)中線(xiàn)交點(diǎn)。

裂紋的起始點(diǎn)通常位于輔助網(wǎng)格單元的內(nèi)部，稱(chēng)作尖端交點(diǎn)，如圖5所示。如果一個(gè)輔助網(wǎng)格單元與相場(chǎng)等高線(xiàn)只有一個(gè)中線(xiàn)交點(diǎn)，搜索等高線(xiàn)在單元內(nèi)部的缺口，那么缺口的中點(diǎn)即為尖端交點(diǎn)。

圖5 情況3示意圖Fig.5 Illustration of Case 3

總之，上述算法可使我們找到每個(gè)輔助單元上的中線(xiàn)交點(diǎn)以及每個(gè)中線(xiàn)交點(diǎn)的方向向量和方向向量的屬性(即“離開(kāi)”“進(jìn)入”或“掠過(guò)”)，方向向量是中線(xiàn)線(xiàn)段局部方向的表征。有了這些信息后，就能夠生成中線(xiàn)線(xiàn)段。

3.2 生成輔助網(wǎng)格單元的中線(xiàn)線(xiàn)段

尋找輔助網(wǎng)格單元內(nèi)中線(xiàn)線(xiàn)段的方法與單元中線(xiàn)交點(diǎn)的數(shù)目有關(guān)。如果某單元只有2個(gè)中線(xiàn)交點(diǎn)，則這2個(gè)點(diǎn)可直接相連以創(chuàng)建表征裂紋的線(xiàn)段。若某輔助單元有3個(gè)中線(xiàn)交點(diǎn)，且其中1個(gè)是進(jìn)入狀態(tài)、另外2個(gè)是離開(kāi)狀態(tài)(簡(jiǎn)稱(chēng)為“一進(jìn)兩出”)，則首先找到3個(gè)中線(xiàn)交點(diǎn)的方向向量的交點(diǎn)。若3個(gè)向量的交點(diǎn)位于此輔助單元內(nèi)部，則此交點(diǎn)即為裂紋分叉點(diǎn)，且裂紋線(xiàn)段可通過(guò)連接分叉點(diǎn)和3個(gè)中線(xiàn)交點(diǎn)來(lái)建立(見(jiàn)圖4)。若3個(gè)向量的交點(diǎn)位于此輔助單元外部，則將進(jìn)入狀態(tài)的中線(xiàn)交點(diǎn)作為分叉點(diǎn)。若3個(gè)中線(xiàn)交點(diǎn)的狀態(tài)并非“一進(jìn)兩出”，則使用下面介紹的對(duì)齊算法。

對(duì)于擁有大于等于3個(gè)中線(xiàn)交點(diǎn)的輔助單元而言，可使用對(duì)齊算法確定兩兩相連的中線(xiàn)線(xiàn)段。對(duì)齊算法執(zhí)行時(shí)，程序遍歷所有中線(xiàn)交點(diǎn)。對(duì)于任一中線(xiàn)交點(diǎn)，嘗試將此點(diǎn)與其他所有的中線(xiàn)交點(diǎn)相連。對(duì)于每一個(gè)連接線(xiàn)段，程序計(jì)算此線(xiàn)段與兩端中線(xiàn)交點(diǎn)處方向向量的夾角，其中較大的角被定義為對(duì)齊指標(biāo)，用于衡量中線(xiàn)線(xiàn)段與兩端方向向量的一致程度。對(duì)齊指標(biāo)的定義為

=max(,)

(11)

=

(12)

式中：=1或2；=-；、分別為中線(xiàn)線(xiàn)段兩端中線(xiàn)交點(diǎn)的坐標(biāo)向量；為中線(xiàn)交點(diǎn)處的方向向量。若某中線(xiàn)交點(diǎn)的某一個(gè)中線(xiàn)線(xiàn)段的對(duì)齊指標(biāo)最小，則此線(xiàn)段為唯一有效的中線(xiàn)線(xiàn)段。為了詳細(xì)解釋?zhuān)疾靾D6中的中線(xiàn)交點(diǎn)3。在所有可能的中線(xiàn)線(xiàn)段3-4、線(xiàn)段3-1和線(xiàn)段3-2中，線(xiàn)段3-4及其2個(gè)端點(diǎn)的方向向量的對(duì)齊程度最好。因此，線(xiàn)段3-4是中線(xiàn)交點(diǎn)3唯一有效的中線(xiàn)線(xiàn)段。

圖6 對(duì)齊算法示意圖Fig.6 Alignment algorithm

當(dāng)某輔助網(wǎng)格單元的中線(xiàn)交點(diǎn)數(shù)量為奇數(shù)時(shí)，上述算法可能導(dǎo)致某中線(xiàn)交點(diǎn)同時(shí)從屬于多條中線(xiàn)線(xiàn)段。若對(duì)這類(lèi)中線(xiàn)交點(diǎn)再次使用對(duì)齊算法，可使僅擁有最佳對(duì)齊指標(biāo)的中線(xiàn)線(xiàn)段被保留。

3.3 添加新的中線(xiàn)線(xiàn)段

隨著仿真的進(jìn)行，相場(chǎng)等高線(xiàn)會(huì)動(dòng)態(tài)變化。當(dāng)有新的中線(xiàn)交點(diǎn)生成時(shí)，程序會(huì)嘗試將新中線(xiàn)交點(diǎn)的方向向量和已有的中線(xiàn)線(xiàn)段相交來(lái)生成新的中線(xiàn)線(xiàn)段，如圖7所示。在圖7中，線(xiàn)段1-2是已經(jīng)存在的中線(xiàn)線(xiàn)段，中線(xiàn)交點(diǎn)3是新的中線(xiàn)交點(diǎn)。程序計(jì)算中線(xiàn)交點(diǎn)3處方向向量與中線(xiàn)線(xiàn)段1-2的交點(diǎn)(即圖7中的點(diǎn)4)，并將其與中線(xiàn)交點(diǎn)3相連，從而生成了新的中線(xiàn)線(xiàn)段。

圖7 生成新中線(xiàn)線(xiàn)段的過(guò)程Fig.7 Process of computing new medial segments

4 數(shù)值算例

4.1 靜態(tài)等高線(xiàn)的中線(xiàn)

圖8所示為一個(gè)給定的靜態(tài)等高線(xiàn)的中線(xiàn)生成結(jié)果。在本例中，輔助網(wǎng)格被巧妙處理以創(chuàng)造更多的棘手情況，從而測(cè)試程序的健壯性。圖8中，中線(xiàn)線(xiàn)段以紅線(xiàn)表示。特別地，單元9左上角中線(xiàn)交點(diǎn)的計(jì)算即為3.1節(jié)中節(jié)點(diǎn)塊搜索的典型例子，而單元22的中線(xiàn)交點(diǎn)符合“一進(jìn)兩出”，其分叉結(jié)構(gòu)由圖4所示的方法生成(詳見(jiàn)3.2 節(jié))。此外，使用對(duì)齊算法生成了單元21的中線(xiàn)線(xiàn)段。

圖8 靜態(tài)等高線(xiàn)的中線(xiàn)生成Fig.8 Medial-axis generation of a static isocontour

4.2 動(dòng)態(tài)等高線(xiàn)的中線(xiàn)

4.2.1 裂紋融合

設(shè)計(jì)1個(gè)二維斷裂數(shù)值算例：1個(gè)1 mm×1 mm 的方形板在=0～0.00 308 203 s的時(shí)間內(nèi)受拉，拉力施加于上部邊界。隨后，上部邊界停止受拉，轉(zhuǎn)而在右側(cè)邊界受拉。

上部邊界的位移邊界條件為

(13)

右側(cè)邊界的位移邊界條件為

(14)

將材料設(shè)置為各向同性，楊氏模量=208 kN/mm，泊松比=0.3，本構(gòu)閾值=10kN/mm，相場(chǎng)特征尺度=0.03 mm，人工黏性=10(kN·s)/mm。使用開(kāi)源平臺(tái)MOOSE(Multi-physics Object Oriented Simulation Environment)進(jìn)行仿真，結(jié)果見(jiàn)圖9。可見(jiàn)，紅色的線(xiàn)段所示的相場(chǎng)等高線(xiàn)中線(xiàn)能夠隨著相場(chǎng)等高線(xiàn)適當(dāng)?shù)財(cái)U(kuò)展，且捕捉到了裂紋融合的過(guò)程。

圖9 二維融合裂紋的相場(chǎng)云圖及其中線(xiàn)Fig.9 Phase-field contour and medial-axes of 2D merging cracks

4.2.2 裂紋分叉

4.3 使用相場(chǎng)中線(xiàn)切割網(wǎng)格

圖10 二維分叉裂紋的相場(chǎng)云圖及其中線(xiàn)Fig.10 Phase-field contour and medial-axis of a 2D branching crack

圖11 使用單條裂紋的相場(chǎng)中線(xiàn)切割網(wǎng)格Fig.11 Mesh-cutting using medial-axis of a single crack

圖12對(duì)比了不切割網(wǎng)格與切割網(wǎng)格時(shí)的位移-載荷曲線(xiàn)，二者幾乎完全重合。結(jié)合圖11、圖12 可知，由X-FEM造成的網(wǎng)格切割不影響相場(chǎng)結(jié)果。此時(shí)，可安全無(wú)虞地在清晰的斷裂面上引入接觸、摩擦等算法。值得注意的是，使用X-FEM切割網(wǎng)格后，位于破壞區(qū)域的單元的長(zhǎng)寬比得到了改善。

圖12 裂紋擴(kuò)展問(wèn)題的位移-載荷曲線(xiàn)Fig.12 Load-displacement curves of single edge-notched tension test

5 結(jié) 論

提出了一種能夠捕捉二維斷裂力學(xué)問(wèn)題相場(chǎng)等高線(xiàn)中線(xiàn)的算法。通過(guò)引入擴(kuò)展有限元方法，使用相場(chǎng)等高線(xiàn)的中線(xiàn)切割網(wǎng)格，從而創(chuàng)建了清晰的裂紋界面。本文算法使得在斷裂面上添加物理現(xiàn)象成為可能，可以緩解裂紋附近區(qū)的域網(wǎng)格畸變，且改善了矩陣的病態(tài)程度。

本文算法的基礎(chǔ)是計(jì)算相場(chǎng)等高線(xiàn)與輔助網(wǎng)格的交點(diǎn)，與材料本構(gòu)無(wú)關(guān)。因此，雖然本文中算例皆使用各項(xiàng)同性材料，但本方法亦適用于各向異性材料的情況。注意，欲使本算法成功運(yùn)行，輔助網(wǎng)格單元的尺寸須足夠大。一般地，輔助網(wǎng)格單元尺寸應(yīng)至少為相場(chǎng)等高線(xiàn)寬度的2倍。使用時(shí)，還應(yīng)盡量避免裂紋在輔助單元內(nèi)部極度彎曲的場(chǎng)景，因?yàn)榇藭r(shí)程序生成的相場(chǎng)中線(xiàn)線(xiàn)段可能部分位于相場(chǎng)等高線(xiàn)外部。數(shù)值算例表明，本文提出的中線(xiàn)生成算法可以成功捕捉到脆性材料典型斷裂形貌的相場(chǎng)中線(xiàn)，且使用相場(chǎng)中線(xiàn)切割網(wǎng)格不會(huì)對(duì)相場(chǎng)本身產(chǎn)生影響。