完全拉格朗日SPH 在沖擊問題中的改進(jìn)和應(yīng)用1)

王璐 徐緋 楊揚(yáng)

*(長安大學(xué)建筑工程學(xué)院,西安 710061)

?(西北工業(yè)大學(xué)航空學(xué)院,西安 710072)

引言

光滑粒子流體動力學(xué)(smoothed particle hydrodynamics,SPH)模擬拉伸狀態(tài)的粒子時,由控制方程得知粒子間會存在吸引力,有可能導(dǎo)致粒子聚集從而產(chǎn)生拉伸不穩(wěn)定現(xiàn)象[1-2].Swegle等[3]指出當(dāng)粒子的應(yīng)力狀態(tài)和核函數(shù)不匹配,即核函數(shù)的二階導(dǎo)數(shù)與應(yīng)力乘積大于零時,SPH 計算會出現(xiàn)不穩(wěn)定現(xiàn)象.由此,不少學(xué)者在消除拉伸不穩(wěn)定方面做了嘗試:理論層面有修正核函數(shù)插值方式[4]和重構(gòu)SPH 核函數(shù)[5].數(shù)值層面有密度均勻化[6]、施加人工黏性[7]、增加人工應(yīng)力項[8-9]等.Belytschko等[10]認(rèn)為SPH 產(chǎn)生拉伸不穩(wěn)定性的根源在于SPH 采用的是歐拉核函數(shù),該種核函數(shù)會扭曲材料的應(yīng)力空間.此外Belytschko等[11]通過穩(wěn)定性分析發(fā)現(xiàn)拉格朗日核函數(shù)下SPH穩(wěn)定性與應(yīng)力狀態(tài)無關(guān),并在此基礎(chǔ)上提出了TLSPH(total Lagrangian SPH)方法.但由于固體邊界支持域的缺失,TL-SPH 仍會面臨精度不足的缺陷[12],由此本文擬在完全拉格朗日框架下推導(dǎo)高階TLSPH 方程以提高整個計算域的計算精度.

SPH 粒子法的特性使得計算模型的邊界條件難以像網(wǎng)格法一樣嚴(yán)格實施,而SPH 界面處理對于沖擊問題尤為關(guān)鍵.目前對SPH 界面處理技術(shù)的研究有不少成果:理論層面上,DSPH(discontinuous SPH)公式最初由Liu等[13]在處理非連續(xù)沖擊波問題時提出,隨后Xu等[14]將DSPH 擴(kuò)展到了多維情形,并進(jìn)一步提出高階形式的DFPM(discontinuous finite particle method)方法[15],DFPM 可以改善傳統(tǒng)SPH 方法在非連續(xù)區(qū)域上的計算精度,若將DFPM移植到完全拉格朗日框架,可提高完全拉格朗日SPH 方法對界面的模擬精度.數(shù)值層面上,由于TLSPH 中粒子的運(yùn)動狀態(tài)只由初始構(gòu)型下鄰域粒子的應(yīng)力狀態(tài)決定,所以需考慮沖擊過程中相互作用的傳遞.目前SPH 方法中,虛粒子法和邊界力法最為常用[16-19],但虛粒子法在復(fù)雜邊界處難以實現(xiàn),邊界力法作用范圍和強(qiáng)度難以控制容易導(dǎo)致數(shù)值振蕩或穿透問題.為了避免上述弊端,本文擬在虛粒子法的基礎(chǔ)上提出一種基于粒子間黎曼模型[20-21]的接觸算法.

SPH 模擬固體的損傷破壞具有其獨特的優(yōu)勢,不像網(wǎng)格法中需要定義損傷單元的刪除準(zhǔn)則,甚至需對網(wǎng)格進(jìn)行重新劃分[22-23].SPH 方法在模擬損傷破壞時,一方面可直接通過材料的本構(gòu)和損傷模型體現(xiàn)粒子自身的損傷程度[24-25].另一方面可通過修正核函數(shù)或其作用范圍實現(xiàn)裂紋的模擬[26].TLSPH 方法中,粒子的鄰域粒子在物理演變過程中保持不變,始終是初始時刻下的鄰域粒子,對破碎、斷裂等間斷現(xiàn)象的模擬能力有限.因此,需根據(jù)當(dāng)前時刻粒子的損傷破壞分布,自適應(yīng)地更新拉格朗日核函數(shù).本文將在完全拉格朗日框架下提出一種損傷斷裂模型以準(zhǔn)確模擬固體沖擊問題中出現(xiàn)的斷裂現(xiàn)象.

