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

1.華中科技大學 數(shù)字制造裝備與技術國家重點實驗室，武漢 430074 2.鄭州大學 機械與動力工程學院，鄭州 450001

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

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

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

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

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

1 相場斷裂模型

1.1 裂紋正則化

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

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

滿足以上斷裂相場條件的標量場可通過求解以下邊值問題得到：

(1)

(2)

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

(3)

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

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

(4)

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

(5)

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

(6)

1.2 控制方程強形式

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

(7)

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

通過引入相場標量(,)對裂紋進行正則化描述，即由全域彌散分布的裂紋近似描述真實裂紋，彈性體自由能可近似表述為

(8)

(9)

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

()=()+()

(10)

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

(11)

=+

(12)

(13)

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

(14)

式中：為引入的一個極小值，用于保證材料完全破壞時單元剛度矩陣的非奇異性；()為彈性應變能退化函數(shù)，描述斷裂相場對彈性應變能拉伸部分的影響，需滿足：

(15)

通常選取

(16)

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

(17)

對式(17)進行位移場和相場的變分可得

(18)

相應邊值問題為

(19)

=-Δ

(20)

(21)

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

(22)

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

1.3 控制方程弱形式

忽略體力的作用，位移場的弱形式為

(23)

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

(24)

式中：為二階單位張量。

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

(25)

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

(26)

1.4 有限元離散

對位移場和相場進行有限元離散，可得位移場平衡方程和相場演化方程

(27)

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

(28)

(29)

2 子結構建模及損傷識別

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

2.1 子結構建模

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

=

(30)

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

(31)

式中：下標i、b分別表示該自由度對應節(jié)點位于子域內(nèi)部和邊界。由式(32)可得

(32)

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

(33)

定義以下矩陣和向量：

(34)

可得邊界節(jié)點響應的有限元列式：

=

(35)

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

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

圖2 周期性結構的子結構建模Fig.2 Substructuring of periodic structure

2.2 損傷識別

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

(36)

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

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

2.3 算法流程

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

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

對特征子域結構的位移和相場的剛度矩陣進行靜態(tài)凝聚，組裝結構整體剛度矩陣。

對特征子域結構進行離線典型工況測試，標定基于應變能的損傷識別閾值。

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

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

3 數(shù)值算例

考慮由方孔單胞周期性排列構成結構的斷裂模擬。首先，對方孔單胞進行離線測試，標定損傷識別的能量閾值；隨后，基于子結構建模和損傷識別對兩類周期性結構進行在線的相場斷裂模擬，并與常規(guī)相場法的模擬結果和計算規(guī)模進行比較。

3.1 方孔單胞離線測試

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

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

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

選擇拉伸和剪切加載過程中方孔單胞相場最大值分別為0.05、0.10、0.20時的彈性應變能作為損傷判定的閾值，并進行編號，如表1所示。對

圖5 單胞A測試工況下彈性應變能-位移曲線和裂紋軌跡Fig.5 Elastic energy-displacement curves and crack paths of cell A under testing loads

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

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

3.2 Ⅱ 型斷裂模擬

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

圖6 周期性結構的Ⅱ型斷裂問題Fig.6 Mode Ⅱ fracture of periodic structure

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

圖7 周期性結構Ⅱ型斷裂問題的載荷-位移曲線及裂紋擴展情況Fig.7 Load-displacement curves and crack propagation of mode Ⅱ fracturing of periodic structure

圖8 周期性結構Ⅱ型斷裂問題的最終裂紋軌跡Fig.8 Ultimate crack paths of mode Ⅱ fracturing of periodic structure

圖9 Ⅱ型斷裂問題過程中含裂紋子結構占比Fig.9 Proportion of substructures with cracks during mode Ⅱ fracturing

3.3 L形梁斷裂

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

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

圖11 單胞A組成的L形梁斷裂模擬載荷-位移曲線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

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

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

圖14 單胞B組成的L形梁斷裂模擬載荷-位移曲線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

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

4 結 論

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