基于高精度Boussinesq方程的三维浅水晃荡数值研究

2021-06-01袁心怡刘祖源

上海交通大学学报 2021年5期

袁心怡，苏焱，刘祖源

(武汉理工大学交通学院，武汉 430063)

晃荡可以描述为部分填充的容器中液体的运动.液化天然气(LNG)船、浮式液化天然气装置(FLNG)等采用液舱形式的液货船在部分装载的情况下，即使受到小幅激励，也有可能发生剧烈的晃荡现象，对液舱结构产生不利影响，威胁船舶航行安全.

研究晃荡问题多采用数值方法，基于黏性流理论在流域内求解Navier-Stokes方程或基于势流理论在流域内求解Laplace方程实现数值模拟.Wang等[1]基于比例边界有限元法开发了一种研究挡板对舱内晃荡影响的新方法，经过收敛性和研究，发现该方法具有计算简单、精度高、速度快等优点.Zhao等[2]基于势流理论建立了LNG液舱在强制激励下的非线性晃荡数值模型，并加入了人工阻尼模型考虑耗能，研究了部分充液的矩形容器在平移和旋转激励下的晃荡特性.Hu等[3]利用边界元法求解控制方程，采用线性自由表面条件研究了多种形状带挡板的二维液舱在小幅激励下晃荡的固有频率和振型，并提出了带挡板的矩形液舱内晃荡的固有频率计算经验公式.Ünal等[4]研究了带T形挡板的液舱内的晃荡现象，考察了挡板高度、充液深度等参数对于舱壁水动力载荷和自由表面高程的影响，计算采用了两种不同的方法，分别使用黏性的层流和紊流求解器.李金龙等[5]对一种新的几何有限体积法isoAdvector进行修正，使之可应用于液舱晃荡的模拟中，并基于不同的有限体积法对各种工况进行了模拟，发现修正后的方法所获得的自由面位置和整体水动力载荷载荷精度更高，且能较好地模拟波面的翻卷和破碎.Fu等[6]提出了一种基于局部径向基函数和二阶显式龙格库塔法的半拉格朗日无网格模型来分析二维舱内的晃荡现象，通过文献对比证明了该方法的有效性和准确性.

对于强非线性的晃荡问题，粒子法可以很好地捕捉自由面的大变形.但对于低充液率的容器，稳定和精确地模拟容器内长时间的晃荡具有难度.Green等[7]提出了一种有效的光滑粒子流体动力学(SPH)方法来克服这些困难，使边界条件、耗散项和算法的误差最小化，与试验数据比较发现，数值模拟的结果非常精确.

布西内斯克(Boussinesq)方程是一类描述波浪传播变形的数学模型.Bingham等[8]建立了数值离散方法，主要用于近岸变化底部处波浪场精确模拟.Antuono等[9]使用Boussinesq型平均水深方程建立了一个二维模型用于分析浅水中液体晃荡现象.Su等[10]采用高精度Boussinesq型速度势方程，对二维矩形舱内的晃荡现象进行了数值模拟，结果表明基于Boussinesq型方程的数值模型在预测小激励幅值晃荡方面具有较好的性能.

本研究旨在探讨浅水、微幅激励条件下三维液舱内的晃荡现象.假定液舱内为无黏无旋不可压流体，流域内存在速度势.将总速度势分解为两部分，一部分为满足拉普拉斯方程和壁面边界条件的特解，另一部分采用基于高精度速度势型Boussinesq方程的数值模型求解.首先给出三维液舱中流体受强迫运动的控制方程和边界条件，然后基于有限差分法对自由面进行空间离散，建立数值模型.通过与文献中二维晃荡的结果进行对比验证了数值模型的有效性，并对三维工况的结果进行分析.

1 数学模型

考虑一个长为l、宽为b、静止水深为h的液舱，如图1所示，固定坐标系O-xyz位于静止水面上，原点位于静水面的几何中心，z轴垂直向上.液舱底部记为z=-h(x,y)，流体自由表面记为z=η(x,y).

图1 晃荡运动模型Fig.1 Sloshing motion model

假定舱内为无黏无旋不可压流体，基于势流理论，流域内存在一个速度势Φ.记水平速度为u，垂向速度为w，则速度势与速度的关系为

(1)

ηt+Φxηx+Φyηy-Φz=0,z=η

(2)

(3)

(4)

∂Φ/∂n=v·n

(5)

其中，式(2)为自由面运动学边界条件；式(3)为自由面动力学边界条件；式(4)为拉普拉斯方程，Xij(i,j=x,y,z)表示物理量X对i方向求一阶导数，对j方向求二阶导数；式(5)为物面不可穿透条件；g为重力加速度；v=[vxvyvz]为液舱的运动速度；n为物面的法向量，由液舱内指向外.对于惯性系中定义的自由面边界条件，利用如下关系将其转换到随液舱运动坐标系中：

则可以给出随舱运动坐标系下，自由面满足运动学和动力学边界条件：

(8)

(9)

