基于ABAQUS的SPH粒子生成及其在高速沖擊計算中應用

肖毅華，董晃晃，平學成

(華東交通大學 機電與車輛工程學院，南昌　330013)

?

基于ABAQUS的SPH粒子生成及其在高速沖擊計算中應用

肖毅華，董晃晃，平學成

(華東交通大學 機電與車輛工程學院，南昌330013)

提出一種新的基于ABAQUS的SPH粒子生成方法。采用Python語言編寫腳本，通過ABAQUS運行腳本實現粒子劃分。在腳本中，利用ABAQUS提供的內部函數，先導入幾何模型，再進行幾何區域均勻正交分割，最后將分割所得的各子區域轉化成粒子。該方法適用于任意區域，生成的粒子以規則的矩形分布為主。通過對半圓環區域進行粒子劃分，對比了基于該方法和傳統方法生成的粒子的近似精度。以截卵形彈體為例，說明了通過該方法結合參數化幾何建模易于實現參數化粒子模型的建立。應用生成的粒子進行了高速沖擊問題的SPH模擬，驗證了該方法在實際應用中的可行性。

SPH；粒子生成；ABAQUS；腳本；高速沖擊

光滑粒子動力學法(SPH)是一種無網格拉格朗日粒子法，在計算中無需網格，能有效地處理自由表面、移動界面、材料大變形、破壞和飛濺等現象。它在流體流動[1]、爆炸[2-3]、高速沖擊[4-5]等領域中得到了廣泛應用。

SPH需在初始時刻對問題域用粒子離散。初始粒子分布對計算結果有一定的影響。研究表明，規則、均勻的粒子分布可以給出更準確的結果[6]、使數值計算更穩定[7]。SPH中最常用的粒子分布形式是矩形粒子分布，有時為了保證計算精度且減少計算量，采用分區變化粒子間距的矩形粒子分布[5]、漸進變化粒子間距的圓形粒子分布[8]等。這類粒子分布對于規則的問題域如矩形區域容易生成，而對于具有曲線邊界的復雜問題域生成很繁瑣。傳統的SPH粒子生成方法一般是根據具體問題域的邊界等設計專門的算法和程序，如熱合買提江·依明等[9]研究了齒輪模型的粒子生成方法。這樣的方法通用性差，對于復雜區域生成困難、效率低。國內外關于通用粒子生成方法的研究很少。倪國喜等[6,10]在這方面進行了一定的研究，提出了基于邊界預處理結合Delaunay三角劃分的任意區域粒子生成方法，能保證粒子的Voronoi面積較均勻。另外，一些學者[2]采取先生成三角形單元再轉化為粒子的單元轉化法，也容易實現復雜區域的粒子離散。然而，這兩種方法的初始粒子分布一般都是不規則的。

本文提出一種基于ABAQUS的粒子生成方法。ABAQUS不僅具備強大的復雜系統分析能力，而且向用戶提供了豐富的內部函數以進行二次開發。用戶利用Python語言編程，可以調用其內部函數，避開圖形界面操作，直接操縱其內核，實現幾何建模、網格劃分、邊界和載荷定義、作業提交、結果后處理等有限元分析全過程的操作。本文方法正是利用ABAQUS的這種二次開發功能，調用其內部函數，先讀取問題域的幾何模型，然后搜索得到問題域的矩形包圍盒，再根據指定的粒子間距在包圍盒范圍內對問題域進行正交分割，最后將分割得到的子區域轉化生成粒子。該方法不僅通用性強，且生成的粒子以規則的矩形分布為主。

1　粒子生成方法

本文方法包括幾何模型導入、區域自動分割、粒子物理量計算3個步驟。下面詳細敘述每個步驟的具體過程以及該方法的實現與運用。

1.1幾何模型導入

