一種新型SPH-FEM 耦合算法及其在沖擊動力學(xué)問題中的應(yīng)用*

張志春，強洪夫，高巍然

(第二炮兵工程學(xué)院，陜西 西安710025)

SPH 最早由L.B.Lucy［1］、R.A.Gingold 等［2］用于解決三維開放空間天體物理學(xué)問題，特別是多變性問題。SPH 的近似函數(shù)建立在一系列離散的粒子上，不需要借助網(wǎng)格，通過使用一系列任意分布的粒子來求解相應(yīng)的偏微分方程組，從而得到精確穩(wěn)定的數(shù)值解。由于不需要借助網(wǎng)格，并且具有Lagrange性質(zhì)，SPH 方法在處理大變形、自由表面等問題時具有明顯的優(yōu)勢，目前已經(jīng)廣泛應(yīng)用于連續(xù)固體力學(xué)和流體力學(xué)。但是在固體力學(xué)領(lǐng)域，SPH 的計算精度和效率較有限元方法低，并且其核函數(shù)不具有Kronecher delta 張量的性質(zhì)，施加邊界條件困難。因此耦合SPH 和FEM，在大變形區(qū)域使用SPH 粒子離散，小變形區(qū)域使用有限元離散，是計算沖擊動力學(xué)問題的一種有效手段。

在同一材料中同時使用SPH 和FEM，耦合界面的計算是問題的關(guān)鍵。在這方面G.R.Johnson［3-4］和S.Ferna＇ndez-Me＇ndez 等［5］分別提出了各自的算法。G.R.Johnson 將SPH 粒子固結(jié)在有限元節(jié)點，通過固結(jié)的粒子傳遞與有限單元之間的相互作用力，但在計算界面粒子應(yīng)變率時沒有考慮鄰近有限元節(jié)點的影響，無法保證耦合界面物理量的連續(xù)性。S.Ferna＇ndez-Me＇ndez 在有限單元和SPH 粒子之間設(shè)置了過渡區(qū)域，在過渡區(qū)域使用混合插值函數(shù)，不需要將有限元節(jié)點替換為粒子，但混合插值函數(shù)的計算過于復(fù)雜。

本文中，為解決SPH 和FEM 耦合界面的計算問題，提出一種新型耦合算法。將SPH 粒子固結(jié)在有限元節(jié)點，在有限元節(jié)點設(shè)置背景粒子，通過背景粒子的方式將有限元節(jié)點納入SPH 粒子的鄰近搜索列表。這種方法可避免SPH 的光滑核函數(shù)被邊界截斷，消除其邊界效應(yīng)，確保耦合界面處物理量的連續(xù)性。使用該耦合算法對圓柱形鋼彈正沖擊鋼板進行三維數(shù)值模擬，其中鋼彈和鋼板中心區(qū)域采用SPH 離散，鋼板其余區(qū)域采用有限元離散，施加固支邊界條件;鋼板參數(shù)的計算采用含損傷的Johnson-Cook 本構(gòu)模型和Grüneisen 狀態(tài)方程。子彈與靶板之間采用SPH 粒子接觸算法。

1 SPH-FEM 耦合算法

1.1 算法流程

如圖1 所示，SPH 粒子固結(jié)在有限元節(jié)點，其中小實線圓代表SPH 粒子，大虛線元代表粒子i 的支持域，小虛線元代表設(shè)置在有限元節(jié)點處的背景粒子。背景粒子具有SPH 粒子的屬性，其質(zhì)量、位置、速度、應(yīng)力與相應(yīng)有限元節(jié)點保持一致。計算流程如圖2 所示，在每個時間步開始前，有限元節(jié)點的質(zhì)量、位置、速度、應(yīng)力傳遞給相應(yīng)的背景粒子。在每個時間步結(jié)束后，界面處SPH 實粒子的位置、速度傳遞給相應(yīng)的有限元節(jié)點。SPH 粒子采用跳蛙格式求解Navier-Stokes 方程，有限元采用中心差分法求解顯示動力學(xué)方程。對SPH 粒子進行數(shù)值積分時，有限元節(jié)點以背景粒子的形式加入SPH 鄰近搜索算法中，即任何位于鄰近搜索范圍內(nèi)的SPH 粒子和有限元節(jié)點都被加入到臨近列表。如對實粒子i 的積分包含SPH 粒子n1、n2、…、n5和有限元節(jié)點n6、n7、n8。背景粒子僅被動地被實SPH 粒子搜索，而自身不進行SPH 數(shù)值積分，數(shù)據(jù)的更新由相應(yīng)的有限元節(jié)點確定。

