基于SPH 模型的2 种压力振荡抑制算法比较分析

孙雷，程聪，郑博学，刘昌凤

1 大连理工大学船舶工程学院，辽宁大连116024

2 高新船舶与深海开发装备协同创新中心，上海200240

3 大连海洋大学海洋与土木工程学院，辽宁大连116023

0 引 言

在船舶流体力学中，液舱晃荡为典型的周期性水动力学问题［1］。海上航行的液化天然气（liquefied natural gas，LNG）船在风、浪、流的联合作用下，其液舱内会产生巨大的砰击压力，这将对舱壁造成直接的损坏［2］。研究如何准确预报舱壁所受的砰击载荷，对于船舶的结构设计具有重要意义。

由于光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）方法在模拟自由表面的流动、翻卷和破碎时无需特别处理［3-5］，对交界面捕捉优势明显［6-7］，因此十分适合求解液舱晃荡问题。然而，该方法仍存在压力振荡的问题，威胁到流场的稳定性。虽然利用不可压缩求解器［4］或对核函数和核梯度进行修正［8］可以在一定程度上缓解压力振荡，但仍未从根本上解决该问题。

国内外学者针对该问题进行了深入研究，并提出了2 种可行的方法：密度正则化和数值耗散技术。

第1 种方法是对粒子进行密度正则化处理。Randles 和Libersky［9］提出了一种零阶精度的密度正则化（Shepard 密度正则化）算法。Shao 等［8］基于该模型改进了离散方法，在模拟二维矩形液舱晃荡问题中获得了稳定的流场压力；Chen 等［6］在该方法的基础上改进，提出了适用于大密度比两相流的密度正则化算法，模拟了溃坝和瑞丽泰勒不稳定性等问题，改善了压力结果的稳定性。Dilts［10-11］提出了一种一阶精度的密度正则化（MLS密度正则化）算法。Colagrossi 和Landrini［12］利用这种方法缓解了流场的压力噪声，稳定了大密度比情况下的界面流动问题；Chen 等［4］将MLS 密度正则化算法应用到二维液舱晃荡问题中，提高了流场压力的计算精度。虽然MLS 密度正则化算法能够一定程度上抑制压力振荡，但是在长时期的模拟过程中，并不能够保持系统内粒子体积的守恒［13］。因此，Shepard 密度正则化算法因其简捷性和经济性，被广泛使用。

第2 种方法是应用数值耗散技术来抑制压力振荡。Monaghan［14］考虑速度的耗散，在动量方程中引入了人工黏性项，人工黏性能够有效减小压力噪声，但是该方法是在整个流场施加速度耗散的，会引入极大的数值耗散。为了减小密度的非物理性振荡，Molteni 和Colagrossi［13］通过在连续性方程中引入密度耗散项建立了δ-SPH 模型，对二维弱可压缩自由表面流动问题进行了模拟，获得了光滑稳定的压力数值结果。Antuono 等［15］在Molteni 和Colagrossi 工作的基础上改进了密度耗散项，解决了流体体积不守恒和自由表面不收敛的问题，有效抑制了流场压力的振荡，使之更适用于静力学问题。Antuono 等［16］对δ-SPH 模型的能量平衡进行了深入分析，发现该模型的耗散机制比传统SPH 模型要好，因此δ-SPH 模型也获得了广泛的应用。

在计算一些船舶水动力学问题时，选取合适有效且耗散性比较小的方法来抑制流场压力的振荡尤为重要。而目前就此类问题进行比较和探讨的资料甚少。Calogrossi 等［13，15］在人工黏性修正的基础上加入密度耗散项综合比较，证实了密度耗散项加入后能获得更光滑的压力场，但是一些低冲击性的水动力学问题没必要施加人工黏性，因为这会带来不必要的人为数值耗散。另外，人工黏性是一种速度耗散算法，而密度耗散项是一种密度耗散算法，为了体现实际的物理意义，本文将在仅考虑物理黏性、不加入人工黏性的情况下，比较2 种密度耗散算法对压力振荡的抑制效果。在实际工程领域，液舱晃荡问题的运动形式复杂，本文将考虑运动速度不同情况下2 种算法的优劣，以便为耗散项在复杂工况下的选取提供指导。此外，先前讨论晃荡问题的流动形态时都是直观比较波面形状，关于定量的波高测量算法没有太多参考信息。基于以上原因，本文将比较2 种密度耗散算法（Shepard 密度正则化算法和δ-SPH 密度耗散项算法）对压力振荡的抑制效果。

1 控制方程

1.1 基本控制方程

