基于SPH与FEM的结构入水分析方法

李上明

摘要：建立基于光滑粒子动力学（smoothed particle hydrodynamics， SPH）、有限元法（finite element method， FEM）和无反射边界耦合的结构入水分析方法，将无限水域利用无反射边界条件截断成有限水域，将有限水域分为流体变形大的SPH区域、流体变形小的FEM区域和声学流体FEM区域，结构用FEM离散。采用通用接触算法模拟SPH与FEM的耦合，采用声固耦合方法处理FEM区域之间的耦合，建立流固耦合的SPH-FEM分析方法。该方法结合SPH模拟大变形的优点和FEM的高效性，可实现含自由液面变形、液体飞溅和无限水域等特点的流固耦合问题的模拟，为结构入水分析缩小离散区域、降低自由度和SPH粒子数等提供一种有效的分析方法。

关键词：入水；流固耦合；SPH；SPH-FEM耦合

中图分类号：O353.4；TB115.1

文献标志码：B

0 引 言

飞机水上迫降、返回舱水上回收等问题均可认为是结构入水问题。结构入水常面临局部损坏、器件失灵、姿态失控等问题。随着现代军事和航空、航海领域的发展，结构入水受到越来越多的关注。

随着计算机性能的提高和数值计算方法的发展，数值方法已成为模拟结构入水过程的主要手段之一。张岳青等[1]利用基于有限元网格和ALE算法的流固耦合算法进行飞船返回舱入水分析，杨衡等[2]利用有限元与势流理论分析结构入水过程及其载荷情况。与有限元网格类方法相比，光滑粒子动力学（smoothed particle hydrodynamics， SPH）方法[3]因其能有效模拟大变形、界面跟踪问题，可避免网格畸变等优点，被广泛应用于结构入水问题中。

OGER等[4]应用弱可压缩SPH方法模拟研究楔形体入水问题，分析结果逼近解析解和试验值。SHAO[5]应用不可压缩SPH方法模拟结构入水问题。CHEN等[6]比较弱可压缩和不可压缩SPH方法模拟不可压缩流体的差异，结果表明弱可压缩SPH方法通过人工状态方程模拟流体压力形式简单、易计算执行，但会遇到压力波动问题。不可压缩SPH方法能获取更光滑的粒子压力[7-8]，但其以牺牲计算效率为代价。随着弱可压缩SPH理论的发展，压力波动可通过高阶SPH离散格式[9]、人工黏性[10]、密度重新初始化[11]、改进的边界处理算法[12]和自由液面处理算法[13]等手段得到一定的改进。基于弱可压缩SPH理论，LIU等[14-15]和HU等[16]模拟结构入水或超弹性结构流固耦合问题。

上述研究结果表明，SPH方法能较好地模拟结构入水流固耦合问题，但在边界处理和无限水域处理方面还需改进，主要表现在以下2个方面。

（1）结构入水过程涉及强烈的流固耦合作用，目前的动态固壁处理算法与流固耦合模型有待改进。LIU等[14]相对系统地论述各类固壁处理算法，主要包括排斥力固壁处理算法、动态固壁处理算法、传统固壁处理算法和耦合动态固壁处理算法等的优缺点。排斥力固壁处理算法形式简单、易计算，适用于复杂的固壁边界形状，但其在边界上流体粒子数量不足，无法保证计算精度。动态固壁处理算法在固壁附近增加虚拟粒子，可在一定程度上提高SPH的计算精度，同样适用于复杂的固壁边界形状，但存在压力波动问题，并且流体粒子会出现穿透固壁的情况。传统固壁处理算法在固壁面上增加边界粒子并在固壁附近增加虚拟粒子提升计算精度，其适合简单构型的固壁边界，并且每个时间步需要通过镜像重新生成虚拟粒子，分析效率低。耦合动态固壁处理算法通过将动态固壁处理算法和排斥力固壁处理算法耦合，采用软排斥力模型，可有效避免固壁穿透和压力大幅波动现象，已应用于模拟刚体和超弹性体流固耦合问题，但其软排斥力模型因结构刚性不同而有所差异，是否适用于其他弹塑性结构还需进一步验证。目前，固壁处理算法是SPH方法的研究热点之一。

