基于SPH與FEM的結構入水分析方法

李上明

摘要：建立基于光滑粒子動力學（smoothed particle hydrodynamics， SPH）、有限元法（finite element method， FEM）和無反射邊界耦合的結構入水分析方法，將無限水域利用無反射邊界條件截斷成有限水域，將有限水域分為流體變形大的SPH區(qū)域、流體變形小的FEM區(qū)域和聲學流體FEM區(qū)域，結構用FEM離散。采用通用接觸算法模擬SPH與FEM的耦合，采用聲固耦合方法處理FEM區(qū)域之間的耦合，建立流固耦合的SPH-FEM分析方法。該方法結合SPH模擬大變形的優(yōu)點和FEM的高效性，可實現含自由液面變形、液體飛濺和無限水域等特點的流固耦合問題的模擬，為結構入水分析縮小離散區(qū)域、降低自由度和SPH粒子數等提供一種有效的分析方法。

關鍵詞：入水；流固耦合；SPH；SPH-FEM耦合

中圖分類號：O353.4；TB115.1

文獻標志碼：B

0 引 言

飛機水上迫降、返回艙水上回收等問題均可認為是結構入水問題。結構入水常面臨局部損壞、器件失靈、姿態(tài)失控等問題。隨著現代軍事和航空、航海領域的發(fā)展，結構入水受到越來越多的關注。

隨著計算機性能的提高和數值計算方法的發(fā)展，數值方法已成為模擬結構入水過程的主要手段之一。張岳青等[1]利用基于有限元網格和ALE算法的流固耦合算法進行飛船返回艙入水分析，楊衡等[2]利用有限元與勢流理論分析結構入水過程及其載荷情況。與有限元網格類方法相比，光滑粒子動力學（smoothed particle hydrodynamics， SPH）方法[3]因其能有效模擬大變形、界面跟蹤問題，可避免網格畸變等優(yōu)點，被廣泛應用于結構入水問題中。

OGER等[4]應用弱可壓縮SPH方法模擬研究楔形體入水問題，分析結果逼近解析解和試驗值。SHAO[5]應用不可壓縮SPH方法模擬結構入水問題。CHEN等[6]比較弱可壓縮和不可壓縮SPH方法模擬不可壓縮流體的差異，結果表明弱可壓縮SPH方法通過人工狀態(tài)方程模擬流體壓力形式簡單、易計算執(zhí)行，但會遇到壓力波動問題。不可壓縮SPH方法能獲取更光滑的粒子壓力[7-8]，但其以犧牲計算效率為代價。隨著弱可壓縮SPH理論的發(fā)展，壓力波動可通過高階SPH離散格式[9]、人工黏性[10]、密度重新初始化[11]、改進的邊界處理算法[12]和自由液面處理算法[13]等手段得到一定的改進?；谌蹩蓧嚎sSPH理論，LIU等[14-15]和HU等[16]模擬結構入水或超彈性結構流固耦合問題。

上述研究結果表明，SPH方法能較好地模擬結構入水流固耦合問題，但在邊界處理和無限水域處理方面還需改進，主要表現在以下2個方面。

（1）結構入水過程涉及強烈的流固耦合作用，目前的動態(tài)固壁處理算法與流固耦合模型有待改進。LIU等[14]相對系統(tǒng)地論述各類固壁處理算法，主要包括排斥力固壁處理算法、動態(tài)固壁處理算法、傳統(tǒng)固壁處理算法和耦合動態(tài)固壁處理算法等的優(yōu)缺點。排斥力固壁處理算法形式簡單、易計算，適用于復雜的固壁邊界形狀，但其在邊界上流體粒子數量不足，無法保證計算精度。動態(tài)固壁處理算法在固壁附近增加虛擬粒子，可在一定程度上提高SPH的計算精度，同樣適用于復雜的固壁邊界形狀，但存在壓力波動問題，并且流體粒子會出現穿透固壁的情況。傳統(tǒng)固壁處理算法在固壁面上增加邊界粒子并在固壁附近增加虛擬粒子提升計算精度，其適合簡單構型的固壁邊界，并且每個時間步需要通過鏡像重新生成虛擬粒子，分析效率低。耦合動態(tài)固壁處理算法通過將動態(tài)固壁處理算法和排斥力固壁處理算法耦合，采用軟排斥力模型，可有效避免固壁穿透和壓力大幅波動現象，已應用于模擬剛體和超彈性體流固耦合問題，但其軟排斥力模型因結構剛性不同而有所差異，是否適用于其他彈塑性結構還需進一步驗證。目前，固壁處理算法是SPH方法的研究熱點之一。

