， ， ，

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，工程力學(xué)系，大連 116024)

1 引 言

裂紋擴(kuò)展是很多工程結(jié)構(gòu)破壞失效的主要形式之一。對(duì)該過程實(shí)施有效的數(shù)值模擬與分析是研究斷裂機(jī)理乃至防止斷裂破壞事故發(fā)生的主要手段之一。目前，Belytschko等[1,2]提出的擴(kuò)展有限元法是模擬裂紋擴(kuò)展的主要方法之一。該方法在有限元法的基礎(chǔ)上，通過向位移插值空間中引入描述裂紋間斷的Heaviside函數(shù)和描述裂尖奇異性的精化(增強(qiáng))函數(shù)，理性地實(shí)現(xiàn)了裂紋的精確描述，日益成為裂紋問題數(shù)值方法的研究熱點(diǎn)[3-5]。

以無單元伽遼金法EFG(Element-free Galerkin)[6]為代表的無網(wǎng)格法[7]也是裂紋問題數(shù)值分析的主要方法之一。與有限元法相比，無網(wǎng)格法的插值(近似)函數(shù)不依賴于網(wǎng)格單元，因而具有易于形成高階近似以及易于實(shí)現(xiàn)局部節(jié)點(diǎn)加密等優(yōu)勢(shì)，給裂紋的數(shù)值分析帶來便利。實(shí)際上，Belytschko等[6]在提出EFG方法時(shí)就已經(jīng)利用其易于局部加密節(jié)點(diǎn)的特性有效改善了應(yīng)力強(qiáng)度因子的計(jì)算精度；隨后，EFG方法用于穩(wěn)態(tài)和動(dòng)態(tài)裂紋擴(kuò)展模擬[8,9]，可在很大程度上克服有限元法在模擬裂紋擴(kuò)展時(shí)需要不斷重新劃分網(wǎng)格的困難；Krysl等[10]將EFG方法應(yīng)用于三維動(dòng)態(tài)裂紋擴(kuò)展模擬；Fleming等[11]針對(duì)裂尖應(yīng)力場的奇異性提出了兩種精化(強(qiáng)化)函數(shù)，不僅減小了應(yīng)力場的虛假數(shù)值振蕩，而且提高了應(yīng)力強(qiáng)度因子的計(jì)算精度和效率；Duflot[12]通過向無網(wǎng)格法的權(quán)函數(shù)中引入精化函數(shù)，有效提高了裂尖應(yīng)力場的解析精度，并將該方法用于三維裂紋擴(kuò)展分析；近期，無網(wǎng)格法還通過耦合有限元法應(yīng)用于韌性斷裂分析[13]中。

上述工作大多采用了低階(線性)無網(wǎng)格近似，未能充分利用無網(wǎng)格法易于形成高階近似，從而有效改善應(yīng)力場求解精度的優(yōu)點(diǎn)；另一方面，無網(wǎng)格形函數(shù)為非多項(xiàng)式的有理函數(shù)，弱形式的數(shù)值積分須采用較多積分點(diǎn)，導(dǎo)致計(jì)算效率較低[14]，積分精度也不能令人滿意。由段慶林等[15,16]提出的基于混合變分原理發(fā)展的一致性無網(wǎng)格伽遼金法CEFG(Consistent EFG)通過修正積分點(diǎn)上形函數(shù)的空間導(dǎo)數(shù)，大幅度減少了所需積分點(diǎn)，并有效改善了計(jì)算精度，因而顯著提高了無網(wǎng)格法的計(jì)算效率，已在三維彈塑性[17]、自適應(yīng)分析[18]以及不可壓縮固體的變形分析[19,20]等問題中獲得成功應(yīng)用，展現(xiàn)出良好的發(fā)展?jié)摿?。?yīng)強(qiáng)調(diào)的是，二階CEFG方法在每個(gè)背景三角形積分子域內(nèi)僅使用三個(gè)積分點(diǎn)，但能精確通過線性和二次分片試驗(yàn)，尤其是能得到精度高、光滑無振蕩的應(yīng)力場?？紤]到應(yīng)力解析精度對(duì)于裂紋問題至關(guān)重要，高階CEFG方法對(duì)于裂紋分析具有顯著優(yōu)勢(shì)。