幾何模型通過調用ABAQUS的內部函數PartFromGeometryFile來導入。模型可以由ABAQUS/CAE或CATIA、 UG、ProE等專業建模軟件建立，保存為ABAQUS能夠讀取的文件格式。以ABAQUS/CAE生成的“.sat”格式的幾何模型文件為例，導入二維幾何模型的命令如下：

acis =mdb.openAcis('myFile.sat',scaleFromFile=OFF)

mdb.models[' Model-1'].PartFromGeometryFile(

name=' myPart',geometryFile=acis,

dimensionality=TWO_D_PLANAR,

type=DEFORMABLE_BODY)

其中：myFile為待進行粒子劃分的幾何模型的文件名，myPart 為由導入的幾何模型創建的部件名稱，TWO_D_PLANAR和DEFORMABLE_BODY為ABAQUS的預定義常數，分別表示二維模型和變形體。

圖1　支架模型及其矩形包圍盒Fig.1　The bracket model and its rectangular bounding box

為后續敘述方便同時說明本文方法的通用性，以具有復雜形狀的支架結構為例，設其幾何模型文件名即為myFile.sat，運行上述命令后，將在ABAQUS中創建圖1所示的支架部件。該支架的主要結構尺寸參見文獻[1]，其中增加的橢圓孔結構的長、短半軸分別為60 mm和28 mm。

1.2區域自動分割

為了對導入幾何模型的區域進行分割，先找到該區域的最小矩形包圍盒(如圖1所示的虛線框)。該包圍盒的邊界通過調用ABAQUS提供的內部函數getBoundingBox來自動尋找。實現邊界尋找的命令如下：

myPart =mdb.models[' Model-1'].parts[' myPart ']

myFace=myPart.faces

myBound=myFace.getBoundingBox()

通過上述命令，myBound返回矩形包圍盒各個坐標方向的下界和上界，記x和y方向的下界分別為xmin和ymin、上界分別為xmax和ymax。根據邊界位置，先沿x方向在邊界線x=xmin和x=xmax之間對幾何區域按從左至右的順序進行等距分割，然后再沿y方向在邊界線y=ymin和y=ymax之間按從下至上進行等距分割。這樣，幾何區域被劃分為許多子區域。每個子區域將在后續步驟中轉化生成1個粒子。分割時，x和y方向的分割份數nx和ny根據指定的粒子間距Δ來計算，分別為

nx=int[(xmax-xmin)/Δ]

(1)

ny=int[(ymax-ymin)/Δ]

(2)

相應地，x和y方向的實際分割間距Δx和Δy分別為

Δx=(xmax-xmin)/nx

(3)

Δy=(ymax-ymin)/ny

(4)

上述分割過程通過循環調用ABAQUS的內部函數PartitionFaceByDatumPlane來實現。該函數通過基準平面來對幾何區域進行分割。以x方向為例，分割幾何區域的命令如下：

axis1 =myPart.DatumAxisByTwoPoint(

point1=(0.0,0.0,0.0),point2=(1.0,0.0,0.0))

xcoor =xmin+deltax

while xcoor

d=myPart.DatumPointByCoordinate(

coords=(xcoor,0.0,0.0))

myFace=myPart.faces[:]

myDatumplane=

myPart.DatumPlaneByPointNormal(

point=datums[d.id],normal=datums[axis1.id])

myPart.PartitionFaceByDatumPlane(

datumPlane=datums[myDatumplane.id],

faces=myFace)

xcoor=xcoor+deltax

在上述命令流中，首先定義了沿x方向的基準軸axis1、確定了開始分割時基準平面在x方向上的位置xcoor，然后通過while循環進行連續分割。在每次循環中，先調用內部函數DatumPlaneByPointNormal根據基準點d和基準軸axis1定義1個垂直于axis1的基準平面，再通過內部函數PartitionFaceByDatumPlane利用定義的基準平面對幾何區域完成1次分割。每完成1次分割，決定基準平面位置的變量xcoor(即基準點d的x坐標)增加1個分割間距deltax。這樣，基準平面在循環過程中不斷右移，從而實現對幾何區域按從左至右的順序等距分割。類似地，y方向的分割也可以完成。圖2給出了分割支架模型時產生的基準點、軸和平面以及分割后的支架模型。