假设流体不可压缩，考虑绝热条件，流体运动应遵循质量守恒定律和动量守恒定律。在拉格朗日坐标系下，控制方程被简化为如下形式的连续性方程和动量守恒方程：

式中：ρ 为流体密度，kg/m3；v 为速度矢量，m/s；P 为压强，Pa ；F 为体积力矢量，m/s2；D*为黏性耗散项，m/s2。

1.2 状态方程

SPH 方法在求解问题时常把不可压缩流体假设为弱可压缩［5，17］。这样，求解压力时可避免求解关于压力的泊松方程，压力可以直接由状态方程显式求解得到。在本文中，状态方程采用如下形式［6］：

式中：cs为数值声速，m/s ；ρ0为流体的参考密度，kg/m3；P0为背景压力，Pa。为了控制流体的弱可压缩性，需保证流体的密度变化范围在1%以内，数值声速通常取为流场内最大速度vmax的10倍，即满足cs＞10vmax。

1.3 2 种SPH 耗散算法的密度求解方程

1.3.1 Shepard 密度正则化

在传统SPH 方法中，临近边界的流体粒子支持域内没有足够的临近粒子来进行插值，导致计算边界时得到的物理量出现了振荡。在实际计算过程中，会出现粒子分布不均的情况，这种非物理现象会严重影响物理量计算的精度。在计算流场压力时，为了恢复支持域的紧支性，保持SPH 粒子近似的一致性，一种常用的办法就是采用Shepard密度正则化［8］，对计算的密度进行重整：

1.3.2 δ-SPH 方法的耗散项

对于弱可压缩SPH 模型，为了过滤高频压力噪声，抑制压力振荡，Molteni 和Colagrossi［13］基于数值耗散的原理，在连续性方程中加入了密度耗散项，其表达式为：

2 SPH 数值计算模型

2.1 SPH 基本原理

在SPH 方法中，计算空间被离散成一系列任意分布的粒子，粒子携带流场计算的基本物理信息，包括位置、速度、密度、压力、质量等，这些粒子之间通过光滑核函数建立联系并相互作用。SPH方法的近似过程包括两部分［3］：核近似和粒子近似。核近似是指流场内某粒子点的场函数及其导数在该粒子支持域内通过光滑核函数来近似积分的过程。场函数及其导数的核近似形式如下：

式中：r 为某 点的位置矢量，r͂=r-r′ ；Ω 为支持域。

为了对控制方程进行离散，将上述积分形式的表达式转换成对支持域内粒子的求和形式，即粒子近似。场函数及其导数的核近似表达式的离散求和形式为：

式中，q= ||ri-rj/h。

2.2 物理黏性

本文仅考虑动量方程的耗散项为物理黏性项，物理黏性的表达式采用Lo 和Shao［18］给出的形式，其离散形式为：

式中：υ 为流体运动黏性系数，m2/s；rij=ri-rj；μi=ρiυ，μj=ρjυ，μi，μj为i，j 粒子的动力黏性系数，Pa ⋅s；η=0.1h，为防止出现奇点的系数。根据以上描述，控制方程可以表达成如下离散形式：

2.3 边界条件

根据边界粒子的布置及其场变量的离散形式，边界条件可以分为斥力边界条件、动力边界条件和耦合动力学边界条件等。斥力边界条件通过改变斥力系数来人工调节斥力的作用范围，但斥力粒子只负责提供斥力。另外，斥力粒子通常设置为一层，这会导致边界粒子支持域的截断，使边界处出现数值振荡。动力边界条件在边界外设置若干层虚粒子，虚粒子通过参与控制方程的运算进而提供物理斥力，但容易出现边界粒子的非物理性穿透。耦合动力学边界由Liu 等［19］提出，Chen 等［6］在此基础上进行了改进，对二维情况下的大密度比多相流动问题进行了模拟。本文采用Chen等［6］改进的耦合动力学边界条件，如图1所示。

图1 耦合动力学边界条件Fig.1 Schematic diagram of coupled dynamic boundary condition

如图，在边界上布置一层斥力粒子，斥力粒子的属性设为流体的属性。另外，根据选取的核函数确定边界以外所布置虚粒子的层数，斥力粒子的物理信息和虚粒子的物理信息都根据其所在支持域的周围粒子插值得到。当流体靠近边界区域时，斥力粒子对流体施加的斥力为：

式中：Rij为流体粒子所受到的边界对它的排斥力，N；mfluid为流体粒子的质量，kg。

2.4 时间积分

本文的时间积分方案采用预测-校正法，该算法共分为2 步。

