擴展無網(wǎng)格法分析功能梯度材料斷裂問題

張保文， 馬文濤， 黃凌霄

（寧夏大學數(shù)學計算機學院，寧夏銀川 750021）

擴展無網(wǎng)格法分析功能梯度材料斷裂問題

張保文， 馬文濤， 黃凌霄

（寧夏大學數(shù)學計算機學院，寧夏銀川 750021）

提出一種擴展無網(wǎng)格法模擬二維功能梯度材料I、II型復合斷裂問題.基于單位分解思想，在EFGM近似函數(shù)中分別添加階躍函數(shù)和奇異函數(shù)，依次表征不連續(xù)位移場和裂尖奇異應力場.結(jié)合一種不需要求解梯度材料參數(shù)的互作用積分的域形式計算了復合型應力強度因子.兩個功能梯度板的數(shù)值算例驗證了單位分解擴展無網(wǎng)格法的可行性和有效性.

功能梯度材料；EFGM；擴展無網(wǎng)格法；應力強度因子

非均勻性是導致材料力學性質(zhì)復雜的主要原因之一.功能梯度材料（FGMs）是一類十分特殊的非均勻材料，其組成、結(jié)構(gòu)和物理力學性質(zhì)在空間上主要沿某一方向變化而在垂直于這個方向的平面或曲面上的變化往往非常小.天然材料中的生物材料和巖土材料往往具有FGMs的性質(zhì).近年來，F(xiàn)GMs在機械、生物、化學、電子等工程領域中也得到了廣泛的應用，其力學性質(zhì)成為研究的熱點.由于材料結(jié)構(gòu)和組成不同，不同梯度材料的細觀、宏觀斷裂機理也不同.在細觀、宏觀尺度上，裂紋和孔洞對梯度材料結(jié)構(gòu)的力學性質(zhì)有著重要的影響.對這些缺陷進行力學分析成為梯度材料研究的一個分支，也成為斷裂力學研究的熱點［1］.

FGMs的早期研究可追溯到1960年，Gibsin將土視為非均勻材料［2］.由于問題的難度和復雜程度都很高，含裂紋的FGMs往往被假定為彈性模量是關于空間坐標變化的線性函數(shù)或指數(shù)函數(shù).Erdogan等［3］和Eischen等［4］分別采用理論分析和有限單元法證明了FGMs和各向同性材料的裂尖應力場都具有的奇異性.Kim和Paulino［5］采用三種不同的方法：路徑無關J積分方法、修正的裂紋閉合積分法和位移修正方法計算了FGMs的應力強度因子.Dolbow和Gosz［6］在FGM斷裂分析中引入相互作用積分計算了二維復合應力強度因子.他們的研究表明，與傳統(tǒng)的J積分相比，相互作用積分顯得更為方便，因為后者無需計算沿自由裂紋面的應變能密度.Kim和Paulio［7］采用FEM和相互作用積分計算了各向同性FGM上的復合應力強度因子.Rao和Raham［8］采用相互作用積分和無網(wǎng)格法求解了FGMs的應力強度因子.Kim和Paulio［9］還總結(jié)了三種不同輔助場的定義，討論了如何獲取復合型應力強度因子和T應力的方法.Yu等［10］提出了一種無需計算材料梯度的相互作用積分，并結(jié)合擴展有限元計算了非各向同性材料的I、II復合應力強度因子.Mojdehi等［11］采用無網(wǎng)格局部－Petrov－Galerkin（MLPG）分析了三維靜態(tài)和動態(tài)功能梯度厚板問題，但未涉及板內(nèi)包含裂紋的情況.Zhang和Kim［12］推導了二維各向異性FGM的裂尖奇異應力場，指出高階項受材料梯度的影響非常大，在奇異裂尖場中扮演著重要的角色.