目前，CEFG方法還僅限于連續(xù)體問題。本文目的是將CEFG方法拓展到以裂紋為特征的非連續(xù)體問題。位移在裂紋處的間斷采用虛擬節(jié)點(diǎn)法(Phantom node method)[21]描述，無需增加額外的節(jié)點(diǎn)自由度。為進(jìn)一步簡化方法，與原始的虛擬節(jié)點(diǎn)法保持一致，本文也忽略了裂尖奇異場的加強(qiáng)函數(shù)，雖然這會(huì)導(dǎo)致裂尖解析精度有所降低，但是會(huì)大大簡化方法的實(shí)現(xiàn)，并便于模擬裂紋擴(kuò)展。這一策略也體現(xiàn)在近期無網(wǎng)格斷裂分析的研究中[22-24]。

2 裂紋的數(shù)值模型及空間離散

如圖1所示，考慮一含裂紋二維彈性體,所占據(jù)的區(qū)域?yàn)棣?則靜力平衡方程可寫為

·σ+b=0inΩ

(1)

式中σ為Cauchy應(yīng)力，b為體力。邊界條件有

(2)

(3)

σ·n=0onΓc+

(4)

σ·n=0onΓc-

(5)

為描述裂紋處位移場的不連續(xù)性,向位移插值空間中引入間斷的Heaviside函數(shù),因而位移近似為

(6)

式中NI(x)為移動(dòng)最小二乘(MLS)形函數(shù)，uI為節(jié)點(diǎn)位移參數(shù)，qI為附加自由度，H為單位階躍函數(shù)，即

(7)

圖1 含裂紋二維彈性體示意圖

Fig.1 Schematic diagram of a two -dimensional elastic body with a crack

式中f(x)=0表示裂紋面。式(6)右端加uINIH-uINIH,有

(8)

式中H=H(f(x))。為便于推導(dǎo)，將H(-f(xI))記為H-I,H(f(xI))記為H+I。由式(7)可知，當(dāng)H+I≠0時(shí)，H-HI=H-1；當(dāng)H-I≠0時(shí)，H-HI=H,則有

uINIHH+I+uINIHH-I+

qINI(H-1)H+I+qINIHH-I]

(9)

進(jìn)一步整理可得

uIH+INIH+(uI+qI)H-INIH]

(10)

根據(jù)虛擬節(jié)點(diǎn)法[21]，式(10)可寫為

(11)

式中

(12)

(13)

式(11)表明,含裂紋的位移插值可近似寫為兩部分的疊加，這一點(diǎn)可結(jié)合圖2說明。圖2的e0為一背景積分單元，實(shí)心圓點(diǎn)表示該單元關(guān)聯(lián)的計(jì)算節(jié)點(diǎn)(即節(jié)點(diǎn)形函數(shù)在該單元上不等于0的節(jié)點(diǎn))。若該單元發(fā)生斷裂，則可通過引入虛擬節(jié)點(diǎn)(圖2的空心圓點(diǎn))，將該單元拆分成e1和e2兩個(gè)獨(dú)立的單元。式(11)的S1與S2分別表示這兩個(gè)單元的關(guān)聯(lián)節(jié)點(diǎn)集合。顯然，S1∩S2=?,即這兩個(gè)單元沒有共享的關(guān)聯(lián)節(jié)點(diǎn)，因而可以描述裂紋處的位移間斷。同時(shí)，這兩個(gè)單元仍然是連續(xù)體，因而其計(jì)算同未斷裂單元完全一致，只是積分域需取為圖中的陰影區(qū)域。

圖2 虛擬節(jié)點(diǎn)法示意圖

Fig.2 Schematic diagram of phantom node method

采用式(11)的位移插值近似，并經(jīng)過標(biāo)準(zhǔn)的伽遼金過程，可得到控制方程及邊界條件式(1～5)的最終離散化方程為

(K-KΓu+βKp)U=f-fΓu+βfp

(14)

式中K為剛度陣，f為等效節(jié)點(diǎn)載荷向量，其他矩陣向量由Nitsche法[15](用以施加位移固定邊界條件)引入，β為相應(yīng)的罰系數(shù)。對(duì)于圖2所示的斷裂單元，各計(jì)算矩陣和向量為

(15)

(16)

(17)

(18)