第1 步：预测步。

第2 步：校正步。

一个时间步以后，粒子最终的物理量为

由于该积分方法是一种显式积分法，为使积分结果稳定，计算时间步长的选取需要满足柯朗-弗里德里希斯-列维（Courant-Friedrichs-Lewy，CFL）约束条件［3］，并满足一定的准则［17，20］，即

式中：Δt1为满足声速与光滑长度的时间步长；Δt2为满足外力条件的时间步长；Δt3为满足黏性耗散的时间步长；f为粒子加速度，N/m2。

3 数值算例及结果讨论

针对液舱静止和晃荡情况，分别采用2 种SPH 耗散模型对其进行数值模拟。本文的计算模型基于Li 等［21］的实验，将三维液舱简化为二维方箱，如图2 所示。水箱的长度L=1 m，高度H=1 m，水深Dw=0.3 m。不考虑空气的影响，水箱内的介质为水，水的密度取为ρ=1 000 kg/m3，水的运动粘度取为υ=1.0×10-6m2/s，背景压力取为P0=0 Pa。关于该模型的验证参见文献［22］，本文不再赘述。

图2 数值模型示意图Fig.2 Schematic diagram of numerical model

3.1 静水情况

在静水算例中，在距离水箱底部0.2 m 的左侧壁面上设置一个压力监测点，用来实时监测该处的静水压力。粒子初始间距为Δx0=0.005 m ，计算步长根据CFL 条件取为Δt=5×10-5s ，模拟的总物理时间为55 s，除耗散项不同以外，2 种SPH算法的其他数值模拟参数完全一致。

对于传统SPH 方法，流场内部粒子的密度由密度正则化最终得到，对流场内部的粒子每隔30个时间步长执行一次密度正则化。在δ-SPH 算法中，流场内部粒子的密度由引入了密度耗散项的连续性方程得到，柯西数δ取为0.05。将监测点处的压力减去该点处的初始压力P*，并做无量纲处理得到2 种SPH 算法的压力耗散率，用Rdf表示，即Rdf=(P-P*)/P*。图3 所示为传统SPH 方法（即Shepard 密度正则化算法）与耗散系数为0.05 的δ-SPH 算法得到的压力耗散程度随时间变化的曲线。

图3 2 种耗散算法的耗散率Fig.3 Dissipation rate of two diffusive algorithms

从图中可以看出：在模拟的初期，因重力施加的初始效应，压力耗散率曲线呈现出振荡的特点。随着模拟时间的增加，2 种算法的耗散率都呈现出持续减小的趋势，至t=55 s 时压力耗散都趋近于理论值3.5%。事实上，Molteni和Colagrossi［13］已经指出：密度耗散项算法是一种一阶精度的耗散算法，精度较高，而密度正则化算法则是一种零阶精度的空间均化方法，精度较低；在初始阶段，压力和速度在达到平衡的过程中会出现强烈的非线性振荡，此时使用一阶精度的密度耗散项算法和零阶精度的密度正则化算法会出现不同的耗散效果，一阶精度的耗散强度较高，速率较快，因此可以看出两条耗散曲线存在差别；在达到稳定阶段后，压力振荡较小，压力分布趋于静水压力分布，此时2 种耗散算法的耗散效果几乎一致，因此2 条曲线趋于吻合。在其他参数相同的情况下，2种算法在经历了相同的耗散程度后趋于一致，证明2 种算法都达到收敛。

为了比较2 种SPH 耗散项算法在抑制压力振荡上的表现，分别选取2 种算法在5 个时间节点（t=0，12，25，40 和55 s）的流场压力分布图进行了对比，如图4 所示。各分图中，左图和右图分别为密度正则化算法和δ-SPH 密度耗散项算法所得的流场压力分布。

通过对比可以发现：在数值模拟的初期，如图4（a）～图4（b）所示，密度正则化算法与δ-SPH密度耗散项算法均能够获得较为光滑的压力场。在数值模拟的中后期，如图4（c）～图4（e）所示，密度正则化算法获得的流场压力分布变得粗糙，不再保持压力梯度分层过渡的趋势，存在高频压力噪声，压力振荡严重，且振荡的程度随着时间的增加而加剧了，最终导致数值模拟失效。与此相反，在相同的时间节点上，δ-SPH 密度耗散项算法获得的流场压力分布均匀、层次清晰、光滑过渡，虽然因为黏性数值耗散的原因，计算的压力值呈现下降的趋势，但是随着时间的推移，依然能够获得稳定的流场压力。

