黃小華， 李雙， 金艷麗， 何小橋， 王家明

(1.廣西大學 土木建筑工程學院，廣西 南寧 530004；2. 工程防災與結構安全教育部重點實驗室，廣西 南寧 530004；3. 香港城市大學 建筑學及土木工程學系，香港 999077)

材料在生產(chǎn)和使用過程中會產(chǎn)生一系列細小裂紋而帶縫工作，這些裂縫極有可能導致材料的破壞，甚至結構的失效[1]。因此在材料中預制裂紋，研究裂紋擴展規(guī)律是解決含缺陷材料破壞問題的一種較為常見的手段。有限元等基于傳統(tǒng)連續(xù)介質(zhì)力學的數(shù)值方法在結構變形分析中取得了巨大的成功，但是它在不連續(xù)處沒有定義，存在奇異性，在解決材料的破壞問題上舉步維艱。新近興起的近場動力學理論[2]采用位移的空間積分形式完全重構了傳統(tǒng)連續(xù)介質(zhì)力學理論，這種積分型方程在不連續(xù)處仍有定義，且形式不變，實現(xiàn)了對連續(xù)介質(zhì)和非連續(xù)介質(zhì)力學行為的統(tǒng)一描述，能模擬材料中裂紋自發(fā)地萌生和擴展，是一種特別適合模擬材料破壞的力學理論和方法。Silling[2]提出了鍵基近場動力學理論，它將物質(zhì)點間的相互作用視為它們中心受力的鍵，這種中心力模型只有一個參數(shù)，即鍵剛度(或微觀模量)，稱為傳統(tǒng)鍵基近場動力學模型(或單參數(shù)鍵基近場動力學模型，bond-based peridynamics，BPD)。

均勻各向同性線彈性材料的力學屬性須由2個獨立的彈性參數(shù)才能完全表述，采用單參數(shù)鍵基PD模型描述其力學行為時，2個獨立彈性參數(shù)將被縮減為1個參數(shù)，必然導致某一個彈性參數(shù)將以恒定的形式固化到中心力模型中。例如泊松比參數(shù)被固化成常數(shù)，在平面應力狀態(tài)下，泊松比只能是1/3，平面應變狀態(tài)和三維應力狀態(tài)下，泊松比只能是1/4[3]。毫無疑問，這種泊松比恒定的要求大大限制了傳統(tǒng)鍵基近場動力學理論的適用范圍。為此，Silling[4]又提出了態(tài)基近場動力學理論，它是通過考慮兩物質(zhì)點鄰域內(nèi)所有鍵的集合作用來增強材料的本構描述。與鍵基近場動力學模型中2物質(zhì)點間的作用僅由它倆自身決定不同，態(tài)基近場動力學模型中2物質(zhì)點間的作用須考慮2物質(zhì)點鄰域內(nèi)所有物質(zhì)點[5]，是一種典型的多體相互作用。因此，態(tài)基近場動力學方法的計算復雜性以及對資源的耗費都較鍵基近場動力學有了更高的要求。于是，為了保留傳統(tǒng)鍵基近場動力學模型的簡單性以及數(shù)值實施的穩(wěn)定性，部分學者仍立足于鍵基近場動力學模型的兩體間相互作用，通過引入切向鍵[6-8]、引入鍵的轉角[9]、力量補償方案[10]或以梁換鍵[11-12]等方式擴展鍵基近場動力學理論，突破了泊松比的限制。其中，通過引入切向鍵或鍵的轉角等方式增加2物質(zhì)點間的切向作用較為便捷，對傳統(tǒng)鍵基近場動力學模型改動很小。Prakash等[9]通過引入法向彈簧，提出一種二維雙參數(shù)鍵基近場動力學模型。Zhou等[6-7]基于超彈性理論提出一種三維兩參數(shù)共軛鍵基近場動力學模型。Zhu[8]等通過增加鍵的轉角，提出雙參數(shù)鍵基近場動力學模型。當泊松比取上述恒定值時(如在平面應力狀態(tài)下，泊松比取1/3)，兩參數(shù)都能退化為單參數(shù)。但一些模型[6-7,9]退化后的參數(shù)與傳統(tǒng)單參數(shù)鍵基近場動力學的參數(shù)并不相同。

本文在BPD模型的基礎上，通過鍵內(nèi)直接添加切向剛度系數(shù)，提出一種單鍵雙參數(shù)鍵基近場動力學模型(TBPD)，推導模型的雙參數(shù)與材料彈性參數(shù)之間的轉換關系，建立TBPD模型單鍵斷裂條件。然后探討不同泊松比下BPD模型和TBPD模型的適用性，最后，采用TBPD模型，研究裂紋間距和隨機分布的細小裂紋對平行雙裂紋薄板破壞的影響。