圖2　支架模型的分割Fig.2　Partition of the bracket

1.3粒子物理量計算

粒子布置在分割后得到的每個子區域的面心處。粒子的位置坐標和質量分別取為對應子區域的面心坐標和質量。任意子區域的面心坐標和面積分別調用ABAQUS的內部函數getCentroid和getSize獲得，從而粒子的坐標可以直接得到，而其質量則需進一步按下式計算

(5)

式中：mi為粒子i的質量；ρi為密度，即區域的材料密度；Ai為粒子i對應的子區域的面積；t為平面問題中的厚度，ri為軸對稱問題中粒子i的徑向坐標。粒子的光滑長度取為

hi=αhΔ

(6)

式中：Δ為指定的粒子間距，αh為光滑長度系數。在接觸計算中，需要粒子的半徑來表征其實際大小以用于接觸檢測。粒子半徑按下式定義

(7)

圖3為本文方法生成的支架粒子模型。可見，大部分的內部粒子嚴格按矩形分布，具有完全相等的面積和質量。區域邊界處會出現少數非矩形的子區域，粒子不再嚴格按矩形分布。

圖3　支架的SPH粒子模型Fig.3　SPH particle model for the bracket

1.4方法實現與運行

采用Python語言編程，依次實現上述3個步驟的操作，將編寫的程序保存為1個“.py”為后綴的腳本文件。對幾何模型劃分粒子時，只需打開該腳本文件，修改文件名參數“myFile”為實際的模型文件名，同時設置好粒子間距Δ，然后在ABAQUS的“Run Script”菜單下運行該腳本，即可完成區域的粒子劃分，得到所需的粒子信息。

2　粒子生成算例分析

本節進一步給出2個粒子生成算例，說明本文方法的可行性和優越性。

2.1半圓環

針對內、外徑分別為0.03和0.04 m的半圓環區域，用傳統方法、單元轉化方法和本文方法生成其粒子。傳統方法按先沿徑向等分、再沿周向等分來設計粒子生成算法，得到圖4(a)所示的均勻粒子分布，其中每個粒子具有基本相等的面積。圖4(b)為單元轉化方法的粒子分布，通過先在ABAQUS中用自由網格技術劃分三角形單元、再取單元中心布置粒子而得到。圖4(c)為本文方法的粒子分布。三種方法生成的粒子數目近似相等，分別為4 389、4 560和4 588。

下面對比3種方法生成的粒子對場函數1階導數的近似精度。在半圓環區域上定義場函數

f(x,y)=x+y

(8)

采用SPH中常用的1階導數粒子近似式

(9)

(10)

圖4　半圓環粒子模型Fig.4　Particle model for a half of a ring

表1為導數近似的相對誤差絕對值的平均值。表中內部粒子指沿原區域邊界向內偏置0.001 25 m(1倍光滑函數影響域半徑)的區域內的粒子。由該表可見，與單元轉化方法相比，不論對于內部粒子還是所有粒子，本文方法的近似精度都要高。與傳統方法相比，本文方法在所有粒子的近似精度上略低，在內部粒子的近似精度上仍要高。傳統方法在本算例中給出的粒子面積非常均勻，故在整個區域上近似精度最高，但其對于本算例的粒子劃分思路難以用于其他形狀區域，通用性不如單元轉化方法和本文方法。另外，對于復雜幾何區域，傳統方法要得到具有本算例一樣的均勻性的粒子分布也是很困難的。

表1　函數導數的近似誤差

2.2截卵形子彈

圖5為某截卵形子彈[12]示意圖。該子彈幾何形狀較復雜。傳統方法生成其粒子較繁瑣。本文方法仍容易實現其粒子劃分。圖6(a)為本文方法得到的該子彈軸對稱面(圖5所示的一半)上的粒子分布。可見，它仍能保證粒子分布具有以規則矩形分布為主的特點。

