基于SBFE虚粒子边界的结构入水数值模拟

2022-01-11李上明赵梓斌

兵器装备工程学报 2021年12期

李上明，赵梓斌

(中国工程物理研究院总体工程研究所，四川绵阳 621000)

1 引言

在海洋和军事工程中，常常会涉及到结构的出水和入水问题，许多学者做了大量的研究。譬如，陆宇等[1]采用LS-DYNA建立了水下垂直发射装置头罩外层内侧内压对头罩开裂过程及冲击力影响的有限元模型。黄鸿鑫等[2]通过结构入水的有限元模型分别研究了射弹的头部形状和质心位置对高速入水稳定性的影响。由于光滑粒子法(SPH)能够模拟结构入水过程中自由液面的破碎和翻转等问题，因此近些年来SPH算法逐步被广泛应用于结构入水的数值模拟。SPH方法是由 Gingold 和 Monaghan[3]和Lucy[4]提出的一种粒子类方法。对于有限水域中结构入水的数值模拟，流域通常采用固壁边界条件，但是当模拟的流域为无限域时，就会涉及到无反射边界的应用。就目前看来，SPH方法的无反射边界条件主要包括三大类：第一类为开放边界条件[5-8]，第二类为海绵层边界条件[9-13]，第三类为时间线插值边界条件[14]。

开放边界条件是指有进口和出口的固定区域，当粒子通过进口进入计算域时，则添加新的粒子，而当粒子通过出口离开计算域时，则将该粒子删除。Morris等[5]采用周期性边界条件模拟低雷诺数条件下的层流，但是周期性边界允许粒子通过出口以后再经过入口进入计算域，这样随着时间的推进，扰动粒子的数量将会持续下降，会降低模拟的精度。为了提高模拟精度，Lastiwka等[6]在出口的上游设置流入区，在出口的下游设置流出区，并且根据实际的流速确定流入流出的粒子数量，这样能有效地避免扰动粒子数量的下降，但是这种方法仅仅在流出域删除了粒子，依然不能完全避免流出域产生的反射波对上游计算域的影响。Tafuni等[8]在计算域的上游或者下游定义了缓冲区，缓冲区内粒子的物理量等于其对应诡粒子的物理量。但是诡粒子的位置信息是缓冲区内粒子经过边界对称映射得到的，这样就大大降低计算效率，尤其是边界形状比较复杂的时候。第二类无反射边界条件是海绵层边界条件。Lind[9]和Liu[10]在模拟波的传播的时候采用了阻尼层的概念，即粒子的速度随着与边界的靠近而呈指数衰减，这样就能消除一部分边界上的反射波。Altomave等[11]在模拟长凤波的产生和吸收的时候，同样采用了阻尼层的概念，只是根据耗散海滩的需要，将粒子速度的衰减由指数型改为二次型。密度变化层的概念是由Gong等[12]提出的，通过在计算域外围增加几十层密度可变的粒子，粒子密度随时间的变化率与其距边界的距离相关，以此来消除固壁边界的反射波。但是海绵层法在模拟无反射边界时，由于需要在计算域的外面设置几十层粒子，因此计算量较大，导致计算效率偏低，特别是对三维问题，计算量将迅速增加。最近， Pingping Wang[14]提出了时间域插值边界条件。时间域插值边界条件是一种基于特征线法的无反射边界，通过在边界设置3～4层的仅仅向下游传递干扰的边界粒子来避免反射波的产生。边界粒子的任一时刻的性质(速度、压力等)，则是根据波传播的速度，通过上游不同时刻的粒子采用拉格朗日插值来获取的。但是时间插值类边界条件需要存储上一个或几个时间步内的粒子信息，对计算内存要求较高。本文联合SBFEM算法和SPH算法，提出了一种新的无反射边界条件，即SBFE虚粒子边界，该边界具有计算精度高、效率快等优点。最后，我们采用2个结构入水算例验证了该SBFE虚粒子边界的准确性。需要特别说明的是，本文中的所有结果均采用FORTRAN 95自编程序计算所得。

2 SBFE虚粒子边界

无反射边界对结构入水、水下爆炸等领域的数值模拟具有很大的意义。如图1所示，SBFE虚粒子边界只需在SPH流体粒子的外层额外布置2～4层的SBFE虚粒子，虚粒子的层数应该满足大于或者等于SPH粒子的光滑半径。

图1 SBFE虚粒子边界示意图Fig.1 SBFE virtual particle boundary

SBFE虚粒子的速度和加速度是通过核函数从附近的SPH粒子中插值获得的，其插值公式为：

(1)

式(1)中：W表示光滑粒子法的核函数；N表示支持域内的粒子总数；下标i，j表示不同的粒子；m表示粒子质量；ρ表示粒子的密度。

进一步地，虚粒子的加速度值可以采用差分的方法来获得，即：

(2)

式(2)中，Δt表示时间步长。

比例边界单元示意图如图2所示。虚粒子之间可以构成一个个比例边界单元，因此每一层边界的无穷域动力刚度矩阵的刚度阵A和阻尼阵B能够按照SBFEM的高频近似方法获得[15]。

需要特别说明的是，位于同一位置的SBFE虚粒子和比例边界单元节点具有相同加速度值，所以比例边界单元节点的加速度值可以通过同一位置的SBFE虚粒子的加速度值获得。

进一步地就能求解每一个比例边界单元节点的等效加速度，其表达式为：

(3)

式(3)中：a表示由各个虚粒子组成的加速度列向量；ns表示单元外法向的单位向量；N表示比例边界单元的形函数。

图2 比例边界单元示意图Fig.2 The diagram of SBFEM

采用高频近似方法就能求解每层比例边界单元各节点的压力值，即：

(4)

其中，

(5)

