基于二次插值重構有限元法的動態裂紋擴展模擬

王玉光 +李忠文 吳圣川 臧曉蕾 王燕

摘要：

基于二次插值重構有限元法（Twiceinterpolation Finite Element Method， TFEM）分析動態斷裂力學問題并進行數值實驗，考察TFEM在裂紋動態擴展模擬中的準確性和可靠性.由于TFEM保證節點處梯度場的連續性，因此裂尖附近的應力場可以得到較好的逼近.把該算法成功移植到自主開發的三維裂紋擴展仿真軟件（ZonCrack）中.利用ZonCrack進行的裂紋擴展，分析結果表明：TFEM得到裂尖應力強度因子（Stress Intensity Factor， SIF）與解析解基本一致；裂紋擴展的模擬結果與實驗值吻合良好.

關鍵詞：

插值重構有限元； 斷裂力學； 疲勞裂紋擴展； 裂紋尖端場； 損傷容限； ZonCrack軟件

中圖分類號： TB115；O242.21

文獻標志碼： B

0 引 言

動態疲勞裂紋的模擬分析對工程結構服役安全性和剩余壽命評價有重要的工程指導意義.[1]為模擬裂紋擴展，一般采用節點松弛方法，即在計算中把節點沿著裂縫路徑一分為二.[2]由于這種技術適用性差，學者們便又提出一些新型計算方法，如無網格[3]、邊界單元法[4]，以及有效的單元分割策略[5].這些方法雖然可在一定程度上提高計算精度，但計算效率卻降低，顯然不利于工程應用.

為此，ZHENG等[6]提出一種二次插值重構有限元法（Twiceinterpolation Finite Elements Method， TFEM）求解線彈性斷裂力學問題.為提高缺陷體的分析精度，通過2次連續插值得到形函數，其主要特點是把標準有限元法（Finite Elements Method， FEM）插值函數作為第二階段的插值參數，或者說把傳統FEM不連續的節點梯度（節點導數）作為插值參數，這樣，新的形函數便具有連續的節點梯度場.結果發現，TFEM比經典FEM有更好的計算精度和收斂性.

該算法植入新一代基于自適應加密的擴展有限元（eXtend FEM，XFEM）的斷裂仿真軟件ZonCrack.ZonCrack基于概率疲勞斷裂力學理論，對含裂紋工程的結構進行裂紋擴展仿真，獲取無損探傷數據，為工程師提供科學的基礎數據和判據.該算法能夠預測靜載荷或疲勞載荷作用下的開裂行為，確定工程結構損傷容限，不僅可以進行傳統有限元強度分析，而且提供疲勞斷裂及擴展壽命的分析功能.

1 基于平面三角形的TFEM

選取對復雜幾何形狀適應性好的三角形單元構建插值函數以求解斷裂力學問題.

1.1 一次插值函數

給定一個典型的常應變的單純形三角形，單元由節點i，j，m逆時針排列構成.插值重構法在二維單純形單元中的示意見圖1.

如圖1所示，Si，Sj和Sm為有公共節點i，j和k的關系單元集合，組成支撐域.對u（x）有貢獻的節點集（無單元法中一般稱為支撐點）正是Si，Sj和Sm所涵蓋的全部節點.設qs為支撐節點的位移向量，

1.2 二次插值函數

根據經典FEM得到節點導數后，再進入二次插值獲取一個新的形函數

1.3 關于C0節點

從式（16）看出，二次插值得到的試函數在節點處具有C1連續性.這一特點具有提高計算精度的優勢，但位于材料界面和位移邊界上的節點給求解帶來困難.為此，設節點i具有C0連續性，則對于插值點x位于單元e內的節點，可令

分析發現，當問題域內所有節點進行C0連續性處理后，

TFEM型函數將退化為經典FEM型函數，即在某些情況下2種算法可以耦合.也就是說，對于含缺陷體需要更高插值精度的局部區域可采用TFEM，而在遠離裂紋前緣區則可采用一般精度FEM，從而實現計算精度和效率的平衡考慮.

2 裂紋擴展計算關鍵技術

發展TFEM方法用于求解斷裂力學問題，并植入ZonCrack軟件中.專業斷裂仿真軟件ZonCrack是一個面向動態斷裂過程的仿真實驗平臺，其以成熟的疲勞斷裂力學為理論基礎，對含缺陷工程部件進行仿真模擬.[7]為減少對經典FEM源代碼的改動，引入節點松弛技術模擬材料的劈開過程.

2.1 應力強度因子

采用區域互作用積分計算混合型裂紋的應力強度因子（Stress Intensity Factor， SIF）解[8]，以最大周向應力判斷裂紋方向[9]，則極坐標形式的裂尖應力場分量可有I型和II型裂紋應力強度因子表示

2.2 節點投影線技術

研究發現，網格質量對插值重構有限元法的計算精度影響較小（網格敏感性低），也就是說調整裂紋面附近的網格對SIF的影響也較小.

當前求解斷裂力學問題的通用處理方法是確保裂紋面/線與單元邊/節點的重合.本文采用節點投影方法準確模擬裂紋的擴展路徑，具體操作過程為：網格整體上不調整，僅對裂紋附近節點投影到裂紋面/線上（見圖2）.一般來說，裂紋與三角元有3種關系：第一，裂紋與2個邊相交，把相交邊的公共節點移至裂紋；第二，裂紋與1個邊相交，把離裂紋近的節點移至裂紋；第三，裂尖位于單元內，則移動單

元上距裂紋尖端最近的節點.

采用上述處理后，問題域內所有單元邊/節點都會位于裂紋面/線上.自此，采用經典節點分離法人為把裂紋面/線上的節點一分為二，并賦予不同節點號.通過這樣處理，裂紋擴展路徑上就出現一個強不連續性位移，即實現裂紋張開模擬.

