去奇异边界元方法在液舱晃荡模拟中的应用

王庆丰，徐 刚，王树齐，朱仁庆

(江苏科技大学 船舶与海洋工程学院，江苏 镇江 212003)

在船舶领域的应用中，频域理论已经被运用的相当成熟，但是频域法通常只适用于周期稳态问题，即它存在一个本质上的限制就是对于瞬变或者强非线性问题往往无法处理。而频域理论的超集实际上是时域理论，它是一种不可分离时间项的方法，理论上可以解决物体任意运动的问题，还能解决完全非线性问题，在自由度方面存在明显的优势。本文采用时域边界元方法进行数值模拟，但是对于边界元方法，奇点的处理一直是关注的重点，传统方法在交界面处的速度精度很难保障，为了消除积分的奇异性，学者们已经研究出诸多成熟可用的方法，其一径向积分边界元法[1-2]，它基于纯数学处理技术，不借助任何特解和微分算子可将任何问题的域积分转化为边界积分，是一种无内部网格的纯边界元算法，目前已应用于断裂力学、热辐射等方面；其二就是无奇异边界元方法，即将传统方法中原本直接布置于流体计算域表面网格中心点上的所有奇点，移至计算域外部，可以减少存储需求并提高计算效率，目前该方法主要应用于船舶与海洋工程水动力学分析方面。为了能够更好地模拟液舱内液体的运动，本文采用基于去奇异边界元方法(Desingularized Boundary Integral Equation Method,DBIEM)研究了液舱晃荡的基本数值解法。

传统的边界元积分方法，其奇点是分布在计算域的表面，在进行相关数值计算的过程中需要对奇点积分进行特殊的处理，表现在时域分析非线性问题的时候这过程中会出现大量繁琐的计算分析过程。当我们采用无奇异边界元方法的时候会得到两个优点[3]：① 可以得到控制面上高精度的计算结果；② 获得离散边界积分问题代数系统的计算时间变少了。同时，该方法简单易实施。Jensen 等利用去奇异方法来分析在稳定入射波情况下求解阻力的问题[4]。Lee[5]将去奇点边界积分方程和时间步进技术进行结合，模拟了简单潜体及浮体运动的完全非线性效应问题；李宗翰等[6]将去奇异积分方程方法运用到任意三维浮体的完全非线性流场的模拟和仿真中；Scorpio等[7-8]对完全非线性数值水池进行了探索；Kara等[9]利用去奇异边界元方法在三维时域内研究了完全非线性波物相互作用的问题；Zhang等[10]基于去奇异方法研究了双体船的耐波性；Xu等[11]进行了随机波浪下三维液舱的全非线性数值分析。

1 去奇异边界元理论

1.1 控制方程及初始条件

根据势流基本理论和图1液舱晃荡系统，速度势φ在流体域D中满足Laplace方程，定解条件如下

(1)

(2)

(3)

(4)

式中:η代表波高;x，y，z代表坐标轴三个方向上的物理量。初始条件条件如下(φ0是初始速度势、ξ是初始波高)

φ(x,y,z=ξ,t=0)=φ0(x,y,z)

(5)

η(x,y,t=0)=ξ(x,y)

(6)

图1 液舱晃荡坐标系统Fig.1 The coordinate system of sloshing tank

1.2 自由面条件积分格式

有限差分法一般被运用在自由表面上的速度势φ的求解，有限差分法的不足是求解多维问题时不易适应不规则几何边界，并且在求解过程中会常常出现数值误差的累积，从而导致计算结果发散，影响数值的稳定性。本文采用一种新的积分格式的自由面条件(Integral Form of Free Surface Boundary Condition[12]，IFBC )，运用此积分格式的自由面条件，能够在线性自由面条件的基础上得到自由面上一阶速度势的表达式，而且这种积分格式的自由面条件有相对较好的稳定性，可以在很大程度上降低总体误差，能防止出现结果发散的现象。参照文献[13]中方法，在式(3)和(4)的基础上，对线性自由面条件采用先积分后离散的处理方法，首先两边对t在(0,τ1)积分，然后再对τ1在(0,τ)中积分，得到

(7)

再由初始条件，上式可以写为

(8)

然后采用如下的离散模型将式(8)进行离散：在空间上把自由面划分成若干小平面单元，并假定每个单元中心点上的物理量为常数分布；在时间上以Δt为步距，等间隔前进。采用梯形公式离散式(8)，可以得到