式中i=1,2分別對(duì)應(yīng)于圖2的單元e1和e2，Ωe i為斷裂單元的有效域(圖中陰影區(qū))，Γi為有效域的邊界。

將擴(kuò)展有限元法中的虛擬節(jié)點(diǎn)技術(shù)引入到無網(wǎng)格法，無需增加節(jié)點(diǎn)額外自由度即可描述裂紋處的位移間斷，很大程度上簡化了無網(wǎng)格法處理裂紋問題，易于進(jìn)行程序?qū)崿F(xiàn)。

3 區(qū)域積分：一致性無網(wǎng)格法

如前所述，為提高應(yīng)力解析精度，本文充分利用無網(wǎng)格法易于形成高階近似的優(yōu)點(diǎn)，對(duì)位移采用二階無網(wǎng)格近似。然而，由于二階MLS形函數(shù)為非多項(xiàng)式的有理函數(shù)，式(15～17)的區(qū)域積分需使用較多積分點(diǎn)，才能保證積分穩(wěn)定性和精度，嚴(yán)重降低了計(jì)算效率[15]。為此，本文采用段慶林等[15,16]提出的一致性無單元伽遼金法(CEFG)，該方法與標(biāo)準(zhǔn)無網(wǎng)格法的不同之處在于區(qū)域積分格式和積分點(diǎn)上形函數(shù)導(dǎo)數(shù)的計(jì)算方法。對(duì)于本文采用的二階無網(wǎng)格近似，CEFG方法采用如圖3所示的積分格式，即在每個(gè)背景三角形單元中使用三個(gè)區(qū)域積分點(diǎn)(圖3的叉點(diǎn))，在每個(gè)三角形邊界上使用兩個(gè)邊界高斯點(diǎn)(圖3的三角形)。

為得到式(14)的剛度陣K,需計(jì)算區(qū)域積分點(diǎn)上形函數(shù)的導(dǎo)數(shù)。區(qū)別于標(biāo)準(zhǔn)的無網(wǎng)格方法，該導(dǎo)數(shù)在CEFG方法中由如下散度一致性的條件確定[16]。

(19)

式中ΩS為背景積分子域，ΓS為其邊界，NI(x)為節(jié)點(diǎn)形函數(shù)，q(x)=p,x(x)∪p,y(x),p(x)為形成MLS形函數(shù)的基底函數(shù)向量。對(duì)于本文采用的二階基底p(x)=[1xyx2xyy2]T,有q(x)=p,x(x)∪p,y(x)=[1xy]T,因而式(19)實(shí)際上含有3個(gè)獨(dú)立方程。對(duì)于每個(gè)積分子域，采用圖3所示的積分格式計(jì)算式(19)的區(qū)域和邊界積分，可得到相應(yīng)的3個(gè)離散方程，剛好可以求解該子域上3個(gè)區(qū)域積分點(diǎn)處的形函數(shù)導(dǎo)數(shù)。為區(qū)別于MLS形函數(shù)的經(jīng)典導(dǎo)數(shù)[6](通過對(duì)形函數(shù)直接求導(dǎo)得到)，以上由式(19)確定的形函數(shù)導(dǎo)數(shù)稱為修正導(dǎo)數(shù)[15]。文獻(xiàn)[15]的工作表明，對(duì)于二階無單元伽遼金法，采用圖3所示的積分格式，并使用修正導(dǎo)數(shù)計(jì)算剛度陣，可大幅提高計(jì)算效率并顯著改善計(jì)算的精度和收斂性，還能得到精確、光滑的應(yīng)力場，這是本文采用該方法的主要原因。

對(duì)于本文研究的裂紋問題，還需進(jìn)一步考慮背景積分單元發(fā)生斷裂后,如何構(gòu)造斷裂單元的積分格式。如圖4所示，當(dāng)積分單元發(fā)生斷裂時(shí)，首先，要根據(jù)該單元的關(guān)聯(lián)節(jié)點(diǎn)生成相應(yīng)的虛擬節(jié)點(diǎn)；然后，將該單元分為兩個(gè)單元，并根據(jù)裂紋位置確定各單元的關(guān)聯(lián)節(jié)點(diǎn)(包含虛擬節(jié)點(diǎn));最后，將這兩個(gè)單元的實(shí)際積分域劃分成三角形單元，并按圖3的積分格式布置積分點(diǎn)。顯然，一方面，由于引入了虛擬節(jié)點(diǎn)，裂紋處的位移間斷可以得到正確描述；另一方面，將實(shí)際積分域劃分成三角形單元，并按圖3的積分格式布置積分點(diǎn)，保證了計(jì)算的高效性和準(zhǔn)確性。