近年來，無網(wǎng)格法被廣泛應用于求解斷裂力學問題.由于無網(wǎng)格法只需節(jié)點信息，而不需要單元之間的連接信息，因此很容易避免傳統(tǒng)FEM方法中面臨的網(wǎng)格劃分和重構(gòu)帶來的巨大困難.目前，無網(wǎng)格法的研究主要集中于均勻材料的斷裂分析，極少有人關注和將其應用于FGMs中［13］.基于此，本文主要的工作是利用單位分解擴展無網(wǎng)格法求解FGMs的復合應力強度因子.該方法與傳統(tǒng)的內(nèi)部基擴展無網(wǎng)格法不同，無需使用修改權函數(shù)的方法，裂紋的幾何位置由符號距離函數(shù)表達，裂紋所引起的不連續(xù)性和裂尖奇異性則由近似函數(shù)中的不連續(xù)項和裂尖奇異項表達.由不需要計算材料參數(shù)導數(shù)的相互作用積分的域形式，計算了不

同梯度材料的復合應力強度因子.

1　單位分解擴展無網(wǎng)格法

1.1移動最小二乘近似（MLS）

考慮區(qū)域中的函數(shù)，根據(jù)移動最小二乘（MLS）法，其近似函數(shù)可表示為

其中n節(jié)點x 影響域內(nèi)包含的節(jié)點數(shù)目；ΦT（x）＝PT（x）A－1（x）C（x）；PT（x）＝｛1，x1，x2｝，A（x）＝PTWPT，C（x）＝PT（x）W；W（x）＝diag［W1（x），W2（x），…，Wn（x）］為與節(jié)點相關的權函數(shù)，本文選取為四次樣條權函數(shù)；dT＝｛d1，d2，…，dn｝，dI為與節(jié)點I相關的參數(shù).

1.2不連續(xù)近似函數(shù)

傳統(tǒng)的MLS只能構(gòu)造連續(xù)的近似函數(shù)，針對含裂紋的不連續(xù)問題，必須對其進行改進.本文采用單位分解法對MLS進行修改，得到不連續(xù)位移場函數(shù).其基本思想是在MLS近似函數(shù)中添加階躍函數(shù)項和裂尖奇異項，分別表征由裂紋導致的不連續(xù)位移場和裂尖奇異應力場.該方法的思想類似XFEM，其近似函數(shù)在裂紋周圍一狹小區(qū)域內(nèi)（大小由節(jié)點影響域半徑確定）是不連續(xù)的，在遠離裂紋的區(qū)域內(nèi)則是光滑連續(xù)的.單位分解擴展無網(wǎng)格法已被成功應用于模擬黏聚裂紋［13］、裂紋開裂［14］和巖體材料接觸摩擦分析［15］等問題，本文將其應用于求解功能梯度材料斷裂問題.

單位分解擴展無網(wǎng)格法近似函數(shù)可表示為

圖1　含裂紋的功能梯度材料Fig.1 A crack in a functionally graded material

其中：N為整個計算區(qū)域內(nèi)的節(jié)點集；Nc為節(jié)點影響域被裂紋完全切割的節(jié)點集；Nt為節(jié)點影響域部分被裂紋切割的節(jié)點集（見圖1）.公式（6）中ΦI（x）為MLS形函數(shù)，第一項為標準的MLS近似函數(shù)；第二項為階躍擴展函數(shù)，第三項為裂尖擴展函數(shù).αI，βkI為引入的附加變量；通常，階躍函數(shù)H（fI（x））取為符號距離函數(shù)，定義為

其中xup為裂紋尖端坐標，n為裂紋面法向.

其中r，θ為以裂紋尖端為坐標原點的極坐標系（見圖1）.

1.3彈性力學方程的弱形式及其離散化

二維彈性力學問題的虛功方程為

式中Ω為計算域，Γ＝Γu∪Γt為域邊界，δε為虛應變向量，u＝［u v］T為位移向量，σ為與u對應的應力向量，b為體力向量，珋t為應力邊界Γt給定的面力，珔u為位移邊界Γu上給定的位移，λ為拉格朗日乘子向量.

