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

王璐 徐緋 楊揚(yáng)

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

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

引言

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

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

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

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

1 TL-SPH法

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

材料變形時(shí)粒子的柯西應(yīng)力σ包括各項(xiàng)同性壓力p和偏應(yīng)力S,壓力p通過式(1)中的壓力-密度狀態(tài)方程得到.沖擊過程中材料不可避免地會(huì)進(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 核近似精度受粒子分布和支持域完整性影響,固體邊界附近粒子的核估計(jì)精度會(huì)降低[5],下面在完全拉格朗日框架下推導(dǎo)SPH 高階插值公式.借鑒基于當(dāng)前構(gòu)型的高階插值公式的思想,以二維情形為例,將f(X)在初始構(gòu)型下的點(diǎn)Xi處泰勒展開,省去二階及以上導(dǎo)數(shù)項(xiàng),分別在等式兩邊乘以Wi,X和Wi,Y并在Ω區(qū)域內(nèi)積分,可得到[5]

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

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

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

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

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

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

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

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

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

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

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

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

2.1 變形梯度求解

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

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

2.2 加速度求解

若直接使用式(6) 求解加速度方程中的第一Piola-Kirchhoff 應(yīng)力P的梯度,此時(shí)會(huì)導(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對(duì)粒子j的作用力與粒子j對(duì)粒子i的作用力一致.若考慮人工黏性,將人工黏性直接加到加速度方程即可

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

3 界面處理

3.1 TL-DFPM法

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

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

圖4 二維情形下的不連續(xù)計(jì)算域示意圖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ù)性近似分析,該修正式可以保證梯度計(jì)算的精確性,也能夠?qū)崿F(xiàn)不連續(xù)界面處粒子信息的準(zhǔn)確估計(jì).若將上述方法引入完全拉格朗日框架中,可得

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

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

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

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

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

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

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

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

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

3.2.1 流固接觸

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

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

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

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

3.2.2 固固接觸

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

4 粒子損傷破壞模型

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

5 算例驗(yàn)證

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

這節(jié)以流固耦合沖擊問題為例驗(yàn)證本文提出的TL-SFPM 方法和界面處理技術(shù)的優(yōu)勢(shì).如圖7 所示,紅色區(qū)域?yàn)閺椥該醢?上端與剛性固壁固連,擋板內(nèi)藍(lán)色區(qū)域?yàn)樗?在重力作用下,水坍塌沖擊彈性擋板,擋板發(fā)生變形.擋板的材料參數(shù)和模型幾何尺寸見表1.水初始密度為1000kg/m3,水近似為弱可壓模型[35-36],用低耗散GSPH 方法計(jì)算.低耗散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)致流動(dòng)失真,將p*中的聲速限制為[38]

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

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

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

圖10 TL-SFPM 模擬結(jié)果和試驗(yàn)結(jié)果[39]對(duì)比圖(自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 對(duì)比了某一時(shí)刻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] 模擬得到某時(shí)刻彈性擋板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八核心處理器下的計(jì)算時(shí)間.可以看出TL-SFPM 相比TL-SPH 可節(jié)省一定的計(jì)算時(shí)間,這是因?yàn)門L-SFPM 中鄰域內(nèi)參與計(jì)算的粒子數(shù)相比TL-SPH 更少.雖然TL-SFPM 方法采用的是高階格式,但是由于初始構(gòu)型下粒子相對(duì)位置確定,高階格式的修正矩陣Bi只需在初始時(shí)刻計(jì)算一次,后續(xù)直接調(diào)用即可.二者計(jì)算時(shí)間相差并不大是因?yàn)楸舅憷袕椥泽w粒子數(shù)所占總粒子數(shù)比重較少,而且兩種方法均需基于當(dāng)前構(gòu)型計(jì)算流固接觸力.

表2 TL-SPH和TL-SFPM 模擬時(shí)間對(duì)比Table 2 Computational time with TL-SPH and TL-SFPM methods

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

這節(jié)以子彈撞擊靶板算例驗(yàn)證本文提出的TLSFPM 方法、界面處理技術(shù)以及損傷破壞模型的合理性和精確性.圓頭的鋼彈撞擊鋁合金的飛機(jī)蒙皮[40],彈丸直徑12.7 mm,長(zhǎng)度51 mm,靶板長(zhǎng)和寬均為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 無(wú)損傷破壞模型

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

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

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

5.2.2 有損傷破壞模型

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

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

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

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

圖15 鋼彈沖擊靶板計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比[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八核心處理器下的計(jì)算時(shí)間.可以看出TL-SFPM 方法相比CSPM 方法計(jì)算時(shí)間略長(zhǎng),雖然TL-SFPM 方法可節(jié)省粒子鄰域搜索的時(shí)間,但是當(dāng)粒子進(jìn)入塑性后,TL-SFPM 方法的計(jì)算過程相比CSPM 更復(fù)雜,需要將柯西應(yīng)力σ轉(zhuǎn)化為第一Piola-Kirchhoff 應(yīng)力P.另外,TL-SFPM 方法斷裂模型中需要及時(shí)對(duì)核函數(shù)更新.當(dāng)不同材料作用,TL-SFPM 也需在每一時(shí)間步內(nèi)計(jì)算當(dāng)前構(gòu)型下的壓力黎曼近似解p*以求解不同材料間的接觸力.

表4 CSPM和TL-SFPM 模擬時(shí)間對(duì)比Table 4 Computational time with CSPM and TL-SFPM methods

6 結(jié)論

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

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