綜上,本文擬在TL-SPH 方法基礎(chǔ)上,從高階格式推導(dǎo)、界面處理以及損傷模型建立三個方面入手,實現(xiàn)對沖擊問題的精確模擬.

1 TL-SPH法

完全拉格朗日框架下SPH 方法離散的固體動力學(xué)控制方程為

材料變形時粒子的柯西應(yīng)力σ包括各項同性壓力p和偏應(yīng)力S,壓力p通過式(1)中的壓力-密度狀態(tài)方程得到.沖擊過程中材料不可避免地會進(jìn)入塑性屈服階段,粒子進(jìn)入塑性階段后,本文用Von Mises 準(zhǔn)則[14]修正材料偏應(yīng)力.

得到柯西應(yīng)力σ后,粒子i的第一Piola-Kirchhoff應(yīng)力Pi可由柯西應(yīng)力σi經(jīng)過變形梯度及其行列式轉(zhuǎn)化得到[27]

2 TL-SFPM法

SPH 核近似精度受粒子分布和支持域完整性影響,固體邊界附近粒子的核估計精度會降低[5],下面在完全拉格朗日框架下推導(dǎo)SPH 高階插值公式.借鑒基于當(dāng)前構(gòu)型的高階插值公式的思想,以二維情形為例,將f(X)在初始構(gòu)型下的點Xi處泰勒展開,省去二階及以上導(dǎo)數(shù)項,分別在等式兩邊乘以Wi,X和Wi,Y并在Ω區(qū)域內(nèi)積分,可得到[5]

式(5)可進(jìn)一步整理為式(6),可以看出高階插值公式本質(zhì)是利用僅與粒子位置有關(guān)的矩陣Bi對核函數(shù)梯度做以修正.與其他高階方法類似,式(6)會面臨計算繁雜的問題.為了節(jié)省計算量,可簡化參與求和粒子的選取方式.

若借鑒SFPM 方法[15]中的粒子選取方式,并充分考慮粒子間相互作用的一致性,可對稱地選取距離粒子i最近的四個粒子.若粒子i位于邊界,此時粒子選取方式與SFPM 方法一致.如圖1 所示為規(guī)則幾何構(gòu)型下不同區(qū)域粒子的選取方式.

圖1 規(guī)則構(gòu)型下不同區(qū)域粒子選取方式Fig.1 Particles selection modes in different regions under regular configuration

下面給出不同粒子選取模式下函數(shù)梯度計算公式.假設(shè)初始時刻粒子均勻分布且粒子間距為d,光滑長度為1.2d,核函數(shù)取為常用的Wendland 核函數(shù).此時內(nèi)部區(qū)域、邊界區(qū)域和角點區(qū)域的粒子對應(yīng)的函數(shù)梯度計算結(jié)果如下

該種基于初始構(gòu)型的選取四個粒子的函數(shù)梯度計算式所對應(yīng)的粒子近似精度、一致性已在文獻(xiàn)[28]中驗證.

為了說明式(7)～式(9)的收斂性,以估計函數(shù)f(x,y)=x3+y3+2xy的梯度為例,計算域為[1,2]×[1,2].圖2為不同粒子數(shù)(M)下,式(7)～式(9)計算函數(shù)?f/?X的均方誤差MSE,為了清楚顯示,坐標(biāo)取為對數(shù)值.圖中顯示橫坐標(biāo)和縱坐標(biāo)呈現(xiàn)反比例關(guān)系,當(dāng)粒子間距減小,計算誤差呈指數(shù)減小,說明函數(shù)梯度計算式(7)～式(9)具有收斂性.

圖2 不同粒子間距下式(7)～式(9)計算三次函數(shù)?f/?X 的均方誤差Fig.2 Mean square error of cubic function ?f/?X calculated by Eqs.(7)～(9) under different particle spacings