（2）通用SPH公式未較好處理無限水域截斷邊界的能量吸收特性。研究結果表明：無限水域截斷邊界因反射能量未完全吸收，對SPH分析精度和效率有較大影響[17]；雖然已有工作利用海綿層盡可能克服界面反射現象，但該方法中反射系數及相應參數的選擇缺乏理論基礎。目前，SPH方法無限水域截斷邊界無反射特性研究比較薄弱。

針對上述SPH方法模擬結構入水流固耦合問題的不足，采用SPH方法、有限元法（finite element method， FEM）和無反射邊界聯合模擬結構入水過程，建立SPH與FEM的耦合分析模型，實現含無限水域的結構入水分析。

1 結構入水系統(tǒng)離散形式

為有效捕捉結構入水過程的自由液面和流固耦合界面，模擬液體飛濺現象和無限水域無反射特性等，縮小無限水域離散區(qū)域、降低分析自由度，將整個無限水域利用無反射邊界條件截斷成有限水域，將有限水域又分成入水結構附近的近場域和相對較遠的遠場域。因入水結構的沖擊，近場域會產生大位移和液體飛濺現象，所以采用SPH方法離散和模擬；遠場域位移較小，采用FEM離散和模擬。整個系統(tǒng)的離散區(qū)域見圖1，其中：遠場域又分為遠場1區(qū)域和遠場2區(qū)域；結構采用FEM模擬。

對于近場域（SPH離散區(qū)域）與遠場2區(qū)域（聲學流體FEM離散區(qū)域）之間的遠場1區(qū)域，采用以位移為自由度的結構有限元模擬，本文稱其為結構流體FEM；遠場2區(qū)域采用以壓力為自由度的聲學流體有限元模擬；無限水域的無反射特性直接通過在遠場2區(qū)域的截斷邊界上施加無反射邊界條件描述。近場域滿足流體動力學控制方程，用SPH方法離散模擬，其流體壓力采用狀態(tài)方程計算。

上述離散過程中涉及的結構有限元、聲學流體有限元和無反射邊界條件等理論較為成熟，本文不再贅述。

2 SPH方程

3 狀態(tài)方程

對于圖1中的近場域和遠場1區(qū)域，采用Mie-Grüneisen線性Us-Up的Hugoniot狀態(tài)方程表示其本構關系。流體抗拉、抗壓本構（流體壓力p與流體密度ρ）關系為

4 耦合面處理

根據圖1可知，入水過程的耦合面有近場SPH離散區(qū)域與入水結構的耦合面、SPH離散區(qū)域與遠場1區(qū)域的耦合面、遠場1區(qū)域與遠場2區(qū)域的耦合面，以及截斷邊界。

SPH離散區(qū)域與入水結構和遠場1區(qū)域的耦合作用采用接觸算法處理，即SPH離散區(qū)域粒子與結構表面和遠場1區(qū)域的相互作用通過Abaqus中的通用接觸算法傳遞。

遠場1區(qū)域與遠場2區(qū)域的耦合作用采用聲固耦合算法描述，即遠場1區(qū)域（結構流體FEM離散區(qū)域）和遠場2區(qū)域（聲學流體FEM離散區(qū)域）分別采用結構有限元方程和聲學有限元方程模擬，并在耦合面上基于壓力與加速度傳遞關系建立耦合關系，從而實現遠場1區(qū)域與遠場2區(qū)域的耦合。

