一種基于顆粒接觸的二維無網(wǎng)格方法 及其在高速?zèng)_擊模擬中的應(yīng)用*

馮 春,李世海,劉曉宇

(中國科學(xué)院力學(xué)研究所，北京 100190)

基于連續(xù)介質(zhì)力學(xué)的離散元方法(continumm-based discrete element method, CDEM)通過有限元與離散元的有機(jī)結(jié)合，在分析材料的變形斷裂及斷裂后的運(yùn)動(dòng)方面具有一定的優(yōu)勢[1-2],但該方法在模擬高速?zèng)_擊等大變形問題時(shí)，會(huì)出現(xiàn)網(wǎng)格畸變導(dǎo)致的數(shù)值系統(tǒng)能量發(fā)散等現(xiàn)象。為解決上述問題，本文中提出了一種基于顆粒接觸的二維無網(wǎng)格方法(particle contact-based meshfree method, PCMM)。

由于網(wǎng)格畸變導(dǎo)致的數(shù)值系統(tǒng)發(fā)散是所有拉格朗日架構(gòu)下，基于網(wǎng)格的數(shù)值計(jì)算方法所面臨的共同問題。目前的解決方法主要包括生死單元法及自適應(yīng)網(wǎng)格等兩大類。生死單元法是將嚴(yán)重畸變的單元進(jìn)行鈍化處理，從而保證了每一計(jì)算時(shí)步的網(wǎng)格質(zhì)量，該方法雖然簡單易行，但單元挖空后導(dǎo)致的應(yīng)力場不連續(xù)問題值得探討。自適應(yīng)網(wǎng)格[3-4]根據(jù)數(shù)值系統(tǒng)內(nèi)局部區(qū)域的網(wǎng)格畸變程度，對(duì)該區(qū)域網(wǎng)格進(jìn)行調(diào)整或加密，從而保證了網(wǎng)格質(zhì)量，但網(wǎng)格重劃分的耗時(shí)及新生單元各類場信息的插值問題仍然值得關(guān)注。

無網(wǎng)格方法是近年來的研究熱點(diǎn)，該方法的核心是通過空間域內(nèi)一系列離散節(jié)點(diǎn)構(gòu)建連續(xù)介質(zhì)插值函數(shù)，由于擺脫了網(wǎng)格的限制，因此可以解決因網(wǎng)格畸變導(dǎo)致的系統(tǒng)能量發(fā)散等問題。目前發(fā)展的無網(wǎng)格方法有幾十種，這些方法的主要區(qū)別在于使用不同的插值函數(shù)對(duì)微分方程進(jìn)行等效，如光滑粒子法(SPH)、物質(zhì)點(diǎn)法(MPM)等[5]。

SPH最早用于天體物理方面的研究[6-7]，核心思想是將空間場量的偏導(dǎo)數(shù)轉(zhuǎn)化為對(duì)設(shè)定核函數(shù)的偏導(dǎo)數(shù)，從而簡化了計(jì)算。近年來，SPH及其擴(kuò)展形式在自由表面流動(dòng)、高速?zèng)_擊、材料斷裂等方面得到了廣泛應(yīng)用[8-10]。盡管如此，該方法在邊界處的處理方式、張力不穩(wěn)定的改進(jìn)方法及支持域內(nèi)粒子的快速搜索[11-13]等方面仍然需要進(jìn)行深入的研究。

MPM的主要思想是空間離散質(zhì)點(diǎn)系在背景網(wǎng)格中的流動(dòng)計(jì)算。該方法將物質(zhì)點(diǎn)信息映射到背景網(wǎng)格并進(jìn)行動(dòng)量方程求解，獲取網(wǎng)格點(diǎn)的速度、位移等信息，而后將網(wǎng)格點(diǎn)運(yùn)動(dòng)量映射回物質(zhì)點(diǎn)，更新物質(zhì)點(diǎn)的空間坐標(biāo)[14-15]。該方法發(fā)揮了拉格朗日構(gòu)型(物質(zhì)點(diǎn))及歐拉構(gòu)型(背景網(wǎng)格)的雙重優(yōu)勢，但在背景網(wǎng)格尺寸的選取原則及質(zhì)點(diǎn)穿過背景網(wǎng)格產(chǎn)生的數(shù)值振蕩等方面仍然需要仔細(xì)研究。

