基于改進PIB搜索法的SPH方法在高速沖擊模擬中的應用

肖毅華，張浩鋒，胡德安，平學成

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

基于改進PIB搜索法的SPH方法在高速沖擊模擬中的應用

肖毅華1，張浩鋒1，胡德安2，平學成1

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

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

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

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

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

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

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

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

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

2 條形化PIB搜索法

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

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

圖1 傳統PIB搜索法示意圖Fig.1 Schematic diagram of traditional PIB search algorithm

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

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

該算法的基本原理見圖2。它把任意分布的粒子集劃分為一些條形分布的子集。劃分子集時，先計算出所有粒子各方向的最大、最小坐標值，得到一個恰好包含粒子集的矩形;再以該矩形尺寸較大的方向作為條形方向，并確定一個條形尺寸Δs，將該矩形劃分為若干條狀區域，每一個條狀區域內的粒子為一個子集。這些子集沿坐標方向順序編號，且每個子集內的粒子坐標值按各個方向獨立排序。當搜索某個給定矩形區域內的粒子時，首先根據該矩形域的兩條邊界位置坐標，計算出其所在的兩個子集的編號，它們為可能包含位于該矩形域內的粒子的首、末兩個子集。然后，利用與傳統PIB搜索法完全一樣的粒子列表構造方法和粒子列表交集確定方法，得到這兩個子集中位于該矩形域內的粒子;首、末兩個子集之間的粒子集也包含位于該矩形域內的粒子，對于這些子集，只需利用傳統PIB搜索法中的粒子列表構造方法確定沿條形方向的粒子列表，此粒子列表中的粒子都是位于該矩形域內的粒子。最后，合并在每個子集中搜索到的粒子即得到該矩形域內的所有粒子(見圖2)。在本文中，條形尺寸大小按式(7)確定

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

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

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

3 算例分析

3.1 泰勒桿實驗模擬

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

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

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

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

表1 泰勒桿實驗模擬結果Tab.1 Simulation results of Taylor tests

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

3.2 平頭彈撞擊鋼板模擬

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

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

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

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

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

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

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

4 結論

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

［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］楊剛.光滑粒子法的改進及其若干典型應用［D］.長沙:湖南大學，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－hua1，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

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

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

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

平學成男，博士，教授，1975年生