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

孫 標，廖秋慧，何建萍

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

0 引言

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

圖1 界面波Fig.1 Wavy morphology observed

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

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

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

1 理論算法概述

1.1 ALE 算法

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

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

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

1.2 SPH 光滑粒子法

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

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

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

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

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

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

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

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

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

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

2 數(shù)值模型建立

2.1 幾何模型及材料參數(shù)

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

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

圖3 數(shù)值模型Fig.3 Numerical model

2.2 數(shù)值分析模型

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

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

材料模型參數(shù)見表2。

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

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

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

子彈和基板的狀態(tài)方程參數(shù)見表3。

表3 銅與低碳鋼的狀態(tài)方程參數(shù)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 網(wǎng)格法與SPH 粒子法界面處波形

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

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

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

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

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

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

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

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

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

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

3.3 ALE 與SPH 耦合

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

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

因此，本案例利用Ls－dyna 軟件中的耦合關鍵字?CONTAC _TIED_NODES_TO_SURFACE，將SPH 粒子單元與ALE 網(wǎng)格單元耦合；建模時對碰撞的非觀察區(qū)域進行網(wǎng)格劃分，觀察區(qū)域進行粒子法近似處理，如圖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）中可以發(fā)現(xiàn)，耦合后界面處的波形沒有任何影響，圖10（b）為x面上y方向上的剪應力處于同一色塊，應力傳遞具有一致性，意味著耦合的成功。即對于模型復雜或難以簡化的模型，可采用ALE 與SPH 方法耦合的方法建模計算模型。

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

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

4 結束語

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

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

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