圖5　彈體結構示意圖Fig.5　Schematic diagramfor the projectile’s structure

圖6　不同結構參數下的子彈粒子模型Fig.6　Particle model of the projectile at various structural parameters

另外，基于本文方法還易于實現粒子模型的參數化建模。以本算例的子彈為例，采用Python語言編程，調用ABAQUS的幾何建模函數，建立該子彈的參數化幾何模型。幾何建模程序如下：

CRH=0.66; L=0.060; D1=0.018; D2=0.025; D4=0.005

R3=CRH*D2

m=R3-0.5*D2+0.5*D4

n=sqrt(R3*R3-m*m)

mySketch=mdb.models['Model-1'].ConstrainedSketch(

name='mySketch',sheetSize=0.1)

mySketch.Line(point1=(0.0,0.0),point2=(0.5*D4,0.0))

mySketch.ArcByCenterEnds(center=(0.5*D2-R3,n),

point1=(0.5*D4,0.0),point2=(0.5*D2,n),

direction=COUNTERCLOCKWISE)

mySketch.Line(point1=(0.5*D2,n),point2=(0.5*D2,L))

mySketch.Line(point1=(0.5*D2,L),point2=(0.5*D1,L))

mySketch.Line(point1=(0.5*D1,L),point2=(0.5*D1,n))

mySketch.ArcByCenterEnds(center=(0.0,n),point1=(0.5*D1,n),point2=(0.0,n-0.5*D1),direction=CLOCKWISE)

mySketch.Line(point1=(0.0,n-0.5*D1),point2=(0.0,0.0))

myPart=mdb.models['Model-1'].Part(name='myPart',dimensionality=TWO_D_PLANAR,type=DEFORMABLE_BODY)

myPart.BaseShell(sketch=mySketch)

以上述幾何建模程序代替粒子生成程序中幾何模型導入部分(1.1節對應部分)，得到該子彈粒子模型的參數化建模程序。在該程序中，只需設定5個彈體結構參數，即曲徑比CRH、內徑D1、外徑D2、平頭部分的直徑D4以及彈長L，就可以自動地生成該組參數值對應的彈體粒子模型。圖6(b)～(f)是在圖5基礎上分別改變CRH、D1、D2、D4和L時的彈體粒子模型。彈體的參數化粒子模型在進行彈體結構的優化設計時很有用。

3　高速沖擊模擬應用

3.1彈性球撞擊剛性壁

本算例模擬空心彈性球體垂直撞擊剛性壁面。球體的初始沖擊速度為100 m/s，密度為ρ=1 010 kg/m3，楊氏模量為E=451 MPa，泊松比為υ=0.397 5。采用軸對稱SPH計算。空心球體軸對稱面的幾何形狀為半圓環，設其尺寸與2.1節相同。保持其他條件相同，用圖4中的3種粒子分布計算。圖7為0.5 ms時球體的變形和等效應力分布。可見，在變形方面，采用3種粒子分布的結果差別不大；但在等效應力方面，采用單元轉化法的粒子的結果存在顯著的數值振蕩，而采用傳統方法和本文方法的粒子的結果相近，兩者的等效應力分布更光滑，與ABAQUS中的有限元結果更接近。上述結果說明本文方法的粒子的計算精度與傳統的均布粒子接近，優于單元轉化法的粒子。

圖7　t=0.5 ms時球體的變形和等效應力Fig.7　Deformations and effective stresses of the ball at t=0.5 ms

3.2卵形彈侵徹鋁合金板