在截斷邊界上利用無反射邊界條件通過聲學阻抗邊界定義無限水域與遠場2區(qū)域的耦合作用。因該截斷邊界用于截斷無限聲學流體域，故該處采用聲學流體無反射邊界條件。文獻[19]和文獻[20]論述當前主要無反射邊界條件，其合理性已經得到驗證，并且部分無反射邊界條件已經集成到商業(yè)有限元軟件中，故本文不再對無反射邊界條件反射特性進行驗證。

基于上述處理方法，聯合使用結構流體有限元、聲學流體有限元和SPH，可實現SPH、FEM和無反射邊界條件三者耦合分析，充分發(fā)揮SPH方法模擬大變形和無反射邊界模擬無限水域的優(yōu)勢，有效模擬結構入水問題。

5 算 例

所有單元網格、無反射邊界條件和SPH粒子均在Abaqus CAE中完成，SPH粒子通過有限元單元轉換產生，即利用Abaqus將SPH離散區(qū)域離散成六面體單元，每一個單元根據轉換準則轉換成一個SPH粒子。采用時間轉換準則，即從t=0時刻起，將SPH區(qū)域的六面體單元轉換為SPH粒子。SPH核函數為三次B樣條核函數，粒子的質量、密度、光滑長度由Abaqus內部程序根據其六面體母單元自動計算。在下面2個算例中，由于粒子變形較大，隨著分析時間的增加，粒子分布極不均勻，因此粒子光滑長度采用變長度光滑長度。

圖1中不同區(qū)域的耦合模擬方法如下：粒子與有限元單元的相互作用采用Abaqus缺省設置的通用接觸算法模擬，聲學流體與結構流體的耦合采用綁定約束模擬，無反射邊界采用Abaqus中的平面無反射邊界條件模擬。模擬分析考慮重力加速度（取9.8 m/s2）的影響。

5.1 含彈性障礙物的水柱坍塌模擬

某含彈性障礙物的水柱坍塌初始構型見圖2。水柱在重力作用下發(fā)生坍塌并沖擊彈性障礙物，彈性障礙物底部固定并在水沖擊作用下發(fā)生大變形，水與彈性障礙物發(fā)生流固耦合作用。設定水柱寬L=0.146 m、高為2L，彈性障礙物高w=0.012 m、寬h=0.08 m，水箱寬W=4L、高H=0.365 m。

圖 2 含彈性障礙物的水柱坍塌初始構型

利用SPH與FEM耦合模擬水柱坍塌過程，用SPH方法模擬水柱，用FEM模擬彈性障礙物和剛性水箱，水柱與障礙物和水箱壁面的相互作用通過SPH粒子與有限元單元之間的接觸算法模擬。水柱本構關系采用狀態(tài)方程式（6）和（7）描述，密度ρ0=1 000 kg/m3、聲速c0=1 500 m/s、黏性系數μ=0.001 Pa·s。彈性障礙物彈性模量為1×106 N/m2，密度為2 500 kg/m3，泊松比為0。

采用三維等效模型分析圖2的二維問題，在紙面的垂直方向只有一個粒子，所有粒子只有平面位移，彈性障礙物左上角的水平位移模擬結果見圖3。

圖3中FE-SPH1、FE-SPH2、FE-SPH3、FE-NSPH和FE-NSPH-R方法對應的初始粒子距離l和粒子數見圖4和表1。FE-SPH2、FE-SPH3、FE-NSPH、FE-NSPH-R方法的粒子數分別為FE-SPH1方法粒子數的4.00倍、3.36倍、3.36倍和8.97倍，水柱采用的流體模型和SPH公式見表1，其余結果均來自文獻[16]。

圖3和表1表明：水柱采用流體不抗拉模型的結果比抗拉模型結果更加逼近于文獻[16]的分析結果。FE-SPH3分析在0.55 s左右出現提前終止情況，而FE-NSPH和FE-NSPH-R未出現分析終止情況。這說明NSPH公式比標準SPH公式具有更強的粒子不規(guī)則分布的適應能力，并且隨著粒子數的增多，計算結果逐漸逼近其他分析結果。

數值結果表明所采用的通用接觸算法可有效模擬SPH粒子與有限元單元的之間的相互作用，可替代SPH的固壁處理算法，采用的狀態(tài)方程可有效模擬水的不抗拉現象。

5.2 楔形體入水分析