使用該SPH-FEM 耦合算法處理耦合界面問題，對于SPH 部分，由于將耦合界面附近的有限元節(jié)點納入到粒子的鄰近搜索列表，積分時不會被界面截斷，消除了SPH 的邊界效應(yīng)。對于有限元部分，界面處有限元節(jié)點數(shù)據(jù)由對應(yīng)的SPH 粒子確定，相當(dāng)于在耦合界面處對有限單元施加了邊界條件，由于有限元形函數(shù)滿足Kronecher delta 張量性質(zhì)，該方案是可行的。

圖1 SPH 粒子固結(jié)在有限單元Fig.1 SPH particles attached to finite elements

圖2 SPH-FEM 耦合算法流程Fig.2 A solution procedure for the coupled SPH-FEM algorithm

1.2 時間步長控制

如圖2 所示，SPH 部分采用跳蛙格式求解Navier-Stokes 方程，具體步驟如下:

(1)第1 個時間步結(jié)束后

(2)在隨后的每個時間步開始前

(3)在每個時間步結(jié)束后

式中:si可以表示粒子i 的密度、速度和能量，ri和vi分別表示粒子i 的位置和速度，t 表示時間，Δt 表示時間步長。

對于有限元部分，采用中心差分法求解顯示動力學(xué)方程，節(jié)點加速度、速度和位移分別為

式中:M 是系統(tǒng)集中質(zhì)量矩陣，F(xiàn)int是內(nèi)力，F(xiàn)co是粘性力，F(xiàn)ext是外力。

在Lagrangian 框架下，SPH 與FEM 的耦合要求兩者的積分必須同步，數(shù)據(jù)的傳輸必須是在同一時間點。這就要求在每一步計算中采用相同的時間步長，在此選擇SPH 與FEM 時間步長的較小者

式中:SPH 時間步長ΔtSPH=βh/c，F(xiàn)EM 時間步長ΔtFEM≤L/c，其中h 是SPH 的光滑長度，L 是最小單元尺寸，c 是材料聲速，β 是時間步長比例系數(shù)。

2 數(shù)值算例

為驗證SPH-FEM 耦合算法的準(zhǔn)確性，對桿一端受壓和圓柱形鋼彈正沖擊鋼板的問題進行了計算。采用J.J.Monaghan［6］提出的人工應(yīng)力法來消除拉伸不穩(wěn)定性造成的數(shù)值斷裂，采用強洪夫等［7］提出的完全變光滑長度法來消除變光滑長度帶來的影響。

2.1 桿一端受壓

為驗證SPH 粒子與有限單元耦合界面物理量的連續(xù)性，對方形桿一端受壓問題進行計算。桿件的尺寸為5 mm×5 mm×100 mm，采用線彈性材料模型，密度ρ=7.83 t/m3，彈性模量E=207 GPa，泊松比ν=0.3。如圖3 所示，桿件一半離散為1 835 個SPH粒子，另一半離散為1 250 個有限單元。在有限元一端突加矩形脈沖載荷，加載類型為均布壓強p=40 GPa，加載時間為0 ～1 μs，計算總時間為14 μs。為了進行對比驗證，相同的模型完全用有限元離散進行了計算。

圖4 為耦合界面附近的SPH 粒子750、800、1 050、1 100 和有限元節(jié)點909、910、1 215、1 216 沿z 軸的速度和應(yīng)力曲線，其中虛線表示對比算例中與SPH 粒子相同位置處節(jié)點的速度和應(yīng)力。因為4 個粒子和4 個節(jié)點分別沿桿中心線對稱分布，其速度和應(yīng)力曲線分別重合，表明該算法完全滿足對稱性要求。初始時刻，速度和應(yīng)力均為0，之后桿件產(chǎn)生沿z 軸的應(yīng)力波，在t=10.5 μs 和t=11.0 μs，節(jié)點和粒子最大速度vz分別達到-582 和-539 m/s，應(yīng)力最大值σz分別達到-27.3 和-27.9 GPa。圖4 表明耦合界面兩側(cè)粒子的運動與節(jié)點保持一致，并且與對比算例中相同位置處的節(jié)點也保持一致，表明該耦合算法滿足耦合界面處物理量的連續(xù)性要求。

