比例邊界有限元模擬裂紋和夾雜動力相互作用

施明光， 徐艷杰， 張楚漢， 劉鈞玉

(1.清華大學 水沙科學與水利水電工程國家重點實驗室，北京　100084；2. 大連理工大學 建設工程學部水利學院，遼寧 大連　116024;3.三明市水利局,福建 三明　365000；4.沈陽工業大學 建筑工程學院，沈陽　110870)

?

比例邊界有限元模擬裂紋和夾雜動力相互作用

施明光1,3， 徐艷杰1， 張楚漢1， 劉鈞玉2,4

(1.清華大學 水沙科學與水利水電工程國家重點實驗室，北京100084；2. 大連理工大學 建設工程學部水利學院，遼寧 大連116024;3.三明市水利局,福建 三明365000；4.沈陽工業大學 建筑工程學院，沈陽110870)

摘要：基于三角形背景網格，任意結構可用n(n≥3)邊多邊形比例邊界有限元 (Polygon Scaled Boundary Finite Elements, PSBFE)自動離散。相對以往基于比例邊界有限元(SBFEM)的應用，該多邊形單元不但繼承SBFEM半解析求解裂紋尖端奇異性的特性，而且在模擬復雜結構的網格生成和裂紋擴展上具有更高的通用性。首次用該單元模擬了動荷載下復合材料裂紋和夾雜相互作用。動荷載穩定裂紋情況下，PSBFE計算結果同現有文獻吻合良好，在此基礎上，結合基于拓撲的局部網格重剖分方法，模擬了動荷載下夾雜和擴展裂紋相互作用。結果表明，硬質夾雜和軟質夾雜對結構的動力應力強度因子分別起到抑制和放大的作用。夾雜尺寸，夾雜大小也會在一定范圍內影響動力應力強度因子，尺寸越大距離裂紋越近的夾雜影響越大。

關鍵詞：裂紋擴展；復合材料；多邊形單元；網格重剖分；動力應力強度因子

纖維復合材料，顆粒增強材料在工程中得到廣泛應用。一方面，通過在基體中加入異質材料(夾雜物)可以顯著提升復合材料的性能。另一方面，孔洞，裂紋等生產缺陷的存在又會導致材料性能發生劣化。研究基體中夾雜物和裂紋的相互作用對于理解這類復合材料的強化和破壞機制有重要作用。這類問題的解析方法多限于規則的幾何形狀，實驗研究也大部分限于擬靜力實驗，動力實驗并不多見。Jajam等[1]近年來在這方面做了嘗試。他們用氨基環氧樹脂做為基體樹脂，分別采用硼硅酸鹽玻璃和聚亞安酯作為硬質夾雜和軟質夾雜(彈性模量高于基體的夾雜稱為硬質夾雜，反之稱為軟質夾雜), 并考慮夾雜和基體的黏結強度(強黏結和弱黏結)，研究動力荷載作用下裂紋和夾雜的相互作用,證明夾雜物和基體的彈性模量比，夾雜的位置和夾雜-基體的黏結強度等對試樣的斷裂路徑，斷裂速度和斷裂韌度均有影響。除實驗和解析研究外，有限元(FEM),邊界元(BEM),擴展有限元(XFEM)等數值方法也被廣泛運用到復合材料的研究中。如Lipetzky等[2-3]采用FEM研究靜荷載下，夾雜大小，位置和形狀等因素對穩定單邊裂紋的平板應力強度因子的影響，Willams[4]采用BEM對含圓心夾雜的三點彎梁進行了擬靜力加載的裂紋擴展模擬。同實驗一樣，動力荷載下裂紋夾雜的相互作用數值分析較少。雷鈞等[5-6]近年在這方面做了很多研究工作。他們采用時域邊界元法(TDBEM)先后研究了基體裂紋在多夾雜影響下的動力應力強度因子變化[6]，微裂紋和夾雜對基體中主裂紋擴展路徑，速度，應力強度因子的影響等[7]。