本算例模擬Piekutowski等[13]開展的卵形4340鋼彈正侵徹6061鋁合金靶。采用軸對稱SPH計算。在2.2節的參數化建模程序中，設置參數CRH=3、D1=0 mm、D2=12.9 mm、D4=0 mm和L=88.9 mm，取粒子間距為0.645 mm，按本文方法生成彈體軸對稱面的粒子模型。同時，采用單元轉化法生成彈體的粒子模型以作對比。兩種彈體模型的粒子數相近，分別為1 282和1 289。靶板軸對稱面為矩形，離散為矩形均布粒子，粒子數為9 676。圖8給出了兩種粒子模型。計算中，彈體和靶體材料分別采用理想彈塑性模型和Johnson-Cook模型[14]，材料參數見表2。彈體的初始沖擊速度v0=730 m/s。

圖8　卵形鋼彈侵徹鋁合金靶的粒子模型Fig.8　Particle model for ogive-nose projectile penetrating into aluminum alloy target

圖9給出了不同光滑長度下基于兩種粒子模型得到的t=25 μs時的等效應力分布。可見，采用單元轉化法生成的彈體粒子時，等效應力分布隨光滑長度系數的變化而明顯變化，αh為1時，等效應力分布無規律，數值計算不穩定，αh為1.25和1.5時，計算變得穩定，但等效應力存在一定的數值振蕩，振蕩程度隨αh的增大而減少。采用本文方法生成的彈體粒子時，3種光滑長度系數下的等效應力分布均相近且變化較光滑。這說明本文方法生成的粒子在計算穩定性方面更好，可以用更小的光滑長度以減少計算量。當αh為1、1.25和1.5時，采用本文方法的彈體粒子計算得到的彈體剩余速度分別為613、646和645 m/s。可見，光滑長度對剩余速度計算結果有一定影響。αh=1.25和1.5時的計算結果變化較小，兩者與實測剩余速度結果665 m/s[13]吻合較好。

表2　彈體和靶板的材料參數

圖9　彈體和靶板的等效應力Fig.9　Effective stresses of the projectile and the target

4　結　論

本文探討了一種基于ABAQUS二次開發的SPH粒子生成方法。該方法可以自動對任意幾何區域生成粒子，生成過程無需人工干預，效率高。同時，得到的粒子均勻性較好，區域內部粒子嚴格按矩形規則分布。場函數導數的粒子近似計算表明：該方法的粒子的近似精度較好，特別是內部粒子具有很高的近似精度。同時，結合Python語言和ABAQUS的幾何建模函數，易于實現粒子模型的參數化建模。空心球體和彈丸侵徹算例表明：該方法的粒子在高速碰撞過程模擬中應用是可行的，在計算精度和穩定性上優于單元轉化法。另外，該方法的實施簡單，易于拓展到3維。后續研究將開展這方面工作，并進一步改進邊界區域的粒子離散。

[1]MONAGHAN J J.Simulating free surface flows with SPH[J].Journal of Computer Physics,1994,110(2):399-406.

[2]LIU M B,LIU G R,ZONG Z,et al.Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology[J].Computers & Fluids,2003,32(3):305-322.

[3]楊剛,韓旭,龍述堯.應用SPH方法模擬近水面爆炸[J].工程力學,2008,25(4):204-208.

YANG Gang,HAN Xu,LONG Shuyao.Simulation of underwater explosion near air-water surface by SPH method[J].Engineering Mechanics,2008,25(4):204-208.

[4]徐志宏,湯文輝,羅永.SPH算法在高速侵徹問題中的應用[J].國防科技大學學報,2005,27(4):41-44.

XU Zhihong,TANG Wenhui,LUO Yong.Smoothed particle hydrodynamics algorithm applied in penetration problem[J].Journal of National University of Defense Technology,2005,27(4):41-44.

[5]卞梁，王肖鈞，章杰，等.高速碰撞數值計算中的SPH分區算法[J].計算物理,2011,28(2):207-212.

BIAN Liang,WANG Xiaojun,ZHANG Jie,et al.Numerical simulation of hypervelocity impact with subdomian in SPH computation[J].Chinese Journal of Computational Physics,2011,28(2):207-212.

[6]倪國喜,王瑞利,林忠,等.任意區域上的粒子均勻分布方法[J].計算力學學報,2007,24(4):408-413.

