基于改進(jìn)PIB搜索法的SPH方法在高速沖擊模擬中的應(yīng)用

肖毅華，張浩鋒，胡德安，平學(xué)成

(1.華東交通大學(xué)機(jī)電工程學(xué)院，南昌330013;2.湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長沙410082)

基于改進(jìn)PIB搜索法的SPH方法在高速沖擊模擬中的應(yīng)用

肖毅華1，張浩鋒1，胡德安2，平學(xué)成1

(1.華東交通大學(xué)機(jī)電工程學(xué)院，南昌330013;2.湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長沙410082)

SPH方法在高速沖擊模擬中具有優(yōu)勢并已得到廣泛使用，但其相鄰粒子搜索非常耗時，其實(shí)際應(yīng)用受計算效率制約。為了提高計算效率，采用一種改進(jìn)的Point－In－Box(PIB)搜索法進(jìn)行相鄰粒子搜索。利用基于該相鄰粒子搜索法的SPH方法，模擬了泰勒桿、平頭彈撞擊鋼板等一系列高速沖擊問題。計算結(jié)果表明:SPH方法能有效地模擬沖擊過程中的大變形、絕熱剪切破壞等現(xiàn)象，模擬結(jié)果與實(shí)驗結(jié)果符合良好;同時，與傳統(tǒng)相鄰粒子搜索法相比，應(yīng)用改進(jìn)的PIB搜索法使SPH方法的計算效率得到明顯提高。

高速沖擊;SPH方法;相鄰粒子搜索

SPH(Smoothed Particle Hydrodynamics)方法是一種無網(wǎng)格拉格朗日粒子法，具有良好的大變形計算能力，易于處理斷裂、破碎等復(fù)雜物理現(xiàn)象，因而在高速沖擊模擬中具有優(yōu)勢并得到廣泛應(yīng)用［1－2］。SPH方法無需復(fù)雜的數(shù)值積分，其場函數(shù)的近似可通過相鄰粒子進(jìn)行簡單求和而得到，故其計算效率比一般的無網(wǎng)格法高。但是，它需搜索每個粒子的相鄰粒子，這是一個很耗時的過程，是制約其計算效率的關(guān)鍵。對于實(shí)際高速沖擊問題的SPH模擬，往往所需的粒子數(shù)比較多，采用高效的相鄰粒子搜索技術(shù)對于控制計算成本具有重要意義。目前，已有一些比較有效的搜索算法被提出和應(yīng)用。Monaghan［3］采用鏈表搜索法實(shí)現(xiàn)相鄰粒子搜索，該方法的復(fù)雜度可以達(dá)到O(N)(N為粒子數(shù))，但是它不能適應(yīng)變光滑長度的情況。Hernquist等［4］提出用樹形搜索法進(jìn)行相鄰粒子搜索，該方法對于大規(guī)模和變光滑長度問題都具有較高的效率和良好的穩(wěn)定性，是一種綜合性能較好的搜索法。Attaway等［5］應(yīng)用Point－In－Box(PIB)搜索法完成相鄰粒子搜索，該方法比較高效而且非常節(jié)省存儲空間，但其效率受粒子分布影響，即:當(dāng)存在某一坐標(biāo)方向的粒子數(shù)很少時，其效率很高，而各坐標(biāo)方向粒子數(shù)都較多時，其效率會明顯下降。

為了提高高速沖擊模擬的計算效率，本文采用了一種改進(jìn)的PIB搜索法即條形化PIB搜索法，該方法具有比現(xiàn)有搜索法更優(yōu)越的性能。應(yīng)用基于該搜索法的SPH方法，模擬了兩個典型的高速沖擊算例，通過與實(shí)驗和其他數(shù)值方法的模擬結(jié)果對比，驗證了計算結(jié)果的有效性，同時對比了基于不同相鄰粒子搜索法的SPH方法的計算效率。

1 模擬高速沖擊過程的軸對稱SPH方法

本文僅考慮具有軸對稱特性的高速沖擊問題，采用二維軸對稱SPH方法進(jìn)行模擬以提高計算效率。軸對稱SPH方法的基本方程采用文獻(xiàn)［6］中按直接離散方式導(dǎo)出的方程，并引入人工黏性和人工應(yīng)力來穩(wěn)定計算和抑制數(shù)值斷裂。具體的求解格式如下:

式中:t為時間，ρ為密度，m為質(zhì)量，v為速度，σ為應(yīng)力，u為單位質(zhì)量的內(nèi)能，r、z和θ分別為徑向、軸向和環(huán)向坐標(biāo)，i和j表示粒子編號，η=2πρr，W為三次樣條核函數(shù)，Q為人工黏性，Rfn為人工應(yīng)力項。人工黏性Q和人工應(yīng)力項Rfn的具體計算式分別參考文獻(xiàn)［2］和文獻(xiàn)［7］。

為避免產(chǎn)生虛假的剪切應(yīng)力和拉伸應(yīng)力、引起嚴(yán)重的計算誤差，本文采用了一種粒子－粒子接觸算法顯式地處理物體間的接觸界面。同時，為了提高接觸計算效率，本文的邊界粒子識別采用了Randles等［8］提出的方法，即根據(jù)物體“顏色”的核近似來確定。

2 條形化PIB搜索法

條形化PIB搜索法是通過改進(jìn)傳統(tǒng)PIB搜索法而得到的。下面先簡要介紹傳統(tǒng)PIB搜索法的原理。首先對所有粒子各個坐標(biāo)方向的坐標(biāo)值分別進(jìn)行排序(見圖1)。當(dāng)搜索給定矩形區(qū)域內(nèi)的粒子時，根據(jù)排序的坐標(biāo)采用二分法快速查找得到位于x、y兩個方向的邊界間的粒子，得到兩個粒子列表。最后，求出兩個粒子列表的交集即為給定矩形區(qū)內(nèi)的粒子。求粒子列表交集的方法為:遍歷粒子數(shù)最少方向的粒子列表中的每個粒子，判斷其在其他方向的坐標(biāo)排序序號是否在相應(yīng)方向的最小和最大坐標(biāo)排序序號之間。例如，假設(shè)y方向的粒子列表中的粒子數(shù)最少，則檢查粒子

那么點(diǎn)i在給定的矩形區(qū)域內(nèi)，否則不在其中。在式(5)和式(6)中，I和R是記錄坐標(biāo)排序結(jié)果的數(shù)組，數(shù)組元素Ix(i)記錄x坐標(biāo)排序序號為i的粒子號;數(shù)組元素Rx(i)記錄粒子i的x坐標(biāo)排序的序號，Iy和Ry針對y坐標(biāo)設(shè)置，功能分別與Ix和Rx相同;和分別為構(gòu)造粒子列表時通過二分法得到的第1個和最后1個位于x方向邊界間的粒子的x坐標(biāo)排序序號，和對應(yīng)于y方向。

圖1 傳統(tǒng)PIB搜索法示意圖Fig.1 Schematic diagram of traditional PIB search algorithm

確定粒子列表交集是PIB搜索法中最耗時的部分。對于條狀分布粒子集，始終有一個方向的粒子列表中的粒子數(shù)很少，按式(5)和式(6)確定粒子列表交集的計算量很小，故PIB搜索法對于此類分布的粒子集搜索效率很高。條形化PIB搜索法就是利用這一特點(diǎn)而提出的。

圖2 條形化PIB搜索法示意圖Fig.2 Schematic diagram of traditional PIB search algorithm

該算法的基本原理見圖2。它把任意分布的粒子集劃分為一些條形分布的子集。劃分子集時，先計算出所有粒子各方向的最大、最小坐標(biāo)值，得到一個恰好包含粒子集的矩形;再以該矩形尺寸較大的方向作為條形方向，并確定一個條形尺寸Δs，將該矩形劃分為若干條狀區(qū)域，每一個條狀區(qū)域內(nèi)的粒子為一個子集。這些子集沿坐標(biāo)方向順序編號，且每個子集內(nèi)的粒子坐標(biāo)值按各個方向獨(dú)立排序。當(dāng)搜索某個給定矩形區(qū)域內(nèi)的粒子時，首先根據(jù)該矩形域的兩條邊界位置坐標(biāo)，計算出其所在的兩個子集的編號，它們?yōu)榭赡馨挥谠摼匦斡騼?nèi)的粒子的首、末兩個子集。然后，利用與傳統(tǒng)PIB搜索法完全一樣的粒子列表構(gòu)造方法和粒子列表交集確定方法，得到這兩個子集中位于該矩形域內(nèi)的粒子;首、末兩個子集之間的粒子集也包含位于該矩形域內(nèi)的粒子，對于這些子集，只需利用傳統(tǒng)PIB搜索法中的粒子列表構(gòu)造方法確定沿條形方向的粒子列表，此粒子列表中的粒子都是位于該矩形域內(nèi)的粒子。最后，合并在每個子集中搜索到的粒子即得到該矩形域內(nèi)的所有粒子(見圖2)。在本文中，條形尺寸大小按式(7)確定