當(dāng)計算域為不規(guī)則幾何構(gòu)型時,首先,假設(shè)計算域內(nèi)部初始時刻粒子均勻分布,此時,粒子的選取方式和圖1 中內(nèi)部粒子一致.當(dāng)粒子位于不規(guī)則的邊界時,如圖3 所示,則按照圖1 邊界粒子選取方式:首先選擇粒子本身,其余粒子選為距離粒子i最近的三個粒子.此時,可根據(jù)被選擇粒子與粒子i的相對位置計算式(6)即可得到函數(shù)梯度計算式.

圖3 不規(guī)則計算域邊界粒子選取方式Fig.3 Particles selection mode in boundary region under irregular configuration

為了保證計算的穩(wěn)定性,就要保證距離粒子i最近的三個粒子不在同一條直線上.所以初始時刻離散的粒子模型應(yīng)盡可能均勻分布.對于特別復(fù)雜的模型,可采用在常壓p0的作用下對粒子均勻化[29].在此基礎(chǔ)上,仍需檢測所選粒子是否滿足穩(wěn)定性要求,如果粒子在一條直線上,那么將選擇不在同一條直線的最近粒子.

本文是在二維情況下推導(dǎo)TL-SFPM 方程,高階插值公式(6)可通過增加維度擴(kuò)展至三維,這構(gòu)成了TL-SFPM 方法應(yīng)用于三維問題的理論基礎(chǔ).借鑒SFPM 方法[15]以及充分考慮粒子間相互作用的一致性,三維情況下,當(dāng)粒子在空間均勻分布時,可對稱地選取以粒子i為中心的六面體上距離最近的六個粒子,與二維類似,此時粒子i與被選取粒子的間距均為d.當(dāng)粒子位于邊界或者角點時,首先選擇粒子本身,其余粒子選為距離粒子i最近的五個粒子.上述三維模擬策略仍需進(jìn)一步驗證.

下面介紹TL-SFPM 方法中不同物理量梯度的求解式.

2.1 變形梯度求解

作為初始構(gòu)型和當(dāng)前構(gòu)型之間的橋梁,粒子的變形梯度Fi是關(guān)鍵物理量.使用式(7)～式(9)可分別得到內(nèi)部區(qū)域、邊界區(qū)域和角點區(qū)域粒子變形梯度Fi的計算式為

式中,u,v分別為粒子x方向和y方向的位移,下標(biāo)數(shù)字(1,2,3,4)與圖1 的粒子編號對應(yīng),I為單位矩陣.

2.2 加速度求解

若直接使用式(6) 求解加速度方程中的第一Piola-Kirchhoff 應(yīng)力P的梯度,此時會導(dǎo)致粒子i與j之間相互作用不一致,破壞守恒關(guān)系.為了維持守恒關(guān)系,下面推導(dǎo)加速度方程.

在不考慮人工黏性條件下,TL-SPH 求解加速度方程為[27]

若使用高階插值公式修正核函數(shù)梯度,加速度方程可改為[30]

可以看出,該式可保證mi(dvi/dt)j→i=-mj(dvj/dt)i→j,粒子i對粒子j的作用力與粒子j對粒子i的作用力一致.若考慮人工黏性,將人工黏性直接加到加速度方程即可

上式并未修正人工黏性項,這是由于該項所占比重較小,為了節(jié)省計算時間直接使用傳統(tǒng)核函數(shù)計算.為了與變形梯度計算保持一致,加速度計算仍采用圖1 所示的粒子選取方式,此時,加速度方程寫為

3 界面處理

3.1 TL-DFPM法

DFPM 方法是一種考慮界面不連續(xù)的高階離散方法,它基于界面處物理量的不連續(xù)性,在非連續(xù)區(qū)域兩邊分別進(jìn)行連續(xù)性修正,且可避免不連續(xù)位置的確定[15].為了提高TL-SPH 方法在界面處的精度,本節(jié)擬將DFPM 方法引入完全拉格朗日框架中.

如圖4 所示為二維情形下的不連續(xù)計算域示意圖,Ω為粒子的整個支持域,Ω1和Ω2為其子域,Ω1和Ω2之間不連續(xù)界面為Γ(Ω1∪Ω2∪Γ=Ω).

圖4 二維情形下的不連續(xù)計算域示意圖Fig.4 Discontinuous computational domain in 2-D

二維情形下DFPM 方法的基本方程[15]為