圖3 二階一致性無網(wǎng)格法積分格式示意圖

Fig.3 Schematic diagram of the integration scheme for quadraticconsistent element-free Galerkin method

4 裂紋擴(kuò)展模擬：斷裂準(zhǔn)則

為模擬裂紋擴(kuò)展，需確定裂紋的擴(kuò)展方向。本文采用最大周向應(yīng)力準(zhǔn)則，即裂紋沿著最大周向應(yīng)力的方向擴(kuò)展，因而裂紋擴(kuò)展角度θc由式(20)確定[3]，

KIsinθc+KII(3cosθc-1)=0

(20)

式中KI與KII為裂紋的應(yīng)力強(qiáng)度因子。式(20)的解可寫為

(21)

(22)

I(1,2)=2(K(1)IK(2)I+K(1)IIK(2)II)/Eeff

(23)

式中Eeff定義為

圖4 斷裂單元積分格式示意圖

Fig.4 Schematic diagram of the integration schemefor cracked element

(24)

E為彈性模量。因此，若將輔助場選為Mode I斷裂模式的裂尖場，即當(dāng)K(2)I=1,K(2)II=0時(shí)，有

(25)

同理，若將輔助場選為Mode II斷裂模式的裂尖場，即當(dāng)K(2)I=0,K(2)II=1時(shí)，有

(26)

5 數(shù)值算例

采用4個(gè)算例考察本文針對(duì)斷裂問題發(fā)展的二階一致性無單元伽遼金法Q-CEFG處理裂紋問題的有效性。如無特殊說明,所有算例的物理量均采用國際單位制。

圖6顯示了本文方法得到的σx x應(yīng)力場。可以看出，本文處理裂紋問題的Q-CEFG方法能夠正確反映該問題上下互不影響的常應(yīng)力狀態(tài)[25]，因而精確通過了該分片試驗(yàn),充分說明了本文方法處理裂紋間斷的有效性。

圖5 間斷分片試驗(yàn)的幾何構(gòu)型和加載示意圖

Fig.5 Schematic diagram of the configuration and loading of the discontinuous patch test

圖6 本文方法得到的間斷分片試驗(yàn)的σx x應(yīng)力場

Fig.6σx xfield of the discontinuous patch test obtained by the method of this paper

算例的幾何構(gòu)型如圖7所示，頂部受到均勻剪力作用，底部固定，按平面應(yīng)變問題求解，并取彈性模量E=107,泊松比υ=0.3?？紤]到該問題的應(yīng)力強(qiáng)度因子有準(zhǔn)確解KI=34.0,KII=4.55,可以用于比較計(jì)算精度，因而除本文建議的Q-CEFG方法外，二階標(biāo)準(zhǔn)的EFG方法(Q-EFG)和線性一致性無網(wǎng)格法(L-CEFG)也用于該算例的求解。

表1和表2分別顯示了三種方法KI和KII的計(jì)算結(jié)果。可以看出，與采用低階近似的L-CEFG方法相比，本文Q-CEFG方法得到的應(yīng)力強(qiáng)度因子更為準(zhǔn)確，計(jì)算精度高1～2個(gè)數(shù)量級(jí)，說明位移場的高階近似能有效改善應(yīng)力場的解析精度，從而能顯著提高應(yīng)力強(qiáng)度因子的計(jì)算精度。應(yīng)強(qiáng)調(diào)的是，本文方法沒有引入裂尖強(qiáng)化函數(shù)，但應(yīng)力強(qiáng)度因子的計(jì)算已相當(dāng)準(zhǔn)確。