1 單鍵雙參數(shù)鍵基近場動力學方法

1.1 理論模型

傳統(tǒng)鍵基近場動力學模型認為宏觀連續(xù)體在其所在空間域Ω內(nèi)由大量物質(zhì)點組成，物質(zhì)點以初始構型的位置X作為標記，攜有體積VX和質(zhì)量密度ρ。任一物質(zhì)點X僅在其有限的“域”內(nèi)通過鍵與其他物質(zhì)點X′存在相互作用力f。域是以該物質(zhì)點為中心，以δ為半徑的范圍，記域內(nèi)所有其他物質(zhì)點為Hδ={X′∈Ω:|X′-X|≤δ}，據(jù)牛頓第二定律，物質(zhì)點X在t時刻的運動方程為：

(1)

傳統(tǒng)鍵基近場動力學模型將物質(zhì)點間的相互作用視為它們中心受力的鍵，鍵僅有一個剛度系數(shù)，表示法向作用，故這里稱之為單鍵單參數(shù)鍵基近場動力學模型。

在傳統(tǒng)鍵基近場動力學模型的基礎上，通過鍵內(nèi)直接添加切向剛度系數(shù)，使其具有法向和切向雙剛度，并稱之為單鍵雙參數(shù)鍵基近場動力學模型，如圖1所示。

圖1 TBPD模型中鍵的運動示意Fig.1 Motion of a bond in TBPD model

(2)

此時，X和X′間鍵的微觀彈性應變密度(單位體積的應變能密度)為:

(3)

將每個鍵的應變能都等分給鍵兩端的物質(zhì)點，則物質(zhì)點X的宏觀彈性應變能密度為：

(4)

對點力f是鍵的微觀彈性應變能密度ω對該鍵相對位移矢量η的導數(shù)：

(5)

以示區(qū)別，傳統(tǒng)BPD模型的對點力[14]為：

(6)

(7)

則式(6)可改寫為：

(8)

(9)

式中δij為Kronecker符號。于是，Green應變張量E可表示為：

(10)

小變形時有|?ui/?Xj|?1，略去位移梯度的乘積項，Green應變E變?yōu)镃auchy應變ε，則式(9)變?yōu)椋?/p>

Fij=δij+εij+ωij

(11)

(12)

將方程(3)代入方程(4)，得物質(zhì)點X的宏觀彈性應變能密度：

(13)

將方程(12)代入方程(13)化簡得：

(14)

式中h為材料厚度。在經(jīng)典連續(xù)介質(zhì)力學中，平面應力狀態(tài)的應變能密度為：

(15)

(16)

式(16)中的法向剛度參數(shù)cn與傳統(tǒng)鍵基近場動力學的微模量[14]完全相同。若方程(14)和方程(15)中采用平面應變狀態(tài)下或者三維狀態(tài)下的應變能密度，則可得到對應應力狀態(tài)下TBPD模型的微模量與材料彈性參數(shù)間的轉換關系，如表1所示。

表1 TBPD模型的微模量cn和ctTable 1 Micro-modulus cn and ct of TBPD model

不難驗證，當切向剛度系數(shù)ct=0時，TBPD模型退化后參數(shù)cn與BPD模型參數(shù)完全相同，且有方程(5)右邊表達式將退化成方程(8)右邊表達式，則此時，TBPD模型將近似為BPD模型。微模量參數(shù)cn和ct作為鍵基PD的剛度系數(shù)，必須滿足cn≥0且ct≥0。因此TBPD模型的泊松比適用范圍為：平面應變和三維情形時-1<ν≤1/4，平面應力時-1<ν≤1/3。

值得指出的是，表1轉換關系與Zhu等[8]相同，但兩者引入切向剛度的方法和推導過程并不同。TBPD模型通過鍵內(nèi)直接增加切向剛度系數(shù)，在小變形條件下，以均勻應變工況，得到模型的雙參數(shù)與材料彈性參數(shù)的轉換公式。Zhu等[8]未就雙參數(shù)鍵基PD模型開展破壞分析，本文接下來將討論TBPD模型的斷裂準則及其破壞應用。

1.2 單鍵斷裂條件

這里采用能量準則作為判斷單鍵斷裂的依據(jù)，即鍵中包含的微觀應變能密度ω(η,ξ)超過某個極限值ωc(η,ξ)(簡記為ωc)時，認為該鍵斷裂，將不再傳遞力的作用。