本文中在CDEM方法的基礎(chǔ)上，結(jié)合顆粒離散元中接觸對(duì)的概念，提出了基于顆粒接觸的二維無網(wǎng)格方法(PCMM)，并基于VC++編制了相應(yīng)的計(jì)算程序。

1 PCMM方法的基本原理

1.1 核心思想

采用顯式解法計(jì)算高速?zèng)_擊問題，當(dāng)單元出現(xiàn)畸變后會(huì)嚴(yán)重影響計(jì)算時(shí)步，進(jìn)而導(dǎo)致系統(tǒng)能量的發(fā)散。顆粒離散元中顆粒間的接觸信息及拓?fù)潢P(guān)系復(fù)雜而有序，基于顆粒間的接觸關(guān)系可構(gòu)建一系列的連續(xù)介質(zhì)單元，通過在單元中引入宏觀本構(gòu)關(guān)系即可表述連續(xù)介質(zhì)在沖擊載荷作用下的動(dòng)態(tài)響應(yīng)。當(dāng)上述連續(xù)介質(zhì)單元的變形逐漸增大，顆粒間的接觸關(guān)系也將逐漸發(fā)生改變；當(dāng)單元出現(xiàn)嚴(yán)重畸變(如壓扁、拉長、扭曲等)時(shí)，顆粒間的原有接觸關(guān)系被打破，基于原有接觸關(guān)系構(gòu)建的連續(xù)介質(zhì)單元(畸變單元)自動(dòng)刪除；同時(shí)，基于當(dāng)前顆粒接觸關(guān)系的單元將自動(dòng)創(chuàng)建。通過舊單元的刪除及新單元的創(chuàng)建，即可解決高速?zèng)_擊計(jì)算中因單元畸變導(dǎo)致的數(shù)值系統(tǒng)發(fā)散等問題。基于以上想法，文本中提出了PCMM方法，該方法的基本流程如圖1所示。

圖1 PCMM方法的基本流程Fig.1 Basic steps of PCMM

1.2 三角形單元的創(chuàng)建

三角形單元?jiǎng)?chuàng)建的必備條件包括：

(1)構(gòu)成三角形單元的3個(gè)顆粒必須彼此接觸：

dij≤Ri+Rj+δi,j=1,2,3;i≠j

(1)

(2)所構(gòu)建三角形單元的任何一個(gè)內(nèi)角角度應(yīng)在30°～150°之間(保證單元不會(huì)太奇異)：

30°≤θi≤150°i=1,2,3

(2)

(3)構(gòu)成的三角形單元3條邊的任意一條邊不應(yīng)該小于3個(gè)顆粒平均半徑的0.5倍：

(3)

式中：dij為顆粒i與j之間的距離，Ri、Rj為顆粒i與j的半徑，δ為顆粒i與j之間的接觸容差，θi為三角形單元的某個(gè)內(nèi)角。在隨機(jī)排列的顆粒體系中，為了構(gòu)建穩(wěn)定的連續(xù)介質(zhì)單元系統(tǒng)，接觸容差δ必須足夠大以防止空隙的發(fā)生。

為提升顆粒接觸對(duì)的檢索效率，本文中采用了二維空間分割法(staticcell)，計(jì)算復(fù)雜度為O(N)。該方法將二維空間域按照顆粒的最大半徑進(jìn)行正交格子劃分，并將每個(gè)顆粒按照其質(zhì)心坐標(biāo)映射到對(duì)應(yīng)的格子內(nèi)(見圖2)；接觸尋找時(shí)循環(huán)所有顆粒，根據(jù)當(dāng)前顆粒的質(zhì)心坐標(biāo)確定顆粒所在格子的位置(圖2中顆粒A的位置為x=3,y=3),在顆粒A周圍的9個(gè)格子內(nèi)(包含自身所在格子)尋找潛在接觸顆粒。