與同樣采用二階近似的Q-EFG方法相比，一方面，本文Q-CEFG方法的計(jì)算精度更高，尤其是KI的精度約高1個(gè)數(shù)量級(jí)；另一方面，Q-EFG方法由于使用了更多的積分點(diǎn)，因而計(jì)算速度明顯低于本文的Q-CEFG方法。如采用2232個(gè)計(jì)算節(jié)點(diǎn)時(shí)，Q-EFG方法求解位移場消耗的CPU時(shí)間大約是本文Q-CEFG方法的2倍。這說明本文的Q-CEFG方法不僅改善了計(jì)算精度，而且加快了計(jì)算速度，因而顯著提高了計(jì)算效率。

圖7 單邊裂紋受到剪切作用

Fig.7 Edge -cracked plate under shear

表1 應(yīng)力強(qiáng)度因子KI的數(shù)值結(jié)果

Tab.1 Numerical results of stress intensity factorKI

Number of nodesQ-EFG L-CEFGQ-CEFGKIError/%KIError/%KIError/%61238.241245.463334.591.7119636.94841.232134.300.8223236.11539.121534.190.5

表2 應(yīng)力強(qiáng)度因子KII的數(shù)值結(jié)果

Tab.2 Numerical results of stress intensity factorK II

Number of nodesQ-EFG L-CEFGQ-CEFGKIIError/%KIIError/%KIIError/%6124.8056.70474.85611964.6215.81274.66222324.51-0.85.46204.540.2

算例用于考察本文Q-CEFG方法預(yù)測裂紋擴(kuò)展路徑的能力。如圖8所示，一結(jié)構(gòu)部件(圖中長度單位為mm)固定于一工字梁上，結(jié)構(gòu)部件一側(cè)的倒角處有一初始裂紋，載荷豎直向上，具體的實(shí)驗(yàn)設(shè)置可參考文獻(xiàn)[27]。取圖8的陰影區(qū)為計(jì)算域，并假定為平面應(yīng)力問題，材料參數(shù)取彈性模量E=2×1011,泊松比υ=0.3,載荷P=1.0,初始裂紋長度a0=5 mm。根據(jù)文獻(xiàn)[28]，計(jì)算域的底部采用固定和鉸支約束分別相應(yīng)于厚梁和薄梁的情況。

計(jì)算采用2294個(gè)節(jié)點(diǎn)，本文方法得到的對(duì)應(yīng)于兩種約束情況的裂紋擴(kuò)展路徑如圖9所示?？梢钥闯觯玫降穆窂脚cVentura等[28]的結(jié)果十分吻合，驗(yàn)證了本文方法模擬裂紋擴(kuò)展問題的有效性。

圖8 倒角裂紋問題的實(shí)驗(yàn)構(gòu)型和求解域示意圖

Fig.8 Experimental configuration and simulated region ofthe crack growth from a fillet problem

圖9 倒角裂紋問題的裂紋路徑

Fig.9 Crack trajectory of the crack growth froma fillet problem

圖10 帶圓孔的三點(diǎn)彎曲梁示意圖

Fig.10 Schematic diagram of the beam with 3 holes subjectedto 3-point bending

為進(jìn)一步考察本文方法，考慮如圖10所示的三孔梁的三點(diǎn)彎曲問題。按平面應(yīng)力問題求解，并取彈性模量E=3×1010,泊松比υ=0.3。對(duì)于梁底部初始裂紋的長度和位置，本文考慮兩種情況，情況I為a=1.0,b=4.0；情況II為a=1.5,b=5.0。

采用2602個(gè)節(jié)點(diǎn)進(jìn)行計(jì)算，得到的兩種情況下的裂紋路徑如圖11所示。圖12為本文方法得到的裂紋附近局部放大圖及Bittencourt等[29]給出的裂紋路徑結(jié)果?？梢钥闯?，裂紋路徑與文獻(xiàn)[29]的實(shí)驗(yàn)結(jié)果吻合良好，再次表明了本文方法模擬裂紋擴(kuò)展問題的正確性。

圖11 三孔梁問題的裂紋路徑

Fig.11 Crack trajectory of the beam with 3 holes problem

圖12 裂紋路徑與實(shí)驗(yàn)[28]對(duì)比

Fig.12 Comparison of crack trajectories with the experimental results[28]

6 結(jié) 論