圖3 桿一端受壓模型Fig.3 The model of a bar pressed at one end

圖4 耦合界面速度和應(yīng)力曲線Fig.4 Plots for velocity and stress of the coupled interface

2.2 圓柱形鋼彈正沖擊鋼板

對Arne tool 鋼制圓柱彈丸與Weldox 460 E 鋼板正沖擊發(fā)生沖塞破壞的過程進行三維數(shù)值計算，計算模型如圖5 所示。其中鋼彈及靶板中心區(qū)域采用SPH 粒子離散，粒子間距1 mm，粒子總數(shù)47 133。靶板其余部分采用六面體單元離散，粒子附近網(wǎng)格尺寸為1 mm，其余部分網(wǎng)格尺寸由靶板中心向四周逐漸增大，單元總數(shù)18 816。子彈與靶板初始間距2 mm，靶板四周施加固支邊界條件，計算時間160 μs。子彈與靶板的接觸采用R.Vignjevic 等［8］提出的SPH 粒子接觸算法。

圖5 沖擊模型尺寸及離散方式Fig.5 Impact model size and its discretization mode

為描述靶板材料的屈服應(yīng)力及損傷演化，這里引入T.B?rvik 等［9］提出的修正Johnson-Cook 強度模型，該模型將材料的屈服強度表示為損傷變量、等效應(yīng)變、等效應(yīng)變率和溫度的函數(shù)

式中:A、B、C、n、m 是材料常數(shù);D 為損傷變量，D=0 表示材料沒有損傷，D=1 表示材料完全失效;εd是等效累積損傷塑性應(yīng)變是等效累積塑性應(yīng)變是用戶定義參考等效應(yīng)變率;量綱一溫度T*=(T－T0)/(Tm－T0)，T0是室溫，Tm是材料熔點。實際上材料出現(xiàn)宏裂紋時，損傷變量的臨界值小于1，失效準(zhǔn)則可描述為D=Dc≤1，其中Dc是損傷變量臨界值。損傷變量D 是等效累積塑性應(yīng)變ε 的函數(shù)，可描述如下

式中:εt是損傷閾值，εf是等效斷裂塑性應(yīng)變，與材料的應(yīng)力三軸度、應(yīng)變率和溫度相關(guān)。G.R.Johnson和W.H.Cook［10］提出的剪切損傷演化模型中將εf描述如下

式中:D1～D5為材料常數(shù)，σ*=σm/σeq為應(yīng)力三軸度，σm=(σx+σy+σz)/3 為平均正應(yīng)力。假設(shè)材料處于絕熱狀態(tài)，靶板材料溫度的升高是由于沖擊過程中的塑性功轉(zhuǎn)換為熱量造成的，溫度的變化計算為

式中:ρ 為材料密度，cp為材料比定壓熱容，α 為比例常數(shù)，Wp為塑性功。靶板材料參數(shù)見表1，靶板壓強采用Grüneisen 狀態(tài)方程計算。

子彈沖擊過程中變形較小，采用線彈性模型，材料參數(shù)為:E=200 GPa，ν=0.33，ρ=7.838 t/m3。

表1 靶板材料參數(shù)Table 1 Material parameters of target

當(dāng)子彈初速v0=285.4 m/s 時，其數(shù)值計算結(jié)果如圖6 所示，子彈、靶塞以及靶板沖塞型破壞的實驗結(jié)果如圖7 ～8 所示［9］。從圖中可以看出，計算結(jié)果與實驗結(jié)果吻合很好。在t=7 μs 時，子彈與靶板接觸，彈頭正前方處的靶板開始加速，靶板被擠壓產(chǎn)生壓縮變形。之后在彈頭周圍的靶板出現(xiàn)剪切區(qū)，由于接觸部位較大的塑性應(yīng)變，剪切區(qū)域內(nèi)粒子的損傷變量迅速演化，達到臨界值，部分粒子完全失效，裂紋逐漸向靶板后部擴展。在沖擊后期，失效模式同時包含靶板內(nèi)部的剪切斷裂和靶板背面的拉伸損傷，在t=120 μs，靶塞在子彈的推動下與靶板完全脫離。