式中:hmax和hmin分別為所有粒子的最大和最小光滑長度。

應(yīng)用條形化PIB搜索法在SPH方法中進(jìn)行相鄰粒子搜索時，首先以粒子i為中心構(gòu)造邊長為4hi的正方形，然后按前述過程搜索得到該正方形內(nèi)的粒子。如果搜索到的粒子j滿足相鄰粒子條件則將粒子i和粒子j記錄為一個相鄰粒子對。

條形化PIB搜索法能有效地克服傳統(tǒng)PIB搜索法對粒子分布形式敏感的缺陷，同時具有很高的搜索效率。下節(jié)將通過算例分析與傳統(tǒng)搜索法中綜合性能較好的樹形搜索法進(jìn)行對比，以驗證其優(yōu)越性。

3 算例分析

3.1 泰勒桿實(shí)驗?zāi)M

本算例模擬House等［9］開展的泰勒桿實(shí)驗，即直徑為7.595 mm、長度為37.97 mm的4340鋼圓柱桿撞擊剛性壁面，模擬了三種不同初始沖擊速度的實(shí)驗情況，分別為181 m/s、224 m/s和270 m/s。模擬結(jié)果與實(shí)驗結(jié)果［9］以及Batra和Zhang［10］的模擬結(jié)果進(jìn)行了對比。計算模型中，圓柱桿用均勻分布的粒子離散，徑向粒子數(shù)為20，軸向粒子數(shù)為200，總共有4 000個粒子;剛性壁面通過鏡像虛粒子方法進(jìn)行處理;4340鋼的材料特性采用Johnson－Cook模型描述，模型相關(guān)參數(shù)參考文獻(xiàn)［10］。

圖3為三種不同初始沖擊速度下桿在60μs時的變形。該時刻桿的動能已經(jīng)很小，桿的沖擊端基本不再產(chǎn)生塑形變形，只有桿的尾部仍沿軸向做微小幅度的振蕩。以此時的桿件形狀作為SPH模擬的桿的最終形狀，得到三種情況下桿的最終長度和沖擊面直徑，結(jié)果見表1。由該表可見，三種工況下SPH模擬結(jié)果與實(shí)驗結(jié)果均符合良好，沖擊面直徑的最大絕對誤差僅0.3 mm，最終長度的最大絕對誤差也僅0.8 mm。與MSPH相比，SPH在沖擊面直徑方面的計算精度要略高，而在最終長度方面的計算精度稍低。

圖3 不同沖速度下桿的變形Fig.3 Deformations of bars at different impact velocities

圖4對比了基于三種不同相鄰粒子搜索方法的SPH模擬泰勒桿實(shí)驗時一個計算循環(huán)內(nèi)各部分的平均計算時間。由圖4可知，相鄰粒子搜索是SPH模擬中最耗時的部分。采用樹形搜索法時，相鄰粒子搜索時間占總時間的比例達(dá)66%。采用PIB搜索法時，相鄰粒子搜索時間減少約2/3，其占總時間的比例降為42%。采用條形化PIB搜索法時，相鄰粒子搜索時間略少于基于PIB搜索法的搜索時間，其占總時間的比例為41%。在本算例中，基于條形化PIB搜索法的SPH效率最高，略高于基于PIB搜索法的SPH，遠(yuǎn)高于基于樹形搜索法的SPH，與其相比節(jié)省了大概45%的總計算時間。PIB搜索法在本算例中也表現(xiàn)出很高的計算效率，這主要由于泰勒桿沿徑向的粒子數(shù)目較少(僅為20)。

