SPH 粒子法在高速沖擊連接問題中的應用

孫 標，廖秋慧，何建萍

（上海工程技術大學 材料工程學院，上海 201620）

0 引言

高速沖擊類問題由于速度高，物理過程難以觀測，還具有一定的破壞性。然而高速沖擊卻具有潛在的應用，如在焊接方面，高速沖擊可以實現固態連接。基于材料的狀態，可對焊接技術分類，可分為固態焊接和熔化焊接［1］。熔化焊接是利用金屬熔化原子擴散，經一定時間冷卻后實現冶金結合；固態焊接是低于材料熔點條件下將材料結合在一起。因為熔化焊接會在焊縫處出現金屬間化合物，金屬間化合物屬于硬脆相，其將嚴重影響材料連接后的性能。因此，高速沖擊連接異種材料的固態連接方法具有很好的應用前景。實現高速沖擊連接，關鍵點是要呈現出連接界面處的有規律的界面波［1］，如圖1 所示。

圖1 界面波Fig.1 Wavy morphology observed

數值模擬成為現實中工程問題以及科研問題的重要研究手段，該技術為物理世界的現象理論提供了試驗和檢測，并輔助理解工程中復雜問題，同時可指導工程應用。數值模擬的關鍵點是需要將連續體離散化進行計算。有限元法（FEM）和有限差分法（FDM）是基于網格劃分進而離散，然而傳統Lagrangian 算法對于大變形以及沖擊類問題，由于沖擊使得網格變形，該算法很難處理與原始網格不一致的網格，并造成不連續界面；同時，為保證不連續界面的一致性，通過不同時刻網格重構計算，會增加計算成本，同時網格畸變過大導致計算終止［2］。絕對拉格朗日－歐拉算法（Arbitrary Lagrangian－Eulerian）可模擬沖擊大變形問題，例如：利用ALE算法研究彈塑性入水沖擊問題，水下爆炸沖擊。

無網格方法的關鍵思想是使用一組任意分布的節點或粒子為具有各種可能邊界條件的積分方程或偏微分方程提供準確和穩定的數值解［3］，光滑粒子流體動力學（Smoothed Particle Hydrodynamics，簡稱SPH）方法是無網格法中的一種，該方法最早解決三維天體物理學中行星運動和碰撞問題，后被廣泛應用于流體動力學，沖擊模擬，爆炸焊接，水下爆炸沖擊等問題。本文將以旋轉的子彈作為研究對象，該問題關鍵在于當子彈具有一定的轉速后，簡化為二維問題，難以觀察其旋轉的影響，必須利用三維SPH 粒子法，這也會導致SPH 粒子數增多，進而會增加模擬計算的時間。

綜上，本文將對比網格法FEM 與SPH 粒子法的數值模擬結果，確定合適的模擬方法；為子彈沖擊連接進一步研究做探索，并對比子彈在有無轉速條件下的界面波形貌；最后，將FEM 網格法與SPH 粒子法進行耦合模擬，確定沖擊連接合適的算法模型。

1 理論算法概述

1.1 ALE 算法

ALE 算法是將網格建立在物體上，但是網格既可隨著材料的變形而變形，也可以保持在空間中不動，甚至可以在空間的一個方向保持固定，在另一方向隨物體的運動而運動［4］。在每一個時間間隔的計算中有兩個階段，第一階段計算時，根據質量、動量、能量守恒方程運用拉格朗日算法獲得平衡方程（1）和（2）［5］：

在第二階段中，需要將發生移動的網格重新映射到最初的位置，并計算穿越網格邊界的動量、質量，網格節點根據公式（3）和公式（4）更新速度與位移。

其中，u為網格速度；Fext為外力矢量；Fint為內力矢量；M為質量對角矩陣。

1.2 SPH 光滑粒子法

SPH 的核心是核函數，其可以被理解為一種在一定范圍內其它臨近粒子對研究粒子影響程度的權函數，核函數的近似積分表示任意一點的場函數及其導數，最后在偏微分方程組中運用粒子近似法，得到只與時間相關的常微分方程［6－7］。SPH 算法一般需要兩個步驟，第一步是場函數核近似（kernel approximation），第二步是粒子近似（particle approximation）。

對于任意連續宏觀變量（密度、溫度、壓力等），函數f（x）可用表達式（5）描述［8］：

對場函數導數，其可近似為式（6）：