比例邊界有限元(Scaled Boundary Finite Element Method，SBFEM)是由Song等[8-9]提出的一種半解析數值方法，最初是為了研究無限域的波動問題。該法結合了FEM和BEM的優點，同時具有自身的優點。同FEM采用奇異單元和XFEM采用富集函數模擬裂紋尖端奇異性不同，通過將比例中心放置在裂縫尖端，SBFEM無須事先知道裂紋尖端奇異性，可以直接從半解析的應力或者位移解中提取應力強度因子，因此對于任意材料的裂紋問題均可以采用相同方法求解。Song[10-11]首先將SBFEM運用到各向異性材料的應力強度因子求解，發展至今該方法已經被成功運用到溫度荷載，瞬時動力荷載下的應力強度因子求解[12]，單一材料線性[13]，靜力非線性[14]，線彈性動力裂紋擴展模擬[15]等。然而，如果只用少數幾個子域，SBFEM在應用到復雜模型時候就顯得有點困難。文獻[16]通過三角剖分將模型離散成任意的n(n≥3)多邊形比例邊界單元(Polygon Scaled Boundary Finite Elements ，PSBFE)，每個多邊形單元相當于比例邊界元的子域(Subdomains)。同SBFEM一樣，它可以精確求解裂紋尖端的奇異性，且相對以往的基于SBFEM的文獻，該剖分方法具有更高的通用性，可以應用到任意形態的模型。線彈性斷裂力學中,動應力強度因子(Dynamic Stress Intensity Factors, DSIFs)是重要的斷裂參數。本文首先闡述PSBFE的網格生成方法及其半解析求解裂紋尖端奇異性的原理，并介紹基于PSBFE的局部網格重剖分技術及新舊網格映射方法。在此基礎上，分別考慮動載下裂紋穩定和擴展兩種情況，采用PSBFE模擬夾雜對DSIFs的影響。

1多邊形比例邊界有限元(PSBFE)

1.1PSBFE的網格生成

圖1(a)所示帶裂紋的復合材料由基質和夾雜組成。對于該模型，首先采用三角形單元進行離散，得到背景網格圖1(b)。 然后對于圖1(b)中每個三角形節點：

1.1.1若其不在邊界或材料交界面(如圖1(c)中的R點)，則將和其相連的三角形的形心連接形成多邊形單元(圖1(c)點1～6)。該多邊形任意滿足到所有邊界可視的點均可作為比例中心。一般情況下直接以R點為比例中心。

圖1　多邊形單元生成示意圖Fig.1 Illustration of polygon element mesh generation

1.1.2若其位于裂紋尖端(如圖1(d)的O點)，首先在包含該點且位于裂紋面上的三角形邊的中點生成2個結點(圖1(d)點4，5)，然后將中點和周邊三角形形心一起連接形成多邊形單元(圖1(d)點1～8)。比例中心位于該裂紋尖端。

1.1.3若其位于邊界(如圖1(e)的P點)，首先在包含P點且位于邊界的三角形邊中點處生成節點(圖點1,6和)，然后將該三角形節點同周邊中點，形心點逆時針連接行成P-1-2-3-4-5-6的多邊形單元。

1.1.4若其位于材料交界面(如圖1(f)中的Q點)，類似圖1(c)首先生成中點2,3，然后按照材料劃分，形成Q-3-4-1-2，Q-2-5-6-3的2個多邊形單元。

1.2PSBFE動力應力強度因子求解

圖1(c)的多邊形單元的比例中心位于R點，徑向坐標ξ和法向坐標η稱為比例坐標，(r(θ),θ)為原點在比例中心的極坐標系。該多邊形單元在邊界用6個2節點線單元進行離散。不考慮體力和裂紋面上荷載時，以位移u(ξ)為變量的SBFEM的動力學控制方程為[8]：

[E0]ξ2{u(ξ)},ζξ+([E0]-

[E1]+[E1T])ξ{u(ξ)},ξ-

[E2]{u(ξ)}+ω2[M0]ξ2{u(ξ)}=0

(1)

式中，系數矩陣[E0]，[E1]，[E2]，[M0]只和多邊形單元形態和材料參數有關，詳細表達式見文獻[8]。ω為荷載激勵頻率，ω=0時，式(1)表示靜力控制方程。引入Hamiliton矩陣[Z]并進行分塊Schur分解[17]

(2)

[Z][Ψ]=[Ψ][S]

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

{u(ξ,η)}=[N(η)]{u(ξ)}

(11)

(12)

(13)

(14)

式(7)～式(10)中，[N(η)]為線單元的插值函數，{c}為多邊形單元的積分常數，[Ψσi(η)]為應力模態(Stress Modes)[17]。應力{σ(ξ,η)}可看成由一系列應力分量{σi(ξ,η)}疊加而成。對于裂尖多邊形單元，采用PBFEM求解會得到實特征值介于-1和0的分塊[Ss]，這樣當ξ趨于0(裂尖)時，表達式ξ-[Ss]-[I]值將趨于無窮大。由于應力的奇異性包含在解的特征值里面，所以PSBFE無須像FEM那樣事先假定裂紋尖端奇異性，相比而言適應性不但具有更高精度且適用性更廣。在求得對應分塊[Ss]的應力分量{σs(ξ,η)}后，首先根據靜態應力強度因子定義求解等效靜應力強度因子{Keq}。

