基于FE-SPH耦合算法的球体冲击问题研究

李旦， 赵廷渝

(中国民用航空飞行学院飞行技术学院, 四川广汉618307)

基于FE-SPH耦合算法的球体冲击问题研究

李旦， 赵廷渝

(中国民用航空飞行学院飞行技术学院, 四川广汉618307)

为充分利用SPH方法处理大变形问题与FEM方法处理冲击问题时各自的优势，提出了新型的三维SPH-FEM耦合算法，该算法在大变形区域采用SPH粒子离散，从而克服FEM方法单元畸变存在的问题，其他模拟区域使用FEM单元离散，既可以避免SPH边界效应，也能保证耦合界面物理量具有连续性。使用该算法分别对球体冲击土壤与水域的过程进行了三维数值模拟，选取含损伤的Johnson-Cook模型和Mie-Gruneisen状态方程对土壤与水域进行参数计算，对比分析球体在土壤与水域中的速度、加速度以及有效应力数据结果，分析可知，此次模拟结果与其他已有仿真结果基本吻合，验证了采用FE-SPH 耦合算法处理结构与流体相互作用问题是可行和有效的，为求解流固耦合问题提供了一条崭新的途径。

FE- SPH耦合算法；边界条件；流固耦合；数值模拟

引言

诸如飞机水上迫降、返回舱回收以及波浪对海洋结构的作用等情况是工程实际中普遍存在流固耦合问题。由于该耦合效应具有非线性特性，利用传统解析法进行求解具有一定的局限性，因此，数值模拟方法是解决此问题的重要手段[1]。

FEM方法采用拉格朗日格式求解结构域，ALE格式或者欧拉格式求解流体域，该方法虽然计算效率高，但需要划分网格，尤其在拓扑结构变化很大时进行网格划分困难，且后期计算处理较复杂，极有可能导致网格畸变[2]。因此，上述方法在模拟复杂自由表面流动与结构大变形相互作用问题时存在困难。SPH方法起初用于处理三维开放空间天体物理学问题，由于该方法不受网格束缚，在模拟大变形以及自由表面流动等问题时能有效避免网格畸变，且能在拉格朗日框架下对流体和结构同时进行求解，目前已广泛应用于流体力学与连续固体力学。但是，在固体力学领域，用该方法模拟固体介质存在拉伸不稳定问题，计算精度和效率相比FEM方法较低。由于SPH方法中近似核函数没有Kronecher delta张量性质，处理边界条件比较困难[3]，因此，单独采用SPH方法解决此类问题仍存在不足。

近年来，为解决SPH和FEM耦合界面的计算问题，FE-SPH耦合算法被提出并得到快速发展，其基本思想为：采用SPH方法处理局部大变形区域，FEM计算处理小变形区域，便于综合利用有限元法计算精度好、效率高、易于处理复杂边界与SPH方法模拟自由表面流动的优势[4]。目前，FE-SPH耦合算法在工程领域得到广泛应用，尤其是在数值模拟冲击动力学问题方面，例如在处理高速冲击及爆炸问题上普遍采用该耦合算法，同时，在材料加工处理方面，FE-SPH耦合算法也得到了应用。

1FE-SPH耦合算法基本概念

1.1三维FE-SPH耦合算法基本原理

三维FE-SPH耦合算法基本原理如图1所示。Johnson G R提出将SPH粒子固结在有限元节点，通过固结的方式保证两者之间传递相互作用力，且SPH粒子与对应单元节点加速度相同，从而实现SPH与FE耦合[5]。但是，在计算界面粒子应变率时该耦合算法没有考虑邻近有限元节点的影响，因此，耦合界面物理量的连续性无法得到保证。同时，该方法要求耦合界面节点数与粒子数相同，两者位置需要重合，限制了耦合界面周围粒子和单元划分，导致粒子与单元密度不匹配，对计算精度产生影响。

图1粒子和单元节点的连接示意图

