全空間中粒子和裂紋的對偶邊界積分方程及數(shù)值解

馬杭,潘蒙

(1.上海大學(xué)理學(xué)院,上海　200444；

2.上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上?！?00072)

全空間中粒子和裂紋的對偶邊界積分方程及數(shù)值解

馬杭1,潘蒙2

(1.上海大學(xué)理學(xué)院,上海200444；

2.上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上海200072)

針對彈性固體中同時含有粒子和裂紋的情況,建立了全空間中的粒子和裂紋的位移間斷形式的對偶邊界積分方程計算模型,解決了全空間條件下難以對研究對象進(jìn)行直接加載的問題.采用邊界積分方程的離散形式對含有少量粒子和裂紋的典型情況進(jìn)行了數(shù)值分析,其中對粒子邊界(或界面)和裂紋面分別采用邊界點法和高斯配點法進(jìn)行離散,計算了裂紋的應(yīng)力強度因子,探討了粒子與裂紋的相互作用.通過與已有研究結(jié)果比較,驗證了計算模型與計算機程序的正確性與可靠性.

粒子；裂紋；對偶邊界積分方程；數(shù)值解；應(yīng)力強度因子

彈性固體中不可避免地含有某些缺陷,例如裂紋和粒子.彈性固體中的粒子可以是夾雜物或者局部析出相,也可以是人工加入的增強相,例如粒子增強復(fù)合材料.含有缺陷的彈性固體屬于非均質(zhì)問題.含有裂紋和粒子等缺陷使得彈性固體的性能發(fā)生改變,這使其成為固體力學(xué)研究人員和材料工程師所關(guān)注的焦點之一.另外,對全空間中粒子和裂紋彈性狀態(tài)的研究往往是進(jìn)一步研究的基礎(chǔ)和參照基準(zhǔn),因此具有重要的意義.

在線彈性斷裂力學(xué)的研究中,有兩類相輔相成的研究方法,即解析方法和數(shù)值方法.多年來通過兩類方法都取得了大量的研究成果.解析方法適于裂紋構(gòu)型相對簡單的情況,數(shù)值方法則適于分析比較復(fù)雜的情況.在各種數(shù)值方法中,基于邊界積分方程的數(shù)值方法表現(xiàn)出更大的優(yōu)越性,這已逐漸成為研究者的共識[1-2].在含有裂紋的問題中,由于裂紋的兩個面在幾何上是完全重合的,如果僅僅使用一個邊界積分方程,必須采用特殊的手段,例如采用分域法將裂紋置于問題的邊界,否則將導(dǎo)致問題的退化而不可解[3].具體而言,退化問題經(jīng)過離散后得到的系統(tǒng)方程的系數(shù)矩陣是奇異的,因而不可解.采用對偶邊界積分方程,即通過位移邊界積分方程和面力邊界積分方程來共同描述含有裂紋的問題,就能夠避免問題的退化,因此成為邊界積分方程解決裂紋問題的主要方法.

關(guān)于裂紋問題已有大量研究成果,例如各種條件下裂紋之間的相互作用、裂紋與孔洞的相互作用等[4-7].Dong等[8-13]較早地研究了裂紋與粒子間相互作用,首先提出了積分方程方法,其中粒子用區(qū)域積分描述,需要在粒子上剖分內(nèi)部胞元；繼而提出了邊界積分方程方法,并對原有算法進(jìn)行改進(jìn),研究了各向異性、橫觀各向同性條件下裂紋與粒子的相互作用問題以及材料中剛性線的影響.采用邊界積分方程方法只需在界面上剖分單元,但總體而言,采用該方法對裂紋與粒子間的相互作用進(jìn)行研究的成果相對較少.本工作中的粒子是指廣義的粒子,其彈性模量可以大于或小于基體的彈性模量,如果將粒子的彈性模量設(shè)為零,則粒子退化為孔洞,在這個意義上,本工作的研究具有一般性.本工作的長遠(yuǎn)目標(biāo)是研究大規(guī)模問題的計算方法以作為進(jìn)一步研究的基礎(chǔ),因此需要在全空間中研究粒子和裂紋的彈性狀態(tài).由于全空間不存在外邊界,其加載條件有一定的特殊性,在全空間中無法對所研究的對象按有限物體那樣直接進(jìn)行加載.因此在前人[14]工作的基礎(chǔ)上,首先建立了全空間中粒子和裂紋的位移間斷形式的對偶邊界積分方程計算模型,解決加載問題；然后采用邊界積分方程的離散形式——邊界點法[15]對含有少量粒子和裂紋的典型情況進(jìn)行數(shù)值分析,計算裂紋的應(yīng)力強度因子,探討粒子與裂紋的相互作用；最后通過與已有研究結(jié)果進(jìn)行比較,對計算模型與計算機程序的正確性與可靠性進(jìn)行驗證.

