高速沖擊問題的SPH粒子接觸算法三維數值計算

張志春，卞 強，展文豪，劉桐宇

(中國航天員科研訓練中心人因工程國防科技重點實驗室， 北京 100094)

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

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

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

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

1 考慮材料強度的SPH方程

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

連續性方程：

(1)

動量方程：

(2)

能量方程：

(3)

運動方程：

(4)

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

材料強度體現在應力張量σαβ，由各向同性壓力p和黏性剪切應力ταβ組成

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

(5)

其中偏應力率采用Jaumann比率的形式表示為

(6)

(7)

(8)

其SPH離散格式分別表示為

(9)

(10)

2 SPH粒子接觸算法

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

(11)

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

▽xiW(rij)

(12)

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

(13)

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

3 球頭鋼彈斜沖擊鋼板

本節采用SPH粒子接觸算法對球頭鋼彈斜沖擊鋼板進行了三維數值計算。子彈和靶板的材料參數為：密度ρ=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°。為了進行對比驗證，采用LS-DYNA970對相同問題進行了計算，子彈和靶板均采用六面體單元離散，單元數分別為1 824和2 700。計算時間均為100×10-6s。

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

4 平頭鋼彈正沖擊鋼板

本節采用SPH粒子接觸算法對圓柱形鋼彈正沖擊鋼板發生沖塞破壞的過程進行了三維數值計算，計算模型的尺寸及離散方式如圖5所示。子彈和靶板均由SPH粒子離散，粒子總數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區域內粒子間距10-3m，為減小計算量，粒子間距從中心正碰區向邊緣逐漸增大。子彈與靶板初始間距2×10-3m，靶板四周采用固支邊界條件，計算總時間160×10-6s。采用Monaghan[9]提出的人工應力法來消除拉伸不穩定性造成的數值斷裂，采用強洪夫[13-14]提出的完全變光滑長度法來消除光變光滑長度帶來的影響。

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

(14)

(15)

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

(16)

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

(17)

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

(18)

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

表1 子彈材料參數

表2 靶板材料參數

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

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

表3 數值模擬結果與實驗結果

5 結論

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

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