为了有效地克服原算法中存在的缺陷，在有限元节点处设置了与有限元节点对应的背景粒子，背景粒子具有SPH粒子的属性，其质量、位置、速度、加速度与相应有限元节点保持一致，对SPH粒子进行数值积分时，有限元节点以背景粒子的形式加入SPH邻近搜索算法中，通过背景粒子转化方式将有限元节点纳入SPH粒子的邻近搜索列表，SPH使得耦合界面周围的粒子和单元划分变得相当自由，从而避免SPH方法中光滑核函数被边界截断，也能去除边界效应，保证了耦合界面处物理量的连续性[6]。改进的耦合算法不仅实施简单，也能使计算精度提高。

1.2三维FE-SPH耦合算法基本理论

FE-SPH自适应算法分为四类，本文主要对耦合算法进行分析研究，FE-SPH耦合算法示意图如图2所示。由于耦合算法要确定12个速度，计算流程相比其他算法更为繁琐复杂，因此，确保方程组正确求解相当关键[7]。建立上述方程需要根据线动量守恒，角动量守恒以及速度、位置在耦合点匹配的假设条件，计算步骤分为以下三步：

(1)检测四节点法向速度的增量。

(2)检测主面内速度增量。

(3)检测主面三节点存储角动量部分的速度增量。

图2有限元与光滑粒子耦合算法

步骤(1)中，总穿透量表达式：

δtotal=δadjust+δattach

步骤(2)中，光滑粒子与主面相对移动距离为：

在图2(a)中，s′轴、粒子和主面相对运动方向一致，三个节点的速度变化量ΔU1、ΔU2、ΔU3也沿s′轴方向。四节点在主面内的速度增量形式与法向速度增量形式相同，即：

上述各步骤中，步骤(1)中三个轴向保持角动量守恒；步骤(2)在s′轴方向保持角动量守恒，保证节点在主面上的位置和速度得到匹配；步骤(3)保证x′、y′轴方向保持角动量守恒，其中x′、y′方向角动量变化量与主面内粒子点速度变化的大小相等方向相反，且法向动量守恒。通过上述计算最终得到以下三个方程：

ΔV3=-M1ΔV1+M2ΔV2/M3

其中，ΔV1、ΔV2、ΔV3表示主面内粒子点速度变化量，如图2(b)所示。

综合以上计算，可以得到有限单元与光滑粒子耦合时需要调整的参量值。但是，由于耦合界面处周围节点的影响，计算结果并不能一次调整到位，一般需要迭代1～5次之后才能使方程组最终达到平衡。

2数值建模

本文利用SPH粒子模拟具有大变形的流体，FE单元模拟具有不规则形状的固体，固体模型在软件LS-DYNA中建立，同时完成网格划分，流体部分模型在LS-PREPOST完成。

2.1FE-SPH建模

为了减少计算成本，本文只建立1/2模型，如图3所示。图3(a)所示球体采用默认的Belytshco-Lin-Tsay壳单元进行映射式有限元建模，并在SECTION_SHELL中定义壳单元算法、积分规律和相关属性，球体视为刚体，故其旋转自由度受到约束，通过关键字MAT_RIGID来定义[8]；图3(b)所示为FE单元，通过设置单元密度来模拟土壤与水域；图3(c)所示模型主要由SPH粒子组成，用来模拟具有大变形的流体或固体，控制其粒子疏密程度相当重要，粒子过多可能使后期计算时间过长，出错概率增大，最终导致仿真失败；通过初期建模与网格划分，最终得到需要的模型，如图3(d)所示。模型材料具体密度见表1。

图3FE-SPH耦合模型

表1模型密度

2.2关键参数选择

FE-SPH耦合算法中流体与结构交界面处理是一个难点，因此，在关键字BOUNDARY设置上需要仔细分析，选取合适的参数。例如关键字CONTACT、CONTROL、DEFINE、INITIAL、LOAD、SET等必须根据数值模拟要求赋值，本文中球体初始速度为4.88m·s-1，球体加载曲线比例因子为1，终止时间为0.1s。

水的状态方程采用Mie-Gruneisen状态方程[9]，表示为：

P=Pc+PT=A(μ)+B(μ)·eVO

3仿真结果分析

3.1冲击过程分析