將式（2）的不連續(xù)近似函數(shù)代入式（6）并考慮位移變分的任意性，可得節(jié)點的離散線性方程

其中Bstd由標準的MLS近似得到；Benr由擴展部分形成；當xI∈ΩΓ，則ΨI（x）是階躍函數(shù)H（fI（x））；當xI∈ΩΛ，則ΨI（x）是擴展基函數(shù)Tα（x）.方程中的形函數(shù)N為1階拉格朗日插值函數(shù).D（x）為彈性矩陣，對于平面應力問題，

對于平面應變問題，

其中E（x）和v（x）分別為彈性模量和泊松比.對于各向同性材料，E（x）和v（x）均為常數(shù)；對于梯度材料，E（x）和v（x）則是隨空間坐標變化的函數(shù).

2　互作用積分

目前，在求解復合模式荷載作用下的應力強度因子時，互作用積分的域形式被廣泛采用，也被證明是十分準確的方法.

2.1輔助場

在計算互作用積分時，輔助場起著至關重要的作用.通常有幾種不同的形式，本文采用不完備形式.在裂尖局部極坐標下，其不完備形式定義為

其中KauxI和KauxII分別為I型和II型輔助應力強度因子，為一致張量.對應的函數(shù)和為

這里，κup＝3－4υup（平面應變）或κup＝（3－υup）／（1＋υup（平面應力），μup，υup分別為裂尖附近的剪切模量和泊松比.

2.2互作用積分的域形式

I型與II型復合應力強度因子可以通過互作用積分的區(qū)域形式（等效區(qū)域積分）計算.等效區(qū)域積分的具體形式為：

方程（19）中的第一項與均勻材料的公式是相同的.其中Supijkl為裂尖處的一致張量；A為由邊界Γ＝Γ0＋Γ＋－Γs＋Γ－圍成的積分區(qū)域（見圖2）；q為足夠光滑的權函數(shù)，要求在內(nèi)邊界Γs上為1，在外邊界Γ0上為0.

在Dolbow等［6］、Kim等［7］、Rao等［8］的論文中，等效區(qū)域積分的形式中涉及到了材料參數(shù)的導數(shù)（Dijkl，1或Sijkl，1），公式（19）并不涉及材料參數(shù)的導數(shù).在實際問題中，材料參數(shù)的導數(shù)很難確定，甚至根本不存在.因此公式比傳統(tǒng)的等效區(qū)域積分更適合梯度材料的斷裂分析，相應地，應力強度因子與互作用積分之間的關系為

圖2　互作用積分的積分區(qū)域Fig.2 Integral domain of interaction integral

其中E′tip＝Etip（平面應力）或E′tip＝Etip／（1－υtip）（平面應變）.Etip，υtip分別為裂紋尖端處的楊氏模量和泊松比.當KauxI＝1，KauxII＝0時，KI＝E′tipI／2，而當KauxI＝0，KauxII時，KII＝EtipI／2.

2.3互作用積分的無網(wǎng)格離散形式

為了在擴展無網(wǎng)格法中求解互作用積分，需要將方程離散為

其中eN為積分區(qū)域A內(nèi)的積分網(wǎng)格數(shù)目，gN為每個網(wǎng)格中的高斯積分點數(shù).對于梯度材料，在建立剛度矩陣時采用積分點上的材料屬性.

3　數(shù)值算例

采用MATLAB編寫了擴展無網(wǎng)格法的相應程序.近似函數(shù)建立過程中采用四次樣條權函數(shù).圓形影響域半徑取為1.7倍節(jié)點間距.公式（7）和（22）均采用背景積分網(wǎng)格計算.沒有被裂紋切割的背景網(wǎng)格，采用16個高斯積分點；被裂紋切割的背景網(wǎng)格，將其分解為幾個三角形網(wǎng)格，每個網(wǎng)格中采用36個高斯積分點.

