張子瑜，郝林

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

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

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

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

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

1 擴展有限元基本算法

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

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

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

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

本文采用Miehe等提出的相場斷裂力學模型。通常情況下，可以通過將總勢能最小化，或使虛功率等于零，來獲得控制方程的強形式:

+-=0

(1)

式中：為總應(yīng)變能；為耗散能；為外力功。總勢能為位移和相場的函數(shù):

(2)

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

(3)

(4)

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

(5)

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

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

(6)

且滿足邊界條件

(7)

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

(8)

(9)

(10)

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

3 生成相場中線的方法

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

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

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

注意，相場等高線由常規(guī)的行進正方形(Marching Squares)算法生成。此外，使用排序算法將組成等高線的小線段沿著逆時針方向排列存儲。圖2為相場云圖及=0.9所對應(yīng)的等高線示意圖。值得一提的是，即使在行進正方形算法中將相場值設(shè)置為1，通常情況下也無法得到1條單一的等高線。這是因為，有限元方法的離散性使其無法捕捉到相場值嚴格等于1的尖點。

圖2 相場云圖及其對應(yīng)的等高線Fig.2 A phase-field and its isocontour

3.1 計算輔助網(wǎng)格單元的中線交點

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

輔助網(wǎng)格單元的邊和相場等高線只有一個交點。

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

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

相場中線與輔助網(wǎng)格交點的方向向量可基于相場等高線與輔助網(wǎng)格交點的切向量來計算。如前所述，由于組成相場等高線的小線段以逆時針排列存儲，故可將與輔助網(wǎng)格單元邊相交的相場等高線小線段的方向作為等高線交點的切向量。取第1個等高線交點和最后一個等高線交點的切向量，將二者平均后即可得到相場中線與輔助網(wǎng)格交點的方向向量。在圖3中，相場中線交點的方向向量即為0.5(-)，其中，為等高線交點1的切向量，為等高線交點4的切向量。注意，相場中線交點方向向量的方向與等高線交點1切向量的方向一致，因為后者的方向通常代表了裂紋擴展的方向。

輔助網(wǎng)格單元的邊和相場等高線有多個交點。

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

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

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

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

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

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

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

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

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

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

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

=max(,)

(11)

=

(12)

式中：=1或2；=-；、分別為中線線段兩端中線交點的坐標向量；為中線交點處的方向向量。若某中線交點的某一個中線線段的對齊指標最小，則此線段為唯一有效的中線線段。為了詳細解釋，考察圖6中的中線交點3。在所有可能的中線線段3-4、線段3-1和線段3-2中，線段3-4及其2個端點的方向向量的對齊程度最好。因此，線段3-4是中線交點3唯一有效的中線線段。

圖6 對齊算法示意圖Fig.6 Alignment algorithm

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

3.3 添加新的中線線段

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

圖7 生成新中線線段的過程Fig.7 Process of computing new medial segments

4 數(shù)值算例

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

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

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

4.2 動態(tài)等高線的中線

4.2.1 裂紋融合

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

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

(13)

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

(14)

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

圖9 二維融合裂紋的相場云圖及其中線Fig.9 Phase-field contour and medial-axes of 2D merging cracks

4.2.2 裂紋分叉

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

圖10 二維分叉裂紋的相場云圖及其中線Fig.10 Phase-field contour and medial-axis of a 2D branching crack

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

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

圖12 裂紋擴展問題的位移-載荷曲線Fig.12 Load-displacement curves of single edge-notched tension test

5 結(jié) 論

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

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