由于液舱的运动始终是已知量，在数值计算时，将总速度势分为两部分：Φ=φ+φ，其中φ是一个满足拉普拉斯方程和物面不可穿透条件的特解.仅考虑横荡、纵荡、垂荡运动时，可以直接写出φ的表达式如下：

φ=xvx+yvy+zvz

(10)

然后由叠加原理可导出待求解的分速度势φ所满足的物面边界条件和控制方程：

(11)

将Φ=φ+φ以及式(10)引入到式(8)、(9)中，可以得到φ满足的自由面边界条件如下：

ηt+ηxφx+ηyφy-φz=0

(12)

(13)

2 数值计算方法

(14)

则可用新变量表示自由面条件并整理如下：

(15)

(16)

(17)

式中：参数δ≪1，表示展开面是缓慢变化的.由上述递推关系式可推导速度势和垂向速度的表达式如下：

(18)

(19)

将表达垂向速度的无穷级数截断在3阶时，以上两式中：

由此可得到速度势和垂向速度的表达式，将其代入到底部运动学边界条件a1、a3的表达式中得：

(25)

整理后可得底部方程的表达式为

(26)

数值计算过程中，空间导数计算采用有限差分法，时间步进采用四阶龙格库塔方法.计算程序基于MATLAB语言编写，求解线性方程组使用MATLAB中的直接求解器，该求解器可根据系数矩阵的特征选择高效的求解方法.

3 二维工况验证

为验证数值模型的有效性，将三维液舱的宽度设置为小量并在相同工况下与文献中的二维结果[10]进行对比.三维液舱的宽度b=0.1 m，长度和充液高度都与二维工况相同，分别为l=1.175 m，h=0.06 m.液舱受长度方向的正弦激励做纵荡运动，激励幅值A=0.39 cm，记线性理论预测的共振频率ωr=2.042 5 s-1，激励频率ω′分别为0.98ωr、1.02ωr及1.08ωr.

图2依次展示了激励频率为0.98ωr、1.02ωr及1.08ωr时左壁面处的液面高度变化对比图，图中：t′为时间与激励周期之比，η′为自由面高度与充液高度之比.3个工况下的波高变化都呈现出强非线性现象，波峰陡峭但波谷平坦.当外部激励频率处于共振频率附近时，矩形液舱内的晃荡运动形式对外部激励频率的变化非常敏感，在3种工况下观察到了完全不同的运动形式，稳定状态下的一个变化周期内分别包含3个波峰、2个波峰和单个波峰.相应地，在液舱内可以观察到对应数量的行进波在舱内来回移动.

图2 液面高度时历对比图Fig.2 Comparisons of time history of surface height

可以看出，三维数值计算程序能很好地捕捉自由面高度的变化和晃荡运动形式对于外部激励频率变化的敏感性.

4 三维工况

计算选择长宽相等的方型液舱，将二维液舱的长度缩小一半，三维液舱模型参数为l=b=0.587 5 m，h=0.03 m.仅考虑纵荡、横荡或二者混合的正弦形式运动，A=0.1 cm，线性理论预测的舱内晃荡运动共振频率为ωr=2.888 1 s-1.考虑到二维情况下，外部激励频率在共振频率附近变化时，晃荡运动会出现特定的波高变化形式，三维工况计算也在共振频率附近进行.

数值计算过程中，将静止液面划分成了正方形网格，x和y方向都有25个节点.在二维液舱中，波高观测点为左壁面处第1号节点，而在三维液舱中，波高观测点为混合运动合成矢量方向上的第1号节点.

如图3所示，首先选取了ω′=0.98ωr，将二维工况和三维单方向运动工况进行对比.可以看出，两种情况下自由面的高度的变化幅值和周期都是一致的，液舱的宽度并没有对舱内流体的运动产生影响.接着将三维单方向运动工况和混合运动工况进行了对比，从图中可以看出，两种工况下波高变化的周期是一致的，将单方向运动的波高加倍以后，幅值的大小也有很好的一致性.

如图4所示，计算了其他激励频率下，混合运动工况中波高观测点处的波高变化.为了保证舱内的流体运动达到稳定状态，计算时间超过200个激励周期.在1.02ωr激励频率下，观察到了和二维工况中类似的现象：在一个变化周期内出现两个波峰.但是单波峰的现象却出现在1.06ωr激励频率下，并且在1.08ωr激励频率下，观测点波高变化的非线性减弱，说明三维情况下出现强非线性行进波的激励频率域有所缩短.

图3 0.98ωr下波高对比图Fig.3 Comparisons of wave height at 0.98ωr

图4 各频率下波高变化图Fig.4 Change of wave height at different frequencies

数值计算程序可以实时显示自由面变化，以1.06ωr激励频率为例，图5显示了数值模拟过程中某一时刻整个自由面形状.在各个运动方向上，各自形成了行进波，两条波浪相遇的区域形成了一个波峰点.类似地，在1.02ωr激励频率下，会在各个方向上观察到两条行进波，且自由面上观察到2×2个波峰点.而在0.98ωr激励频率下，会在自由面上观察到9条行进波和3×3个波峰点.