1　計算模型與求解

1.1邊界積分方程

對于區(qū)域為?,邊界為Γ的各向同性彈性體,位移和面力的邊界積分方程[1]分別為

1.2區(qū)域的分解與加載

圖1　全空間中受遠(yuǎn)場單位拉伸載荷的粒子和裂紋Fig.1　Particle and crack under far field unit tension in full space

1.3矩陣形式的邊界積分方程

圖2中分離的粒子可用式(1)描述.離散后的式(1)可寫成以下矩陣形式：

圖2　分離的粒子與全空間中的孔洞和裂紋Fig.2　Separated particle and the hole and crack in full space

圖3　受遠(yuǎn)場拉伸的均質(zhì)全空間與不受遠(yuǎn)場載荷的孔洞和裂紋Fig.3　Full uniform space under far-field tension and the hole and crack without far-field load

對于全空間中的孔洞和裂紋(見圖3),可用位移間斷的對偶邊界積分方程(3)和(4)來描述.當(dāng)源點p位于孔洞邊界(或粒子邊界)時采用式(3)描述,而當(dāng)源點p位于裂紋表面時則采用式(4)描述.離散后的矩陣形式為

式中,R為右向量,是已知量,即

利用式(8)求得粒子界面的位移u和裂紋張開位移δ,再利用式(6)求出粒子邊界的面力τ.分別記基體與粒子的彈性模量為E和EI,定義模量比為假設(shè)基體與粒子的泊松比相同,則存在如下關(guān)系：

利用式(11)和(12)可使計算得到一定程度的簡化.

如果粒子的彈性模量為0,即β=0,則粒子和裂紋的問題退化為孔洞和裂紋的問題,孔洞邊界上的面力為0.這時,式(12)不成立,式(8)退化為

1.4離散方式與應(yīng)力強度因子的計算

采用邊界積分方程的兩種離散形式分別對粒子邊界和裂紋表面進(jìn)行離散,其中粒子邊界的離散采用了邊界點法[15].設(shè)N為粒子邊界的節(jié)點數(shù),對于只有一個粒子的情況,列向量u, τ的大小為d·N,系數(shù)矩陣TII,UII的大小為(d·N)2,其中d為維數(shù).如果粒子總數(shù)為NI,則列向量u,τ的大小為d·NI·N,矩陣TII,UII的大小為(d·NI·N)2.

對裂紋表面采用高斯配點法進(jìn)行離散[6,16].當(dāng)源點和場點均位于同一裂紋的上表面時,式(4)轉(zhuǎn)化為

式中,HFP表示Hadamard主值積分[17].利用積分核τ?ijk的具體表達(dá)式,可將式(14)簡化為

式中,μ為剪切模量,ν為泊松比,r為源點和場點間的距離.式(15)表明在二維條件下,張開型和滑開型位移間斷與相應(yīng)的面力(即裂紋的法向面力和切向面力)具有完全相同的對應(yīng)關(guān)系.記NG為裂紋表面的高斯配點數(shù),對于單位長度的裂紋,即設(shè)裂紋半長為a=1,則式(15)的離散形式可寫為

式中,ξk和wk分別為高斯點的位置和權(quán).式(7)和(8)中的系數(shù)矩陣TCC可用式(16)計算得到.式(16)中含有裂紋張開位移δ(δ1δ2)的導(dǎo)數(shù),為了便于計算,沿裂紋面對其進(jìn)行Lagrange插值：

式中,lk為NG+2階Lagrange函數(shù),滿足關(guān)系式δ(?a)=δ(+a)=0.設(shè)裂紋總數(shù)為NC,則列向量δ的大小為d·NC·NG,系數(shù)矩陣TCC的大小為(d·NC·NG)2.

