基于ABAQUS 平臺考慮T 應力的I 型裂紋擴展模擬開發

楊立云，韋 鵬，王青成，陳美霞，楊登輝

(中國礦業大學(北京)，力學與建筑工程學院，北京 100083)

隨著社會生產的發展，新工藝、新技術、高強度材料得到廣泛應用，工程部件和結構在服役期間或由于制造缺陷會產生微裂紋，之后由于應力集中、疲勞、腐蝕等因素的影響，微裂紋不斷成核、擴展、生長為宏觀裂紋，最終導致材料斷裂失效。斷裂準則是含裂紋材料在各種應力作用下斷裂和破壞的重要判據，解釋了裂紋何時起裂和判斷裂紋擴展方向的問題，因此在預測裂紋擴展路徑方面得到廣泛應用。傳統線彈性斷裂力學中常用的MTS[1]和SED[2]準則僅采用應力強度因子K表征裂紋尖端應力集中的程度[3]，但經過研究發現[4]，當裂紋受到高度約束時，使用單個參數(即臨界應力強度因子)來表征線彈性材料的脆性斷裂韌度是有缺陷的。SUN 等[5]研究發現邊界約束能顯著改變裂紋尖端應力場，尤其是由奇異應力場(K-主導區)主導的區域大小，隨著K-主導區尺寸的減小，奇異應力可能不足以描述裂紋的起裂，需要引入一個新的參數來表征裂紋尖端附近非奇異應力場對斷裂的影響。

國外研究者LIM 等[6]和KHAN 等[7]在分析線彈性材料的脆性斷裂時發現，不同幾何構型的試件通過斷裂實驗得到的材料斷裂韌度差別很大，而且應用傳統的斷裂準則計算的裂紋起裂角和擴展角與試驗結果有一定的偏差，傳統斷裂力學理論無法對這一現象作出解釋。CHAO 等[8]指出I 型裂紋擴展發生偏轉是由于裂紋尖端約束效應的作用，即試件厚度、平面幾何形狀和加載方式等與材料的斷裂行為的相關性，并用裂紋尖端Williams級數展開式的第2 項-平行于裂紋方向的常數項即T應力來表征約束效應，發現其對裂紋起裂和擴展的影響顯著。研究結果表明[8-9]：經過T應力修正后的斷裂準則能更準確地預測裂紋的偏轉角和擴展路徑。

II 型裂紋與I-II 復合型裂紋的擴展相對于I 型裂紋更加復雜，需要確定裂紋起裂角度以及方向。對于II 型裂紋，T應力對裂紋的斷裂韌度有較大影響，SMITH 等[10]基于理論與實驗研究了T應力對線彈性脆性材料II 型斷裂韌度的影響，結果表明：較大的負T應力使II 型斷裂韌度增加，T應力和斷裂韌性之間幾乎呈線性反比關系。

對于I-II 復合型裂紋，隨著T應力的增大，裂紋擴展角逐漸變大。高文和王生楠[11]將理論預測結果和試驗數據進行比較得出結論：負T應力可以提高裂紋擴展阻力，正T應力則可以降低材料的裂紋擴展阻力。

ABAQUS 是目前國際上通用的非線性有限元分析軟件之一，它提供了基于常規有限元法或擴展有限元法的斷裂力學參數計算方法，采用ABAQUS軟件內置的交互積分方法[12]可以計算裂紋尖端的應力強度因子KI、KII，J積分和T應力，在模擬裂紋起裂和擴展問題時具有可靠性。ABAQUS 中模擬裂紋擴展路徑的方法主要有：1) 增量裂紋擴展法(The incremental crack growth)[13]是在利用傳統有限元法求解斷裂力學問題時，將裂紋尖端鄰域網格單元進行細化，并設置為奇異性單元，模擬裂紋擴展時還需要對網格進行重劃分，耗時且效率低；2)擴展有限元法 (Extended finite element method, XFEM)[14]消除了裂紋對網格的依賴，可以自動模擬裂紋擴展路徑，但不能準確得到應力強度因子等斷裂力學參數；3)虛擬裂紋閉合法(Virtual crack closure technique, VCCT)[15]基于線彈性斷裂力學，具有對有限元網格尺寸不敏感和無需對裂紋尖端應力單元特殊處理等優點，分析相對穩定，但同樣無法求取準確的裂尖應力參數；4)離散元法[16](Discrete element method, DEM)作為分析不連續介質問題的數值方法，通過將顆粒相互黏結形成的不連續黏結體來模擬各種材料，可以較為真實地模擬非線性變形和裂紋擴展特征；5)近場動力學法[17](Peridynamics, PD)是一種基于非局部作用思想并以積分形式為運動方程的數值方法，解決了傳統連續介質力學在處理不連續問題時的奇異性問題；6)相場法[18](Phase field method,PDM)相對于描述離散裂紋的計算裂紋方法，無需額外的裂紋起裂準則，可以描述從無裂紋到裂紋起裂擴展的全過程。