式(5)中：p表示边界节点压力的列向量；p(i)(i=1，2，3，…，M)为辅助变量。

位于同一位置的SBFE虚粒子和比例边界单元节点具有相同的压力值，而位于比例边界单元节点间的SBFE虚粒子的压力，则可以采用线性插值的方式来获得。如图3所示，虚粒子m的压力值可以通过比例边界单元节点k和节点h的压力值获得，即：

(6)

图3 SBFE虚粒子和比例边界单元节点示意图Fig.3 SBFE virtual particles and SBFE nodes

3 SBFE虚粒子边界的计算效率和精度

我们验证了SBFE虚粒子边界的准确性和计算效率。算例中，水池的长0.2 m，宽0.1 m。楔形体的一条边长为0.04 m，底升角为30°(如图4所示)，以初速度3.15 m/s竖直落入水中，重力加速度为9.8 m/s2。另外设置一个同样大小的水池，该水池的边界采用固壁边界条件。为了验证SBFE虚粒子边界的准确性，还需要一个长为0.4 m，宽为0.2 m的水池(见图4(a))，以保证其在计算过程中，没有反射波进入感兴趣的区域。

图4 二维楔形体入水模型示意图Fig.4 The model of wedge water entry

表1表示不同边界条件SPH方法的计算耗时。从表1中可以明显的看出，SBFE虚粒子边界能够极大地减少计算域内的粒子数量，从而提高SPH方法模拟无穷域流体的计算效率。更重要的是， SBFE虚粒子边界的计算耗时和固壁边界的计算耗时几乎相同。这表明，SBFE虚粒子边界不会显著增加计算耗时，具有计算效率高的优点。

表1 不同边界条件的计算耗时

图5表示水池底边中点的压力曲线。在固壁边界条件下，我们能很明显地发现压力值存在2个波峰点，第1个波峰点是入射波和反射波叠加的峰值，其大小约为SBFE虚粒子边界结果的2倍，而第2个波峰点则是反射波与底边反射波，侧边反射波的多重叠加。从图5中可以发现，无穷远边界中点的压力值与SBFE虚粒子边界的中点压力值几乎相同。这说明了SBFE虚粒子边界能够准确的模拟边界的无反射特征。

图5 水池底边中点的压力曲线Fig.5 The pressure of point 1

4 无穷水域楔形体入水数值模拟

楔形体几何结构如图6所示，楔形体的高度h为0.01 m，底升角α为20°，结构的密度为600.3 kg/m3，结构的质量为0.18 kg。结构分别以初速度100 m/s、80 m/s、60 m/s、40 m/s、20 m/s、5 m/s落入水中。水池的大小为0.2 m×0.2 m。其中一个水池的边界采用SBFE虚粒子边界来模拟无穷水域的情形，另一个水池的边界则采用固壁边界条件来模拟有限域的结构入水问题。SPH粒子的初始间距为1e-3 m，整个流体域被离散为40 401个粒子。固体结构的初始粒子间距同样采用1e-3 m，整个固体域被离散为297个粒子。在数值模拟过程中，时间步长均为5e-7 s，总共模拟5e-3 s，模拟总步数为10 000步。

图7表示楔形体结构以100 m/s速度入水时流体的压力云图。图7中，左侧压力云图表示SBFE虚粒子边界的模拟结果，右侧压力云图表示固壁边界的模拟结果。

图6 楔形体几何结构图Fig.6 The geometry of two-dimensional wedge

图7 楔形体结构入水压力云图Fig.7 The pressure distribution of wedge water entry

从图7中可以看出，当压力波到达边界时，能穿过SBFE虚粒子边界传向无穷远，而固壁边界则会将压力波反射重新进入计算域，从而影响整个流域的压力场。图7表明SBFE虚粒子边界能够模拟压力波的透射过程，从而实现准确高效模拟无穷域中结构入水问题。

进一步的，我们模拟了刚体结构以100 m/s的速度入水的速度衰减曲线，如图8，结构入水深度曲线如图9所示。从图9中可以看出，在结构入水的初期，即入水砰击阶段，SBFE虚粒子边界的模拟结果和固壁边界的模拟结果几乎相同。这是因为此时反射波仍然没有进入计算域的缘故。但是，在结构入水模拟的后期，固壁边界的速度衰减比SBFE虚粒子边界的速度衰减更快，这是因为固壁边界的反射波抵达了结构底部，对结构施加了一个额外阻力的缘故。

图8 楔形体结构入水的速度衰减曲线Fig.8 The velocity of the wedge in water entry

图9 楔形体结构的入水深度曲线Fig.9 The depth of the wedge in water entry

不同初速度的结构入水深度曲线如图10。从图10可以看出，随着初速度的提升，结构入水深度迅速增加。表2展示了2种不同边界条件下结构在t=5 ms时的入水深度。

图10 不同初速度的结构入水深度曲线Fig.10 The depth of the wedge with different velocities in water entry

表2 结构不同入水初速度的入水深度Table 2 The depth error of the wedge with different velocities in water entry

在表2中，入水初速度越大，2种边界条件下模拟结果的相对误差越大。这表明对于无穷水域中结构高速入水问题，固壁边界反射波对模拟结果的影响已经不能再忽略。这是因为入水初速度越大，流体域中压力波的最大值越大，固壁边界的反射波也会越强烈，这样就会导致流体域内的压力场变动越剧烈，从而对无穷域内结构入水的数值模拟结果产生不可忽略的影响。因此，在SPH方法中模拟无穷域中结构低速入水时，采用固壁边界可以得到理想的结果，但是当结构入水速度提高时，则需要进一步考虑边界反射波的影响。SBFE虚粒子边界能很好地模拟边界波的透射过程，从而能够准确地模拟无穷域中结构入水的问题。