其中，Ω為求解域；x為待求粒子在空間中的坐標；x′為待求粒子在求解域內的相鄰粒子在空間中的坐標，經過函數核近似計算，上述場函數及其導數均可轉化為場函數的值及核函數的積分表達式；W（x－x′，h）為核函數，滿足歸一化條件和緊支性條件［9－10］。

核函數W的選取決定粒子緊支域大小以及核近似和粒子近似的精度［3］。常用核函數包括鐘形、高斯型、分段三次樣條核函數。

在SPH 算法中，粒子近似法作為SPH 的第二個關鍵步驟，在核近似過程中將場函數轉化為近似連續性積分方程后，通過粒子近似法，得到相鄰粒子的疊加求和的離散形式。粒子近似法如圖2 所示。

圖2 粒子近似法Fig.2 Method of particle approximation

若將積分近似式中表示粒子j處的無窮小體積dx′用離子體積ΔVj代替，則粒子質量mj表達式（7）：

其中，ρj為粒子j的密度。

第i個粒子的核近似函數f（xj）最終轉化為粒子近似［3］，式（8）：

2 數值模型建立

2.1 幾何模型及材料參數

采用圓柱型子彈沖擊平面的基板建立數值模擬模型。子彈材料為銅，基板材料為低碳鋼材料，性能參數見表1。模型的尺寸分別是φ9×6 mm，10×10×5 mm，網格單元尺寸為0.02 mm，如圖3 所示。相比網格法，SPH 粒子法，沒有單元信息，取而代之的是節點信息。SPH 粒子建模過程是將建立好的有限元網格模型，在軟件LS－dyna 前處理軟件中SPH generation 功能將單元信息生成節點信息即可。

表1 模擬材料銅與低碳鋼的性能［11］Tab.1 Property of material copper and low carbon steel［11］

圖3 數值模型Fig.3 Numerical model

2.2 數值分析模型

由于子彈沖擊基板屬于大變形和高應變率問題，Johnson－Cook 方程的材料行為適用于強沖擊載荷，式（9）。因此，選擇了Johnson－Cook 模型作為模擬中子彈和基板的材料模型。

其中，εp是材料等效塑性應變；ε′p為材料實際應變率；為參考應變率；為無量綱的等效應變速 率；為無量綱的等效應變速率；T?m ＝；T是無量綱溫度；Tm是材料熔點；Tr是室溫，A、B、C、n、m為材料參數。

材料模型參數見表2。

表2 銅與低碳鋼J－C 本構模型參數Tab.2 J－C Constitutive model parameters of Copper and lowcarbon steel

當利用JC 材料模型時，必須有狀態方程才能進行模擬計算。材料的狀態方程與密度、壓力及一些熱力學參數有關，能夠反應材料的體積特性。Grüneisen 方程適用于固體中的波的傳播，能夠與波陣面中的沖擊跳躍存在直接聯系，且方程中參數也可以通過實驗確定。因此，子彈和基板材料模型的狀態方程為Grüneisen 方程，其表達式（10）為：

式中，ρ0為參考密度；C1為uS－ uP曲線的截距；γ0、μ、α為Grüneisen方程系數；S1、S2、S3為uS－uP曲線斜率的系數；，V為相對體積；α為對γ0的一階體積修正；E為材料內能。

子彈和基板的狀態方程參數見表3。

表3 銅與低碳鋼的狀態方程參數Tab.3 Parameters of the equation f state for copper and low carbon steel

沖擊連接可成功的速度范圍一般在150～1 500 m／s［6］。因此，選取模擬的初始條件為：沖擊速度為150 m／s，角速度1 177 rad／s；為更好對比有無旋轉影響的影響，將沖擊角度設置為0°，目的是有利于子彈施加旋轉速度。

3 結果與討論

3.1 ALE 網格法與SPH 粒子法界面處波形

分別對SPH 數值模型以及ALE 數值模型計算，獲得如圖4 所示的界面波。

圖4 SPH 與ALE 界面波形圖Fig.4 SPH and ALE interface waveform

并對比高速沖擊實驗的界面波圖1 與圖4 模擬結果可知，采用SPH 粒子法可模擬沖擊界面處波浪狀的界面波形如圖4（a）；而采取ALE 網格法其沖擊界面處出現較大的扭曲和網格畸變如圖4（b）。因此，對于沖擊類問題，SPH 粒子法呈現沖擊過程的細節比ALE 方法合適。