考慮二維楔形體跌入深1 m的水池，楔形體橫截面見圖5，離散區(qū)域見圖6，參數見表2。楔形體考慮為剛體，采用有限元離散，密度為466.3 kg/m3，初始速度方向向下、大小為5.05 m/s，彈性模量為210 GPa，泊松比為0.3。水池用SPH、結構流體FEM和聲學流體FEM離散，水池底部固定，無吸收性，在聲學流體FEM離散區(qū)域左側邊界施加平面無反射邊界條件。聲學流體參數設置密度為1 000 kg/m3，聲速為1 438.00 m/s；SPH離散區(qū)域和結構流體FEM離散區(qū)域采用流體不抗拉狀態(tài)方程式（7），設置密度為1 000 kg/m3、聲速為1 438.00 m/s，黏性系數μ=0.001 Pa·s。

不同計算方法和不同網格楔形體速度計算結果比較分別見圖7和8。圖中mesh 1、mesh 2和mesh 3的結果是基于表2參數，利用SPH、FEM和無反射邊界條件計算得到的，mesh 4為未考慮結構流體和無反射邊界條件（除自由液面外，SPH離散區(qū)域其余邊界均與剛性壁面接觸）的結果，其余結果均來自文獻[22]。

與其他方法相比，本文結果逼近ZHAO模型的結果和試驗結果。圖7表明，mesh 2結果優(yōu)于mesh 1結果，mesh 2的網格比mesh 1密，說明隨著網格加密，計算收斂性更強。圖8中mesh 3和mesh 2的結果幾乎一致，表明縮小SPH離散區(qū)域幾乎不會對結果產生影響，SPH、FEM和無反射邊界耦合可有效模擬含無限水域和流體大變形的入水問題。mesh 4的速度計算結果偏小，表明在較小區(qū)域僅采用SPH離散水域、不考慮邊界能量吸收時，隨著分析終止時間的加長，其結果會無法接受。

6 結 論

建立基于SPH、FEM和無反射邊界條件耦合的結構入水分析方法，為結構入水提供一種分析思路。數值算例表明該方法正確，并有如下特點。

（1）利用Mie-Grüneisen狀態(tài)方程模擬流體特性，可有效模擬流體對結構的作用和流體不抗拉的特點。

（2）利用通用接觸算法代替虛擬粒子法，可有效模擬SPH固體邊界。

（3）利用FEM，可有效耦合SPH和無反射邊界條件。

（4）充分發(fā)揮SPH模擬大變形、FEM模擬小變形、無反射邊界條件模擬無限水域的優(yōu)點，為結構入水分析提供一種能縮小離散區(qū)域、降低自由度和SPH粒子數的分析方法，同時也為含液體飛濺或結構大變形流固耦合問題提供一種分析思路。

參考文獻：

[1] 張岳青， 徐緋， 金思雅， 等. 飛船返回艙水上回收的沖擊響應和入水姿態(tài)分析[J]. 振動與沖擊， 2014， 33（18）： 204-208.

[2] 楊衡， 孫龍泉， 龔小超， 等. 彈性結構入水砰擊載荷特性三維數值模擬研究[J]. 振動與沖擊， 2014， 33（19）： 28-34.

[3] MONAGHAN J J. Smoothed particle hydrodynamics and its diverse applications[J]. Annual Review of Fluid Mechanics， 2012， 44（1）： 323-346. DOI： 10.1146/annurev-fluid-120710-101220.

[4] OGER G， DORING M， ALESSANDRINI B， et al. Two-dimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics， 2006， 213（2）： 803-822. DOI： 10.1016/j.jcp.2005.09.004.

[5] SHAO S. Incompressible SPH simulation of water entry of a free-falling object[J]. International Journal for Numerical Methods in Fluids， 2009， 59（1）： 91-115. DOI： 10.1002/fld.1813.

[6] CHEN Z， ZONG Z， LIU M B， et al. A comparative study of truly incompressible and weakly compressible SPH methods for free surface incompressible flows[J]. International Journal for Numerical Methods in Fluids， 2013， 73（9）： 813-829. DOI： 10.1002/fld.3824.