表1 泰勒桿實(shí)驗?zāi)M結(jié)果Tab.1 Simulation results of Taylor tests

圖4 泰勒桿模擬中一個計算循環(huán)內(nèi)各部分的平均計算時間Fig.4 Average computational time for each part of a calculation cycle inTaylor test simulation

3.2 平頭彈撞擊鋼板模擬

本算例模擬B?rvik等［11］的平頭彈撞擊金屬板沖塞破壞實(shí)驗，即直徑為20 mm、長度為80mm的Arne工具鋼平頭彈正撞擊直徑為500 mm、厚度為8 mm的Weldox 460E鋼板。模擬工況的初始沖擊速度為173.7m/s。圖5為模擬該實(shí)驗的計算模型。模型共使用104 047個粒子，彈體和靶板中心2倍彈體半徑以內(nèi)的區(qū)域采用均勻粒子離散，粒子間距為0.1 mm;靶板外圍區(qū)域采用沿徑向向外粒子間距按等比例逐漸增大的粒子離散，當(dāng)粒子間距達(dá)到約1 mm時再采用恒定粒子間距離散，保證厚度方向有足夠多的粒子。彈體材料采用雙線性硬化模型，靶板材料采用耦合黏塑性和延性損傷的本構(gòu)模型［12］。

圖5 平頭彈撞擊金屬板的計算模型Fig.5 Computationalmodel for steel plate impacted by blunt projectile

圖6為SPH模擬得到的3個典型時刻的彈體和靶板的變形圖像。在整個過程中，彈體變形很小。靶板發(fā)生絕熱剪切變形，產(chǎn)生沖塞破壞。它在彈體半徑附近發(fā)生了大的塑性變形并產(chǎn)生裂紋，裂紋沿厚度方向發(fā)展，最終穿透靶板的整個厚度，形成一個半徑與彈體半徑大致相等的圓柱塞塊而被彈體沖出。上述模擬得到的現(xiàn)象與實(shí)驗觀察結(jié)果基本相符。圖7給出了彈體的速度時間歷程曲線。計算的彈體速度變化與實(shí)驗結(jié)果趨勢大體相同，在沖擊過程前期計算的彈體速度比實(shí)驗結(jié)果衰減稍快，但最終的彈體剩余速度相差很小。上述結(jié)果說明:SPH方法有效地模擬絕熱剪切失效過程、預(yù)測彈體的剩余速度。

圖6 不同時刻彈體和靶板的變形Fig.6 Deformations of projectile and target plate at different time

圖8對比了基于三種不同相鄰粒子搜索法的SPH模擬平頭彈撞擊鋼板時一個計算循環(huán)中各部分的平均計算時間。在本算例中，粒子總數(shù)多且沿徑向和軸向粒子數(shù)都較多，同時光滑長度在空間上變化也大(1個數(shù)量級)。由圖8可知，相鄰粒子搜索仍是最耗時的計算部分。采用樹形搜索法時，相鄰粒子搜索時間占總時間的比例為60%。采用PIB搜索法時，相鄰粒子搜索時間比樹形搜索法減少45%，占總時間的比例為44%，雖然搜索效率仍有所提高，但提高幅度比上一算例有大幅下降，這主要是受到粒子分布和總數(shù)的影響。而采用條形化PIB搜索法時，相鄰粒子搜索時間仍比樹形搜索法減少約2/3，且比PIB搜索法減少了41%，占總時間的比例為32%。在本算例中，基于條形化PIB搜索法的SPH的計算效率明顯高于基于PIB搜索法和樹形搜索法的SPH，與其相比分別節(jié)省了20%和40%左右的總計算時間。

圖7 彈體的速度變化歷程Fig.7 Velocity history of projectile

圖8 平頭彈撞擊鋼板模擬中一個計算循環(huán)內(nèi)各部分的平均計算時間Fig.8 Average computational time for each part of a calculation cycle in simulation of plate impacted by blunt projectile

4 結(jié)論