ABAQUS 為用戶提供了兩種二次開發接口：1)使用Fortran 語言二次開發本構模型或單元[19-20]；2)使用Python 語言二次開發ABAQUS 腳本接口。很多國內外學者利用Python 語言對ABAQUS 進行二次開發，實現裂紋自動擴展數值模擬，SONG等[21]采用子模型技術二次開發ABAQUS 軟件，成功模擬了混凝土裂縫擴展過程；李強等[22]通過ABAQUS 腳本接口進行二次開發，對閉口斜裂紋的分支裂縫的擴展過程進行了數值模擬；WANG等[23]基于Python 程序進行ABAQUS 預處理二次開發，建立了中尺度輕骨料混凝土有限元模型，模擬結構中離散裂紋的起裂和擴展行為；王海鵬和吳永東[24]利用Python 語言編寫了裂紋自動擴展的程序模塊；竹宇波[25]利用圍道積分及網格重剖分技術，實現了裂紋自由擴展數值模擬功能，并對SCB 模型進行了應力強度因子及裂紋擴展路徑數值仿真。關于ABAQUS 裂紋擴展模塊的開發主要是采用最大周向應力準則(MTS)或最小應變能密度因子準則(SED)判斷裂紋的擴展方向，卻沒有考慮T應力的影響。

基于此，本文采用Python 語言編寫基于修正最大周向應力和最小應變能密度因子準則的程序模塊，將T應力引入裂紋起裂角的計算中，研究T應力對起裂角度的影響，并利用編寫的程序進行裂紋擴展路徑模擬，與相關文獻中含I 型裂紋的實驗結果進行對比，驗證程序的有效性。

1 理論分析

1.1 修正的最大周向應力準則

ERDOGAN 等[1]根據純II 型裂紋擴展實驗結果建立了最大周向應力準則，即在脆性材料中，裂紋將沿著裂尖周圍的最大周向應力方向θ0開裂。CHAO 等[8-9]提出了經過T應力修正的最大周向應力準則，修正后的準則考慮了K和常數項T應力的影響。GMTS(Generalized maximum tangential stress)準則可以表示如下：

求解式(1)和式(2)，兩個可能的結果可以寫為：

式中：Tc、Kc分別為斷裂時的臨界T應力和應力強度因子；rc為裂紋尖端的臨界距離或斷裂過程區域的半徑。

裂紋一旦發生彎曲，裂紋模式將由I 型轉變為I-II 復合型，此時應力強度因子KI、KII為應力場的主項，因此不考慮高階項T應力的影響。根據MTS 準則并且僅考慮Williams 級數展開的第一奇異項時，裂紋擴展角θ0可以通過以下公式求得：

1.2 修正的最小應變能密度準則

SIH[2]于1973 年提出了基于應變能密度場的斷裂概念，得到了應變能密度因子準則，該準則在預測含裂紋結構的疲勞和斷裂方面得到廣泛應用。SED 準則認為當以裂尖為圓心、臨界距離rc為半徑的小圓周上的應變能密度因子達到臨界值Sc時裂紋失穩擴展，裂紋沿最小的應變能密度因子方向擴展，記為θm。rc和Sc均為材料常數。SED 準則可表示為：

式中，S是應變能密度因子，表達式為：

式中： dW/dV為單位面積的應變能；rc為臨界半徑；ν、E分別為材料的泊松比和彈性模量；σx、σy、σz為空間直角坐標系中的正應力分量；τxy、τxz、τyz為空間直角坐標系中的剪切應力分量。將包含常數項T應力的裂紋尖端應力場表達式代入式(5)～式(7)中，得到GSED(Generalized Minimum Strain Energy Density)準則，經求導化簡后其表達式為：

其中：

2 裂紋擴展模擬方法

ABAQUS 作為常用的非線性有限元分析軟件之一，其內置程序為二次開發用戶提供了豐富的數據庫，通過Python 語言來調用這些數據庫可以實現從建模到后處理等過程的全程控制，還可以開發用戶自定義界面以實現交互式操作等功能。

2.1 裂紋自動擴展程序流程

利用Python 語言對ABAQUS 進行二次開發，編寫裂紋自動擴展程序。設置裂紋擴展增量，自動提取新裂尖坐標，刪除上一步模型的網格，根據ABAQUS 軟件內置的交互積分方法[12]得到的裂尖應力強度因子和T應力計算裂紋起裂和擴展角度，自動生成新的模型以及裂尖網格細化，重復上述步驟，直至裂紋擴展到模型邊界，完成自動擴展程序模塊的編寫。程序運行流程圖如圖1 所示，圖中Cpd 為內置程序計算的裂紋起裂角和擴展角，Seam 為指定裂縫，Remeshing 為重新劃分網格。