(p∈SF)

(9)

1.3 去奇异边界元积分方法

采用边界元方法的关键一步是积分方程的离散，本文采用分布源法，则流场中的速度势可以表示为

(10)

因此式(2)和(9)中提到的边界条件可改写成如下形式

φ(p,t)=φ0(p,t)

(11)

(12)

将式(11)和(12)代入式(10)，那么关于未知量σ0(q,t)的积分方程可表示为

(13)

(14)

对于式(13)、式(14)的积分方程，传统的边界元方法都是将源强直接分布在物体及自由表面上，也就是说积分表面S1即为流域的边界。这样往往会出现由于格林函数1/rpq的分母趋于零所带来的奇异积分的问题。本文采用的方法为去奇异边界元方法，即是将传统作法中原本直接布置于流体计算域表面节点上的所有奇点，移至计算域外面，也就是将传统作法中叠在一起的节点和奇点分开了。在积分边界SF和SW每个单元上分别布置强度为σF和σW的点源，SF和SW是实际流体域外距离流体域很近的积分表面，则：

(15)

将式(13)式(14)代入式(15)，则：

(p∈ΓF)

(16)

(17)

假设研究的问题中总共有N个孤立的点源构成，则式(16)和(17)可简化成如下的N个孤立点源的代数和的形式

(18)

(19)

当采用上述离散方法对总共N个配置点列写方程后，可最终得到如下N×N的线性代数方程组

Aijσj=bi

(20)

式中：Aij为影响系数矩阵；bi为边界条件所形成的确定矩阵；σj为待求的未知源强。

线性代数方程组(20)求解后即得出未知源强，由式(9)求得速度势，可以通过下式求得水动压力p、波高η和水动力F

(21)

1.4 基于DBIEM晃荡程序的流程

基于DBIEM理论，以Fortran语言为工具，本文编写了一套能够模拟任意尺寸任意形状液舱晃荡的程序，图2是晃荡程序的基本流程图。

2 数值方法有效性验证

2.1 水平激励下液舱晃荡解析解

对于周期性激励ue=Acosωt，其中ue是液舱激励速度，A=bω速度幅值，b为位移幅值，ω为激励频率，Faltinsen给出了求解速度势φ的线性解析方法，通过下式可以很容易的将求得的φ值转化为自由液面的位置分布[14]。(液舱水深h、长度2a、宽度2w，且坐标原点位置在自由面上液面的中心处。)

图2 晃荡程序基本流程Fig.2 The basic process of sloshing program

(22)

2.2 单向水平激励下晃荡数值解法验证

液舱的长(L) 和宽(B)分别为1 m和0.5 m；液舱内液体高度(h)为0.5 m。固有频率ω0x和ω0y分别为5.316 rad/s和7.835 rad/s。为验证程序是否准确，选取合适的网格大小、时间步长以及去奇异距离(如表1所示)，模拟不同激励频率作用下液舱晃荡问题。

表1 主要计算参数Tab.1 The main calculation parameters

表1中，Nx为液舱长度方向的划分网格个数；Ny为液舱宽度方向的划分网格个数；Nz为液舱高度方向的划分网格个数；dt为时间步长；T为激励周期；ld为去奇异距离；A为激励幅值。

考虑到液面条件为线性的，单向激励幅值取为0.01h，h为液舱内液面高度，入射频率分别为0.5ω0(ω0为液舱在激励方向上的固有频率，此处即为液舱的纵向固有频率)、0.9ω0、0.95ω0、0.999 9ω0、1.05ω0、1.1ω0。计算结果如图3所示，图中实线为数值解，虚线为解析解。

(a) ωP=0.5ω0(b) ωP=0.9ω0

(c) ωP=0.95ω0(d) ωP=0.999 9ω0

(e) ωP=1.05ω0(f) ωP=1.1ω0

图3 不同激励频率下计算点波高时历曲线
Fig.3 Time history of free surface elevation for different frequencies