式中,Vj為鄰域粒子j的體積,f,x(xi)=?f(xi)/?x,f,y(xi)=?f(xi)/?y,Wij,x=?Wij/?x,Wij,y=?Wij/?y,x,y為當(dāng)前構(gòu)型下的坐標(biāo)分量.上式只考慮了鄰域內(nèi)與粒子i同側(cè)的粒子信息.通過非連續(xù)性近似分析,該修正式可以保證梯度計算的精確性,也能夠?qū)崿F(xiàn)不連續(xù)界面處粒子信息的準(zhǔn)確估計.若將上述方法引入完全拉格朗日框架中,可得

可以看出上式將式(6)方程中的所有鄰域粒子求和修正為子域Ω1內(nèi)的粒子求和,這可以避免異側(cè)子域Ω2粒子物理信息間斷帶來的誤差.基于此,界面處TL-SFPM 方法中四個粒子的選取方式可修正為圖5 所示方式,界面附近粒子被當(dāng)作邊界粒子處理.

圖5 界面處TL-SFPM 方法中四個粒子的選取方式Fig.5 Particles selection mode at interface in TL-SFPM method

3.2 接觸模型——黎曼接觸法

為了保證界面處的粒子計算精度,式(19)將子域Ω1粒子的計算與子域Ω2區(qū)域隔離開,這時需引入不同區(qū)域間的接觸模型,定義粒子i和j間的相互作用.本節(jié)將提出基于黎曼解的接觸算法.

如圖6 所示為SPH 粒子間黎曼模型示意圖,假設(shè)粒子間存在間斷面,將其相互作用按照連續(xù)介質(zhì)中的間斷問題處理,間斷界面位于中心點(ri+rj)/2處,此時,間斷面左側(cè)密度、速度、壓力和聲速為

圖6 SPH 粒子間黎曼模型Fig.6 Riemann model with SPH particles

間斷面右側(cè)密度、速度、壓力和聲速為

式中,eij=為粒子i指向粒子j的單位向量.若通過求解該黎曼問題得到中間應(yīng)力p*,將其代入粒子i動量方程中,即可實現(xiàn)異側(cè)粒子j對粒子i的作用力fij,fij的具體表達(dá)式為式(22).很明顯fij=-fji,滿足相互作用的一致性

式中,p*為壓力黎曼近似解,一種近似的線性解為[31-32]

3.2.1 流固接觸

SPH 在利用虛粒子方法施加流體的邊界條件時,需在初始階段提前布置虛粒子,其位置固定或者由流體粒子鏡像決定.這種策略往往比較費(fèi)時且對復(fù)雜邊界受限[33].然而該種方法能夠保證守恒性,且滑移與無滑移邊界條件易實現(xiàn)[34],為了保留該方法的優(yōu)勢,本節(jié)將修正虛粒子方法提出動虛粒子黎曼法.

動虛粒子黎曼法的具體措施:將固體粒子i看作提前布置的虛粒子,然而該種虛粒子位置并不固定,根據(jù)自身的運(yùn)動方程更新其位置.可以賦予該種動虛粒子新的物理量——虛壓力,虛壓力通過流體壓力插值得到[34],即

式中,f為固體粒子i鄰域內(nèi)的流體粒子,Wif為核函數(shù),ai為固體粒子自身的加速度,g為重力加速度.得到虛壓力后,流固間相互作用可通過式(22)實現(xiàn),此時,pL為固體粒子虛壓力,pR為流體粒子壓力pf,ρL和cL為固體粒子密度ρi、聲速ci,ρR和cR為流體粒子密度ρf、聲速cf.vL為動虛粒子即固體粒子速度vi,vR為流體速度vf.此時,固體粒子i所受到來自流體的外力fif為

式中,pif*為固體粒子i和流體粒子f之間的近似壓力黎曼解.

3.2.2 固固接觸

該壓力不同于式(1)狀態(tài)方程求解的壓力,僅用來求解固體間接觸力.當(dāng)不同固體材料相互作用時,黎曼問題的左狀態(tài)物理信息為固體粒子i的密度pi、速度vi、壓力和聲速ci,右狀態(tài)物理信息為固體粒子j的密度pj、速度vj、壓力和聲速cj.通過式(22)即可得到粒子i和j間的相互作用力fij.

4 粒子損傷破壞模型