（2）通用SPH公式未较好处理无限水域截断边界的能量吸收特性。研究结果表明：无限水域截断边界因反射能量未完全吸收，对SPH分析精度和效率有较大影响[17]；虽然已有工作利用海綿层尽可能克服界面反射现象，但该方法中反射系数及相应参数的选择缺乏理论基础。目前，SPH方法无限水域截断边界无反射特性研究比较薄弱。

针对上述SPH方法模拟结构入水流固耦合问题的不足，采用SPH方法、有限元法（finite element method， FEM）和无反射边界联合模拟结构入水过程，建立SPH与FEM的耦合分析模型，实现含无限水域的结构入水分析。

1 结构入水系统离散形式

为有效捕捉结构入水过程的自由液面和流固耦合界面，模拟液体飞溅现象和无限水域无反射特性等，缩小无限水域离散区域、降低分析自由度，将整个无限水域利用无反射边界条件截断成有限水域，将有限水域又分成入水结构附近的近场域和相对较远的远场域。因入水结构的冲击，近场域会产生大位移和液体飞溅现象，所以采用SPH方法离散和模拟；远场域位移较小，采用FEM离散和模拟。整个系统的离散区域见图1，其中：远场域又分为远场1区域和远场2区域；结构采用FEM模拟。

对于近场域（SPH离散区域）与远场2区域（声学流体FEM离散区域）之间的远场1区域，采用以位移为自由度的结构有限元模拟，本文称其为结构流体FEM；远场2区域采用以压力为自由度的声学流体有限元模拟；无限水域的无反射特性直接通过在远场2区域的截断边界上施加无反射边界条件描述。近场域满足流体动力学控制方程，用SPH方法离散模拟，其流体压力采用状态方程计算。

上述离散过程中涉及的结构有限元、声学流体有限元和无反射边界条件等理论较为成熟，本文不再赘述。

2 SPH方程

3 状态方程

对于图1中的近场域和远场1区域，采用Mie-Grüneisen线性Us-Up的Hugoniot状态方程表示其本构关系。流体抗拉、抗压本构（流体压力p与流体密度ρ）关系为

4 耦合面处理

根据图1可知，入水过程的耦合面有近场SPH离散区域与入水结构的耦合面、SPH离散区域与远场1区域的耦合面、远场1区域与远场2区域的耦合面，以及截断边界。

SPH离散区域与入水结构和远场1区域的耦合作用采用接触算法处理，即SPH离散区域粒子与结构表面和远场1区域的相互作用通过Abaqus中的通用接触算法传递。

远场1区域与远场2区域的耦合作用采用声固耦合算法描述，即远场1区域（结构流体FEM离散区域）和远场2区域（声学流体FEM离散区域）分别采用结构有限元方程和声学有限元方程模拟，并在耦合面上基于压力与加速度传递关系建立耦合关系，从而实现远场1区域与远场2区域的耦合。

在截断边界上利用无反射边界条件通过声学阻抗边界定义无限水域与远场2区域的耦合作用。因该截断边界用于截断无限声学流体域，故该处采用声学流体无反射边界条件。文献[19]和文献[20]论述当前主要无反射边界条件，其合理性已经得到验证，并且部分无反射边界条件已经集成到商业有限元软件中，故本文不再对无反射边界条件反射特性进行验证。

基于上述处理方法，联合使用结构流体有限元、声学流体有限元和SPH，可实现SPH、FEM和无反射边界条件三者耦合分析，充分发挥SPH方法模拟大变形和无反射边界模拟无限水域的优势，有效模拟结构入水问题。

5 算 例