(15)

(16)

(17)

(18)

式中，c為裂紋尖端速度，cd，cs，cR分別代表縱波，剪切波和Rayleigh波波速。

1.3基于拓撲的單元局部網格重剖分技術

文獻[18]提出了基于拓撲的局部網格重剖分計算方法。對于圖2(a)的多邊形單元，假定裂紋尖端從O移動到O′點，通過更新位于O-B-C-D-E-F-G-H-I-J-K-L-O′內的背景三角形網格，再依照1.1節的方法，就可以生成擴展后新的多邊形網格(圖2(b))。具體步驟參見文獻[18]。

圖2　局部重剖分過程Fig.2 Local remeshing scheme

1.4網格映射

進行動力裂紋擴展模擬時，由于網格形態的變換，需要將舊網格的狀態映射到新網格，由于任意一點的位移，速度，加速度在PSBFE中均是半解析表達，所以只需根據式(11)～(13)就可得到新網格的狀態，且保持很高的精度。

2算例

2.1動載下單夾雜對穩定裂紋DSIFs的影響

2.1.1不同材料組合影響

當d=10 mm，r=4 mm時，該問題同采用TDBEM的文獻[6]參數一致。選取3種材料組合(剪切模量比μ01=0.5,1,2)分析硬質夾雜和軟質夾雜對結構動力應力強度因子的影響。多邊形網格見圖4。圖5給出了PSBFE和BEM[6]計算結果對比，可以看出兩種方法計算結果吻合較好。加載初段，3種材料組合動力應力強度因子很小，當縱波第一次從荷載邊界到達裂紋尖端時候(W/cs=2.73 μs)，動力應力強度因子開始增加，但是幅度不一樣，硬質材料相當于起到一個抑制DSIFs的效果(shielding effect)，導致DSIFs峰值降低，而軟質夾雜則對DSIFs有放大作用(amplication effect)。

圖3 帶邊裂紋和圓心夾雜的方形平板Fig.3Asquareplatewithanedgecrackandacircularinclusion圖4 多邊形網格Fig.4Polygonalmesh

2.1.2夾雜半徑影響(硬質夾雜)

以硬質夾雜μ=0.5為例，考慮4種半徑不同大小的夾雜(r=4,3.5,3,1 mm)對結構DSIFs的影響。圖6的結果進一步確定硬質夾雜對DSIFs起抑制效果，夾雜半徑越小，該效果越小，當夾雜半徑接近1 mm時，幾乎和均質材料響應一樣。

2.1.3夾雜位置不同的影響(硬質夾雜)

圖7給出了3種不同位置的硬質夾雜的DSIFs對比。可以看出，夾雜和裂紋越接近，夾雜帶來的抑制效果越明顯，根據圖中趨勢可以推斷出，當夾雜和裂紋距離大到一定程度時，夾雜的抑制效果可以忽略不計。

圖5 穩定裂紋下材料組合影響Fig.5Effectofthematerialcombinationforstationarycrack圖6 穩定裂紋下夾雜半徑影響Fig.6Effectoftheradiusofinclusionforstationarycrack圖7 穩定裂紋下夾雜位置影響Fig.7Effectoftheinclusionpositionforstationarycrack

2.2動載下單夾雜對擴展裂紋DSIFs的影響

為了研究硬質和軟質夾雜在裂紋擴展中的作用，將2.1節中平板夾雜取半徑r=4 mm的并將其移至d=15 mm。同樣選取不同材料組合(剪切模量比μ01=0.5,1,2)進行分析。參照文獻[19]的算例，假定裂紋在3 μs前處于靜止狀態，然后按照1 000 m/s速度沿著水平方向擴展，這相當于動力裂紋擴展模擬中的Mixed-Phase問題[20]。圖8給出了時刻t=0,4.25,8.01 μs的局部裂紋擴展形態，可以看出，基于拓撲的局部網格重剖分方法僅改變裂紋擴展路徑周邊的單元形態，大部分單元形態得以保留，這節省了裂紋擴展的計算量。圖9的結果表明，和動力荷載下穩定裂紋類似，硬質夾雜和軟質夾雜對動力荷載擴展裂紋的DSIFs分別起到抑制和放大作用。

圖8　裂紋擴展過程局部放大圖Fig.8 Close up of crack propagation process of square plate with edge crack under impact load