圖6 沖擊過程的數(shù)值模擬結(jié)果Fig.6 Numerical simulations of the impact process

圖7 沖擊后的子彈和靶塞Fig.7 The projectile and plug after impact

圖8 沖擊后的靶板橫截面Fig.8 The cross section of the target after impact

表2 為不同初速時數(shù)值模擬值與實驗值的對比，其中vpr是子彈余速，vtr是靶塞速度，df是靶板前部穿孔直徑，dr是靶板后部穿孔直徑。計算和實驗存在的靶板穿孔直徑差異，主要是因為當(dāng)子彈和靶板間距達到2 倍光滑長度時，產(chǎn)生接觸力，這些與實際的物理過程不能完全相符。通過加密SPH 粒子設(shè)置，可以在一定程度上減少該誤差。

表2 計算值與實驗值對比Table 2 Comparisons between computed results and experimental data

3 結(jié) 論

將SPH 粒子固結(jié)在有限單元節(jié)點，提出了一種新型SPH-FEM 耦合算法，為SPH 和FEM 的耦合計算提供了一種新的途徑。該方法可以在FEM 計算困難的區(qū)域用SPH 粒子離散，能夠充分發(fā)揮SPH 和FEM 的特點:SPH 具有處理大變形問題的優(yōu)勢;FEM 計算精度和效率更高。應(yīng)用該耦合算法，在SPH粒子的邊界區(qū)域使用有限單元離散，以有限元方的式施加邊界條件，可以解決SPH 施加邊界條件的難題。該方法可以方便地推廣到FEM 與其他無網(wǎng)格方法的耦合計算，如EFG 和RKPM。

提出的SPH-FEM 耦合算法要求在耦合界面處有限單元的尺寸和SPH 粒子間距保持一致，當(dāng)計算模型幾何形狀復(fù)雜時，前處理建模比較麻煩。在工程實際中使用該SPH-FEM 耦合算法時，為了降低對前處理建模要求，需要根據(jù)耦合界面附近有限元單元的尺寸自動調(diào)整耦合界面附近SPH 粒子的光滑長度，這一技術(shù)還有待于在今后的工作中進一步深入研究。

［1］ Lucy L B.Numerical approach to testing the fission hypothesis［J］.Astronomical Journal，1977，82(12):1013-1024.

［2］ 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.

［3］ Johnson G R.Linking of Lagrangian particle methods to standard finite element methods for high velocity impact computations［J］.Nuclear Engineering and Design，1994，150(2-3):265-274.

［4］ Johnson G R，Stryk R A，Beissel S R.SPH for high velocity impact computations［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1-4):347-373.

［5］ Ferna＇ndez-Me＇ndez S，Bonet J，Huerta A.Continuous blending of SPH with finite elements［J］.Computers＆Structures，2005，83(17-18):1448-1458.

［6］ Monaghan J J.SPH without a tensile instability［J］.Journal of Computational Physics，2000，159(2):290-311.

［7］ 強洪夫，高巍然.完全變光滑長度SPH 法及其實現(xiàn)［J］.計算物理，2008，25(5):569-575.QIANG Hong-fu，GAO Wei-ran.SPH method with fully variable smoothing lengths and implementation［J］.Chinese Journal of Computational Physics，2008，25(5):569-575.

［8］ Vignjevic R，De Vuyst T，Campbell J C.A frictionless contact algorithm for meshless methods［C］∥ICCES，2007，3(2):107-112.

［9］ B?rvik T，Langseth M，Hopperstad O S，et al.Ballistic penetration of steel plates［J］.International Journal of Impact Engineering，1999，22(9-10):855-886.

［10］ Johnson G R，Cook W H.Fracture characteristics of three metals subjected to various strains，strain rates，temperatures and pressures［J］.Engineering Fracture Mechanics，1985，21(1):31-48.