介紹了一種基于條形化PIB搜索法的軸對稱SPH方法，應(yīng)用其成功模擬了泰勒桿和平頭彈撞擊鋼板兩個問題。結(jié)果表明，SPH方法可以較準(zhǔn)確地模擬泰勒桿變形和鋼板的絕熱剪切破壞，能有效地預(yù)測彈體的剩余速度。文中對計算效率進(jìn)行了對比分析，結(jié)果表明:基于條形化PIB搜索法的SPH有很高的計算效率，與基于樹形搜索法的SPH相比可以大大減少相鄰粒子搜索時間，使整體計算效率得到明顯提高，與基于PIB搜索法的SPH相比能克服其效率受粒子分布影響的缺陷，在計算效率方面具有良好的穩(wěn)定性。此外，基于條形化PIB搜索法的SPH對于大規(guī)模和光滑長度空間不均的問題具有良好的適應(yīng)性。

［1］Libersky L D，Petschek A G，Carney TC，etal.High strain Lagrangian hydrodynamics:a three－dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109(1):67－75.

［2］Johnson G R，Beissel S R.SPH for high velocity impact computations［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1/2/3/4):347－373.

［3］Monaghan J J.Particle methods for hydrodynamics［J］.Computer Physics Reports，1985，3(2):71－124.

［4］Hernquist L，Katz N.TREESPH－A unification of SPH with the hierarchical tree method［J］.The Astrophysical Journal Supplement Series，1989，70:419－446.

［5］Attaway SW，Heinstein M W，Mello F J，et al.Coupling of smooth particle hydrodynamics with PRONTO［R］.SAND93－1781C.Albuquerque:Sandia National Laboratories，1993.

［6］楊剛.光滑粒子法的改進(jìn)及其若干典型應(yīng)用［D］.長沙:湖南大學(xué)，2011.

［7］Xiao Y H，Hu D A，Han X，et al.Simulation of normal perforation of aluminum plates using axisymmetric smoothed particle hydrodynamics with contact algorithm［J］.International Journal of Computational Methods，2013，10(3)1350039 (21pages).

［8］Randles PW，Libersky L D.Smoothed particle hydrodynamics:some recent improvements and applications［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1/2/3/4):375－408.

［9］House JW，Lewis JC，Gillis P P，et al.Estimation of flow stress under high rate plastic deformation［J］.International Journal of Impact Engineering，1995，16:189－200.

［10］Batra R C，Zhang G M.Modified smoothed particle hydrodynamics(MSPH)basis functions for meshless methods，and their application to axisymmetric Taylor impact test［J］.Journal of Computational Physics，2008，227:1962－1981.

［11］B?rvik T，Hopperstad O S，Berstad T，et al.Numerical simulation of plugging failure in ballistic penetration［J］.International Journal of solids and structures，2001，38(34/35): 6241－6264.

［12］B?rvik T，Hopperstad O S，Berstad T，et al.A computationalmodel of viscoplasticity and ductile damage for impact and penetration［J］.European Journal of Mechanics－A/Solids，2001，20(5):685－712.

Application of SPH w ith modified PIB search algorithm to high-velocity impact simulation

XIAO Yi－h(huán)ua1，ZHANG Hao－feng1，HU De－an2，PING Xue－cheng1
(1.College of Mechanical and Electrical Engineering，East China Jiaotong University，Nanchang 330013，China; 2.College of Mechanical and Vehicle Engineering，Hunan University，Changsha 410082，China)

SPH method hasmany advantages and iswidely used in high－velocity impact simulation，but its practical applications are still restricted by its low computational efficiency due to the time－consuming neighboring particle search (NPS).To improve its efficiency，a modified point－in－box(PIB)search algorithm was employed to determine neighboring particles.By using the SPH with modified PIB search algorithm，a series of impact problems，including Taylor tests and steel plate impacted by blunt projectile，were simulated.The calculated results indicated that SPH can effectively simulate large deformations and adiabatic shear failure during the impact processes，and give results in good agreementwith experiments.The results also showed that the modified PIB search algorithm can make the SPH method achieve amuch better efficiency than a traditional search algorithm.

high－velocity impact;SPH method;neighboring particle search

O353.4;O242.1

A

10.13465/j.cnki.jvs.2015.23.016

國家自然科學(xué)基金項目(11302077);江西省教育廳科學(xué)技術(shù)研究項目(GJJ14398)

2014－09－12修改稿收到日期:2014－11－06

肖毅華男，博士，講師，1984年生

平學(xué)成男，博士，教授，1975年生