裂紋自動擴展程序功能主要包含以下幾個部分：應力分析及裂紋擴展檢測、ODB 文件數據收集、更新模塊、網格重劃分、將結果寫入指定文件。建立模型，定義材料，施加荷載和邊界條件等裂紋體前處理操作均由ABAQUS/CAE 和ABAQUS/Standard 完成。運行腳本程序并讀取模型信息，自動獲取裂紋及其尖端位置，建立分析步、劃分網格后提交計算，根據計算結果和斷裂準則判斷裂紋是否擴展，如果不擴展則終止程序；如果擴展則根據斷裂準則模塊計算的角度和給定的擴展增量得到新裂紋尖端的坐標，返回定義裂紋的步驟并重新劃分網格，往復循環，直至裂紋擴展到模型邊界或止裂時程序終止。通過導入CAE 文件建立初始模型進行分析，是為了拓寬程序模塊的使用范圍，例如對含中心斜裂紋、多條裂紋等試件的模擬，只需提前在ABAQUS 中建立前處理模型，再運行腳本程序即可實現對其裂紋擴展路徑的分析。

2.2 粒子群算法

考慮T應力對裂尖應力場的影響，對最大MTS或最小SED 準則進行修正，需要利用Python 語言將修正后的準則編入裂紋自動擴展程序腳本，判斷裂紋是否起裂。這涉及利用Python 語言求解三角函數方程的問題，由于式(8)所示的方程沒有解析解，Python 中的科學計算庫Sympy 模塊無法解決，本文通過Python 調用Matlab 中的函數解決方程問題，再將計算結果返回原Python 程序，此種數據轉換方法的運算效率非常低。

1997 年KENNEDY[26]根據對鳥群覓食行為的研究和模擬，建立了基于群體智能的全局隨機搜索算法——粒子群算法(PSO)，主要用于解決最優化問題，PSO 算法可以描述為n個粒子在N 維搜索空間里，通過不斷更新速度和位置來尋找最優解的過程。在N 維空間中有n個粒子，用vi為粒子的速度，xi為粒子的位置，速度和位置更新公式如下：

速度更新公式為：

位置更新公式為：

式中：w為慣性權重(有速度就有運動慣性)；pbest為本粒子歷史上最好的位置；gbest 為種群中所有粒子當前最好的位置；c1和c2為學習因子，分別向本粒子歷史最好位置和種群當前最好位置進行學習；r1和r2為兩個隨機函數，取值范圍[0, 1]，以增加更新速度的隨機性。

PSO 算法具有較強的函數極值求解能力，其算法流程如圖2 所示，首先定義粒子速度和位置的取值范圍并隨機初始化每個粒子，根據適應度函數(需要求解極值的函數方程)計算粒子適應度值，通過粒子間的相互搜索、學習得到個體極值和群體極值；利用式(10)～式(11)更新粒子的速度和位置，并相應地更新個體極值和群體極值，反復迭代，直至輸出全局最優結果并結束程序。

圖2 粒子群算法流程圖Fig.2 Flow chart of particle swarm algorithm

為了求解式(8)所示的三角函數方程，本文適應度函數應為其微分后得到的原函數，適應度值為其函數值，種群粒子數為20，每個粒子的維數為2，算法迭代進化次數為100，輸入相關參數后利用Python 編寫的PSO 算法進行計算，得到的關于原函數的最優解(極值點)即為三角函數方程的數值解，經過換算可得裂紋的起裂角或擴展角。

2.3 程序編制

利用Python 將修正后的MTS 或SED 準則編入裂紋自動擴展程序腳本，基于智能群體全局隨機搜索的粒子群算法，計算裂紋起裂角度和擴展角度，部分程序如下：

3 應用驗證

AYATOLLAHI 等[13]在室溫靜態荷載下利用含I 型裂紋的雙懸臂梁(Double cantilever beam, DCB)試件進行了斷裂試驗并給出了擴展路徑，試件模型如圖3 所示，其材質為聚甲基丙烯酸甲酯(PMMA 或有機玻璃)，該試件模型的斷裂韌度為1.3 MPa·m0.5。本節將利用ABAQUS 模擬相同條件下的裂紋擴展路徑，使用增量裂紋擴展法自動迭代處理，并將模擬結果與文獻中兩種尺寸的DCB 試樣實驗結果進行對比，以此驗證程序的有效性。

圖3 試樣幾何形狀Fig.3 Specimen geometry