需要指出的是,通常裂紋面上的實際面力為0,式(14)～(16)中的面力可理解為虛擬面力[6],由于粒子對裂紋存在影響,裂紋面上的虛擬面力并不等于τ0C(見圖3).式(14)～(16)表明裂紋張開位移與虛擬面力存在特定的關(guān)系,利用這個關(guān)系可更方便地計算應(yīng)力強度因子,避免使用裂紋尖端附近某個位置上的位移或應(yīng)力值.具體而言,當(dāng)問題求解之后,首先利用式(16)由裂紋張開位移求高斯節(jié)點上的虛擬面力,然后將求得的面力用下列多項式擬合,即

式中,ci為多項式待定系數(shù),NP為多項式的階數(shù).求得待定系數(shù)ci后,代入下式：

便可得到計算應(yīng)力強度因子的公式如下：

式中,KL和KR分別表示裂紋左端和右端的應(yīng)力強度因子.事實上,這種計算應(yīng)力強度因子的方法充分利用了裂紋的整體信息,即裂紋張開位移和虛擬面力,從而使得應(yīng)力強度因子的計算結(jié)果更為精確和穩(wěn)定.

2　算例

如前所述,本工作的算例中采用邊界積分方程的兩種離散形式分別對粒子邊界和裂紋表面進(jìn)行離散,其中對粒子邊界的離散采用了邊界點法,N為粒子邊界的節(jié)點數(shù)；對裂紋表面采用高斯配點法進(jìn)行離散,NG為裂紋表面的高斯配點數(shù).

2.1裂紋與雙孔

首先考慮裂紋與雙孔問題(見圖4),目的是校核計算機程序.由于該算例的幾何構(gòu)型與加載條件具有對稱性,裂紋兩端的應(yīng)力強度因子(stress intensity factor,SIF)的值相同,而且滑開型(Ⅱ型)的應(yīng)力強度因子為0,因此只需計算張開型(Ⅰ型)的應(yīng)力強度因子.

對于圖4(a)中的情形,令

其中裂紋的無量綱應(yīng)力強度因子

圖4　裂紋與雙孔Fig.4　Crack and two holes

表1　縱向加載情況下裂紋的無量綱應(yīng)力強度因子Table 1　Dimensionless SIF of the crack for longitudinal load

表2　水平加載情況下裂紋的無量綱應(yīng)力強度因子Table 2　Dimensionless SIF of the crack for horizontal load

2.2裂紋與粒子的方位

考察兩種情況下裂紋與各種不同硬度粒子的方位變化對應(yīng)力強度因子的影響,結(jié)果如圖5所示.第一種情況為粒子方位變化時裂紋保持水平位置不變(見圖5(a)),第二種情況為粒子方位變化時裂紋的方向同步變化(見圖5(b)).粒子的硬度,即模量比的變化范圍從0到100.對裂紋面和粒子界面分別采用NG=15和N=40個節(jié)點進(jìn)行了離散.計算結(jié)果如圖6～9所示.

圖5　裂紋與粒子的方位Fig.5　Positions between crack and particle

圖6　裂紋左端和右端Ⅰ型應(yīng)力強度因子與粒子方位角的關(guān)系(水平裂紋)Fig.6　TypeⅠSIF at left tip and right tip of crack as a function of particle bearing angle (horizontal crack)

圖7　裂紋左端和右端Ⅱ型應(yīng)力強度因子與粒子方位角的關(guān)系(水平裂紋)Fig.7　TypeⅡSIF at left tip and right tip of crack as a function of particle bearing angle (horizontal crack)

圖6為裂紋保持水平位置不變時,左端和右端Ⅰ型應(yīng)力強度因子與粒子方位角的關(guān)系.圖7為Ⅱ型應(yīng)力強度因子與粒子方位角的關(guān)系.由圖6和7可知,由于裂紋右端與粒子的距離較近,Ⅰ型和Ⅱ型應(yīng)力強度因子的變化幅度都比裂紋左端的變化幅度要大.值得注意的是,軟粒子與硬粒子對應(yīng)力強度因子變化幅度的影響恰恰相反.

圖8為裂紋方向同步變化時裂紋左端和右端Ⅰ型應(yīng)力強度因子與粒子方位角的關(guān)系.由圖可知,無論在裂紋的左端還是右端,Ⅰ型應(yīng)力強度因子隨粒子方位角均呈單調(diào)變化,但裂紋右端的Ⅰ型應(yīng)力強度因子對模量比的變化比左端更加敏感.圖9為裂紋左端和右端Ⅱ型應(yīng)力強度因子與粒子方位角的關(guān)系.由圖可知,無論在裂紋的左端還是右端,Ⅱ型應(yīng)力強度因子隨粒子方位角均呈先增后減的變化,裂紋右端的Ⅱ型應(yīng)力強度因子對模量比的變化同樣比左端更加敏感.