粒子損傷累積到一定程度后會發(fā)生破壞,定義損傷因子D描述損傷程度[26],D的范圍為0～ 1,D=0時,粒子未發(fā)生破壞,D=1 時,粒子完全破壞,0＜D＜ 1 時,粒子發(fā)生部分破壞,表明與周圍粒子的相互作用減弱.不同材料根據(jù)材料特性有不同形式的損傷因子式.SPH 方法是利用核函數(shù)衡量粒子間的相互作用,因此可通過下列修正體現(xiàn)粒子損傷后與鄰域粒子相互作用的減弱

5 算例驗證

5.1 流固沖擊——帶彈性擋板潰壩

這節(jié)以流固耦合沖擊問題為例驗證本文提出的TL-SFPM 方法和界面處理技術(shù)的優(yōu)勢.如圖7 所示,紅色區(qū)域為彈性擋板,上端與剛性固壁固連,擋板內(nèi)藍(lán)色區(qū)域為水,在重力作用下,水坍塌沖擊彈性擋板,擋板發(fā)生變形.擋板的材料參數(shù)和模型幾何尺寸見表1.水初始密度為1000kg/m3,水近似為弱可壓模型[35-36],用低耗散GSPH 方法計算.低耗散GSPH 方法離散流體的控制方程如下[37]

圖7 帶彈性擋板潰壩模型Fig.7 Model of the dam with elastic baffle

表1 彈性擋板材料參數(shù)Table 1 Material parameters of elastic baffle

式中v*和p*是速度和壓力的近似黎曼解[37].低耗散GSPH 中,為了降低耗散導(dǎo)致流動失真,將p*中的聲速限制為[38]

式中,假設(shè)流體粒子i和j的聲速保持一致為c0,ω為常數(shù),該算例中取1.0.該算例流固接觸模型使用第3 節(jié)的黎曼接觸算法.

圖8為不同粒子間距下,水采用低耗散GSPH方法,彈性體采用TL-SFPM 方法的策略下,彈性擋板自由端豎直位移變化曲線,可以看出粒子間距減小到一定程度時,位移變化趨于收斂.為了保證一定的計算精度和效率,本算例選取粒子間距為0.5 mm.

圖8 不同粒子間距彈性擋板自由端豎直位移變化曲線Fig.8 Vertical displacement of the free end of elastic baffle with different particle spacings

為了說明TL-SFPM 方法相比傳統(tǒng)TL-SPH 方法的優(yōu)勢,圖9(a)比較了相同分辨率下,TL-SFPM方法和傳統(tǒng)TL-SPH 方法模擬的彈性擋板自由端水平和豎直位移變化曲線,可以看出TL-SFPM 的結(jié)果與試驗[39]的結(jié)果更吻合.對比二者模擬結(jié)果,在位移增大的過程中,TL-SFPM 的結(jié)果較大,在達(dá)到最大位移后,反而位移較小,這是由于TL-SFPM 方法修正了邊界區(qū)域的一致性,在計算變形梯度或壓力梯度時更精確.圖9(b)也對比了擋板附近水面高度隨時間的變化曲線,TL-SFPM 仍然與試驗更吻合.圖10給出不同時刻模擬得到的粒子分布與試驗形態(tài)對比圖,結(jié)構(gòu)變形和流體坍塌形狀均與試驗保持一致.但是試驗中可以觀測到有部分水從擋板兩端噴濺出,而模擬的二維模型并未出現(xiàn)該現(xiàn)象.這是由于:(1) 試驗是三維模型,水可從擋板兩端逃逸;(2)為了計算穩(wěn)定,低耗散GSPH 方法會稍微增大水的耗散,導(dǎo)致水飛濺的現(xiàn)象并不明顯.

圖9 彈性擋板和液面位置隨時間變化曲線Fig.9 Curves of the elastic baffle and the water surface positions

圖10 TL-SFPM 模擬結(jié)果和試驗結(jié)果[39]對比圖(自0～0.32 s,每隔0.08 s)Fig.10 Comparison of TL-SFPM simulation and experiment results[39]from0s to0.32 s(at0.08 s intervals)