在文獻[13]中初始裂紋長度與模型長度比值a/W設定為0.5，彈性模量E=2.9 GPa，泊松比ν=0.35，臨界半徑rc=0.1 mm，模型其他參數見表1。研究發現[27]裂紋偏轉具有隨機性，可能會向上或者向下偏轉，如圖3 所示。由于參考的文獻實驗結果是向下擴展，所以本文模擬選取向下的偏轉角度。

表1 模型相關信息[13]Table 1 Model related information[13]

由于試樣的厚度(10 mm)相對于其他尺寸比例較大，應采用平面應變條件進行建模。模型圓孔中心采用集中力形式進行約束，在ABAQUS中實現的具體操作：在圓孔中心設置參考點，在相互作用模塊中將參考點和圓孔內邊界進行耦合(耦合形式為Kinematic Coupling)，在加載模塊中在參考點處施加集中力荷載以及邊界條件(其中在上圓孔的參考點處約束U1 方向的位移/轉角，在下圓孔的參考點處約束U1、U2 方向的位移/轉角)。裂尖網格采用Isoparametric 8-node quadrilateral elements(CPE8R)單元，并將裂尖周圍單元設置為奇異性單元，以滿足裂紋尖端應力與應變的奇異性要求[13]，模型有限元網格如圖4 所示。最后進行了收斂性研究，以確認在有限元建模中使用了適當數量的元素。

圖4 模型有限元網格Fig.4 Model Finite Element Mesh

裂紋以設置的1 mm 的增量長度進行擴展，為了判斷裂紋的擴展方向，采用上述修正準則計算每一增量步的偏轉角度：① 對于GMTS 準則，將KI、KII和T應力的數值代入式(3)和式(4)即可求得下一增量步的裂紋擴展角；② 對于GSED 準則，將數值結果代入式(8)即可計算裂紋擴展角。

裂紋在不同擴展階段下(選取擴展步1、6、11、15)的Mises 應力云圖如圖5 所示(以DCB2 為例)。圖6～圖7 分別給出了兩種模型的GMTS 和GSED 準則預測的理論裂紋擴展路徑與實驗結果之間的對比，具體結果見表2。由對比結果可知，兩種模型的對比結果相似，其中GMTS 準則未能預測實驗中發生的裂紋偏轉，而是沿初始直線路徑擴展，與文獻[13]中的計算結果吻合?；贕SED準則的實驗結果與理論預測結果之間可以觀察到良好的一致性，裂紋的整個擴展過程與試驗數據基本符合，證明了程序的有效性。

表2 理論結果與實驗結果對比Table 2 Comparison of theoretical and experimental results

圖6 DCB1 裂紋擴展路徑對比圖Fig.6 DCB1 Comparison of crack propagation paths

圖7 DCB2 裂紋擴展路徑對比圖Fig.7 DCB2 Comparison of crack propagation paths

進一步研究發現上述不同準則預測裂紋擴展路徑存在差異的主要原因是SED 準則考慮了裂紋尖端附近6 個應力分量的作用，計算裂紋尖端附近局部的應變能密度，建立裂紋失穩判據。GSED準則相對于GMTS 準則，還考慮泊松比的影響。AYATOLLAHI 等[28]基于式(8)給出了斷裂起始角θ0相對于Bα 的變化規律(ν=0.35)，如圖8 所示。對于GMTS 準則，SMITH 等[29]在考慮Bα 對裂紋起始角θ0的影響時得出了類似的結論并繪制于圖8，當Bα＞0.375 時(Bα 見式(8))，I 型裂紋彎曲；當0.22＜Bα＜0.375 時，GMTS 準則無法預測裂紋的偏轉；而GSED 準則可以預測在此Bα 范圍內裂紋的擴展路徑。

圖8 斷裂起始角θ0 相對于Bα 的變化[[28-29]]Fig.8 Variation of fracture initiation angle θ0 relative to Bα[[28-29]]

4 結論

本文采用Python 語言對ABAQUS 進行二次開發，將考慮T應力的GMTS、GSED 準則引入裂紋起裂角的計算中并研究T應力對起裂角度的影響，將模擬結果與相關文獻中含I 型裂紋的實驗結果進行對比，得到以下結論：

(1)通過ABAQUS 二次開發接口，利用Python語言編寫裂紋自動擴展程序，實現了模擬裂紋擴展過程中的自動化功能，并將考慮T應力的GSED和GMTS 準則編入程序腳本中，實現多種斷裂準則對裂紋擴展路徑的模擬。

(2)分別利用GSED 和GMTS 斷裂準則研究了I 型載荷下裂紋擴展路徑的差異，并將數值模擬結果與相關文獻結果進行對比，其中GSED 準則能較好地預測該模式下的裂紋擴展路徑，表明了程序的有效性。對于GMTS 準則，當Bα＞0.375 時才能預測裂紋的彎曲現象。