图4 静水下几个典型时刻2种算法的流场压力分布比较Fig.4 Comparison of the pressure distribution in flow field at some typical moments in still water

以上现象表明：在静水工况数值模拟中，密度正则化算法不能长期有效过滤流场的压力噪声，抑制压力振荡的效果较差；而δ-SPH 密度耗散项算法能够长期有效地光滑流场的压力，抑制压力振荡的效果较好。

从能量守恒的角度分别计算了2 种密度耗散算法中流体的动能变化情况（对动能进行无量纲处理：将动能与初始时刻流场的总势能作比），得到如图5 所示的无量纲动能时历曲线。

图5 无量纲动能曲线Fig.5 Non-dimensional kinetic energy curves

分析图5 可以发现：含有密度耗散项的δ-SPH算法与含有密度正则化的传统SPH 算法相比，能够更快地达到收敛，并且δ-SPH 算法能够在初始效应过后的整个数值模拟过程中基本保持动能为0，即满足动能守恒；而传统SPH 方法在数值模拟的中后期出现了速度波动，导致动能出现0.02%～0.05%的增长，动能不稳定性增加。根据这个趋势可以预测：在长时期的数值模拟中，密度正则化方法的动能守恒性会越来越差；而δ-SPH 算法的动能守恒性更好，证明其具有更好的鲁棒性，更适用于长时期的数值模拟。

3.2 液舱晃荡

为了比较2 种SPH 耗散项算法在周期性运动条件下对流场压力振荡的抑制效果，对二维矩形水箱在横摇状态下的液体晃荡过程进行了数值模拟。在本算例中，密度耗散算法参数的选取均同静水算例一致。水箱以底部中心为旋转点做横摇运动，横摇角度为正弦形式：

式中：θ 为横摇角，rad；A 为横摇角幅值，rad；ω为横摇激励频率，rad/s。

液舱晃荡流体固有频率计算采用Faltinsen等［1］的公式：

为探讨不同转动角速度情况下2 种算法对压力振荡的抑制程度，根据数值实验结果，分别选取并讨论2 组工况：工况a——角速度较低情况；工况b——角速度较高情况。具体工况设置如下：

式中，T 为横摇激励周期，s。

3.2.1 角速度较低情况

对于工况a，由于边界运动的角速度较小，所以流体运动基本处于低速流动状态。为观察在长期模拟过程中2 种耗散算法对压力振荡的抑制效果，模拟的总物理时间设为200 s，图6（a）～图6（d）分别给出了该工况下对应于t = 0.25T，1.75T，2.25T 和2.75T 时刻2 种算法的流场压力分布情况。各分图中，左图和右图分别为密度正则化算法和δ-SPH 密度耗散项算法所得的流场压力分布。

图6 工况a 下几个典型时刻2 种算法的流场压力分布比较Fig.6 Comparison of the pressure distribution in flow field at some typical moments for case a

通过对比发现：2 种耗散算法在流动模拟的初期能够获得稳定光滑的压力场。如图6（b）左图所示，在t=1.75T 时刻，密度正则化算法流场的内部产生了高频压力噪声，流场压力分布的分层状态被破坏。随着时间的增加，如图6（c）～图6（d）左图所示，这种压力振荡的现象一直持续，其压力分布规律与密度正则化算法计算静水问题的规律一致。而在同样时刻，利用δ-SPH 密度耗散项算法得到的流场压力分布却分层清晰，压力分布连续。虽然由于数值耗散的原因，2 种方法获得的压力幅值都在衰减，但总体看来，δ-SPH 密度耗散项算法得到的流场压力的光滑程度要优于传统SPH 密度正则化算法的结果。

2 种算法结果存在差别的原因在于：在转动角速度较低的情况下，速度变化量不大，物理黏性还没有完全起效，而密度正则化算法每隔一定的时间才会对流场压力关于空间均化，一定程度上会造成压力振荡的积累，而该算法又不能每步施加，因为这将给数值结果带来极大的耗散，因此，密度正则化算法得到的压力场存在振荡。而密度耗散项是加在连续性方程中，通过每一个时间步调节密度的振荡来控制压力的振荡，因此密度耗散项算法的压力场更为稳定。

3.2.2 角速度较高情况

对于工况b，流动情况为中、高速状态，横摇周期较短。为便于比较长期流场压力的变化情况，总物理时间取为50 s。另外，在距离液面以下0.1 m 的左侧箱体壁面位置处设置一个压力监测点，用来监测该点处水压随时间的变化；在距离水箱左侧壁面0.05 m 的位置设置一个液面高度监测点，用来实时记录该点处水面高度随时间的变化。