圖11 對比了某一時刻TL-SFPM 方法和文獻(xiàn)[39]得到的彈性板應(yīng)力分布云圖,可以看出本文得到的應(yīng)力分布更加穩(wěn)定均勻,文獻(xiàn)中的結(jié)果仍存在應(yīng)力震蕩現(xiàn)象.其原因可能是:一方面,文獻(xiàn)[39]中采用人工應(yīng)力法改善拉伸不穩(wěn)定,而本文使用的是TLSFPM 方法.另一方面,本文采用的是基于黎曼解的接觸算法,文獻(xiàn)[39]中使用的是傳統(tǒng)虛粒子算法.

圖11 TL-SFPM和文獻(xiàn)[39] 模擬得到某時刻彈性擋板y 方向應(yīng)力分布云圖Fig.11 Stress distribution at y-direction of the baffle obtained by the TL-SFPM and Ref.[39]

表2 比較了TL-SFPM 方法和TL-SPH 方法在Intel(R) Core(TM) i7-3770八核心處理器下的計算時間.可以看出TL-SFPM 相比TL-SPH 可節(jié)省一定的計算時間,這是因為TL-SFPM 中鄰域內(nèi)參與計算的粒子數(shù)相比TL-SPH 更少.雖然TL-SFPM 方法采用的是高階格式,但是由于初始構(gòu)型下粒子相對位置確定,高階格式的修正矩陣Bi只需在初始時刻計算一次,后續(xù)直接調(diào)用即可.二者計算時間相差并不大是因為本算例中彈性體粒子數(shù)所占總粒子數(shù)比重較少,而且兩種方法均需基于當(dāng)前構(gòu)型計算流固接觸力.

表2 TL-SPH和TL-SFPM 模擬時間對比Table 2 Computational time with TL-SPH and TL-SFPM methods

5.2 固固沖擊——子彈撞擊靶板

這節(jié)以子彈撞擊靶板算例驗證本文提出的TLSFPM 方法、界面處理技術(shù)以及損傷破壞模型的合理性和精確性.圓頭的鋼彈撞擊鋁合金的飛機(jī)蒙皮[40],彈丸直徑12.7 mm,長度51 mm,靶板長和寬均為300mm,厚3 mm,子彈以150m/s 速度向右沖擊鋁板,子彈和靶板均采用Johnson-Cook 本構(gòu)模型,具體材料參數(shù)見表3,表中參數(shù)含義與Johnson-Cook 本構(gòu)模型一致.

表3 鋼和鋁材料參數(shù)Table 3 Material parameters of steel and aluminum

5.2.1 無損傷破壞模型

首先在該模型中不加入材料的損傷破壞模型以驗證TL-SFPM 方法改善拉伸不穩(wěn)定的能力.可以分別用傳統(tǒng)高階SPH 方法——CSPM(corrective smoothed particle method)[3]和TL-SFPM 計算該模型.圖12給出某時刻鋼彈以150m/s 的速度撞擊鋁板后子彈和鋁板的速度云圖及局部放大圖,從圖可以看出:(1) CSPM 結(jié)果中,在子彈撞擊鋁板后,鋁板出現(xiàn)了一定程度的破壞,并且在鋁板的變形兩端也出現(xiàn)了不同程度的斷裂,出現(xiàn)了粒子聚集現(xiàn)象,這體現(xiàn)了傳統(tǒng)SPH 拉伸不穩(wěn)定的缺陷;(2)同一時刻TL-SFPM方法模擬的結(jié)果并未出現(xiàn)斷裂現(xiàn)象,粒子分布及速度分布整齊且均勻,改善了拉伸不穩(wěn)定缺陷.圖13為TL-SFPM和有限元(FEM)模擬的不同時刻子彈和鋁板的速度云圖,看出兩種方法得到的結(jié)果吻合,速度分布均勻,這驗證了TL-SFPM 方法具有較好的精度.

圖12 CSPM和TL-SFPM 得到鋼彈撞擊鋁合金板同一時刻速度分布云圖Fig.12 Velocity distribution obtained by the traditional SPH and high-order TL-SPH

圖13 TL-SFPM和有限元模擬的不同時刻子彈和板的速度(m/s)云圖Fig.13 Velocity(m/s) distribution at different times obtained by TL-SFPM and FEM

5.2.2 有損傷破壞模型

若在計算中加入以下?lián)p傷因子