[7] RAFIEE A， THIAGARAJAN K P. An SPH projection method for simulating fluid-hypoelastic structure interaction[J]. Computer Methods in Applied Mechanics and Engineering， 2009， 198（33）： 2785-2795. DOI： 10.1016/j.cma.2009.04.001.

[8] LIU X， XU H H， SHAO S D， et al. An improved incompressible SPH model for simulation of wave-structure interaction[J]. Computers & Fluids， 2013， 71： 113-123. DOI： 10.1016/j.compfluid.2012.09.024.

[9] CHEN J K， BERAUN J E. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems[J]. Computer Methods in Applied Mechanics and Engineering， 2000， 190（1）： 225-239. DOI： 10.1016/S0045-7825（99）00422-3.

[10] COLAGROSSI A， LANDRINI M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics[J]. Journal of Computational Physics， 2003， 191（2）： 448-475. DOI： 10.1016/S0021-9991（03）00324-3.

[11] DILTS G A. Moving least-squares particle hydrodynamics II： Conservation and boundaries[J]. International Journal for Numerical Methods in Engineering， 2000， 48（10）： 1503-1524. DOI： 10.1002/1097-0207（20000810）48：10<1503：：AID-NME832>3.0.CO；2-D.

[12] LIU M B， SHAO J R， CHANG J Z. On treatment of solid boundary in smoothed particle hydrodynamics[J]. Science China：Technological Sciences， 2012， 55（1）： 244-254. DOI： 10.1007/s11431-011-4663-y.

[13] ZHENG X， DUAN W Y， MA Q W. A new scheme for identifying free surface particles in improved SPH[J]. Science China： Physics， Mechanics and Astronomy， 2012， 55（8）： 1454-1463. DOI： 10.1007/s11433-012-4809-3.

[14] LIU M B， SHAO J R， LI H Q. Numerical simulation of hydro-elastic problems with smoothed particle hydrodynamic method[J]. Journal of Hydrodynamics： B， 2013， 25（5）： 840-847. DOI： 10.1016/S1001-6058（13）60412-6.

[15] LIU M B， SHAO J R， LI H Q. An SPH model for free surface flows with moving rigid objects[J]. International Journal for Numerical Methods in Fluids， 2014， 74（9）： 684-697. DOI： 10.1002/fld.3868.

[16] HU D A， LONG T， XIAO Y H， et al. Fluid-structure interaction analysis by coupled FE-SPH model based on a novel searching algorithm[J]. Computer Methods in Applied Mechanics and Engineering， 2014， 276（7）： 266-286. DOI： 10.1016/j.cma.2014.04.001.

[17] GONG K， LIU H， WANG B. Water entry of a wedge based on SPH model with an improved boundary treatment[J]. Journal of Hydrodynamics： B， 2009， 21（6）： 750-757. DOI： 10.1016/S1001-6058（08）60209-7.

[18] LIU G R， LIU M B. Smoothed particle hydrodynamics： A meshfree particle method[M]. Beijing： World Scientific Publishing Company， 2003.

[19] ESPINOZA H， RAMON C A， SANTIAGO B. A Sommerfeld non-reflecting boundary condition for wave equation in mixed form[J]. Computer Methods in Applied Mechanics and Engineering， 2014， 276（6）： 122-148. DOI： 10.1016/j.cma.2014.03.015.

[20] GIVOLI D. High-order local non-reflecting boundary conditions： A review[J]. Wave Motion， 2004， 39（4）： 319-326. DOI： 10.1016/j.wavemoti.2003.12.004.

[21] LIBERSKY L D， RANDLES P W， CARNEY T C， et al. Recent improvements in SPH modelling of hypervelocity impact[J]. International Journal of Impact Engineering， 1997， 20（6）： 525-532. DOI： 10.1016/S0734-743X（97）87441-6.

[22] YETTOU E M， DESROCHERS A， CHAMPOUX Y. Experimental study on water impact of a symmetrical wedge[J]. Fluid Dynamics Research， 2006， 38（1）： 47-66. DOI： 10.1016/j.fluiddyn.2005.09.003.

（編輯 武曉英）