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

馬杭,潘蒙

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

2.上海大學上海市應用數學和力學研究所,上海　200072)

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

馬杭1,潘蒙2

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

2.上海大學上海市應用數學和力學研究所,上海200072)

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

粒子；裂紋；對偶邊界積分方程；數值解；應力強度因子

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

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

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

1　計算模型與求解

1.1邊界積分方程

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

1.2區域的分解與加載

圖1　全空間中受遠場單位拉伸載荷的粒子和裂紋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　受遠場拉伸的均質全空間與不受遠場載荷的孔洞和裂紋Fig.3　Full uniform space under far-field tension and the hole and crack without far-field load

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

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

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

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

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

1.4離散方式與應力強度因子的計算

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

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

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

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

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

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

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

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

便可得到計算應力強度因子的公式如下：

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

2　算例

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

2.1裂紋與雙孔

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

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

其中裂紋的無量綱應力強度因子

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

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

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

2.2裂紋與粒子的方位

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

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

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

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

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

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

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

圖9　裂紋左端和右端Ⅱ型應力強度因子與粒子方位角的關系(傾斜裂紋)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個節點進行離散,計算粒子長寬比對附近裂紋的應力強度因子的影響,結果如圖11所示.由圖11可知,當粒子長寬比變化時,裂紋左端和右端應力強度因子的變化趨勢是相同的,只是右端應力強度因子變化的幅度較大.在圖10所示構型的情況下,具有較小長寬比的軟粒子使得應力強度因子增大,而具有較大長寬比的硬粒子使得應力強度因子減小.

圖11　裂紋左端和右端Ⅰ型應力強度因子與粒子長寬比的關系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　裂紋應力強度因子與圓形孔洞或圓形粒子水平距離的關系Fig.13　SIF as a function of horizontal distance of circular hole or particle

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

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

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

3　結束語

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

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

[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]黎在良,王元漢,李廷芥.斷裂力學中的邊界數值方法[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.

本文彩色版可登陸本刊網站查詢：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

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

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