劉 寧,王 帥,胡夢凡

(南京理工大學(xué) 機械工程學(xué)院,江蘇 南京 210094)

用傳統(tǒng)數(shù)值方法研究陶瓷裝甲、混凝土工事等抗彈材料侵徹毀傷問題時,還不能提供完整的沖擊破壞圖像對沖擊破壞過程給予動力學(xué)一致的分析,限制了結(jié)構(gòu)防護的性能優(yōu)化和評價。材料沖擊破壞是一個裂紋形成、擴展、分叉直到最終破碎的復(fù)雜過程,經(jīng)典連續(xù)介質(zhì)力學(xué)理論在處理這種高度非線性和不連續(xù)性問題時具有一定的局限性,主要難點在于經(jīng)典連續(xù)理論的數(shù)學(xué)模型是基于連續(xù)性假設(shè)下的包含空間導(dǎo)數(shù)的偏微分方程,一旦所研究的問題出現(xiàn)損傷、裂紋擴展、分層、斷裂和穿透等不連續(xù)問題時則不存在位移偏導(dǎo)數(shù),從而出現(xiàn)求解困難。

近年來一種新興的基于非局部思想的近場動力學(xué)方法(Peridynamics,PD)從根本上解決了傳統(tǒng)數(shù)值方法面臨的求解困難,該方法已在模擬裂紋萌生、擴展、分離等不連續(xù)問題時表現(xiàn)出了獨特的優(yōu)勢[1-2]。PD理論與經(jīng)典連續(xù)介質(zhì)力學(xué)最大的區(qū)別在于,PD方法使用空間積分方程來描述物質(zhì)力學(xué)行為,而經(jīng)典連續(xù)介質(zhì)力學(xué)是基于位移的導(dǎo)數(shù)建立運動方程。PD理論的這一特征表明,材料域內(nèi)任意位置的損傷的萌發(fā)和擴展是自發(fā)產(chǎn)生的,并且可沿任意路徑進(jìn)行,而不需要特別引入任何裂紋擴展準(zhǔn)則[3]。

近場動力學(xué)理論體系中鍵基模型(bond-based PD)得到了廣泛應(yīng)用[4-5],但常規(guī)鍵基模型只考慮了物質(zhì)點間的軸向力作用,僅有一個“材料常數(shù)”,導(dǎo)致了對材料泊松比的限制,如平面應(yīng)力問題材料泊松比為1/3,三維問題泊松比限制為1/4,制約了鍵基近場動力學(xué)模型的應(yīng)用范圍。本文基于有限元Euler-Bernoulli梁單元模型思想[6-9],提出一種新型Beam-based近場動力學(xué)本構(gòu)關(guān)系,突破常規(guī)鍵基近場動力學(xué)模型對材料泊松比的限制,通過算例驗證新型PD模型對沖擊動力學(xué)問題的適用性。

1 Beam-based PD模型

1.1 Bond-based PD基本理論

PD理論將物質(zhì)體系看作由眾多物質(zhì)點構(gòu)成,每個物質(zhì)點具有確定的體積和密度,空間域Ω內(nèi)的物質(zhì)點x與其周圍δ范圍內(nèi)的物質(zhì)點x′∈Ω:‖x′-x‖≤δ存在相互作用力f,如圖1所示,在任意時刻t,物質(zhì)點x的運動方程為

(1)

式中:ρ為物質(zhì)密度,u為物質(zhì)點x的位移,e為體力密度,H代表了物質(zhì)點x的近場鄰域范圍。

定義參考構(gòu)形中物質(zhì)點x與x′之間相對位置矢量ξ和相對位移矢量η分別為

ξ=x′-x,η=u′-u

(2)

圖1 PD理論模型

對于Bond-based PD模型,式(1)中物質(zhì)點間的相互作用力f可以寫為

(3)

式中:c為材料微彈性模量,可看作近場動力學(xué)理論中的“材料常數(shù)”;s為物質(zhì)點之間鍵相對伸長率:

(4)

通過PD理論中應(yīng)變能密度與經(jīng)典彈性力學(xué)中應(yīng)變能密度的等價關(guān)系[5],可推導(dǎo)出PD微彈性模量c,如平面應(yīng)力問題微彈性模量:

(5)

式中:E為彈性模量;ν為泊松比,且為1/3,三維問題的泊松比限制為1/4,固定的泊松比限制了Bond-based PD模型的應(yīng)用范圍和求解精度。

在近場動力學(xué)中,某物質(zhì)點的破壞程度用近場域內(nèi)與該粒子有相互作用的鍵的斷裂比例表征:

(6)

式中:μ為物質(zhì)點之間鍵是否斷裂的判據(jù),即

(7)

式中:s0為鍵的臨界伸長率。μ=1,表示物質(zhì)點之間鍵未發(fā)生斷裂;μ=0,表示鍵已斷裂。

