基于二次插值重構(gòu)有限元法的動(dòng)態(tài)裂紋擴(kuò)展模擬

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

摘要：

基于二次插值重構(gòu)有限元法（Twiceinterpolation Finite Element Method， TFEM）分析動(dòng)態(tài)斷裂力學(xué)問(wèn)題并進(jìn)行數(shù)值實(shí)驗(yàn)，考察TFEM在裂紋動(dòng)態(tài)擴(kuò)展模擬中的準(zhǔn)確性和可靠性.由于TFEM保證節(jié)點(diǎn)處梯度場(chǎng)的連續(xù)性，因此裂尖附近的應(yīng)力場(chǎng)可以得到較好的逼近.把該算法成功移植到自主開(kāi)發(fā)的三維裂紋擴(kuò)展仿真軟件（ZonCrack）中.利用ZonCrack進(jìn)行的裂紋擴(kuò)展，分析結(jié)果表明：TFEM得到裂尖應(yīng)力強(qiáng)度因子（Stress Intensity Factor， SIF）與解析解基本一致；裂紋擴(kuò)展的模擬結(jié)果與實(shí)驗(yàn)值吻合良好.

關(guān)鍵詞：

插值重構(gòu)有限元； 斷裂力學(xué)； 疲勞裂紋擴(kuò)展； 裂紋尖端場(chǎng)； 損傷容限； ZonCrack軟件

中圖分類號(hào)： TB115；O242.21

文獻(xiàn)標(biāo)志碼： B

0 引 言

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

為此，ZHENG等[6]提出一種二次插值重構(gòu)有限元法（Twiceinterpolation Finite Elements Method， TFEM）求解線彈性斷裂力學(xué)問(wèn)題.為提高缺陷體的分析精度，通過(guò)2次連續(xù)插值得到形函數(shù)，其主要特點(diǎn)是把標(biāo)準(zhǔn)有限元法（Finite Elements Method， FEM）插值函數(shù)作為第二階段的插值參數(shù)，或者說(shuō)把傳統(tǒng)FEM不連續(xù)的節(jié)點(diǎn)梯度（節(jié)點(diǎn)導(dǎo)數(shù)）作為插值參數(shù)，這樣，新的形函數(shù)便具有連續(xù)的節(jié)點(diǎn)梯度場(chǎng).結(jié)果發(fā)現(xiàn)，TFEM比經(jīng)典FEM有更好的計(jì)算精度和收斂性.

該算法植入新一代基于自適應(yīng)加密的擴(kuò)展有限元（eXtend FEM，XFEM）的斷裂仿真軟件ZonCrack.ZonCrack基于概率疲勞斷裂力學(xué)理論，對(duì)含裂紋工程的結(jié)構(gòu)進(jìn)行裂紋擴(kuò)展仿真，獲取無(wú)損探傷數(shù)據(jù)，為工程師提供科學(xué)的基礎(chǔ)數(shù)據(jù)和判據(jù).該算法能夠預(yù)測(cè)靜載荷或疲勞載荷作用下的開(kāi)裂行為，確定工程結(jié)構(gòu)損傷容限，不僅可以進(jìn)行傳統(tǒng)有限元強(qiáng)度分析，而且提供疲勞斷裂及擴(kuò)展壽命的分析功能.

1 基于平面三角形的TFEM

選取對(duì)復(fù)雜幾何形狀適應(yīng)性好的三角形單元構(gòu)建插值函數(shù)以求解斷裂力學(xué)問(wèn)題.

1.1 一次插值函數(shù)

給定一個(gè)典型的常應(yīng)變的單純形三角形，單元由節(jié)點(diǎn)i，j，m逆時(shí)針排列構(gòu)成.插值重構(gòu)法在二維單純形單元中的示意見(jiàn)圖1.

如圖1所示，Si，Sj和Sm為有公共節(jié)點(diǎn)i，j和k的關(guān)系單元集合，組成支撐域.對(duì)u（x）有貢獻(xiàn)的節(jié)點(diǎn)集（無(wú)單元法中一般稱為支撐點(diǎn)）正是Si，Sj和Sm所涵蓋的全部節(jié)點(diǎn).設(shè)qs為支撐節(jié)點(diǎn)的位移向量，