球体冲击土壤与水域过程基本相似，因此选取冲击土壤过程分析研究。图4所示为球体冲击土壤历程图，从中可以看出FE单元与SPH粒子模拟土壤受冲击后作用后变形效果明显不同。由于球体与SPH粒子直接碰撞，SPH粒子相比FE单元变形量明显要大，产生粒子溅射现象，效果更为直观，可以看出SPH方法在模拟自由表面流动问题时极具优势[10]。但是，同样可以看出FEM单元在模拟土壤固体结构动力问题上有其独特优势，体现出有限元法求解结构动力问题时计算精度好、效率高、易于处理复杂边界的优点。

图4球体冲击土壤历程图

图5所示为SPH粒子受力云图，从俯视图中可以看出，SPH粒子在冲击力作用下受力区域向四周延伸，其中心区域受力最大，变形量最大，粒子溅射现象最为明显。

图5粒子受力云图

图6所示为FE单元模拟土壤受力云图，可以清晰体现土壤各个部分受力大小，充分体现了有限元法计算精度高，效率高的特点，数据结果对后期分析处理土壤结构动力问题相当重要[11]。

图6FEM单元受力云图

3.2数据结果分析

首先对球体撞击土壤与水域时速度、加速度数据进行滤波处理，得到最终的球体速度、加速度时间历程曲线，再分别选取FE单元与SPH粒子的一个单元得到有效应力随时间变化曲线。通过对比球体撞击土壤与水域计算结果，可以发现FE-SPH耦合算法在模拟冲击动力学问题上的优越性，球体速度时间历程曲线如图7与图8所示。

图7球体冲击土壤速度时间历程曲线

图8球体冲击水域速度时间历程曲线

图7与图8所示分别为球体冲击土壤与水域速度时间历程曲线，对比可以看出，在初始速度、能量、加载方式等因素相同的条件下，0.1s时间内球体撞击水域条件下速度首先减到最小值(图中取绝对值)，并逐渐趋于零，受其密度、粘性等因素影响，冲击力比较均匀，速度变化比较平滑[12]。相比较而言，从图7曲线变化趋势可知，由于土壤密度较大，冲击过程中阻力明显高于水域，导致速度减小率更大，这与其他仿真模拟结果基本吻合。

球体加速度时间历程曲线如图9与图10所示。

图9球体冲击土壤加速度时间历程曲线

图10球体冲击水域加速度时间历程曲线

比较图9与图10加速度时间历程曲线，发现球体冲击土壤过程中加速度峰值明显高于水域条件下加速度峰值，最大值达到540.41m/s2，这与土壤密度较大有直接关系。同时，球体撞击土壤与水域都出现二次加速度峰值，但水域条件下速度明显要快，这是由于材料变形涉及冲击波、复杂的材料特性，材料属性不同，变形效应会有差别[13]。SPH粒子所产生的冲击波间接作用于FE单元，使其发生形变，并储存一部分能量，当FE单元吸收的能量重新释放时，会使球体产生反向加速度，而这种结构动力学效应在撞击土壤时更为明显，本次模拟结果与其他已有数值仿真结果基本吻合。

SPH粒子与FE单元有效应力随时间变化历程曲线如图11所示。

图11有效应力时间历程曲线

分别选取耦合界面处一个SPH粒子与对应的单元节点，经过数据处理得到如图11所示的曲线，可以看出SPH粒子承受的有效应力高于FE单元，在0.05s之后，对FE单元施加的有效应力趋于恒定值，这与单元材料属性、结构动力特性有关，仿真结果与已有数值结果基本吻合[14]。

4结束语

对于诸如民机水上迫降、返回舱回收、高速冲击及爆炸问题，由于材料变形极大，因此是数值模拟中的一类难题。同时，材料变形涉及冲击波、复杂的材料特性等因素，想要成功实现数值模拟比较困难[13]。综合FEM法与SPH方法在模拟材料大变形问题时各自的优点，即FEM方法计算效率高和SPH方法模拟大变形能力强的优点，发展为比较成熟的FE-SPH耦合算法，该算法为高速冲击及爆炸问题进行数值模拟提供了一种高效、准确的途径。在对FE-SPH耦合算法不断完善的基础上，开发实用的可市场化的数值仿真软件，进一步将其应用于工程实际问题，这是该研究方向后续发展的重点[14]。