圖8　裂紋左端和右端Ⅰ型應(yīng)力強度因子與粒子方位角的關(guān)系(傾斜裂紋)Fig.8　TypeⅠSIF at left tip and right tip of crack as a function of particle bearing angle (inclined crack)

圖9　裂紋左端和右端Ⅱ型應(yīng)力強度因子與粒子方位角的關(guān)系(傾斜裂紋)Fig.9　TypeⅡSIF at left tip and right tip of crack as a function of particle bearing angle (inclined crack)

2.3裂紋與粒子的長寬比

圖10為裂紋和粒子的長寬比.對裂紋面和粒子界面分別采用NG=15和N=40個節(jié)點進(jìn)行離散,計算粒子長寬比對附近裂紋的應(yīng)力強度因子的影響,結(jié)果如圖11所示.由圖11可知,當(dāng)粒子長寬比變化時,裂紋左端和右端應(yīng)力強度因子的變化趨勢是相同的,只是右端應(yīng)力強度因子變化的幅度較大.在圖10所示構(gòu)型的情況下,具有較小長寬比的軟粒子使得應(yīng)力強度因子增大,而具有較大長寬比的硬粒子使得應(yīng)力強度因子減小.

圖11　裂紋左端和右端Ⅰ型應(yīng)力強度因子與粒子長寬比的關(guān)系Fig.11　TypeⅠSIF at left tip and right tip of crack as a function of particle aspect ratio

圖12　裂紋與粒子的水平距離Fig.12　Horizontal distance between crack and particle

圖13　裂紋應(yīng)力強度因子與圓形孔洞或圓形粒子水平距離的關(guān)系Fig.13　SIF as a function of horizontal distance of circular hole or particle

2.4裂紋與粒子(或孔洞)的水平距離

裂紋與粒子(或孔洞)水平排列,如圖12所示.對裂紋面和粒子界面分別采用NG= 15和N=40個節(jié)點進(jìn)行離散,考察裂紋與粒子(或孔洞)的水平距離對應(yīng)力強度因子的影響,計算結(jié)果如圖13和14所示.由圖13可知,圓形孔洞對應(yīng)力強度因子的影響比圓形粒子要大,特別是當(dāng)水平距離較小時.圖14的結(jié)果與圖13相似,方形孔洞對應(yīng)力強度因子的影響比方形粒子要大.比較圖13和14的結(jié)果可知,同樣條件下方形孔洞對應(yīng)力強度因子的影響比圓形孔洞要大,因為方形孔洞對周圍基體應(yīng)力場的干擾更大.

圖14　裂紋應(yīng)力強度因子與方形孔洞或方形粒子水平距離的關(guān)系Fig.14　SIF as a function of horizontal distance of square hole or particle

3　結(jié)束語

針對含有缺陷的非均質(zhì)彈性固體,即同時含有粒子(或孔洞)和裂紋的情況,建立了全空間中粒子和裂紋位移間斷形式的對偶邊界積分方程計算模型,解決了全空間條件下難以對所研究對象進(jìn)行直接加載的問題.采用邊界積分方程的離散形式對含有少量粒子和裂紋的典型情況進(jìn)行了數(shù)值分析,其中對粒子邊界(或界面)采用邊界點法進(jìn)行離散,對裂紋面采用高斯配點法進(jìn)行離散.采用裂紋整體信息,即通過與裂紋張開位移對應(yīng)的裂紋虛擬面力來計算應(yīng)力強度因子,可以獲得穩(wěn)定和精確的計算結(jié)果.

對含有少量粒子和裂紋的典型情況進(jìn)行了數(shù)值分析,計算了裂紋的應(yīng)力強度因子,通過與已有研究結(jié)果的比較,對計算模型與計算機程序的正確性與可靠性進(jìn)行了驗證,并探討了粒子與裂紋的相互作用.

[1]ALIABADI M H,ROOKE D P.Numerical fracture mechanics[M].Southampton：Computational Mechanics Publications,1991：140-192.

[2]CRUSET A.Boundary element analysis in computational fracture mechanics[M].Dordrecht：Kluwer Academic Publishers,1988：17-60.

[3]黎在良,王元漢,李廷芥.斷裂力學(xué)中的邊界數(shù)值方法[M].北京：地震出版社,1996：237-265.

