上海交通大学学报(自然版), 2021, 55(5): 521-526 doi: 10.16183/j.cnki.jsjtu.2020.053

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

袁心怡, 苏焱,, 刘祖源

武汉理工大学 交通学院, 武汉 430063

Numerical Investigation of Three-Dimensional Shallow-Water Sloshing Based on High Accuracy Boussinesq Equations

YUAN Xinyi, SU Yan,, LIU Zuyuan

School of Transportation, Wuhan University of Technology, Wuhan 430063, China

通讯作者: 苏 焱,男,副教授;E-mail:suyan@whut.edu.cn.

责任编辑: 陈晓燕

收稿日期: 2020-02-28  

基金资助: 国家自然科学基金资助项目(51609187)

Received: 2020-02-28  

作者简介 About authors

袁心怡(1995-),男,湖北省荆州市人,硕士生,研究方向为液舱晃荡. 。

摘要

基于高精度速度势型布西内斯克(Boussinesq)方程对三维液舱内的浅水晃荡现象进行模拟.研究在势流理论框架下进行,总的速度势被分成了两个部分:一部分是在流域内满足拉普拉斯方程且在边界处满足不可穿透条件的特解,而另一部分将由Boussinesq模型求解.数值计算过程中,空间导数离散采用有限差分法,时间步进采用四阶龙格库塔法.首先将三维液舱的长宽比设置为远小于1以模拟二维液舱情况,并与文献的结果对比,验证了数值模型的有效性.三维工况中,各个外部激励频率下,观察到了4种不同的晃荡运动形式,同时在自由面上会观察到相应数量的行进波.此外,讨论了外部激励频率和耦合激励形式对于液舱内晃荡运动形式的影响.

关键词: 布西内斯克方程; 有限差分法; 势流理论; 晃荡; 浅水

Abstract

Highly accurate Boussinesq-type equations in terms of velocity potential are used for the simulation of shallow-water sloshing in a three-dimensional tank under the framework of the potential flow theory. The total velocity potential is separated into two parts: one part is a particular solution which satisfies the Laplace equation in the fluid domain and the no-flow condition on the walls while the other part is solved by the Boussinesq-type model. In the process of numerical calculation, the finite difference method is used for spatial derivative discretization and the 4th Runge-kutta method is used for time iteration. To verify the numerical model, the aspect ratio of the tank is set to be much less than 1 for simulation of 2D cases and is compared with the results published. In the 3D cases, four different sloshing motion forms are observed at each external excitation frequency, and a corresponding number of traveling waves are observed on the free surface. Moreover, the effects of external excitation frequency and coupling excitation on the sloshing motion in the tank are discussed.

Keywords: Boussinesq equations; finite difference method; potential flow theory; sloshing; shallow-water