本文將針對(duì)連續(xù)體的一致性無網(wǎng)格法拓展到以裂紋為特征的非連續(xù)體,實(shí)現(xiàn)了裂紋擴(kuò)展過程的數(shù)值模擬和裂紋路徑的準(zhǔn)確預(yù)測。一致性無網(wǎng)格法的積分格式和導(dǎo)數(shù)修正技術(shù)在很大程度上改善了計(jì)算精度和效率。其中一個(gè)重要發(fā)現(xiàn)是高階無網(wǎng)格法在不引入裂尖強(qiáng)化函數(shù)的前提下也能實(shí)現(xiàn)應(yīng)力強(qiáng)度因子的準(zhǔn)確計(jì)算，而低階無網(wǎng)格法則導(dǎo)致了相當(dāng)大的誤差。應(yīng)強(qiáng)調(diào)的是，引入裂尖強(qiáng)化函數(shù)對(duì)于裂紋擴(kuò)展模擬極其不便，因?yàn)橛?jì)算節(jié)點(diǎn)的強(qiáng)化函數(shù)需隨著裂尖的移動(dòng)不斷地刪除和增加。因此，本文的一致性高階無網(wǎng)格法既能保證應(yīng)力強(qiáng)度因子的計(jì)算精度，又能十分方便地實(shí)現(xiàn)裂紋擴(kuò)展過程的數(shù)值模擬，具有很好的發(fā)展?jié)摿Α?/p>

:

[1] Belytschko T，Black T.Elastic crack growth in finite elements with minimal remeshing[J].InternationalJournalforNumericalMethodsinEngineering，1999,45(5):601-620.

[2] Mo?s N，Dolbow J，Belytschko T.A finite element method for crack growth without remeshing[J].InternationalJournalforNumericalMethodsinEngineering，1999,46(1)：131-150.

[3] 莊 茁，柳占立，成斌斌,等.擴(kuò)展有限單元法[M].北京：清華大學(xué)出版社，2012.(ZHUANG Zhuo，LIU Zhan-li，CHENG Bin-bin,et al.TheExtendedFiniteElementMethod[M].Beijing:Tsinghua University Press，2012.(in Chinese))

[4] Tian R，Wen L F.Improved XFEM—An extra-DOF free，well-conditioning，and interpolating XFEM[J].ComputerMethodsinAppliedMechanicsandEn-gineering，2015，285:639-658.

[5] 王 振，余天堂.模擬三維裂紋問題的自適應(yīng)多尺度擴(kuò)展有限元法[J].工程力學(xué)，2016，33(1)：32-38.(WANG Zhen,YU Tian-tang.Adaptive multiscale extended finite element method for modeling three -dimensional crack problems[J].EngineeringMecha-nics,2016,33(1):32-38.(in Chinese))

[6] Belytschko T，Lu Y Y，Gu L.Element-free Galerkin methods[J].InternationalJournalforNumericalMethodsinEngineering，1994,37(2):229-256.

[7] 張 雄，劉 巖.無網(wǎng)格法[M].北京：清華大學(xué)出版社，2004.(ZHANG Xiong，LIU Yan.MeshlessMe-thod[M].Beijing：Tisinghua Universtiy Press，2004.(in Chinese))

[8] Lu Y Y，Belytschko T，Tabbara M.Element-free Galerkin method for wave propagation and dynamic fracture[J].ComputerMethodsinAppliedMecha-nicsandEngineering，1995,126(1):131-153.

[9] Belytschko T，Lu Y Y，Gu L，et al.Element-free Galerkin methods for static and dynamic fracture[J].InternationalJournalofSolidsandStructures，1995，32(17/18):2547-2570.

[10] Krysl P，Belytschko T.The Element Free Galerkin method for dynamic propagation of arbitrary 3-D cracks[J].InternationalJournalforNumericalMethodsinEngineering,1999,44(6):767-800.

[11] Fleming M，Chu Y A，Moran B，et al.Enriched element-free Galerkin methods for crack tip fields[J].InternationalJournalforNumericalMethodsinEn-gineering，1997，40(8):1483-1504.

[12] Duflot M.A meshless method with enriched weight functions for three -dimensional crack propagation[J].InternationalJournalforNumericalMethodsinEngineering，2006,65(12):1970-2006.

[13] Shedbale A S，Singh I V，Mishra B K，et al.Ductile failure modeling and simulations using coupled FE-EFG approach[J].InternationalJournalofFracture,2017，203(1-2):183-209.