1.2 Beam-based PD模型的構(gòu)建

基于有限元Euler-Bernoulli梁單元模型,在Bond-based PD模型的基礎(chǔ)上加入了物質(zhì)點間相對轉(zhuǎn)動效應(yīng),建立了新型Beam-based PD模型,物質(zhì)點間作用力由軸向力密度f1與切向力密度f2兩部分構(gòu)成,如圖2所示,(x,x′)為初始構(gòu)形,(y,y′)為現(xiàn)時構(gòu)形。

Beam-based PD本構(gòu)方程可以表示為

f=CS

(8)

式中:C為微彈性模量矩陣,可通過應(yīng)變能密度等效方法推導(dǎo);S為物質(zhì)點相對伸長率矢量。在小變形情況下,參考有限元Euler-Bernoulli梁單元模型,將物質(zhì)點相對伸長率矢量S進(jìn)行泰勒展開和線性化處理:

(9)

式中:u1,u2,u′1,u′2分別為物質(zhì)點變形前后軸向和切向位移矢量,r為物質(zhì)點間相對位置矢量。

圖2 梁模型物質(zhì)點運動

1.2.1 二維Beam-based PD模型

二維情況下Beam-based PD本構(gòu)函數(shù)可寫為

(10)

式中:c1,c2分別為軸向和切向微彈性模量。

由式(9)可得物質(zhì)點間軸向力與切向力表達(dá)式:

(11)

在二維笛卡爾坐標(biāo)系下,物質(zhì)點間軸向相對位移與切向相對位移可通過物質(zhì)點坐標(biāo)相對位移分量表示,如圖3所示。

圖3 二維梁模型物質(zhì)點位移分量示意圖

(12)

將式(12)帶入式(11)可得:

(13)

物質(zhì)點間“梁單元”具有的微勢能為

(14)

物質(zhì)點x的應(yīng)變能密度為

(15)

同時,在彈性力學(xué)中,平面應(yīng)力問題的應(yīng)變能密度為

(16)

根據(jù)應(yīng)變能密度等效法,由式(15)、式(16)可得Beam-based PD模型軸向微彈性模量c1與切向微彈性模量c2:

(17)

寫成矩陣形式:

(18)

式(17)中c1與c2為2個相互獨立的PD常數(shù),其中c1>0,c2≥0。可見對于平面應(yīng)力問題,Beam-based PD模型的泊松比取值范圍為(-1,1/3],大大擴展了傳統(tǒng)鍵基PD模型的應(yīng)用范圍。

同理,可推導(dǎo)平面應(yīng)變問題中Beam-based PD模型微彈性模量矩陣:

(19)

式中的泊松比范圍為(-1,1/4]。

1.2.2 三維Beam-based PD模型

三維情況下Beam-based PD本構(gòu)方程可寫為

(20)

寫成分量形式:

(21)

三維條件下物質(zhì)點間“梁單元”微勢能為

(22)

積分可得三維PD應(yīng)變能密度:

(23)

而經(jīng)典彈性力學(xué)中三維應(yīng)變能密度為

(24)

通過應(yīng)變能密度等效方法,可得三維Beam-based PD模型微彈性模量:

(25)

寫成矩陣形式:

(26)

式中:c1>0,c2≥0,三維Beam-based PD模型泊松比范圍為(-1,1/4]。

2 數(shù)值方法

2.1 離散方程

將計算域進(jìn)行離散,運動方程式(1)的離散形式可表示為

(27)

式中:i表示當(dāng)前物質(zhì)點,n為時間迭代步數(shù),Vj為當(dāng)前物質(zhì)點i的鄰域內(nèi)物質(zhì)點j的體積。

2.2 PD沖擊接觸算法

將沖擊彈丸視為剛體,在沖擊過程中彈丸不斷把被沖擊材料物質(zhì)點從其運動軌跡上排擠開,如圖4所示,受排擠物質(zhì)點速度為

(28)

(29)

圖4 剛性體沖擊接觸模型

3 算例研究

3.1 平板沖擊

圖5為剛性圓盤沖擊矩形板模型[10],圓盤質(zhì)量m=1.57 kg,直徑D=50 mm,以速度v0=32 m/s沖擊無約束矩形板上邊緣中間位置,矩形板長l=200 mm,寬b=100 mm,厚度d=9 mm,密度ρ=8 000 kg/m3,彈性模量E=191 GPa,泊松比ν=1/3。采用Beam-based PD模型計算圓盤沖擊過程,離散物質(zhì)點間距Δx=1 mm,近場鄰域范圍δ=3.015dx,時間步長Δt=0.1 μs。