可以看出，各激励频率下，数值解都能够很好的与解析解吻合，误差仅在1%以内，这有效说明了去奇异边界元法在求解液舱晃荡问题的有效性。由以上各图还可以看出计算点处波高时历曲线随着激励频率的不同其变化规律也发生了变化，外部激励频率远离液舱固有频率时如图3(a)所示，曲线变化较快，有一定规律性，波高最大值在0.006 m左右；当入射频率接近液舱固有频率时，计算点处波峰首先随时间不断增加达到最大值后又开始下降，并以此规律不断循环，有很强的规律性，图3(b)、(c)、(e)、(f)都显示了这一特征，且波高最大值分别为0.06 m、0.11 m、0.15 m、0.1 m。不难看出越是接近液舱固有频率波高最大值就越大，这是与共振密切相关的，理论上当激励频率无限接近共振频率时，如图3(d)所示频率为0.999 9ω0，随着计算时间的推移该点处的波高将趋于无穷大，但由于程序算法和网格的限制，本文仅模拟了35 s，此时的最大幅值已经约是图3(a)的80倍。可见，为了保证所设计液舱的安全性，应尽量避免出现共振现象，基于去奇异边界元方法模拟液舱晃荡问题对实际工程问题具有一定的指导作用。

3 不规则水平激励作用下液舱晃荡分析

在单向水平激励下液舱晃荡问题研究的基础上，本文也研究了不规则激励作用下的液舱晃荡，如图4所示。

图4 液舱晃荡模型图Fig.4 Model of sloshing tank

3.1 单向4成分叠加不规则激励

首先采用四个不同于液舱固有频率的单向激励进行研究，激励参数选取如表2所示。

表2 四成分不规则激励参数Tab.2 Parameters for irregular excitation with four

图5为单向不规则激励下点(-L/2,0)处波高时历和相对体积误差曲线，此处不规则激励由四种不同单色激励叠加而成。与单色激励相比，不规则激励下波高时历曲线比较复杂，波高变化范围为-0.004～0.004 m。相对体积误差几乎为0，符合质量守恒物理规律。

(a) 波高时历曲线(b) 相对体积误差

图5 不规则激励下波高时历和相对体积误差曲线
Fig.5 Wave elevation history and relative volume error

for irregular excitation

图6液舱波面图，研究时间区间为13.252 0～16.345 9 s，约为一个波面运动周期，每隔Δt≈0.52的6个不同时刻波面图，从图中可以看到一个周期内波面的基本变化情况。

3.2 单向512成分叠加不规则激励

通常基于波浪谱，一般200个左右规则波浪成分就已经足够模拟不规则波。为了更好的验证本文的数值方法，文中将采用512个成分的不规则波。本文将选用ITTC波浪谱研究随机激励作用下的液舱晃荡情况，波浪谱如下

图6 波面图Fig.6 Free surface profiles

(23)

式中：Hs为有义波高；ωp为谱峰频率。与随机波浪组成相似，液舱运动速度和幅值写成

(24)

(25)

式中：Nω为激励成分的个数，本文取512；φi为0～2π间随机相位角，在程序中设定；ωi为激励频率。

选取有义波高Hs=0.005 h，谱峰频率为ωp=0.9ω0，截断区域上限选取ωi≤5ω0波浪谱部分。

图7为单向512成分不规则激励下点(-L/2，0)处波高时历和相对体积误差曲线。其中，每个成分的波幅、频率都不一样，初始相位采用随机数。模拟结果显示，波高变化范围为-0.015～0.015 m。由于模拟中使用了随机相位，因此无法得到准确的解析解，但根据相对体积误差图可以看出512成分不规则波激励下模拟准确度仍然很高，体积误差接近于0。

(a) 波高时历曲线(b) 相对体积误差

图7 不规则激励下波高时历和相对体积误差曲线
Fig.7 Wave elevation history and relative volume error for irregular excitation

图8是时间区间为13.252 0～16.345 9 s，每隔Δt≈0.52的6个不同时刻波面图。

图8 波面图Fig.8 Free surface profiles

4 结 论

针对传统边界元方法求解波浪中结构物运动响应时会遇到奇点及交界面处速度分布精度较低的问题，本文引入去奇异边界元理论，利用Fortran语言开发了一套可以针对不同形状液舱晃荡模拟的程序，程序中加入了各自由度运动耦合的模块，可模拟液舱在六个自由度激励联合作用下内部液体晃荡各要素随时间变化情况。本文主要对单自由度4成分不规则激励、单自由度512成分不规则激励下液舱内部流体运动进行了模拟，研究结果表明，文中方法可有效模拟随机激励下的液舱晃荡问题，且具有较高的精度。