[4]YAN X Q.Stress intensity factors for cracks emanating from a triangular or square hole in an infinite plateby boundary elements[J].Engineering Failure Analysis,2005,12(3)：362-375.

[5]YAN X Q.A brief evaluation of approximation methods for microcrack shielding problems[J]. ASME Journal of Applied Mechanics,2006,73(4)：694-696.

[6]MA H,GUO Z G,DHANASEKAR M,et al.Efficient solution of multiple cracks in great number using eigen COD boundary integral equations with iteration procedure[J].Engineering Analysis with Boundary Elements,2013,37(3)：487-500.

[7]GUO Z,LIU Y,MA H,et al.A fast multipole boundary element method for modeling 2-D multiple crack problems with constant elements[J].Engineering Analysis with Boundary Elements,2014, 47(1)：1-9.

[8]DONG C Y,CHEUNG Y K,LO S H.An integral equation approach to the inclusion-crack interactions in three-dimensional infinite elastic domain[J].Computational Mechanics,2002,29(4/5)：313-321.

[9]DONG C Y,LO S H,CHEUNG Y K.Numerical analysis of the inclusion-crack interactions using an integral equation[J].Computational Mechanics,2003,30(2)：119-130.

[10]DONG C Y,LEEK Y.A new integral equation formulation of two-dimensional inclusion-crack problems[J].International Journal of Solids and Structures,2005,42(18/19)：5010-5020.

[11]DONG C Y,LEE K Y.Boundary element analysis of infinite anisotropic elastic medium containing inclusions and cracks[J].Engineering Analysis with Boundary Elements,2005,29(6)：562-569.

[12]DONG C Y.The integral equation formulations of an infinite elastic medium containing inclusions,cracks and rigid lines[J].Engineering Fracture Mechanics,2008,75(13)：3952-3965.

[13]DONG C Y,YANG X,PAN E.Analysis of cracked transversely isotropic and inhomogeneous solids by a special BIE formulation[J].Engineering Analysis with Boundary Elements,2011, 35(2)：200-206.

[14]CHEN Y Z.Boundary integral equation method for two dissimilar elastic inclusions in an infinite plate[J].Engineering Analysis with Boundary Elements,2012,36(1)：137-146.

[15]MA H,ZHOU J,QIN Q H.Boundary point method for linear elasticity using constant and quadratic moving elements[J].Advances in Engineering Software,2010,41(3)：480-488.

[16]TELLES J C F,CASTOR G S,GUIMARAES S.A numerical Green’s function approach for boundary elements applied to fracture mechanics[J].International Journal for Numerical Methods in Engineering,1995,38(19)：3259-3274.

[17]HUI C Y,SHIA D.Evaluations of hypersingular integrals using Gaussian quadrature[J].International Journal for Numerical Methods in Engineering,1999,44(2)：205-214.

[18]MURAKAMI Y,KEER L M.Stress intensity factors handbook[J].Journal of Applied Mechanics, 1993,60(4)：1063.

本文彩色版可登陸本刊網(wǎng)站查詢：http://www.journal.shu.edu.cn

Dual boundary integral equations and numerical solutions for particles and cracks in full space

MA Hang1,PAN Meng2
(1.College of Sciences,Shanghai University,Shanghai 200444,China；
2.Shanghai Institute of Applied Mathematics and Mechanics,Shanghai University, Shanghai 200072,China)

As elastic solids contain particles and cracks,a computational model is proposed for the analysis of particles and cracks in full space using dual boundary integral equations in the displacement discontinuity formulation.The method avoids the difficulty of loading the objects under study in full space.Numerical analysis is carried out for some typical cases with a few particles and cracks using a discrete form of boundary integral equations.A boundary point method and Gauss collocation are used for discretization of the particle boundary or interface,and the crack surface,respectively.The stress intensity factors of cracks are computed.Mutual effects between particles and cracks are investigated and compared with those in the literature,verifying correctness and reliability of the proposed computational model and the program.

particle；crack；dual boundary integral equation；numerical solution；stress intensity factor

O 241

A

1007-2861(2016)05-0533-12

10.3969/j.issn.1007-2861.2015.02.020

2015-05-11

國家自然科學(xué)基金資助項目(11272195,11332005)

馬杭(1951—),男,教授,博士生導(dǎo)師,博士,研究方向為計算固體力學(xué). E-mail:hangma@staff.shu.edu.cn