1.2 二次插值函數(shù)

根據(jù)經(jīng)典FEM得到節(jié)點(diǎn)導(dǎo)數(shù)后，再進(jìn)入二次插值獲取一個(gè)新的形函數(shù)

1.3 關(guān)于C0節(jié)點(diǎn)

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

分析發(fā)現(xiàn)，當(dāng)問(wèn)題域內(nèi)所有節(jié)點(diǎn)進(jìn)行C0連續(xù)性處理后，

TFEM型函數(shù)將退化為經(jīng)典FEM型函數(shù)，即在某些情況下2種算法可以耦合.也就是說(shuō)，對(duì)于含缺陷體需要更高插值精度的局部區(qū)域可采用TFEM，而在遠(yuǎn)離裂紋前緣區(qū)則可采用一般精度FEM，從而實(shí)現(xiàn)計(jì)算精度和效率的平衡考慮.

2 裂紋擴(kuò)展計(jì)算關(guān)鍵技術(shù)

發(fā)展TFEM方法用于求解斷裂力學(xué)問(wèn)題，并植入ZonCrack軟件中.專業(yè)斷裂仿真軟件ZonCrack是一個(gè)面向動(dòng)態(tài)斷裂過(guò)程的仿真實(shí)驗(yàn)平臺(tái)，其以成熟的疲勞斷裂力學(xué)為理論基礎(chǔ)，對(duì)含缺陷工程部件進(jìn)行仿真模擬.[7]為減少對(duì)經(jīng)典FEM源代碼的改動(dòng)，引入節(jié)點(diǎn)松弛技術(shù)模擬材料的劈開(kāi)過(guò)程.

2.1 應(yīng)力強(qiáng)度因子

采用區(qū)域互作用積分計(jì)算混合型裂紋的應(yīng)力強(qiáng)度因子（Stress Intensity Factor， SIF）解[8]，以最大周向應(yīng)力判斷裂紋方向[9]，則極坐標(biāo)形式的裂尖應(yīng)力場(chǎng)分量可有I型和II型裂紋應(yīng)力強(qiáng)度因子表示

2.2 節(jié)點(diǎn)投影線技術(shù)

研究發(fā)現(xiàn)，網(wǎng)格質(zhì)量對(duì)插值重構(gòu)有限元法的計(jì)算精度影響較小（網(wǎng)格敏感性低），也就是說(shuō)調(diào)整裂紋面附近的網(wǎng)格對(duì)SIF的影響也較小.

當(dāng)前求解斷裂力學(xué)問(wèn)題的通用處理方法是確保裂紋面/線與單元邊/節(jié)點(diǎn)的重合.本文采用節(jié)點(diǎn)投影方法準(zhǔn)確模擬裂紋的擴(kuò)展路徑，具體操作過(guò)程為：網(wǎng)格整體上不調(diào)整，僅對(duì)裂紋附近節(jié)點(diǎn)投影到裂紋面/線上（見(jiàn)圖2）.一般來(lái)說(shuō)，裂紋與三角元有3種關(guān)系：第一，裂紋與2個(gè)邊相交，把相交邊的公共節(jié)點(diǎn)移至裂紋；第二，裂紋與1個(gè)邊相交，把離裂紋近的節(jié)點(diǎn)移至裂紋；第三，裂尖位于單元內(nèi)，則移動(dòng)單

元上距裂紋尖端最近的節(jié)點(diǎn).

采用上述處理后，問(wèn)題域內(nèi)所有單元邊/節(jié)點(diǎn)都會(huì)位于裂紋面/線上.自此，采用經(jīng)典節(jié)點(diǎn)分離法人為把裂紋面/線上的節(jié)點(diǎn)一分為二，并賦予不同節(jié)點(diǎn)號(hào).通過(guò)這樣處理，裂紋擴(kuò)展路徑上就出現(xiàn)一個(gè)強(qiáng)不連續(xù)性位移，即實(shí)現(xiàn)裂紋張開(kāi)模擬.

3 算例驗(yàn)證

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