考慮含裂紋CD的無限大板，板厚h，裂紋長度為2a。當裂紋CD穿過某一物質(zhì)點的鄰域時，所有被裂紋穿過的鍵都將斷裂，存儲在這些鍵中的應變能也被耗散。如圖2所示，鍵AB將發(fā)生斷裂，存儲在該鍵中的應變能ωcVAVB將被耗散，其中VA和VB分別為物質(zhì)點A和物質(zhì)點B攜帶的體積，當VA和VB趨于無限小時，所有因裂紋CD而耗散的鍵的總應變能用積分表示為：

圖2 斷裂面未完全穿過物質(zhì)點鄰域的微勢能積分域Fig.2 Integration domain of the micropotentials not completely crossing a fracture surface

(17)

式中：ΩA和ΩB分別表示位于裂紋CD兩側的物質(zhì)點A和物質(zhì)點B的取值范圍，它們須滿足一定的斷鍵條件，即線段AB的距離rAB≤δ，且線段AB和線段CD有交點。

ΩB1+ΩB2=ΩB

(18)

則對于所有距裂紋CD垂直距離為z的物質(zhì)點，因裂紋CD而耗散的鍵的總應變能為：

(19)

式中：dz為物質(zhì)點A攜帶體積在z軸的長度，S′A為裂紋CD在直線AA2所在水平面上的投影區(qū)。

記Gc為產(chǎn)生每單位裂紋表面所需的能量，則：

(20)

與傳統(tǒng)BPD模型相類似，將ωc假定為ξ的線性函數(shù)，即:

ωc(η,ξ)=ωc0ξ

(21)

將式(21)代入式(20)，得:

(22)

故

(23)

由式(3)、式(21)和式(23)知，TBPD模型的單鍵斷裂條件為：

(24)

在采用TBPD方法計算過程中，一旦鍵滿足式(24)，則該鍵斷裂且永久失效。不難驗證，當切向剛度系數(shù)ct=0時，TBPD模型的單鍵斷裂條件可退化后為BPD模型的單鍵斷裂條件，即：

(25)

式中c和sc分別為BPD模型中鍵的微觀模量以及鍵的臨界伸長率。

在TBPD模型中，鍵斷裂與否采用狀態(tài)函數(shù)μ(ξ,t)來記錄。物質(zhì)點的損傷通過物質(zhì)點鄰域內(nèi)鍵的斷裂程度，用函數(shù)φ(X,t)來計算，μ(ξ,t)和φ(X,t)具體的表達式可參考文獻[15]。

2 模型驗證

算例1尺寸為1 m×0.5 m的長方形薄板，板厚0.01 m，其彈性模量200 GPa，密度7 850 kg/m3，沿兩短邊作用p=200 MPa的均勻法向拉力，采用BPD和TBPD模型分別計算，物質(zhì)點間距Δ=0.01 m，鄰域半徑δ=3.015Δ，拉伸荷載通過轉變成體力荷載p/Δ，施加在兩短邊最外一層的真實物質(zhì)點上[14]，如圖3所示。

圖3 單軸拉伸荷載作用下的薄板近場動力學離散示意Fig.3 Geometry of a plate under uniaxial tension and its PD discretization

1)當泊松比ν=1/3時，如圖4所示，采用BPD模型或TBPD模型計算的位移結果都與解析解很好地吻合。采用TBPD方法計算的位移ux的最大相對誤差為3.25%，分布在薄板4個角點附近，這是由近場動力學的邊界效應導致的。

圖4 2種模型計算的位移ux的相對誤差Fig.4 Relative error on ux of two models

2)當泊松比ν∈[0,1/3]時，采用反算泊松比方法[8-9]和計算物質(zhì)點的位移ux的相對誤差來驗證BPD和TBPD模型的合理性，結果列于表2和表3。

表2 BPD模型與TBPD模型的泊松比ν反算結果 Table 2 Poisson′s ratio inverse calculation result of BPD and TBPD model

表3 BPD模型與TBPD模型計算的ux最大相對誤差Table 3 Maximum relative error on ux calculated by BPD model and TBPD model %

由表2可知，BPD模型無論材料泊松比怎么取值，反算泊松比總是0.332 8，當ν≠1/3時相對誤差大，這證實了它僅適合用來描述前面所說的特定泊松比材料。而采用TBPD模型時，反算泊松比值與真實泊松比非常接近，相對誤差為0.16%～0.70%。故它較適合用來描述泊松比在一定范圍內(nèi)的材料力學行為。此外，采用TBPD模型計算，在ν∈[0,1/3]時，位移ux的誤差全部小于5.0%，且隨著泊松比的減小，誤差逐漸減??；而采用BPD模型計算，在ν∈[0,0.20]時，位移ux的誤差全部大于5.0%，且泊松比越小，相對誤差越大。