如潛在顆粒i到顆粒A的距離小于設(shè)定容差(見圖3)，將顆粒i增加至顆粒A的接觸鏈表中。完成顆粒A的鄰居檢索后，計(jì)算鄰居顆粒到顆粒A的方位角，并將鄰居顆粒按照方位角進(jìn)行排序。排序過程中，如果相鄰2個(gè)鄰居顆粒之間的方位角之差小于30°，將距離顆粒A較遠(yuǎn)的那個(gè)接觸顆粒從顆粒A的接觸鏈表中刪除(防止產(chǎn)生低質(zhì)量單元)。完成了顆粒A鄰居顆粒的排序及局部刪除后，重新循環(huán)顆粒A的接觸鏈表，判斷顆粒A，第i個(gè)及i+ 1個(gè)鄰居顆粒組成的三角形是否滿足式(1)～(3)所示條件，如果滿足條件且該三角形單元不在顆粒A的單元鏈表中(新單元)，創(chuàng)建三角形單元并增加至顆粒A的單元鏈表中。

圖2 用于鄰居搜索的二維空間盒Fig.2 2D bin array for neighbour searching

圖3 基于方位角的鄰居排序Fig.3 Neighbour sequence according to position angle

1.3 三角形單元的刪除

經(jīng)過一段時(shí)間的迭代計(jì)算后，單元將出現(xiàn)大變形，如果不進(jìn)行任何修正繼續(xù)計(jì)算，將出現(xiàn)由于單元畸變導(dǎo)致的系統(tǒng)發(fā)散等問題。考慮到所有場量(如密度、速度、溫度、應(yīng)力、應(yīng)變等)均位于顆粒上，且單元畸變后某一個(gè)方向尺寸非常小(單元體積可忽略)，將畸變后的單元?jiǎng)h除不會(huì)引起系統(tǒng)總質(zhì)量及總能量的改變。因此，提出如下的三角形單元?jiǎng)h除原則(滿足以下3個(gè)條件的任何一個(gè)，該三角形單元即被刪除)：

(1)組成該三角形單元的3個(gè)顆粒不再彼此接觸；

(2)該三角形單元的任何一個(gè)內(nèi)角小于30°或者大于150°；

(3)該三角形單元的任何一條邊長小于3個(gè)顆粒平均半徑的0.5倍。

1.4 單元內(nèi)力及節(jié)點(diǎn)運(yùn)動(dòng)的計(jì)算

三角形單元?jiǎng)?chuàng)建完畢后，即可采用經(jīng)典連續(xù)介質(zhì)理論進(jìn)行單元變形力的求解。為了模擬高速?zèng)_擊問題，引入流體彈塑性模型，并采用增量的計(jì)算方式，基本計(jì)算流程為：

(4)

(6)根據(jù)物態(tài)方程計(jì)算壓力(體應(yīng)力)p，本文中采用Mie-Grüeisen物態(tài)方程。

(8)計(jì)算數(shù)值阻尼應(yīng)力：采用人工q阻尼消除強(qiáng)間斷帶來的計(jì)算穩(wěn)定性問題

(5)

(10)根據(jù)彈性力學(xué)中斜截面公式，將應(yīng)力轉(zhuǎn)換為節(jié)點(diǎn)力

(6)

式中：N表示第N個(gè)節(jié)點(diǎn)，k為節(jié)點(diǎn)N所在棱的序號(hào)(三角形中每個(gè)節(jié)點(diǎn)隸屬于2個(gè)棱)。

(11)計(jì)算顆粒所受合外力Fi，根據(jù)牛頓第二運(yùn)動(dòng)定律計(jì)算顆粒的加速度、速度、位移：

(7)

(12)回到步驟(1)進(jìn)行下一時(shí)步速度梯度的求解。

2 流體彈塑性模型的引入

流體彈塑性模型是描述材料在爆炸沖擊作用下動(dòng)態(tài)響應(yīng)的經(jīng)典模型，該模型中體應(yīng)力由物態(tài)方程控制，偏應(yīng)力由胡克定律控制并受塑性準(zhǔn)則修正。本文中采用Mie-Grüeisen物態(tài)方程表征材料的流體特征，采用胡克定律表征材料的彈性特征，采用Johnson-Cook模型表征材料的塑性特征。