NI Guoxi,WANG Ruili,LIN Zhong,et al.Equi-distribution of particles in arbitrary domain[J].Chinese Journal of Computational Mechanics,2007,24(4):408-413.

[7]劉謀斌,常建忠.光滑粒子動力學方法中粒子分布與數值穩定性分析[J].物理學報,2010,59(6):3654-3662.

LIU Moubin,CHANG Jianzhong.Particle distribution and numerical stability in smoothed particle hydrodynamics method[J].Acta Physica Sinica,2010,59(6):3654-3662.

[8]OGER G,DORING M,ALESSANDRINI B,et al.Two-dimensional SPH simulations of wedge water entries[J].Journal of Computational Physics,2006,213(2):803-822.

[9]熱合買提江·依明江，買買提明·艾尼.SPH分析中的齒輪建模方法研究[J].新疆大學學報(自然科學版),2008,25(1):122-126.

RAHMATJAN Imin,MAMTIMIN Geni.Study on the gear modeling in SPH analysis[J].Journal of of Xinjiang University(Natural Science Edition),2008,25(1):122-126.

[10]倪國喜，王瑞利，林忠，等.無網格方法中粒子分布與自適應研究[J].計算物理,2006,23(4):419-424.

NI Guoxi,WANG Ruili,LIN Zhong,et al.Adaptive distribution of particles in a meshfree method[J].Chinese Journal of Computational Physics,2006,23(4):419-424.[11]張強,馬永,李四超.基于Python的ABAQUS二次開發方法與應用[J].艦船電子工程,2011,31(2):131-134.

ZHANG Qiang,MA Yong,LI Sichao.Method and application of second-developed ABAQUS based on Python[J].Ship Elec tronic Engineering,2011,31(2):131-134.

[12]楊剛,梁超,劉平,等.基于三維FE_SPH自適應耦合算法的子彈侵徹混凝土靶跳飛問題模擬[J].工程力學,2013,30(9):276-282.

YANG Gang,LIANG Chao,LIU Ping,et al.Numerical simulation of ricochet problem of projectile penetrating into concrete target based on 3D FE-SPH adaptive coupling method[J].Engineering Mechanics,2013,30(9):276-282.

[13]PIEKUTOWSKI A J,FORRESTAL M J,POORMON K L,et al.Perforation of aluminum plates with ogive-nose steel rods at normal and oblique impacts[J].International Journal of Impact Engineering,1996,18(7/8):877-887.

[14]JOHNSON G R,COOK W H.A constitutive model and data for metals subjected to large strains,high strain rates and high temperatures[C]// Proceedings of the Seventh International Symposium on Ballistics.Hague,1983:541-547.

A SPH particle generation method based on abaqus and its application in high velocity impact calculation

XIAO Yihua,DONG Huanghuang,PING Xuecheng

(School of Mechatronics and Vehicle Engineering,East China Jiaotong University,Nanchang 330013,China)

A new method for SPH particle generation was proposed based on ABAQUS.For this method,a script was written using the programming language Python,and runs with ABAQUS to generate particles.In the script,built-in functions of ABAQUS were firstly used to import a geometry model,and then to search for its boundaries and partition the geometric domain.A particle model was finally established by converting sub-domains obtained with partition into particles.The presented method was suitable for arbitrary domains to produce particles with a basically rectangular distribution.An example of particle generation for a half of ring is used to compare the approximate accuracies of particles generated by the present method and the traditional one.The case of a truncated-ogive-nose projectile was given to illustrate the ablity to establish parameterized particle model by combining the presented method with the parameterized geometric modeling of ABAQUS.Finally,particles generated with the new method were employed to SPH simulation of high velocity impact,and the feasibility of the proposed method in practical application was verified.

SPH; particle generation; ABAQUS; script; high velocity impact

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

2015-05-18修改稿收到日期：2015-08-09

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

O242; O347

A DOI:10.13465/j.cnki.jvs.2016.17.024