3.1單邊傾斜裂紋

考慮含有一條單邊傾斜裂紋的二維平板問題，如圖3所示.平板的長L＝2，寬W＝1.邊裂紋傾角γ＝

E（x1）＝E珚exp［η（x1－0.5）］，0≤x1≤W

（22）

其中珚E和η是兩個材料參數(shù).在計算過程中，布置21 ×41個節(jié)點，20×40背景積分網(wǎng)格；珚E＝1和η依次取0，0.1，0.25，0.5和1.泊松比υ＝0.3.平板上部邊界施加荷載σ22＝珋ε珚Eexp［η（x1－0.5）］（珋ε＝1），底部施加位移邊界條件.表1比較了本文方法與Kim等［5］采用有限元法計算不同值對應的正則化應力強度因子及其相對誤差.可以看出，本文方法與FEM的計算結(jié)果吻合得非常好.

3.2中心傾斜裂紋

表1　傾斜裂紋的正則化應力強度因子Table 1 Normalized SIFs for a slanted crack in a plate

如圖4所示，平板尺寸為2L＝2W＝20，中心裂紋長2a＝2，裂紋傾角為γ.計算過程中，泊松比υ＝0.3，彈性模量取為指數(shù)函數(shù)，具體形式為

圖3　復合荷載作用下的單邊傾斜裂紋板Fig.3 Slanted crack in plate under mixed－mode loading

圖4　復合荷載作用下的中心傾斜裂紋板Fig.4 Central inclined rack in plate under mixed－mode loading

其中材料參數(shù)珚E＝1，η依次取0.25和0.5.在平板頂部施加荷載σ22＝珋ε珚Eexp（ηx1）（珋ε＝1）.布置21×41個節(jié)點，20×40背景積分網(wǎng)格；γ／π依次取為0，0.1，0.2，0.3，0.4和0.5.Konda等［16］將該問題視為無限大平板問題得到相應的解析解.無網(wǎng)格法不能直接求解無限大區(qū)域問題，但當取a／W＝a／L≤0.01時，近似結(jié)果是可以接受的.表2和表3分別比較了當η＝0.25和η＝0.5時兩個裂尖在不同γ／π值下對應的正則化應力強度因子的解析解和數(shù)值解.顯然，本文方法與解析解十分接近.

表2　中心裂紋正則化應力強度因子（η＝0.25）Table 2 Normalized SIFs for a plate with an interior inclined crack（η＝0.25）

4　結(jié)論

本文提出單位分解擴展無網(wǎng)格法求解了二維功能梯度材料斷裂力學問題.為了準確描述功能梯度材料裂紋的不連續(xù)位移場和裂尖奇異應力場，根據(jù)單位分解思想，在EFGM近似函數(shù)中添加了階躍函數(shù)項和裂尖奇異函數(shù)項.采用一種修正的互作用積分域形式計算I型、II型復合應力強度因子.復合荷載作用下的邊界斜裂紋算例和中心斜裂紋算例，研究了本文方法的有效性和計算精度.計算結(jié)果表明，本文方法非常適合求解功能梯度材料斷裂問題，具有廣闊的發(fā)展空間.

表3　中心裂紋正則化應力強度因子（η＝0.5）Table 3 Normalized SIFs for a plate with an interior inclined crack（η＝0.5）

［1］ XIZO H T，YUE Z.New Boundary element analysis of fracture mechanics in functionally graded materials［M］.Beijing：Higher Education Press，2011：25.

［2］ GIBSON R E.Some results concerning displacements and stresses in a non－h(huán)omogenous elastic layer［J］.Geotechnique，1967，（17）：58－67.

［3］ DELALE F，ERDOGAN F.The crack problem for a nonhomogeneous plane［J］.International Journal of Applied Mechanics，1983，（50）：609－613.