2.1 Mie-Grüeisen物態(tài)方程的引入

以沖擊Hugoniot態(tài)為參考態(tài)的Grüeisen物態(tài)方程可表述為：

(8)

式中：μ=ρ/ρ0-1，ρ0為初始密度，ρ為當(dāng)前密度(ρ=(V0/V)ρ0)，c0、λ為高壓下沖擊波速度與粒子速度擬合函數(shù)的系數(shù)(D=c0+λu，D沖擊波速度，u粒子速度)，γ0、a為量綱一材料參數(shù)。

上述狀態(tài)方程中還有一個(gè)變量E，它是系統(tǒng)內(nèi)能，可用下式求得，

(9)

2.2 Johnson-Cook模型的引入

Johnson-Cook模型(JC模型)是Johnson和Cook于1983年提出的用于高應(yīng)變率和高溫下材料大變形的本構(gòu)模型：

(10)

當(dāng)達(dá)到上述屈服應(yīng)力后，需對(duì)6個(gè)偏應(yīng)力分量進(jìn)行修正。由于JC模型只定義了屈服條件，沒有定義流動(dòng)法則，因此實(shí)際計(jì)算時(shí)，當(dāng)應(yīng)力超過上述屈服強(qiáng)度，讓其等于此屈服強(qiáng)度。

對(duì)于當(dāng)前溫度T的處理有2種方式：一種方式認(rèn)為數(shù)值計(jì)算過程中當(dāng)前溫度T保持不變，這種方式在參數(shù)輸入時(shí)設(shè)定初值即可。還有一種認(rèn)為當(dāng)前溫度T在數(shù)值計(jì)算過程中會(huì)因?yàn)樗苄怨Χ淖儭?/p>

對(duì)于第2種方式，需根據(jù)塑性功與溫度的等效關(guān)系進(jìn)行當(dāng)前溫度T的求解，每時(shí)步單位體積內(nèi)的塑性功增量可由下式得出，

(11)

每時(shí)步的溫度增量可由下式得出，

(12)

式中：c為材料的比熱，α為功熱轉(zhuǎn)化率。當(dāng)前時(shí)步的溫度值可表示為T=Troom+∑ΔT。

3 高速?zèng)_擊算例

3.1 彈性桿撞擊

根據(jù)文獻(xiàn)[9]，建立直徑1 mm、高20 mm的桿件，并用2 211個(gè)規(guī)則排布的顆粒進(jìn)行離散，顆粒半徑為0.05 mm(見圖4)。彈性桿密度7 850 kg/m3，彈性模量206 GPa，泊松比0.3。彈性桿以50 m/s的速度向基板撞擊，觀察桿件正中間測點(diǎn)軸向速度隨時(shí)間的變化(見圖5)。由圖5可得，數(shù)值解與理論解基本一致，由此表明了PCMM方法在求解彈性動(dòng)力問題方面的準(zhǔn)確性。

圖4 彈性桿撞擊數(shù)值模型Fig.4 Numerical model of elastic bar impacting rigid wall

圖5 桿件中部測點(diǎn)軸向速度曲線Fig.5 Axial velocity curve on the middle point of bar

3.2 泰勒桿

根據(jù)文獻(xiàn)[16]，建立直徑10 mm、高70 mm的OFHC銅桿，撞擊速度165 m/s。桿件采用17 901個(gè)規(guī)則排布的顆粒進(jìn)行離散，顆粒半徑為0.1 mm，并采用流體彈塑性模型進(jìn)行撞擊分析。取A=89.63 MPa，B=45MPa，n=0.31，C=0.025，m=1.09；Mie-Grüeisen物態(tài)方程中取ρ=8 940 kg/m3，c0=3 940m/s，λ=1.49，γ=2.02，a=1.5；胡克定律中剪切模量G取60.74 GPa。不同時(shí)刻桿件的變形如圖6所示，最終狀態(tài)下桿件各位置的半徑如圖7所示。由圖6～7可得，PCMM計(jì)算獲得的桿件運(yùn)動(dòng)過程及最終形態(tài)與文獻(xiàn)[16]的實(shí)驗(yàn)結(jié)果基本吻合，表明了PCMM方法在模擬高速?zèng)_擊破壞效應(yīng)時(shí)的準(zhǔn)確性及合理性。