图7 示出了2 种SPH 算法计算得到的监测点压力随时间变化的曲线。由于初始时刻横摇激励的突然施加，压力的数值结果呈现振荡不稳定性。为了截取稳定有效的压力结果，压力值从t=2 s 开始记录。从图中可以发现，在整个数值模拟过程中，2 种耗散项算法得到的两条曲线的相位分布一致，压力的幅值近似一致，曲线的趋势基本吻合，这表明2 种算法具有几乎同等的有效获得监测点压力的能力。

图7 监测点压力时历曲线Fig.7 Time histories of pressure at the monitoring point

图8 所示为监测点在t=20～30 s 的压力曲线。从中可以很明显地看出2 种耗散项算法在计算压力中的稳定性和振荡程度，2 种算法计算的压力均有一定程度的振荡，但结果相似，表明在该组工况下2 种算法具有同等抑制压力振荡的效果。

图8 监测点在t=20～30 s 的压力时历曲线Fig.8 Time histories of pressure at the monitoring point（t=20～30 s）

为方便比较2 种耗散项算法在周期性运动情况下对流场压力振荡的抑制效果，选取了t=2.75T，12.25T，17.5T 和24.25T 这几个时间节点的流场压力分布图，如图9（a）～图9（d）所示。各分图中，左图和右图分别为密度正则化算法和δ-SPH 密度耗散项算法所得的流场压力分布。从流场压力的对比中发现：2 种算法获得的流场压力分布均匀，随着时间的增长，流场压力的光滑程度没有发生太大的变化，这表明2 种算法都能够长时间起到抑制压力振荡、稳定数值结果的作用。

2 种算法结果趋于一致的原因在于：在转动角速度较高的情况下，速度的变化量相对较大，黏性开始起效，虽然2 种算法的精度不同，但获得的结果相似。

此外，在数值模拟过程中数值耗散的一个直观表现是自由液面高度的变化。为了定量比较2 种方法在自由液面高度上的耗散程度，需要事先标记处在自由液面上的粒子，然后基于线性插值的方式计算监测点位置的液面高度。自由液面粒子可用式（31）计算并标记。

由式（31）不难发现，当粒子处于自由液面时，粒子的支持域在自由液面处被截断，因此该式的结果是一个比1 小的数，这里取为0.9。

在标记的自由液面粒子中，基于线性插值的方式，计算监测点处的液面高度，其计算方法［22］的表达式为

图10（a）记录了距离左侧壁面0.05 m 处自由液面高度随时间的变化情况。分析2 条曲线发现：因选取的角频率远离水箱共振频率，液面高度曲线呈现出周期性的特点，中后期基本上呈现出平稳波动的规律，2 条曲线在整个计算过程中能够达到90%以上的拟合程度，这在一定程度上表明了2 种算法在长时期的周期性水动力学问题中耗散程度基本持平。

图10（b）所示为t=20～30 s 的自由液面高度随时间变化的曲线，可见2 种算法得到的曲线相位和趋势一致，这表明2 种算法得到的流动状态同步，2 条曲线的幅值基本一致，表明2 种算法得到的自由液面形态和高度基本吻合，进一步证明了2 种算法在较高速晃荡运动中能够达到同样的模拟效果。

图10 自由液面高度时历曲线Fig.10 Time histories of the height of free surface

4 结 论

本文不考虑人工黏性，仅考虑密度耗散，比较并分析了2 种SPH 密度耗散算法对流场压力振荡的抑制效果。对二维液舱的静水和晃荡问题进行数值模拟，探究了不同角速度情况下2 种SPH 耗散项算法在数值结果上的差别。得到如下主要结论：

1）当液舱处于静止状态时，δ-SPH 密度耗散项算法比Shepard 密度正则化算法更能够在长时间的模拟中稳定流场压力，保证动能的严格守恒。

2）当液舱在较低角速度下晃荡时，Shepard 密度正则化算法仅能在有限的时间内过滤压力噪声，而δ-SPH 密度耗散项算法则能够在长时间的模拟中光滑压力场，后者更适用于长时间的数值计算。

3）当液舱在较高角速度下晃荡时，Shepard 密度正则化算法和δ-SPH 密度耗散项算法在数值模拟中均能够持续发挥过滤压力噪声、抑制压力振荡的作用，两者抑制压力振荡的程度相当。

总之，当液舱在角速度较低和角速度较高情况下晃荡时，δ-SPH 密度耗散项算法均可有效抑制压力振荡，建议在SPH 模型中采用该密度耗散项算法。