3 算例驗證

為考察新型二次插值有限元的求解精度和穩定性，定義2范數位移誤差和能量誤差為

式中：下標e表示精確解或者理論解；下標n表示數值方法參考解或者TFEM解.在TFEM實現中，單元剛度矩陣積分采用4點Hammer積分方法.

在ZonCrack中調用TFEM求解器分析一些典型的斷裂力學問題，以測試TFEM法準確性與合理性.文中所用變量均采用國際單位制.

3.1 單邊直裂紋

一個單邊直裂紋的無限大平板，兩端受均勻拉力，假設平面應變狀態.模型的幾何及加載見圖3，且有H=L=P=1.力學性能參數E=1010，ν=0.25.解析解為

應力強度因子解對比曲線見圖4.由此可知，對于不同擴展長度的裂紋，TFEM數值解與理論解和經典FEM解吻合良好，且TFEM法的計算精度高于傳統FEM法.TFEM與傳統FEM的裂紋尖端應力場分布下見圖5.由此可見，TFEM結果與理論解給出的裂紋簡單應力場吻合更好，而傳統有限元法卻存在著明顯的應力跳躍現象（即不連續）.

3.2 中心斜裂紋

為進一步考察TFEM求解SIF的準確性與可靠性，以中心斜裂紋板為例研究，見圖6.在無限大平板x方向和y方向分別搜到σ1和σ2大小的拉應力，β為斜裂紋與x軸的夾角，顯然，KI和KII為關于裂紋傾角β的函數，理論解為

有限元計算選取的計算域為邊長L=10的方板，所含裂紋的長度a=0.25，假設軸向拉力P=1.考慮到裂紋長度遠小于矩形板的長度，所以可以視為無限大平板的中心斜裂紋問題.

3.3 4點剪切含裂紋矩形板

4點剪切受載問題的裂紋擴展研究較多，幾何模型及受載情況見圖7，L=4，b=1，P=1.材料參數為E=200 MPa，ν=0.25，假設平面應變條件.第12擴展步后的裂紋路徑見圖8.由此可見，插值重構法的計算結果與實驗基本吻合.

4 結 論

將新型TFEM植入高級三維裂紋擴展仿真與分析系統ZonCrack中用于求解動態斷裂問題.通過標準算例分析，得出以下結論.

（1）在ZonCrack中調用TFEM求解器，能夠獲得與理論解可比的I型、II型及混合型SIF解，且在同一單元網格下，計算精度高于傳統FEM.

（2）與傳統FEM解相比較，TFEM可在裂尖區域得到更加準確和光滑連續的應力分布.

（3）TFEM法預測的裂紋路徑與實測解基本一致，表明TFEM算法的可靠性.

必須指出，與傳統有限元比較，雖然插值重構有限元法使用經典FEM的導數進行二次插值，但并未因此顯著增加問題的總求解自由度.

參考文獻：

[1] 楊小彬， 莊茁， 莊傳晶. 富氣輸送管道裂紋動態擴展的數值模擬[J]. 清華大學學報（自然科學版）， 2008， 48（8）： 13551358. DOI： 10.3321/j.issn：10000054.2008.08.032.

YANG X B， ZHUANG Z， ZHUANG C J. Numerical study of crack propagation in rich gas transporting pipelines[J]. Journal of Tsinghua University （Science and Technology）， 2008， 48（8）： 13551358. DOI： 10.3321/j.issn：10000054.2008.08.032.

[2] 解德， 錢勤， 李長安. 斷裂力學中的數值計算方法及工程應用[M]. 北京： 科學出版社， 2009.

[3] BELYTSCHKO T， LU Y Y， GU L. Crack propagation by elementfree Galerkin methods[J]. Engineering Fracture Mechanics， 1995， 51（2）： 295315. DOI： 10.1016/00137944（94）001539.

[4] 李俊， 馮偉哲， 高效偉. 一種基于直接計算高階奇異積分的斷裂力學雙邊界積分方程分析法[J]. 力學學報， 2016， 48（2）： 387398.DOI： 10.6052/0459187915342.

LI J， FENG W Z， GAO X W. A dual boundary integral equation method based on direct evaluation of higher order singular integral for crack problems[J]. Chinese Journal of Theoretical and Applied Mechanics， 2006， 48（2）： 387398.DOI： 10.6052/0459187915342.

[5] MELENK J M， BABUS K. The partition of unity finite element method： basic theory and applications[J]. Computer Methods in Applied Mechanics and Engineering， 1996， 139（14）： 289314. DOI： 10.1016/S00457825（96）010870.

[6] ZHENG C， WU S C， TANG X H， et al. A novel twice interpolation finite element method for solid mechanics problems[J]. Acta Mechanica Sinica， 2010， 26（2）： 265278. DOI： 10.1007/s1040900902653.

[7] 吳圣川， 吳玉程. 斷裂分析及CAE軟件的現狀與發展[J]. 計算機輔助工程， 2011， 20（1）： 12.

WU S C， WU Y C. Current status and development of fracture mechanics based CAE software[J]. Computer Aided Engineering， 2011， 20（1）： 12.

[8] BUDIANSKY B， RICE J R. Conservation laws and energyrelease rates[J]. Journal of Applied Mechanics， 1973， 40（2）： 201203. DOI： 10.1115/1.3422926.

[9] FLEMING M， CHU Y A， MORAN B， et al. Enriched elementfree Galerkin methods for cracktip fields[J]. International Journal for Numerical Methods in Engineering， 1997， 40（8）： 14831504. DOI： 10.1002/（SICI）10970207（19970430）40：8<1483：：AIDNME123>3.0.CO；26.

（編輯 武曉英）