[1] 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.

[2] 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.

[3] MONAGHAN J J.SPH without a tensile instability[J].Journal of Computational Physics,2000,159(2):290-311.

[4] VIGNJEVIC R,DE VUYST T,CAMPBELL J C.A frictionless contact algorithm for meshless methods[C]∥Proceedings of International Conference on Computational & Experimental Engineering and Sciences(ICCES 2007),Miami,USA,November 27-29,2007,3(2):107-112.

[5] JOHNSON G R,STRYK R A.Conversion of 3D distorted elements into meshless particles during dynamic deformation[J].International Journal of Impact Engineering,2003,28(9):947-966.

[6] BONET J,KULASEGARAM S.Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[J].International Journal for Numerical Method in Engineering,2000,47(6):1189-1214.

[7] 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.

[8] DE VUYST T,VIGNJEVIC R,CAMPBELL J C.Coupling between meshless and finite element methods[J].International Journal of Impact Engineering,2005,31(8):1054-1064.

[9] BEISSEL S R,GERLACH C A,JOHNSON G R.Hypervelocity impact computations with finite elements and meshfree particles[J].International Journal of Impact Engineering,2006,33(1-12):80-90.

[10] JOHNSON G R,COOK W H.A constitutive model and data for metals subjected to large strains,high strain rates and high temperatures[C]∥Proceedings of the 7th International Symposium on Ballistics,Hague,Netherlands,April 19-21,1983:541-547.

[11] HU D A,HAN X,XIAO Y H,et al.Research developments of smoothed particle hydrodynamics method and its coupling with finite element method[J].Chinese Journal of Theoretical and Applied Mechanics,2013,45(5):639-652.

[12] XIAO Y H,HU D A,HAN X.A coupling algorithm of finite element and smoothed particle hydrodynamics[J].Chinese Journal of Computational Physics,2011,28(2):219-225.

[13] YANG G,LIANG C,LIU P,et al.Numerical simulation of ricochet problem of projectile penetrating into concrete target based on 3d FE-SPH adaptive coupling algorithm[J].Engineering Mechanics,2013,30(9):276-282.

[14] SCHWER L.Optional strain-rate forms for the Johnson-cook constitutive model and the role of the parameter epsion_0[C]∥Proceedings of the 6th European LS-DYNA Users’ Conference,Gothenburg,June 15-17,2007:1-17.

Study on Sphere Shock Problem Based on FE-SPH Coupling Algorithm

LIDan,ZHAOTingyu

(Civil Aviation Flight University of China, Guanghan 618307, China)

In order to make full use of SPH method to deal with the problem of large deformation and FEM, the new SPH-FEM coupling algorithm is proposed in this paper. In this algorithm, SPH particle discretization is used in the large deformation area to overcome the FEM distortion problem, and FEM unit discretization are used in other simulation area, both to avoid SPH boundary effects, but also to ensure the physical continuity of the coupling interface. The three-dimensional numerical simulation of the impact of the sphere on the soil and water is carried out. The Johnson-Cook model and the Mie-Gruneisen equation of state are used to calculate the soil and water parameters, and the velocity of the sphere in the soil and water , acceleration and effective stress data are compared. It can be seen that the simulation results are in good agreement with other existing simulation results, which proves that it is feasible and effective to deal with the structure-fluid interaction problem by using FE-SPH coupling algorithm. Coupling problem provides a new way.

FE-SPH coupling algorithm; boundary condition; fluid-solid coupling; numerical simulation

2016-12-02

李 旦(1992-),男,甘肃天水人,硕士生,主要从事飞机水上追降数值方法方面的研究,(E-mail)1171085449@qq.com; 赵廷渝(1965-),男,重庆人，教授,硕士,主要从事民用航空发动机性能及控制技术方面的研究,(E-mail)2494238988@qq.com

1673-1549(2017)02-0053-06

10.11863/j.suse.2017.02.11

O351

A