［4］ EISCHEN J W.Fracture of nonhomogeneous materials［J］.International Journal of Fracture，1987，（34）：3－22.

［5］ KIM J H，PAULINO G H.Finite element evaluation of mixed mode stress intensity factors in functionally graded materials［J］.International Journal of Numerical Method in Engineering，2002，53（8）：1903－1938.

［6］ DOLBOW J E，GOSZ M.On the computation of mixed stress intensity factors in functionally graded materials［J］.International Journal of solids and structure，2002，39（9）：2557－2631.

［7］ KIM J H，PAULINO G H.An accurate scheme for mixed－mode fracture analysis of functionally graded materials using the interaction integral and micromechanics models［J］.International Journal of Numerical Method in Engineering，2003，（58）：1457－1479.

［8］ RAO B N，RAHMAN S.Mesh－free analysis of cracks in isotropic functionally graded materials［J］.Engineering Fracture Mechanics，2003，（5）：1－27.

［9］ KIM J H，PAULINO G H.T－stress，mixe－mode stress intensity factors，and crack initiation angles in functionally graded materials：a unified approach using the interaction integral method［J］.Computer methods in applied mechanics and engineering，2003，（192）：1463－1494.

［10］ YU H J，WU L Z，LI C G.Investigation of mixed－mode stress intensity factors for nonhomogeneous materials using an interaction integral method［J］.International journal of solids and structures，2009，（46）：3710－3724.

［11］ MOJDEHI A R，DARVIZEH A，BASTI A.Three dimensional static and dynamic analysis of thick functionally graded plates by the meshless local Petrov－Galerkin（MLPG）method［J］.Engineering analysis with boundary elements，2011，（35）：1168－1180.

［12］ ZHANG L H，KIM J H.Mixed－mode cracked－tip fields in an anisotropic functionally graded material［J］.Journal of applied mechanics，2012，（79）：1－10.

［13］ RABCZUK T，ZI G.A meshfree based on the local partition of unity for cohesive cracks［J］.Computational Mechanics，2006，39（6）：743－760.

［14］ 馬文濤，師俊平，李寧.水平集與無網(wǎng)格耦合法在裂紋擴展中應用［J］.巖土力學，2012，33（11）：3447－3453.

［15］ 馬文濤，師俊平，李寧.模擬摩擦接觸問題的新型無網(wǎng)格數(shù)值方法［J］.巖土力學，2012，33（10）：3145－3150.

［16］ KONDA N，ERDOGAN F.The mixed mode crack problem in a nonhomogeneous elastic［J］.Engineering Fracture Mechanics，1994，47 （4）：533－545.

Enriched Meshless Analysis of Cracks in Functionally Graded Materials

ZHANG Bao－wen， MA Wen－tao， HUANG Ling－xiao
（School of Mathematics and Computer Science，Ningxia University，Yinchuan 750021，China）

An enriched meshless method was presented for analysing the mixed mode I and mode II fracture problem in two－dimension functionally graded materials.Both enriched functions，including jump function and crack tip singularity function，were added in the approximation of EFGM based on partition of unity.The mixed mode stress intensity factors for cracks in functionally graded materials were numerically evaluated using the modified domain form of interaction integral.The integrand does not involve any derivatives of materials properties.The numerical results were obtained for edge and center cracks，and were found to be in good with the reference solutions for the functionally graded material crack problems.

functionally graded materials；EFGM；enriched meshless method；stress intensity factors

TB301

A

1001－2443（2016）03－0230－06

10.14182／J.cnki.1001－2443.2016.03.005

2015－08－13

國家自然科學基金（51269024、41204041）.

張保文（1975－），男，碩士，講師，復分析及其在力學中的應用.

引用格式：張保文，馬文濤，黃凌霄.擴展無網(wǎng)格法分析功能梯度材料斷裂問題［J］.安徽師范大學學報：自然科學版，2016，39（3）：230－236.