高速沖擊問題的SPH粒子接觸算法三維數(shù)值計(jì)算

張志春，卞 強(qiáng)，展文豪，劉桐宇

(中國航天員科研訓(xùn)練中心人因工程國防科技重點(diǎn)實(shí)驗(yàn)室， 北京 100094)

由于在計(jì)算大變形問題方面的優(yōu)勢，光滑粒子流體動力學(xué)方法(Smoothed Particle Hydrodynamics, SPH)在高速沖擊問題的計(jì)算中，得到了廣泛的應(yīng)用。肖毅華[1]基于改進(jìn)PIB搜索法的SPH方法模擬了泰勒桿的沖擊問題；于永強(qiáng)[2]使用SPH方法計(jì)算了鳥撞復(fù)合材料層合板問題；呂東喜[3]使用SPH方法對磨粒沖擊工件表面過程進(jìn)行了數(shù)值模擬。

使用SPH方法模擬高速沖擊問題，子彈和靶板的接觸是計(jì)算的關(guān)鍵。傳統(tǒng)上采用SPH連續(xù)速度算法模擬彈體和靶板的接觸。連續(xù)速度算法基于SPH粒子求和思想，在求解SPH動量方程時(shí)，將支持域內(nèi)所有的靶板和彈體粒子均考慮在內(nèi)，實(shí)現(xiàn)子彈和靶板的相互作用。由于SPH連續(xù)速度算法不能區(qū)分不同材料的SPH粒子，當(dāng)子彈和靶板靠近時(shí)，該方法能夠描述兩者的“排斥”作用；但當(dāng)子彈和靶板分離時(shí)，該方法卻使兩者相互“吸引”，明顯與實(shí)際不符合。同時(shí)連續(xù)速度算法很容易使接觸面兩側(cè)粒子相互“融合”，不同材料的粒子侵入彼此“領(lǐng)域”，造成接觸界面模糊，難以追蹤。

為了描述SPH粒子之間的接觸作用，Monaghan[4]提出了XSPH算法，包含了相鄰粒子的影響作用，使粒子的速度與相鄰粒子的平均速度相近。該方法存在如下缺陷：不能計(jì)算拉伸力，確保兩物體分離；不能計(jì)算剪切力。Campbell等[5]提出了一種粒子-粒子罰接觸算法，使用Randles和Libersky[6]提出的方法來確定邊界，在三維計(jì)算中較為復(fù)雜。Vignjevic等[7-8]借鑒Monaghan[9]SPH粒子排斥力算法的思想，提出了一種SPH粒子接觸算法，通過在SPH粒子上施加接觸力的方式實(shí)現(xiàn)SPH粒子間的接觸作用，取得了較好的效果，文獻(xiàn)[10]對該算法進(jìn)行了應(yīng)用。

本文基于SPH粒子接觸算法，采用C++編寫的SPH計(jì)算程序，分別對球頭鋼彈斜沖擊鋼板以及Arne tool鋼制圓柱子彈正沖擊Weldox 460 E鋼板發(fā)生沖塞破壞的過程進(jìn)行了三維數(shù)值計(jì)算，通過在子彈和靶板相關(guān)粒子上施加接觸力實(shí)現(xiàn)兩者之間的接觸。將SPH粒子接觸算法的計(jì)算結(jié)果、LS-DYNA970的計(jì)算結(jié)果、SPH連續(xù)速度算法的計(jì)算結(jié)果、實(shí)驗(yàn)結(jié)果進(jìn)行了對比，驗(yàn)證了SPH粒子接觸算法在模擬高速沖擊問題中的有效性。

1 考慮材料強(qiáng)度的SPH方程

在具有材料強(qiáng)度的高速沖擊問題中，沖擊波沿著碰撞體內(nèi)傳播，此時(shí)碰撞體具有流體特性。其中運(yùn)動方程和高壓狀態(tài)方程是控制材料力學(xué)行為的關(guān)鍵，并且材料強(qiáng)度只有在這種具有能量傳動問題的后期才能顯示出其重要性。由SPH離散的具有材料強(qiáng)度的流體動力學(xué)控制方程如下[11-12]：