在數值計算過程中不同算法會出現能量損失，算法是否合適，能量守恒可作為判斷標準：工程上一般在10%以內的能量波動是可接受的。圖5（a）和（b）分別為SPH 粒子法與ALE 網格法在Ls－dyna求解器中求得碰撞過程中的總能量變化；根據圖5（a）可計算得SPH 粒子法總能量損失為9.8%，該值在可接受的范圍內；圖5（b）ALE 網格法能量損失1.3%，相比之下ALE 網格法能量損失較少。這是SPH 粒子法本身特點決定的，因為每個粒子是單獨個體，碰撞后期部分粒子被沖擊散落而失效。綜合SPH 粒子法具有2 個特點：可呈現出界面處的界面波形以及能量損失在可接受范圍內；因此，采用SPH 粒子法模擬沖擊類問題更合適。

圖5 碰撞過程的能量變化Fig.5 Energy change curve during collision

3.2 有、無轉速條件下的界面波

SPH 算法明顯的呈現出高速沖擊連接的界面波形，為探索轉速對界面處形成波形的影響，數值模擬了在沖擊條件一致時模擬有、無轉速條件下界面波差異如圖6 所示。

圖6 高速沖擊條件下界面波形成過程Fig.6 Formation process of interface wave under high speed impact

圖6 中可觀測到在有轉速沖擊條件下，界面處波幅較小，其原因是沖擊過程中的賦予子彈旋轉產生的切應力所致，具體地，主要受到τxy（X面上沿Y方向的切應力即圖7 中xy－應力）的影響，圖7（a）為在無轉速的條件下產生的τxy切應力值為6 MPa，圖7（b）為有轉速的條件下界面處的τxy切應力值達到25 MPa。所以，利用SPH 數值模擬方法，研究旋轉因素對界面波形貌影響具有可行性。

圖7 xy_應力對比圖Fig.7 Xy_stress contrast chart

3.3 ALE 與SPH 耦合

實際上，當三維問題的結構復雜后，SPH 算法建模所需的粒子數會增加，從而增加了計算的時間成本；為減少計算成本，通常是進行并行計算處理，簡化模型，減少粒子數。并行計算可以減少需要計算設備的支持，不是所有模型都可以被簡化，例如本案例，當有轉速時候，很難獲得轉速對界面的影響。如果減少粒子數，無法呈現出界面波形，減少粒子數會增加計算誤差，粒子數與計算精度關系如圖8 所示。

圖8 粒子數與計算精度的關系Fig.8 The relationship between the number of particles and calculation accuracy

因此，本案例利用Ls－dyna 軟件中的耦合關鍵字?CONTAC _TIED_NODES_TO_SURFACE，將SPH 粒子單元與ALE 網格單元耦合；建模時對碰撞的非觀察區域進行網格劃分，觀察區域進行粒子法近似處理，如圖9 所示。最終模擬計算出ALE 與SPH 耦合的結果如圖10 所示。

圖9 ALE 與SPH 耦合模型Fig.9 Model of ALE and SPH coupling

圖10 ALE 與SPH 耦合結果Fig.10 Result of ALE and SPH coupling

圖10（a）中可以發現，耦合后界面處的波形沒有任何影響，圖10（b）為x面上y方向上的剪應力處于同一色塊，應力傳遞具有一致性，意味著耦合的成功。即對于模型復雜或難以簡化的模型，可采用ALE 與SPH 方法耦合的方法建模計算模型。

通過對比3 種不同算法模型單位時間計算的粒子數見表4，可確定當粒子數增加后，ALE 與SPH耦合的方式可減少直接利用SPH 粒子法的時間成本。這對尺寸較大，結構復雜的三維模型作沖擊連接數值計算具有一定的借鑒作用。

表4 3 種不同算法模型單位時間計算的粒子數Tab.4 Number of particles calculated per unit time by three different algorithm models

4 結束語

本文對比了ALE 法和SPH 粒子法旋轉子彈沖擊基板的案例，分析了SPH 粒子法更適合應用于沖擊類問題，并針對SPH 粒子法存在的粒子數過多導致的計算成本增加，進行ALE 與SPH 粒子法耦合優化。綜上，可獲得如下結論：

（1）高速沖擊連接問題，相比于ALE 算法，SPH粒子法能夠很好的呈現連接界面的波形；同時，SPH粒子法也可呈現有轉速與無轉速時界面波波形差異，可研究子彈轉速對沖擊界面波形貌的影響；

（2）隨著模型尺寸變大導致使用SPH 粒子法粒子數增加，耗時過長；ALE 算法與SPH 算法耦合的方式建模，可減少直接利用SPH 粒子法的計算時間。