所有单元网格、无反射边界条件和SPH粒子均在Abaqus CAE中完成，SPH粒子通過有限元单元转换产生，即利用Abaqus将SPH离散区域离散成六面体单元，每一个单元根据转换准则转换成一个SPH粒子。采用时间转换准则，即从t=0时刻起，将SPH区域的六面体单元转换为SPH粒子。SPH核函数为三次B样条核函数，粒子的质量、密度、光滑长度由Abaqus内部程序根据其六面体母单元自动计算。在下面2个算例中，由于粒子变形较大，随着分析时间的增加，粒子分布极不均匀，因此粒子光滑长度采用变长度光滑长度。

图1中不同区域的耦合模拟方法如下：粒子与有限元单元的相互作用采用Abaqus缺省设置的通用接触算法模拟，声学流体与结构流体的耦合采用绑定约束模拟，无反射边界采用Abaqus中的平面无反射边界条件模拟。模拟分析考虑重力加速度（取9.8 m/s2）的影响。

5.1 含弹性障碍物的水柱坍塌模拟

某含弹性障碍物的水柱坍塌初始构型见图2。水柱在重力作用下发生坍塌并冲击弹性障碍物，弹性障碍物底部固定并在水冲击作用下发生大变形，水与弹性障碍物发生流固耦合作用。设定水柱宽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粒子与有限元单元之间的接触算法模拟。水柱本构关系采用状态方程式（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公式具有更强的粒子不规则分布的适应能力，并且随着粒子数的增多，计算结果逐渐逼近其他分析结果。

数值结果表明所采用的通用接触算法可有效模拟SPH粒子与有限元单元的之间的相互作用，可替代SPH的固壁处理算法，采用的状态方程可有效模拟水的不抗拉现象。

5.2 楔形体入水分析

考虑二维楔形体跌入深1 m的水池，楔形体横截面见图5，离散区域见图6，参数见表2。楔形体考虑为刚体，采用有限元离散，密度为466.3 kg/m3，初始速度方向向下、大小为5.05 m/s，弹性模量为210 GPa，泊松比为0.3。水池用SPH、结构流体FEM和声学流体FEM离散，水池底部固定，无吸收性，在声学流体FEM离散区域左侧边界施加平面无反射边界条件。声学流体参数设置密度为1 000 kg/m3，声速为1 438.00 m/s；SPH离散区域和结构流体FEM离散区域采用流体不抗拉状态方程式（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离散区域其余边界均与刚性壁面接触）的结果，其余结果均来自文献[22]。

与其他方法相比，本文结果逼近ZHAO模型的结果和试验结果。图7表明，mesh 2结果优于mesh 1结果，mesh 2的网格比mesh 1密，说明随着网格加密，计算收敛性更强。图8中mesh 3和mesh 2的结果几乎一致，表明缩小SPH离散区域几乎不会对结果产生影响，SPH、FEM和无反射边界耦合可有效模拟含无限水域和流体大变形的入水问题。mesh 4的速度计算结果偏小，表明在较小区域仅采用SPH离散水域、不考虑边界能量吸收时，随着分析终止时间的加长，其结果会无法接受。

6 結 论

建立基于SPH、FEM和无反射边界条件耦合的结构入水分析方法，为结构入水提供一种分析思路。数值算例表明该方法正确，并有如下特点。

（1）利用Mie-Grüneisen状态方程模拟流体特性，可有效模拟流体对结构的作用和流体不抗拉的特点。

（2）利用通用接触算法代替虚拟粒子法，可有效模拟SPH固体边界。

（3）利用FEM，可有效耦合SPH和无反射边界条件。

（4）充分发挥SPH模拟大变形、FEM模拟小变形、无反射边界条件模拟无限水域的优点，为结构入水分析提供一种能缩小离散区域、降低自由度和SPH粒子数的分析方法，同时也为含液体飞溅或结构大变形流固耦合问题提供一种分析思路。

参考文献：

[1] 张岳青， 徐绯， 金思雅， 等. 飞船返回舱水上回收的冲击响应和入水姿态分析[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.

（編辑 武晓英）