連續(xù)性方程：

(1)

動量方程：

(2)

能量方程：

(3)

運(yùn)動方程：

(4)

其中mi，ρi，vi，ei分別表示SPH粒子i的質(zhì)量、密度、速度和內(nèi)能，x表示位置。Wij=W(xi-xj,h)表示SPH光滑核函數(shù)，h是光滑長度。Πij表示人工黏度，用于消除非物理振蕩。

材料強(qiáng)度體現(xiàn)在應(yīng)力張量σαβ，由各向同性壓力p和黏性剪切應(yīng)力ταβ組成

σαβ=-pδαβ+ταβ

(5)

其中偏應(yīng)力率采用Jaumann比率的形式表示為

(6)

(7)

(8)

其SPH離散格式分別表示為

(9)

(10)

2 SPH粒子接觸算法

傳統(tǒng)的處理SPH接觸問題的算法如圖1所示，其中小實(shí)線圓表示SPH粒子，大虛線圓表示粒子的支持域。當(dāng)物體A、B靠近時(shí)，采用SPH動量方程對A、B內(nèi)粒子速度進(jìn)行更新。進(jìn)行SPH連續(xù)速度算法時(shí)，粒子i的支持域包含A內(nèi)的粒子n4,n5,…,n8和B內(nèi)的粒子n1,n2,n3。Vignjevic等[7-8]提出的SPH無摩擦接觸算法通過在SPH粒子上施加接觸力的方式實(shí)現(xiàn)SPH粒子間的接觸作用。如圖2所示，當(dāng)A、B間距達(dá)到兩倍光滑長度時(shí)，在相關(guān)SPH粒子上產(chǎn)生接觸力(SPH粒子支持域的半徑選為兩倍光滑長度)；當(dāng)A、B間距超過2倍光滑長度時(shí)，沒有接觸力產(chǎn)生。首先定義一個(gè)接觸勢函數(shù)

(11)

其中ncont表示SPH粒子i支持域內(nèi)屬于不同體的粒子數(shù)，如圖2中，對于粒子i，ncont=3，即計(jì)算粒子i的接觸力時(shí)其支持域僅包含物體B內(nèi)的粒子n1,n2,n3。當(dāng)xA和xB屬于同一體時(shí)，SPH核函數(shù)W(xA-xB,h)=0。Δdavg表示粒子間距的平均值。K和n表示用戶定義參數(shù)，K的選取與有限元接觸中的罰剛度類似，與材料性質(zhì)和沖擊速度相關(guān)。該接觸勢函數(shù)具有如下特點(diǎn)：在物體內(nèi)部為0；通常為正值；隨粒子間距的減小而增大。對該接觸勢函數(shù)的梯度進(jìn)行SPH格式離散，得到施加在粒子上的接觸力

▽xiW(rij)

(12)

接觸力的方向由接觸勢函數(shù)的梯度確定。將式(12)代入式(2)，得到含有接觸力的SPH動量方程

(13)

采用連續(xù)速度算法處理SPH的接觸問題，當(dāng)物體A、B分離時(shí)，由于粒子求和的作用，A、B粒子會產(chǎn)生相互“吸引”，與實(shí)際情況不符合。而接觸力算法能夠客觀的反應(yīng)A、B之間的接觸作用，較為真實(shí)的描述A、B的靠近和分離。該粒子接觸算法與SPH的求解格式保持一致，可以充分利用SPH的臨近搜索列表。只要當(dāng)兩物體接觸界面粒子間距達(dá)到兩倍光滑長度，就在相關(guān)粒子上施加接觸力，當(dāng)接觸界面粒子間距超過兩倍光滑長度時(shí)，沒有接觸力，因此不需要定義接觸粒子。該算法和有限元接觸算法不同，在每個(gè)時(shí)間步不需要定義接觸界面和界面法線方向，便于進(jìn)行三維程序設(shè)計(jì)。

3 球頭鋼彈斜沖擊鋼板