PDF (1430KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

袁心怡, 苏焱, 刘祖源. 基于高精度Boussinesq方程的三维浅水晃荡数值研究[J]. 上海交通大学学报(自然版), 2021, 55(5): 521-526 doi:10.16183/j.cnki.jsjtu.2020.053

YUAN Xinyi, SU Yan, LIU Zuyuan. Numerical Investigation of Three-Dimensional Shallow-Water Sloshing Based on High Accuracy Boussinesq Equations[J]. Journal of shanghai Jiaotong University, 2021, 55(5): 521-526 doi:10.16183/j.cnki.jsjtu.2020.053

晃荡可以描述为部分填充的容器中液体的运动.液化天然气(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

图1   晃荡运动模型

Fig.1   Sloshing motion model


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

u=ΔΦ, w=Φz

式中:Δ≡[∂/x/y];Xi(i=x,y,z)表示物理量Xi方向的导数.流体满足的控制方程和边界条件如下:

ηt+Φxηx+Φyηy-Φz=0, z=η
Φt+ 12(ΔΦ)2+ 12(Φz)2+gz=0, z=η
Δ2Φ+Φzz=0, -h<z<η
Φ/n=v·n

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

η/t:η/t-v·Δη
Φ/t:Φ/t-v·ΔΦ-vzΦz

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

ηt-v·Δη+Φxηx+Φyηy-Φz=0
Φt-v·ΔΦ-vzΦz+ 12(ΔΦ)2+ 12(Φz)2+gz=0

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

φ=xvx+yvy+zvz

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

Δ2ϕ+ϕzz=0,-h<z<ηϕx=0,x=±l/2ϕy=0,y=±b/2ϕz=0,z=-h

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

ηt+ηxϕx+ηyϕy-ϕz=0
ϕt+ xv·x+ yv·y+ zv·z++ 12(ϕx2+ ϕy2+ ϕz2)=0

式中: v·表示该速度分量对时间的导数.上述自由面方程中独立于空间坐标系的项被忽略了,由此得到ϕ所满足的所有条件,下文将由Boussinesq模型求解ϕ.

2 数值计算方法

根据Zakharov[11]的研究,引入直接定义在自由表面的速度势 Φ~=(x,y,η,t)和垂向速度 w~=(Φz)z=η.基于链式求导法则, Φ~Φ的空间、时间导数存在如下关系:

Φ~t=(Φt)z=η+(Φz)z=ηηtΔΦ~=(ΔΦ)z=η+Δη(Φz)z=η

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

ηt=vxηx+vyηy-ηxΦ~x-ηyΦ~y+ (1+ ηx2+ ηy2) w~
Φ~t=vxΦ~x+vyΦ~y+vzw~--0.5(Φ~x2+ Φ~y2)+0.5(1+ ηx2+ ηy2) w~2

以上两个式子中可用于时间步进求解 Φ~η,但还需要一个计算未知项 w~的方法.由速度势分解的关系可知 w~=ϕz+φz,而φ是已知的,若求出了ϕ则可以求得 w~.

根据Bingham等[8]的研究,将分速度势在z= z^处进行泰勒级数展开并代入拉普拉斯方程,略去高阶小量后进行化简可得到如下递推关系:

ϕ^(n)=-(1-δ2Δ z^2)Δ2ϕ^(n-2)+2δΔ z^·Δ ϕ^(n-1)+δ2Δ2z^ϕ^(n-1)

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

ϕ= ϕ^*+J01Δ2ϕ^*+a1z^(Δ z^)·(Δ ϕ^*)+J11p[z^x(ϕ^xxx*+ ϕ^yyx*)+ z^y(ϕ^xxy*+ ϕ^yyy*)]+ψw^*+J02Δ2w^*+(ψ2+b1ψz^)(Δ z^)·(Δ w^*)+J12p[z^x(w^xxx*+ w^yyx*)+ z^y(w^xxy*+ w^yyy*)]
w=-ψΔ2ϕ^*-J02(ϕ^xxxx*+2 ϕ^xxyy*+ ϕ^yyyy*)-J12w[z^x(ϕ^xxx*+ ϕ^yyx*)+ z^y(ϕ^xxy*+ ϕ^yyy*)]+ w^*+J01Δ2w^*+(2ψ+ z^b1)(Δ z^)·(Δ w^*)+J11w[z^x(w^xxx*+ w^yyx*)+ z^y(w^xxy*+ w^yyy*)]

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

J01=-ψ2/2+ z^2/10J02=-ψ3/6+ z^2ψ/10
J11p=-ψ3/3-ψ2z^(1/5+a1/2)+ z^3a3
J12p=-ψ4/6+ψ2z^2/10-ψ3z^(1/15+b1/6)+ψz^3b3
J12w=ψ2+ψz^(2/5+a1)
J11w=-2ψ3/3+ψz^2/5-ψ2z^(1/5+b1/2)+ z^3b3

a1=-0.465737, a3=-0.0653131

b1=0.138568, b3=-0.0138568

式中:ψ=z- z^, ϕ^*w^*为展开面处的速度势和垂向速度.

由此可得到速度势和垂向速度的表达式,将其代入到底部运动学边界条件a1a3的表达式中得:

w+(Δh)·(Δϕ)=0

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

w+(Δh)·(Δ ϕ^*)+J01[hx(ϕ^xxx*+ ϕ^yyx*)+hy(ϕ^xxy*+ ϕ^yyy*)]+ψh)·(Δ w^*)+J02[hx(w^xxx*+ w^yyx*)+hy(w^xxy*+ w^yyy*)]=0

将式(18)、(26)联立即可形成线性方程组, 由于自由面分速度势 ϕ~在初始条件中是已知的,因此上述方程组可由初始条件以及边界条件求解未知项 ϕ^*w^*,进而求解φ.

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

3 二维工况验证

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

图2   液面高度时历对比图

Fig.2   Comparisons of time history of surface height


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

4 三维工况

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

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

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

图3

图3   0.98ωr下波高对比图

Fig.3   Comparisons of wave height at 0.98ωr


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

图4

图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个波峰点.

图5

图5   自由面形状图

Fig.5   Shape of free surface


5 结论

本研究基于高精度Boussinesq型速度势方程对二维和三维液舱中在小幅激励下的浅水晃荡现象进行了数值模拟.主要结论有:

(1) 首次采用Boussinesq方程实现了三维浅水晃荡模拟,建立的数值模型能很好地捕捉三维液舱内的强非线性晃荡运动以及晃荡运动形式对于外部激励频率变化的敏感性,且具有较高的计算效率.

(2) 三维数值模拟中,观察到了4种不同的晃荡运动形式.在各个工况下,会观察到不同数量的行进波,行进波数量取决于外部激励的频率.在行进波十字交汇的区域,会叠加产生波峰.

(3) 在靠近共振频率的工况下,即使是微幅激励的耦合,都在舱内产生了大幅的波浪.对于远离共振频率的工况,相同外部激励频率下,与二维相比,三维情况下的非线性有所减弱.

(4) 目前数值模型仅考虑单方向的平动或两方向耦合运动,未来将对其他运动形式(单方向转动或耦合)作进一步的研究.

参考文献

WANG W Y, PENG Y, ZHANG Q, et al.

Sloshing of liquid in partially liquid filled toroidal tank with various baffles under lateral excitation

[J]. Ocean Engineering, 2017, 146:434-456.

DOI:10.1016/j.oceaneng.2017.09.032      URL     [本文引用: 1]

ZHAO D, HU Z, CHEN G, et al.

Nonlinear slo-shing in rectangular tanks under forced excitation

[J]. International Journal of Naval Architecture and Ocean Engineering, 2018, 10(5):545-565.

DOI:10.1016/j.ijnaoe.2017.10.005      URL     [本文引用: 1]

HU Z, ZHANG X Y, LI X W, et al.

On natural frequencies of liquid sloshing in 2-D tanks using Boundary Element Method

[J]. Ocean Engineering, 2018, 153:88-103.

DOI:10.1016/j.oceaneng.2018.01.062      URL     [本文引用: 1]

ÜNAL U O, BILICI G, AKYILDIZ H.

Liquid slo-shing in a two-dimensional rectangular tank: A numerical investigation with a T-shaped baffle

[J]. Ocean Engineering, 2019, 187:106183.

DOI:10.1016/j.oceaneng.2019.106183      URL     [本文引用: 1]

李金龙, 尤云祥, 陈科.

一种几何VOF方法在液舱晃荡流动模拟中的应用

[J]. 上海交通大学学报, 2019, 53(8):943-951.

[本文引用: 1]

LI Jinlong, YOU Yunxiang, CHEN Ke.

Application of a geometric VOF method in the simulation of sloshing flow

[J]. Journal of Shanghai Jiao Tong University, 2019, 53(8):943-951.

[本文引用: 1]

FU Z J, ZHANG J, LI P W, et al.

A semi-Lagrangian meshless framework for numerical solutions of two-dimensional sloshing phenomenon

[J]. Engineering Analysis with Boundary Elements, 2020, 112:58-67.

DOI:10.1016/j.enganabound.2019.12.003      URL     [本文引用: 1]

GREEN M D, PEIRÓ J.

Long duration SPH simulations of sloshing in tanks with a low fill ratio and high stretching

[J]. Computers & Fluids, 2018, 174:179-199.

DOI:10.1016/j.compfluid.2018.07.006      URL     [本文引用: 1]

BINGHAM H B, MADSEN P A, FUHRMAN D R.

Velocity potential formulations of highly accurate Boussinesq-type models

[J]. Coastal Engineering, 2009, 56(4):467-478.

DOI:10.1016/j.coastaleng.2008.10.012      URL     [本文引用: 2]

ANTUONO M, BOUSCASSE B, COLAGROSSI A, et al.

Two-dimensional modal method for shallow-water sloshing in rectangular basins

[J]. Journal of Fluid Mechanics, 2012, 700:419-440.

DOI:10.1017/jfm.2012.140      URL     [本文引用: 1]

SU Y, LIU Z Y.

Numerical model of sloshing in rectangular tank based on Boussinesq-type equations

[J]. Ocean Engineering, 2016, 121:166-173.

DOI:10.1016/j.oceaneng.2016.05.033      URL     [本文引用: 2]

ZAKHAROV V E.

Stability of periodic waves of finite amplitude on the surface of a deep fluid

[J]. Journal of Applied Mechanics and Technical Physics, 1968, 9(2):190-194.

DOI:10.1007/BF00913182      URL     [本文引用: 1]

/