圖6 不同時(shí)刻桿件的變形Fig.6 Deformation of Taylor bar at different times

圖7 最終狀態(tài)下桿件各位置的半徑Fig.7 Bar radius at different position in final state

3.3 碎片云

碎片云數(shù)值計(jì)算模型中，彈丸直徑2 mm，沖擊速度5 km/s，靶板厚1 mm，靶板高16 mm。彈丸及靶板采用8 491個(gè)隨機(jī)排布的顆粒進(jìn)行離散，顆粒半徑0.015～0.040 mm。彈丸及靶板的材料均為鋁，Johnson-Cook模型中，取A=324 MPa，B=114 MPa，n=0.42，C=0.016，m=1.34，初始溫度298 K，融化溫度877 K，比熱容875 J/(K·kg)，功熱轉(zhuǎn)化率0.9；Mie-Grüeisen物態(tài)方程中取ρ=2 703 kg/m3，c0=5 350m/s，λ=1.34，γ=1.97，a=1.5；胡克定律中剪切模量G取26.7 GPa。彈丸撞擊靶板后2.68 μs時(shí)刻碎片云的空間形態(tài)如圖8所示，該結(jié)果與文獻(xiàn)[17-19]給出的數(shù)值及實(shí)驗(yàn)結(jié)果基本一致，表明了PCMM方法在模擬超高速碰撞問題中的準(zhǔn)確性及合理性。

圖8 2.68 μs 時(shí)碎片云的形態(tài)Fig.8 Shape of debris clouds at 2.68 μs

圖9 不同時(shí)刻子彈、靶板系統(tǒng)的水平速度云圖Fig.9 Horizontal velocity contour of bullet target system at different times

3.4 子彈入射靶板

子彈入射靶板模型中，子彈材料為鎢，入射速度1 140 m/s，長2.84 cm，直徑8 mm；靶板材料為鋁，板高20 cm，寬5 cm。對(duì)子彈及靶板采用11 428個(gè)隨機(jī)排布的顆粒進(jìn)行離散，顆粒半徑0.3～0.7 mm。靶板及子彈均采用流體彈塑性模型，靶板的材料參數(shù)與3.3節(jié)碎片云算例的參數(shù)一致。子彈的材料參數(shù)如下：Johnson-Cook模型中，取A=1.2 GPa，B=1.03 GPa，n=0.019，C=0.034，m=0.4，初始溫度298 K，融化溫度1 723 K，比熱容134 J/(K·kg)，功熱轉(zhuǎn)化率0.9；Mie-Grüeisen物態(tài)方程中取ρ=19.22 g/cm3，c0=4.02km/s，λ=1.24，γ=1.67，a=1.3；胡克定律中剪切模量G取134.8GPa。子彈撞擊靶板后10、32及65μs時(shí),子彈、靶板系統(tǒng)水平速度云圖(按照顆粒組成的連續(xù)介質(zhì)單元顯示)如圖9所示。圖9給出的子彈對(duì)靶板的侵徹作用與實(shí)際較為接近，證明了PCMM方法的合理性。

4 結(jié) 論

PCMM方法通過顆粒的接觸關(guān)系實(shí)現(xiàn)了連續(xù)介質(zhì)單元的快速創(chuàng)建及刪除，解決了傳統(tǒng)CDEM方法中因網(wǎng)格畸變導(dǎo)致的系統(tǒng)能量發(fā)散等問題。在PCMM中引入流體彈塑性模型后，可精確計(jì)算高速?zèng)_擊下材料的動(dòng)態(tài)響應(yīng)及變形破壞特征。彈性桿撞擊、泰勒桿、碎片云及子彈入射靶板等算例，從不同的角度證明了PCMM方法計(jì)算高速?zèng)_擊問題的精確性及合理性。

當(dāng)然，PCMM方法也存在許多值得進(jìn)一步探討的地方，如顆粒質(zhì)量的等效方式、單元?jiǎng)h除與重建的條件、單元重建中力學(xué)信息的繼承與更新、顆粒接觸的快速檢索方法、三維PCMM的實(shí)現(xiàn)等。