本節(jié)采用SPH粒子接觸算法對球頭鋼彈斜沖擊鋼板進(jìn)行了三維數(shù)值計(jì)算。子彈和靶板的材料參數(shù)為：密度ρ=7 830 kg/m3，彈性模量E=2.07×1011Pa，泊松比v=0.3。子彈長0.078 m，直徑0.018 m，靶板尺寸為0.09 m×0.09 m×0.009 m，子彈和靶板的SPH粒子分別為2 740和2 700，粒子間距均為0.003 m。子彈初速v0=300 m/s，沿Z軸正向，子彈與靶板的初始夾角為60°。為了進(jìn)行對比驗(yàn)證，采用LS-DYNA970對相同問題進(jìn)行了計(jì)算，子彈和靶板均采用六面體單元離散，單元數(shù)分別為1 824和2 700。計(jì)算時(shí)間均為100×10-6s。

在t=70×10-6s時(shí)，SPH接觸算法和LS-DYNA970數(shù)值計(jì)算結(jié)果如圖3所示。在沖擊載荷作用下，靶板發(fā)生較大的變形，而子彈變形相對較小。子彈和靶板沿Z軸的速度時(shí)程曲線如圖4所示，表明兩種算法的計(jì)算結(jié)果吻合較好。兩種方法的誤差主要有兩方面的原因：由于離散方式的不同，在材料密度相同的前提下，兩種方法中子彈和靶板的質(zhì)量存在一定差異；SPH粒子接觸算法使子彈和靶板發(fā)生接觸的時(shí)間比LS-DYNA970要早。通過減小SPH粒子的初始間距，可以一定程度上減小這兩種原因造成的誤差。通過球頭鋼彈斜沖擊鋼板的數(shù)值計(jì)算表明：SPH粒子接觸算法適用于斜沖擊和接觸界面形狀復(fù)雜的問題。

4 平頭鋼彈正沖擊鋼板

本節(jié)采用SPH粒子接觸算法對圓柱形鋼彈正沖擊鋼板發(fā)生沖塞破壞的過程進(jìn)行了三維數(shù)值計(jì)算，計(jì)算模型的尺寸及離散方式如圖5所示。子彈和靶板均由SPH粒子離散，粒子總數(shù)57 728。子彈直徑0.02 m，高0.08 m，粒子間距為10-3m。靶板尺寸0.28 m×0.28 m×0.012 m，中心0.04 m×0.04 m×0.012 m區(qū)域內(nèi)粒子間距10-3m，為減小計(jì)算量，粒子間距從中心正碰區(qū)向邊緣逐漸增大。子彈與靶板初始間距2×10-3m，靶板四周采用固支邊界條件，計(jì)算總時(shí)間160×10-6s。采用Monaghan[9]提出的人工應(yīng)力法來消除拉伸不穩(wěn)定性造成的數(shù)值斷裂，采用強(qiáng)洪夫[13-14]提出的完全變光滑長度法來消除光變光滑長度帶來的影響。

子彈材料為Arne tool鋼，沖擊過程中變形較小，采用線彈性模型。靶板材料為Weldox 460 E鋼，沖擊中發(fā)生沖塞型破壞，變形較大，為描述靶板材料的屈服應(yīng)力及損傷演化，采用B?rvik等[15]提出的修正Johnson-Cook強(qiáng)度模型和Gruneisen狀態(tài)方程。該模型中將材料的屈服強(qiáng)度表示為損傷變量、等效應(yīng)變、等效應(yīng)變率和溫度的函數(shù)

(14)

(15)

T0是室溫，Tm是材料熔點(diǎn)。假設(shè)為絕熱條件，靶板材料溫度的升高是由于沖擊過程中的塑性功轉(zhuǎn)換為熱量造成的，溫度的變化計(jì)算公式為

(16)

Cp是材料比熱，α是比例常數(shù)，Wp是塑性功。D為損傷變量，D=0表示材料沒有損傷，D=1表示材料完全失效。實(shí)際上材料出現(xiàn)宏觀裂紋時(shí)，損傷變量的臨界值小于1，失效準(zhǔn)則可描述為D=DC≤1，其中DC是損傷變量臨界值。損傷變量D是累積塑性應(yīng)變p的函數(shù)，可描述如下：