圖5 剪切邊界條件幾何模型示意Fig.5 Geometric model of shear boundary conditions

如圖6所示，采用TBPD模型計算的位移結果與解析解高度吻合，位移ux和uy的最大相對誤差都為0.202%。

圖6 TBPD方法計算的位移ux與解析解間的相對誤差Fig.6 Relative error on ux for TBPD method

算例3邊長為50 mm的正方形薄板，中心有一長度2a=10 mm的水平裂紋，彈性模量192 GPa，密度8 000 kg/m3，泊松比ν=1/3，臨界斷裂能83 kJ/m2，在薄板上、下邊界作用大小20 m/s的均勻法向速度荷載。分別采用BPD和TBPD模型進行計算，物質(zhì)點間距Δ=0.25 mm，鄰域半徑δ=3.015Δ，計算時步取1.3367E-8 s，上、下兩邊的速度邊界條件是施加在薄板上下邊以外寬度為δ的3層虛擬物質(zhì)點上[14]。

采用BPD和TBPD模型計算時，薄板都是在裂紋左右兩端同時啟裂，然后各自向左右兩側沿水平方向擴展，最終貫通薄板。薄板在啟裂和貫通時刻稍有差異，但薄板破壞耗時基本相同，相關數(shù)據(jù)列于表4中。如圖7所示，取t=6.68 μs時刻的位移進行對比，兩模型計算的位移uy差異甚微。如圖8所示，取t=24.09 μs時刻的薄板破壞形態(tài)進行對比，兩模型計算的結果幾乎相同。

圖7 泊松比ν=1/3時2種模型在t=6.68 μs時刻的uy云圖Fig.7 Comparison of the uyat ν=1/3, t=6.68 μs of two models

圖8 2種模型在t=24.09 μs時刻的破壞形式對比Fig.8 Comparison of the failure modes at t=24.32 μs of two models

表4 薄板啟裂和貫通時刻Table 4 Initiating and through time of plate μs

算例4邊長200 mm的正方形薄板，其彈性模量為30 GPa，密度2 265 kg/m3，Ⅰ型能量釋放率GΙc=110 J/m2，Ⅱ型能量釋放率GIIc=10GIc，泊松比0.2。在試件中部的左右兩側各有一條長25 mm，寬5 mm的預制裂紋[16]， 如圖9所示。首先在薄板左側上半部分和右側下部分都水平施加合力為5 kN的均布力p，然后在薄板上下邊界施加v=10 mm/s的速度荷載。采用TBPD模型計算，物質(zhì)點間距Δ=4 mm，鄰域半徑δ=3.015Δ，計算步長1.0×10-7s，左右兩側的均布力p通過轉變成體力荷載p/Δ后，施加在薄板該部位的最外一層，寬度為Δ的真實物質(zhì)點上，速度荷載v作用在薄板以外寬度為δ的3層虛擬物質(zhì)點上。

圖9 混合破壞模式斷裂試樣的幾何形狀Fig.9 The geometry of the mixed mode fracture specimen

圖10(a)給出了采用TBPD模型計算得到的薄板的最終破壞形態(tài)，它與O?BOLT等采用基于微平面模型的有限元代碼模擬的結果[16]和試驗結果[17]基本相同，如圖10(b)、(c)所示。

圖10 混合破壞模式斷裂試樣的模擬結果Fig.10 Simulation results of the mixed mode fracture specimen

上述算例分別從薄板連續(xù)變形的定量計算和薄板破壞形態(tài)模擬2個方面都驗證了TBPD模型的合理性和有效性。

3 平行雙裂紋薄板的受拉破壞

如圖11所示，一邊長為50 mm的正方形薄板，中部預制有長2a=10 mm的2條水平裂紋，記為主裂紋①和主裂紋②，兩裂紋縱向間距為d。薄板的彈性模量為203 GPa，泊松比0.3，密度為7 850 kg/m3，在薄板上、下邊界作用大小20 m/s的均勻法向速度荷載。所用參數(shù)與趙金海等[18]相同，采用TBPD模型進行計算，物質(zhì)點間距Δ=0.5 mm，鄰域半徑δ=3.015Δ，計算步長取1.336 7×10-8s。按以下2種情形分別計算和討論：

圖11 平行雙裂紋薄板的裂紋布置Fig.11 Micro-crack distribution of plates with parallel double cracks