式中：下標(biāo)e表示精確解或者理論解；下標(biāo)n表示數(shù)值方法參考解或者TFEM解.在TFEM實(shí)現(xiàn)中，單元?jiǎng)偠染仃嚪e分采用4點(diǎn)Hammer積分方法.

在ZonCrack中調(diào)用TFEM求解器分析一些典型的斷裂力學(xué)問(wèn)題，以測(cè)試TFEM法準(zhǔn)確性與合理性.文中所用變量均采用國(guó)際單位制.

3.1 單邊直裂紋

一個(gè)單邊直裂紋的無(wú)限大平板，兩端受均勻拉力，假設(shè)平面應(yīng)變狀態(tài).模型的幾何及加載見(jiàn)圖3，且有H=L=P=1.力學(xué)性能參數(shù)E=1010，ν=0.25.解析解為

應(yīng)力強(qiáng)度因子解對(duì)比曲線見(jiàn)圖4.由此可知，對(duì)于不同擴(kuò)展長(zhǎng)度的裂紋，TFEM數(shù)值解與理論解和經(jīng)典FEM解吻合良好，且TFEM法的計(jì)算精度高于傳統(tǒng)FEM法.TFEM與傳統(tǒng)FEM的裂紋尖端應(yīng)力場(chǎng)分布下見(jiàn)圖5.由此可見(jiàn)，TFEM結(jié)果與理論解給出的裂紋簡(jiǎn)單應(yīng)力場(chǎng)吻合更好，而傳統(tǒng)有限元法卻存在著明顯的應(yīng)力跳躍現(xiàn)象（即不連續(xù)）.

3.2 中心斜裂紋

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

有限元計(jì)算選取的計(jì)算域?yàn)檫呴L(zhǎng)L=10的方板，所含裂紋的長(zhǎng)度a=0.25，假設(shè)軸向拉力P=1.考慮到裂紋長(zhǎng)度遠(yuǎn)小于矩形板的長(zhǎng)度，所以可以視為無(wú)限大平板的中心斜裂紋問(wèn)題.

3.3 4點(diǎn)剪切含裂紋矩形板

4點(diǎn)剪切受載問(wèn)題的裂紋擴(kuò)展研究較多，幾何模型及受載情況見(jiàn)圖7，L=4，b=1，P=1.材料參數(shù)為E=200 MPa，ν=0.25，假設(shè)平面應(yīng)變條件.第12擴(kuò)展步后的裂紋路徑見(jiàn)圖8.由此可見(jiàn)，插值重構(gòu)法的計(jì)算結(jié)果與實(shí)驗(yàn)基本吻合.

4 結(jié) 論

將新型TFEM植入高級(jí)三維裂紋擴(kuò)展仿真與分析系統(tǒng)ZonCrack中用于求解動(dòng)態(tài)斷裂問(wèn)題.通過(guò)標(biāo)準(zhǔn)算例分析，得出以下結(jié)論.

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

（2）與傳統(tǒng)FEM解相比較，TFEM可在裂尖區(qū)域得到更加準(zhǔn)確和光滑連續(xù)的應(yīng)力分布.

（3）TFEM法預(yù)測(cè)的裂紋路徑與實(shí)測(cè)解基本一致，表明TFEM算法的可靠性.

必須指出，與傳統(tǒng)有限元比較，雖然插值重構(gòu)有限元法使用經(jīng)典FEM的導(dǎo)數(shù)進(jìn)行二次插值，但并未因此顯著增加問(wèn)題的總求解自由度.

參考文獻(xiàn)：

[1] 楊小彬， 莊茁， 莊傳晶. 富氣輸送管道裂紋動(dòng)態(tài)擴(kuò)展的數(shù)值模擬[J]. 清華大學(xué)學(xué)報(bào)（自然科學(xué)版）， 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] 解德， 錢勤， 李長(zhǎng)安. 斷裂力學(xué)中的數(shù)值計(jì)算方法及工程應(yīng)用[M]. 北京： 科學(xué)出版社， 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ì)算高階奇異積分的斷裂力學(xué)雙邊界積分方程分析法[J]. 力學(xué)學(xué)報(bào)， 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軟件的現(xiàn)狀與發(fā)展[J]. 計(jì)算機(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.

（編輯 武曉英）