(17)

其中pd是損傷閾值，pf是斷裂等效塑性應(yīng)變，與材料的應(yīng)力三軸度、應(yīng)變率和溫度相關(guān)。Johnson和Cook在文獻(xiàn)[16]提出的剪切損傷演化模型中將pf描述如下

(18)

其中D1～D5為材料常數(shù)，σ*=σm/σeq為應(yīng)力三軸度，σm=(σx+σy+σz)/3為平均正應(yīng)力。子彈和靶板的材料參數(shù)分別見表1、表2。

表1 子彈材料參數(shù)

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

圖6顯示了子彈初速為285.4 m/s時(shí)的數(shù)值計(jì)算結(jié)果，而相關(guān)的實(shí)驗(yàn)結(jié)果如圖7、圖8所示[15]。在t=7×10-6s，子彈與靶板開始接觸，子彈正前方區(qū)域靶板開始加速，子彈頭部周圍的靶板出現(xiàn)剪切區(qū)。靶板被擠壓變形，背部出現(xiàn)凸起。由于子彈-靶板接觸界面的塑性應(yīng)變，剪切區(qū)內(nèi)粒子出現(xiàn)損傷。靶板變形繼續(xù)增大，粒子損傷迅速演化，部分粒子損傷達(dá)到臨界值，開始失效，形成裂紋并逐漸向靶板背部擴(kuò)展。在沖擊后期，靶板的失效模式包含靶板內(nèi)部的剪切斷裂和靶板背部的拉伸損傷。最終形成靶塞，完全從靶板脫離。從中可以看出子彈、靶塞和靶板沖塞型破壞的計(jì)算結(jié)果均與實(shí)驗(yàn)吻合較好。子彈變形很小，表明數(shù)值計(jì)算中子彈選擇線彈性模型可行。

表3給出了SPH粒子接觸算法、SPH連續(xù)速度算法和實(shí)驗(yàn)結(jié)果[15]，其中v0表示子彈初速，vpr表示子彈余速，vtr表示靶塞速度，Δvpr=(vpr計(jì)算值-vpr實(shí)驗(yàn)值)/vpr實(shí)驗(yàn)值，Δvtr=(vtr計(jì)算值-vtr實(shí)驗(yàn)值)/vtr實(shí)驗(yàn)值。對比顯示SPH粒子接觸算法的計(jì)算精度要明顯高于SPH連續(xù)速度算法，表明該方法能夠更為真實(shí)的描述子彈和靶板的接觸問題。并且計(jì)算精度隨子彈初速的增大而提高，表明SPH粒子接觸算法對于高速沖擊領(lǐng)域具有較好的適用性。

表3 數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果

5 結(jié)論

由于在處理大變形方面具有明顯的優(yōu)勢，SPH方法被成功地應(yīng)用于高速沖擊問題的數(shù)值計(jì)算。為了計(jì)算不同物體之間SPH粒子的接觸作用，本文將SPH粒子接觸算法應(yīng)用于高速沖擊。該算法是通過在SPH粒子上施加接觸力的方式實(shí)現(xiàn)，不需要定義接觸界面以及接觸界面的法線方向，便于實(shí)現(xiàn)三維編程。接觸力的計(jì)算與接觸物體的相對運(yùn)動方向、接觸面的幾何形狀無關(guān)，因此該接觸算法對正沖擊、斜沖擊、復(fù)雜接觸界面都有較好的適用性。算例中與其他數(shù)值方法以及實(shí)驗(yàn)結(jié)果的對比顯示了該SPH粒子接觸算法的有效性。

將SPH粒子接觸算法應(yīng)用于高速沖擊問題，當(dāng)子彈和靶板間距達(dá)到兩倍光滑長度時(shí)，在相關(guān)SPH粒子上施加接觸力。因此接觸作用產(chǎn)生的時(shí)機(jī)比實(shí)際情況早，產(chǎn)生一定的誤差。為了減小該誤差，在前處理建模時(shí)，要加密接觸部位的SPH粒子設(shè)置，使接觸作用的計(jì)算與實(shí)際情況更加相符。