圖9　裂紋擴展下材料組合影響Fig.9 Effect of the material combination for crack growth

2.3動載下多夾雜對穩定和擴展裂紋DSIFs影響

圖11為一個4個圓形夾雜的方形平板，邊長2W=40 m。平板中心有長2a=8 mm的直裂紋，4個圓心夾雜圓心到邊界的距離均為d=10 mm。半徑均為r=4 mm。基質和夾雜的材料屬性，外部荷載，應力強度因子歸一化同算例2.1。考慮平面應變問題。由問題的對稱性，取結構右半部分進行分析，并在原結構對稱軸上施加水平向位移約束。圖11給出背景三角形網格和相應的多邊形網格。圖12給出了動荷載下穩定裂紋相應材料組合下的應力強度因子變化規律。可以看出，PSBFE和TDBEM[6]的計算結果基本一致。在此基礎上，采用和算例2.2一致的裂紋擴展速度研究多夾雜下動荷載下裂紋擴展過程中應力強度因子變化規律。圖13和圖14分別給出了裂縫擴展過程中的網格和DSIFs變化情況。從圖13可以看出，整個裂縫擴展過程，只有局部網格發生變化。同算例2.1和2.2一樣，從圖12和圖14可以觀察到硬質夾雜和軟質夾雜對DSIFs分別起到抑制和放大作用這一現象，圖14還可以看出當裂縫尖端離夾雜較遠時候(t≥18 μs)時候，這種抑制作用減弱比較明顯。

圖10　帶中心裂紋和4個圓心夾雜的方形平板Fig.10 A square plate with an central crack and four circular inclusions

圖11　中心裂紋平板背景三角形網格和相應多邊形網格(745個多邊形單元)Fig.11 Background triangular mesh and its polygonal mesh of a square plate with an central crack (745 polygon element)

圖12　穩定裂紋下材料組合影響Fig.12 Effect of the material combination for stationary crack

圖13　裂紋擴展過程局部放大圖Fig.13 Close up of crack propagation process

圖14　裂紋擴展下材料組合影響Fig.14 Effect of the material combination for crack growth

3結論

本文采用多邊形比例邊界有限元模擬動力荷載作用下復合材料夾雜和穩定或者移動裂紋的相互作用。結果表明：

(1) PSBFE可以半解析求解裂紋尖端的奇異性以及單元內任意一點的狀態變量(位移，加速度等)，在裂紋模擬中具有高精度。

(2) 基于拓撲的PSBFE網格重剖分辦法在裂縫擴展中只改變裂縫附近的局部網格形態，使得動力裂縫擴展過程中的網格映射更為方便，精度也更高。

(3) PSBFE可以有效模擬復合材料的裂紋和夾雜的動力相互作用，算例的結果同已有文獻結果吻合良好。

(4) 一般說來，硬質夾雜在基體中對動力應力強度因子起抑制作用，而軟質夾雜則起放大作用。

(5) 夾雜越大，越靠近裂紋，對動力應力強度因子影響越大，反之，影響越小。

致謝

本文研究過程中，得到澳大利亞新南威爾士大學宋崇民教授和Ean Tat博士的大量幫助和指導，特此表示感謝。

參 考 文 獻

[ 1 ] Jajam K C, Tippur H V. Role of inclusion stiffness and interfacial strength on dynamic matrix crack growth: an experimental study[J]. International Journal of Solids and Structures,2012, 49: 1127-1146.

[ 2 ] Lipetzky P, Knesl Z. Crack-particle interaction in a two-phase composite Part II: crack deflection[J]. International Journal of Fracture,1995, 73(1): 81-92.

[ 3 ] Lipetzky P, Schmauder S. Crack-particle interaction in two-phase composites Part I: particle shape effects[J]. International Journal of Fracture,1994, 65(4): 345-358.

[ 4 ] Williams R C, Phan A V, Tippur H V, et al. SGBEM analysis of crack-particle (s) interactions due to elastic constants mismatch[J]. Engineering Fracture Mechanics,2007, 74(3): 314-331.

[ 5 ] Lei J, Wang Y S, Huang Y, et al. Dynamic crack propagation in matrix involving inclusions by a time-domain BEM[J]. Engineering Analysis with Boundary Elements,2012, 36(5): 651-657.

[ 6 ] Lei J, Yang Q, Wang Y S, et al. An investigation of dynamic interaction between multiple cracks and inclusions by TDBEM[J]. Composites Science and Technology,2009,69(7): 1279-1285.