[1] Li S H, Zhao M H, Wang Y N, et al.A continuum-based discrete element method for continuous deformation and failure process[C]∥WCCM VI in Conjunction with APCOM’04.Beijing, 2004:77.

[2] Wang Y N, Zhao M H, Li S H, et al.Stochastic structural model of rock and soil aggregates by continumm-based discrete element method[J].Science in China Series E-Engineering & Materials Science, 2005,48(suppl):95-106.

[3] Littlefield D L.The use ofr-adaptivity with local, intermittent remesh for modeling hypervelocity impact and penetration[J].International Journal of Impact Engineering, 2001,26(1):433-442.

[5] 張雄,劉巖,馬上.無網(wǎng)格法的理論及應(yīng)用[J].力學(xué)進(jìn)展,2009,39(1):1-36.Zhang Xiong, Liu Yan, Ma Shang.Meshfree methods and their applications[J].Advances in Mechanics, 2009,39(1):1-36.

[6] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].The Astronomical Journal, 1977,82:1013-1024.

[7] Gingold R A, Monaghan J J.Smoothed particle hydrodynamics-theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society, 1977,181:375-389.

[8] Antuono M, Colagrossi A, Marrone S, et al.Free-surface flows solved by means of SPH schemes with numerical diffusive terms[J].Computer Physics Communications, 2010,181(3):532-549.

[9] 肖毅華,胡德安,韓旭,等.一種自適應(yīng)軸對(duì)稱FEM-SPH耦合算法及其在高速?zèng)_擊模擬中的應(yīng)用[J].爆炸與沖擊,2012,32(4):384-392.Xiao Yi-hua, Hu De-an, Han Xu, et al.An adaptive axisymmetric FEM-SPH coupling algorithm and its application to high velocity impact simulation[J].Explosion and Shock Waves, 2012,32(4):384-392.

[10] Ma G W, Wang X J, Ren F.Numerical simulation of compressive failure of heterogeneous rock-like materials using SPH method[J].International Journal of Rock Mechanics and Mining Sciences, 2011,48(3):353-363.

[11] Chen J K, Beraun J E, Carney T C.A corrective smoothed particle method for boundary value problems in heat conduction[J].International Journal for Numerical Methods in Engineering, 1999,46(2):231-252.

[12] Sigalotti L D G, López H.Adaptive kernel estimation and SPH tensile instability[J].Computers & Mathematics with Applications, 2008,55(1):23-50.

[13] Ma S, Zhang X, Qiu X M.Comparison study of MPM and SPH in modeling hypervelocity impact problems[J].International Journal of Impact Engineering, 2009,36(2):272-282.

[14] Zhang X, Sze K Y, Ma S.An explicit material point finite element method for hyper-velocity impact[J].International Journal for Numerical Methods in Engineering, 2006,66(4):689-706.

[15] Ambati R, Pan X F, Yuan H, et al.Application of material point methods for cutting process simulations[J].Computational Materials Science, 2012,57:102-110.

[16] 呂劍,何穎波,田常津,等.泰勒桿實(shí)驗(yàn)對(duì)材料動(dòng)態(tài)本構(gòu)參數(shù)的確認(rèn)和優(yōu)化確定[J].爆炸與沖擊,2006,26(4):339-344.Lü Jian, He Ying-bo, Tian Chang-jing, et al.Validation and optimization of dynamic constitutive model constants with Taylor test[J].Explosion and Shock Waves, 2006,26(4):339-344.

[17] Beissel S R, Gerlach C A, Johnson G R.A quantitative analysis of computed hypervelocity debris clouds[J].International Journal of Impact Engineering, 2008,35(12):1410-1418.

[18] Huang J, Ma Z, Ren L S, et al.A new engineering model of debris cloud produced by hypervelocity impact[J].International Journal of Impact Engineering, 2013,56:32-39.

[19] Chi R Q, Pang B J, Guan G S, et al.Analysis of debris clouds produced by impact of aluminum spheres with aluminum sheets[J].International Journal of Impact Engineering, 2008,35(12):1465-1472.