[14] 張 雄，宋康祖，陸明萬.無網(wǎng)格法研究進(jìn)展及其應(yīng)用[J].計(jì)算力學(xué)學(xué)報(bào)，2003,20(6)：730-742.(ZHANG Xiong，SONG Kang-zu，LU Ming-wan.Research progress and application of meshless method[J].ChineseJournalofComputationalMechanics，2003,20(6):730-742.(in Chinese)).

[15] Duan Q L，Li X K，Zhang H W，et al.Second-order accurate derivatives and integration schemes for meshfree methods[J].InternationalJournalforNumericalMethodsinEngineering，2012，92(4):399-424.

[16] Duan Q L,Gao X，Wang B B，et al.Consistent element-free Galerkin method[J].InternationalJournalforNumericalMethodsinEngineering，2014，99(2):79-101.

[17] Duan Q L，Gao X，Wang B B，et al.A four-point integration scheme with quadratic exactness for three -dimensional element-free Galerkin method based on variationally consistent formulation[J].ComputerMethodsinAppliedMechanicsandEngineering，2014，280(10):84-116.

[18] 邵玉龍，段慶林，高 欣,等.自適應(yīng)一致性無單元伽遼金法1)[J].力學(xué)學(xué)報(bào)，2017，49(1)：105-116.(SHAO Yu-long，DUAN Qing-lin，GAO Xin,et al.Adaptive consistent high order element-free Galerkin method1)[J].ChineseJournalofTheoreticalandAppliedMechanics，2017，49(1)：105-116.(in Chinese))

[19] Ortiz-Bernardin A，Hale J S,Cyron C J.Volume -ave -raged nodal projection method for nearly-incompressible elasticity using meshfree and bubble basis functions[J].ComputerMethodsinAppliedMechanicsandEngineering，2015，285:427-451.

[20] Ortiz-Bernardin A，Puso M A，Sukumar N.Improved robustness for nearly-incompressible large deformation meshfree simulations on Delaunay tessellations[J].ComputerMethodsinAppliedMechanicsandEngineering，2015，293:348-374.

[21] Song J H，Areias P M A,Belytschko T.A method for dynamic crack and shear band propagation with phantom nodes[J].InternationalJournalforNumericalMethodsinEngineering，2006,67(6):868-893.

[22] Ai W，Augarde C E.An adaptive cracking particle method for 2D crack propagation[J].InternationalJournalforNumericalMethodsinEngineering，2016,108(13):1626-1648.

[23] Kumar S，Singh I V，Mishra B K，et al.Modeling and simulation of kinked cracks by virtual node XFEM[J].ComputerMethodsinAppliedMechanicsandEngineering,2015,283:1425-1466.

[24] Lua J，Zhang T T，F(xiàn)ang E，et al.Explicit phantom paired shell element approach for crack branching and impact damage prediction of aluminum structures[J].InternationalJournalofImpactEngineering，2016,87:28-43.

[25] Dolbow J E，Devan A.Enrichment of enhanced assumed strain approximations for representing strong discontinuities:addressing volumetric incompressibi-lity and the discontinuous patch test[J].InternationalJournalforNumericalMethodsinEngineering，2004,59(1):47-67.

[26] 劉 豐，鄭 宏，李春光.基于NMM的EFG方法及其裂紋擴(kuò)展模擬1)[J].力學(xué)學(xué)報(bào)，2014，46(4):582-590.(LIU Feng，ZHENG Hong,LI Chun-guang.The NMM-based EFG method and simulation of crack propagation1)[J].ChineseJournalofTheoreticalandAppliedMechanics,2014，46(4):582-590.(in Chinese))

[27] Sumi Y，Yang C，Wang Z N.Morphological aspects of fatigue crack propagation Part II—effects of stress biaxiality and welding residual stress[J].Internatio-nalJournalofFracture,1996，82(3):221-235.

[28] Ventura G，Xu J X，Belytschko T.A vector level set method and new discontinuity approximations for crack growth by EFG[J].InternationalJournalforNumericalMethodsinEngineering，2002,54(6):923-944.

[29] Bittencourt T N，Wawrzynek P A，Ingraffea A R，et al.Quasi-automatic simulation of crack propagation for 2D LEFM problems[J].EngineeringFractureMechanics,1996,55(2):321-334.