[ 7 ] Lei J, Zhang C, Yang Q, et al. Dynamic effects of inclusions and microcracks on a main crack[J]. International Journal of Fracture,2010, 164(2): 271-283.

[ 8 ] Song C, Wolf J P. The scaled boundary finite-element method-alias consistent infinitesimal finite-element cell method-for elastodynamics[J]. Computer Methods in Applied Mechanics and Engineering,1997, 147(3/4): 329-355.

[ 9 ] Song C M, Wolf J P. Body loads in scaled boundary finite-element method[J]. Computer Methods in Applied Mechanics and Engineering,1999, 180(1/2): 117-135.

[10] Song C W J. The scaled boundary finite-element method for anisotropic multimaterial plate with crack.[Z]. Palaiseau, France: 1998.

[11] Song C. Analysis of singular stress fields at multi-material corners under thermal loading[J]. International Journal for Numerical Methods in Engineering,2006, 65(5): 620-652.

[12] Yang Z J, Deeks A J, Hao H. Transient dynamic fracture analysis using scaled boundary finite element method: a frequency-domain approach[J]. Engineering Fracture Mech-anics,2007,74(5): 669-687.

[13] Yang Z J. Fully automatic modelling of mixed-mode crack propagation using scaled boundary finite element method[J]. Engineering Fracture Mechanics,2006, 73(12):1711-1731.

[14] Yang Z J, Deeks A J. Modelling cohesive crack growth using a two-step finite element-scaled boundary finite element coupled method[J]. International Journal of Fracture,2007,143(4): 333-354.

[15] Ooi E T, Yang Z J. Modelling dynamic crack propagation using the scaled boundary finite element method[J]. International Journal for Numerical Methods in Engineering,2011,88(4): 329-349.

[16] Ooi E T, Song C, Tin Loi F, et al. Polygon scaled boundary finite elements for crack propagation modelling[J]. International Journal for Numerical Methods in Engineering,2012,91(3): 319-342.

[17] Song C, Tin-Loi F, Gao W. A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges[J]. Engineering Fracture Mechanics,2010, 77(12): 2316-2336.

[18] Ooi E T, Mingguang S, Song C, et al. Dynamic crack propagation simulation with scaled boundary polygon elements and automatic remeshing technique[J]. Engineering Fracture Mechanics,2013, 106: 1-21.

[19] Fedelinski P, Aliabadi M H, Rooke D P. The time-domain DBEM for rapidly growing cracks[J]. International Journal for Numerical Methods in Engineering,1997, 40(9): 1555-1572.

[20] Nishioka T. Computational dynamic fracture mechanics[J]. International Journal of Fracture,1997, 86(1): 127-159.

Simulation of dynamic interactions between a crack and inclusions with scaled boundary finite element method

SHIMing-guang1,3,XUYan-jie1,ZHANGChu-han1,LIUJun-yu2,4

(1.State Key Laboratory of Hydro Science and Hydraulic Engineering, Tsinghua University, Beijing 100084, China;2.Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, China;3.Water conservancy bureau of Sanming, Sanming 365000, China;4.School of Architecture & Civil Engineering, Shenyang University of Technology, Shenyang 110870, China)

Abstract:Any structural domain can be discretized automatically with a mesh of arbitrary n-sided (n≥3) polygon scaled boundary finite elements (PSBFE) based on Delaunay triangulation background mesh. Compared with previous literatures based on SBFEM, PSBFE retains the characteristics of SBFEM’s accurately representing orders of singularities at crack tips it is more general and flexible in modeling complicated structures and their crack propagation. Here, PSBFE was for the first time applied to simulate the dynamic interactions between a crack and inclusions in composite material. The numerical results of stationary cracks under dynamic load were consistent with available data in literatures. Next, a local remeshing scheme was employed to simulate the dynamic crack propagation. The numerical results demonstrated that stiff and soft inclusions have the restraining and amplification effects on the dynamic stress intensity factor of a structure; the sizes and positions of inclusions also affect the dynamic stress intensity factor, the larger the size and the closer the inclusion, the more the effects.

Key words:crack propagation; composite material; polygon elements; grid remeshing; dynamic stress intensity factor

中圖分類號:O346.1

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.04.003

通信作者張楚漢 男，教授，中科院院士，1933年生

收稿日期：2014-12-05修改稿收到日期：2015-02-06

基金項目：國家自然科學基金(40974063;41274106;51109134);水沙科學與水利水電工程國家重點實驗室科研課題(2011-KY-3)

第一作者 施明光 男，博士，1984年生

E-mail：zch-dhh@tsinghua.edu.cn