即當(dāng)粒子的塑性應(yīng)變εpl累積到一定值時,粒子發(fā)生破壞,該算例εpl,max取為0.12.圖14 給出TLSFPM 計算得到的不同速度下靶板的破壞形態(tài),綠色為破壞粒子.從圖中看出不同速度下鋁板的破壞模式不同.速度為150m/s 時,鋁板在撞擊部位出現(xiàn)張開的裂紋,隨著鋼彈向右運(yùn)動,裂紋不斷張開直到鋼彈穿過板.當(dāng)鋼彈速度為300m/s 時,首先鋁板的內(nèi)部在沖擊部位兩側(cè)發(fā)生破壞,接著裂紋不斷擴(kuò)展,形成了碎片.當(dāng)鋼彈速度增大到500m/s 時,鋁板首先在與鋼彈接觸部位發(fā)生破壞,最終形成了更小的碎片.

圖14 加入斷裂模型后不同沖擊速度下的失效粒子分布Fig.14 Distribution of failure particles at different impact velocities with the fracture model

圖15為鋼彈剩余速度,鋁板的變形直徑和彈坑深度隨沖擊速度的對比曲線,CSPM 由于會出現(xiàn)數(shù)值斷裂導(dǎo)致和試驗[40]相差較大,而TL-SFPM 與試驗吻合較好,表明只有改善數(shù)值斷裂現(xiàn)象,才能模擬出精確的結(jié)果.

圖15 鋼彈沖擊靶板計算結(jié)果與試驗結(jié)果對比[40]Fig.15 Comparison of the steel projectile impacting the aluminum alloy plate results between simulation and experiment[40]

表4 比較了TL-SFPM 方法和CSPM 方法在Intel(R) Core(TM) i7-3770八核心處理器下的計算時間.可以看出TL-SFPM 方法相比CSPM 方法計算時間略長,雖然TL-SFPM 方法可節(jié)省粒子鄰域搜索的時間,但是當(dāng)粒子進(jìn)入塑性后,TL-SFPM 方法的計算過程相比CSPM 更復(fù)雜,需要將柯西應(yīng)力σ轉(zhuǎn)化為第一Piola-Kirchhoff 應(yīng)力P.另外,TL-SFPM 方法斷裂模型中需要及時對核函數(shù)更新.當(dāng)不同材料作用,TL-SFPM 也需在每一時間步內(nèi)計算當(dāng)前構(gòu)型下的壓力黎曼近似解p*以求解不同材料間的接觸力.

表4 CSPM和TL-SFPM 模擬時間對比Table 4 Computational time with CSPM and TL-SFPM methods

6 結(jié)論

本文針對傳統(tǒng)SPH 模擬固體時出現(xiàn)的拉伸不穩(wěn)定缺陷提出了高階格式的完全拉格朗日SPH 方法——TL-SFPM 方法,并分別推導(dǎo)了變形梯度和加速度方程的離散格式.此外,在完全拉格朗日框架下,分別從理論和數(shù)值方面提出了TL-DFPM和黎曼接觸法,形成了雙重耦合的界面處理技術(shù).在此基礎(chǔ)上,建立完全拉格朗日框架下的損傷破壞模型以準(zhǔn)確獲得固體粒子損傷破壞程度.通過帶彈性擋板潰壩算例驗證了TL-SFPM 方法相比已有的消除拉伸不穩(wěn)定策略的優(yōu)勢,也說明了基于黎曼解接觸算法的適用性和精度.通過子彈撞擊靶板算例驗證了TLSFPM 方法可以較好地避免拉伸不穩(wěn)定導(dǎo)致的數(shù)值斷裂.在此基礎(chǔ)上,加入物理斷裂模型分析靶板撞擊后的損傷模式,不同撞擊速度下靶板的破壞模式差異較大,低速時靶板出現(xiàn)斷口較平齊的裂縫,而高速下靶板會產(chǎn)生碎片.通過與試驗中子彈的剩余速度、靶板變形直徑和彈坑深度對比驗證了本文提出算法的有效性和精度.

本文提出的TL-SFPM 方法,在具備可觀計算效率的前提下提高了沖擊荷載下固體動態(tài)響應(yīng)的計算精度,同時具有損傷破壞模型易施加、三維易擴(kuò)展等優(yōu)勢,對于流固耦合、固體沖擊等問題的高精度數(shù)值模擬具有廣闊的應(yīng)用前景.