圖6為不同時刻矩形板內(nèi)部應(yīng)力波云圖,由圖可見,在圓盤沖擊矩形板后,應(yīng)力波以沖擊點為圓心,以半圓形輻射狀向周圍傳播,應(yīng)力波傳播到左、右邊界和底部邊界時發(fā)生反射,多次反射后形成了復(fù)雜的應(yīng)力分布特征,但在整個沖擊過程中矩形板僅發(fā)生彈性形變,沒有發(fā)生破壞。

圖5 圓盤沖擊矩形板模型

圖6 矩形平板應(yīng)力波傳播云圖

為了驗證Beam-based PD模型對不同泊松比材料的適用性,進(jìn)一步計算了泊松比ν=1/4時矩形板中心的沖擊響應(yīng),作為對比同時給出了有限元計算結(jié)果,如圖7所示。

圖7 泊松比為0.25時中心點位移變化曲線

由圖7可見,本文模型計算結(jié)果與有限元結(jié)果吻合較好,驗證了本文模型的正確性,新型PD模型有效拓展了Bond-based模型的應(yīng)用范圍。

3.2 三維Kalthoff-Winkler沖擊模型

Kalthoff-Winkler沖擊試驗是用剛性圓柱彈丸沖擊帶有對稱預(yù)置裂紋的馬氏體時效鋼靶板,觀測裂紋的擴展過程,包括裂紋擴展方向和擴展速度等[11]。試驗結(jié)果表明,圓柱體撞擊靶板時產(chǎn)生沖擊壓縮波,使得裂紋從裂尖處開始擴展,產(chǎn)生Ⅱ型裂紋,最終形成的裂紋與預(yù)制裂紋夾角約為70°。Kalthoff-Winkler沖擊試驗已成為研究材料動態(tài)斷裂本構(gòu)模型和數(shù)值方法的經(jīng)典驗證算例。Kalthoff-Winkler沖擊模型如圖8所示,材料參數(shù)如表1所示。圖中,圓柱直徑D=50 mm,沖擊速度v0=32 m/s,靶板長度l=200 mm,厚度d=9 mm,寬度b=100 mm,預(yù)制裂紋長l0=50 mm,間距h=50 mm,裂紋寬度a=1.5 mm。

圖8 Kalthoff-Winkler試驗幾何模型

表1 靶板材料參數(shù)

E/GPaρ/(kg·m-3)ν1907 8300.3

采用Beam-based PD模型模擬Kalthoff-Winkler沖擊試驗,取物質(zhì)點間距Δx=1 mm,近場鄰域δ=3.015dx,材料臨界伸長率s0=0.01,時間步長Δt=0.58 μs。

圖9為Beam-based PD模擬的裂紋擴展過程,裂紋從預(yù)制裂紋尖端開始擴展,幾乎沿一條直線對稱擴展至靶板左、右兩側(cè)邊緣,裂紋擴展方向與預(yù)制裂紋夾角大約為68.3°。圖10為文獻(xiàn)[12]給出的Kalthoff-Winkler沖擊試驗結(jié)果和擴展有限元仿真結(jié)果。由圖可見,本文模擬的裂紋擴展過程與試驗值吻合較好。

眾所周知,傳統(tǒng)有限元法難以直接模擬材料破壞過程,而改進(jìn)的擴展有限元法需要人為設(shè)置合適的破壞準(zhǔn)則才能模擬裂紋擴展。本文建立的Beam-based PD模型不僅擴展了傳統(tǒng)鍵基模型的適用范圍,還能夠模擬裂紋的自然萌生和自發(fā)擴展過程,無需任何外部干預(yù),在處理材料沖擊破壞問題時顯示出巨大的優(yōu)勢,為沖擊動力學(xué)研究提供了新的研究手段。

圖9 Beam-based PD模擬的裂紋擴展過程

圖10 Kalthoff-Winkler試驗及擴展有限元結(jié)果[12]

4 結(jié)束語

本文針對常規(guī)Bond-based近場動力學(xué)模型對材料泊松比的限制,借鑒有限元Euler-Bernoulli梁單元模型,在Bond-based模型的基礎(chǔ)上加入了物質(zhì)點間相對轉(zhuǎn)動效應(yīng),建立了新型Beam-based PD本構(gòu)模型,有效拓展了泊松比范圍。圓盤沖擊矩形板模擬結(jié)果表明,泊松比對動態(tài)沖擊響應(yīng)具有重要影響,本文模型能夠正確反映材料泊松比效應(yīng)。模擬了Kalthoff-Winkler靶板沖擊過程,裂紋擴展角度和裂紋擴展速度與試驗結(jié)果吻合較好,新型Beam-based PD模型能以材料真實泊松比準(zhǔn)確計算沖擊載荷作用下材料的變形和破壞過程。本文模型對認(rèn)識脆性材料沖擊破壞的復(fù)雜物理機制具有明顯優(yōu)勢,為沖擊動力學(xué)問題研究提供了一條新的技術(shù)途徑。