孫 標(biāo)，廖秋慧，何建萍

（上海工程技術(shù)大學(xué) 材料工程學(xué)院，上海 201620）

0 引言

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

圖1 界面波Fig.1 Wavy morphology observed

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

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

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

1 理論算法概述

1.1 ALE 算法

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

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

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

1.2 SPH 光滑粒子法

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

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

對(duì)場(chǎng)函數(shù)導(dǎo)數(shù)，其可近似為式（6）：

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

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

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

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

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

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

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

2 數(shù)值模型建立

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

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

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

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

2.2 數(shù)值分析模型

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

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

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

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

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

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

子彈和基板的狀態(tài)方程參數(shù)見(jiàn)表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；為更好對(duì)比有無(wú)旋轉(zhuǎn)影響的影響，將沖擊角度設(shè)置為0°，目的是有利于子彈施加旋轉(zhuǎn)速度。

3 結(jié)果與討論

3.1 ALE 網(wǎng)格法與SPH 粒子法界面處波形

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

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

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

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

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

3.2 有、無(wú)轉(zhuǎn)速條件下的界面波

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

圖6 高速?zèng)_擊條件下界面波形成過(guò)程Fig.6 Formation process of interface wave under high speed impact

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

圖7 xy_應(yīng)力對(duì)比圖Fig.7 Xy_stress contrast chart

3.3 ALE 與SPH 耦合

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

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

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

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

圖10 ALE 與SPH 耦合結(jié)果Fig.10 Result of ALE and SPH coupling

圖10（a）中可以發(fā)現(xiàn)，耦合后界面處的波形沒(méi)有任何影響，圖10（b）為x面上y方向上的剪應(yīng)力處于同一色塊，應(yīng)力傳遞具有一致性，意味著耦合的成功。即對(duì)于模型復(fù)雜或難以簡(jiǎn)化的模型，可采用ALE 與SPH 方法耦合的方法建模計(jì)算模型。

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

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

4 結(jié)束語(yǔ)

本文對(duì)比了ALE 法和SPH 粒子法旋轉(zhuǎn)子彈沖擊基板的案例，分析了SPH 粒子法更適合應(yīng)用于沖擊類(lèi)問(wèn)題，并針對(duì)SPH 粒子法存在的粒子數(shù)過(guò)多導(dǎo)致的計(jì)算成本增加，進(jìn)行ALE 與SPH 粒子法耦合優(yōu)化。綜上，可獲得如下結(jié)論：

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

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