1)平行雙裂紋間距d取3 mm、5mm和7 mm時，對薄板破壞形態(tài)的影響；

2)平行雙裂紋間距d=3 mm時，考慮3種隨機分布的細小裂紋對薄板破壞形態(tài)的影響。

3.1 裂紋間距的影響

d取3 mm、5mm和7 mm時，裂紋擴展路徑如圖12(b)、(c)、(d)所示。薄板破壞過程都是在裂紋左右兩端同時啟裂，水平延伸一定距離后，裂紋①斜向上方、裂紋②斜向下方繼續(xù)擴展，兩裂紋縱向間距越來越大，直至貫通薄板。裂紋間距對薄板破壞模式影響不大，但對于裂紋水平延伸距離、擴展角度、裂紋貫通后縱向間距、起裂和貫通時間都有一定影響，如表5所示。隨著裂紋間距的增大，裂紋起裂時刻提前，水平延伸長度逐漸減小，但裂紋偏移角度、貫通間距和貫通時刻卻隨之增大。

圖12 TBPD模型計算的不同間距d下平行雙裂紋擴展Fig.12 Crack propagation paths of parallel double cracks calculated by TBPD model with longitudinal spacing d

表5 不同縱向間距d的平行雙裂紋擴展情況Table 5 Propagation of parallel double cracks with different longitudinal spacing d

值得注意的是，圖12(b)破壞模式與趙金海等基于傳統(tǒng)BPD模型所得結果[18]有很大差異，文獻[18]給出的破壞模式如圖13所示，其雙裂紋水平延伸一定距離后，都斜向薄板水平中心交匯成一條裂紋，然后水平延伸，直至貫通薄板。鑒于傳統(tǒng)BPD模型只能模擬ν=1/3的材料的力學行為[3]，所以將它用于模擬其他泊松比材料時，裂紋擴展路徑將可能出現(xiàn)失真或失效的情形。

圖13 d=3 mm時BPD模型計算泊松比0.3的薄板破壞Fig.13 Crack propagation path of plate with Poisson′s ratio 0.3 calculated by BPD model at d=3 mm

3.2 3種隨機細小裂紋的影響

在平行雙裂紋薄板內(nèi)，裂紋間距d=3 mm時考慮3種隨機細小裂紋的分布，細小裂紋長度取0～2 mm，細小裂紋布置方式和破壞形態(tài)見圖14。

如圖14(a)～(d)所示，當細小裂紋主要分布在某一主裂紋附近時，該主裂紋將沿著垂直于加載方向擴展直至貫通，而另一主裂紋幾乎沒有擴展。如圖14(e)、(f)所示，當細小裂紋分布在主裂紋的某一端附近，則主裂紋沿著分布有細小裂紋的這一端在荷載作用下擴展貫通，而該主裂紋的另一端只發(fā)生較小變化。對比圖12(b)可以看出，隨機細小裂紋的分布在一定程度上會抑制某一主裂紋在單側或兩側的擴展，導致薄板的破壞形態(tài)發(fā)生顯著的改變。

圖14 含隨機細小裂紋的平行雙裂紋板及其擴展路徑Fig.14 Parallel double cracks plate with small cracks and its′ propagation paths

4 結論

1)在不計剛體轉動和小變形條件下，以均勻應變工況，一次得到了模型的雙參數(shù)與材料彈性參數(shù)的轉換公式。且切向剛度系數(shù)的引入，使雙參數(shù)鍵的性能較單參數(shù)鍵更加完備，解決問題的能力更強。BPD模型只能模擬特定泊松比材料的力學行為，否則材料的反算泊松比、位移結果存在較大誤差，破壞模式也可能失真；TBPD模型突破了特定泊松比的限制，能夠模擬泊松比在一定范圍內(nèi)的材料力學行為，退化后的模型參數(shù)與傳統(tǒng)BPD的參數(shù)完全相同。

2) 以能量準則作為判斷單鍵斷裂的依據(jù)，詳細給出了TBPD模型單鍵斷裂條件，在特定泊松比情形下，它可退化為BPD模型的單鍵斷裂條件。

3)在一定條件下，平行雙裂紋的縱向間距對薄板的破壞過程產(chǎn)生一定的影響，隨著雙裂紋縱向間距的增大，裂紋起裂時刻提前，水平延伸長度逐漸減小，但裂紋偏移角度、貫通間距和貫通時刻卻隨之增大。隨機細小裂紋的分布在一定程度上會抑制某一主裂紋在單側或兩側的擴展，導致薄板的破壞形態(tài)發(fā)生顯著的改變。