部分充满液体的容器在外部激励下会引起自由液面的运动，这种现象就是晃荡。晃荡现象通常可以被理解为液面的自由波动，可分为驻波、行波、水跃和破波，甚至是几种运动的组合^[1]。在外部激励频率接近液体的自振频率时，会产生共振现象，此时流体的运动极其复杂，呈现出非常强烈的非线性和随机性，与此同时，剧烈的液体运动会产生极大的局部抨击力，结构物将在瞬时承载巨大的冲击载荷，甚至造成结构的损坏以及船体稳性的丧失^[1]。为了保障工程设计的安全性和船舶航行的稳定性，对于晃荡问题的研究就变得十分必要。

最早进行相关研究主要是在航空航天领域，Abramson等^[2]研究了该问题并对前人的工作进行了梳理，主要采用势流理论分别对几种不同形状的液舱进行了研究，针对线性问题，不仅推导出了经典的解析解，还做了大量的实验予以验证；Faltinsen等^[3-4]同样采用了不可压缩势流理论对该问题进行了研究，并得到了非线性的速度势公式，他的研究成果被许多人广泛使用，这其中就包括Nakayama等^[5]，以及Waterhouse等^[6]，但是其理论依然受到很大的限制，比如容器形状简单，内部光滑无其他结构物，且水深较高等；当水跃现象出现时，浅水波理论是有用的数学模型，这其中包括Verhagen等^[7]以及Armenio等^[8]对该理论和雷诺平均NS方程的比较；随着计算机的广泛运用，越来越多的人试图寻求计算机来求解流体力学问题，CFD的出现让晃荡问题的研究有了新的出路^[9]，事实上很多种数值方法都可以用来研究晃荡现象，这其中VOF方法因为其处理自由界面的优越性而得到了广泛的应用^[10-11]。

本文采用一种将树形自适应网格与VOF方法相结合的并行算法^[12]对几种情况下的晃荡问题进行了数值模拟，并将计算得到的几种工况下的自由表面升高，力和力矩同解析解，实验值进行了对比，验证了该数值方法对于此问题的有效性，并对计算结果进行了相应的分析。

1 液体晃荡数学模型 1.1 控制方程

采用不可压变密度Navier-Stokes方程组作为控制方程:

$ \rho \left( {\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \mathit{\boldsymbol{u}} \cdot \nabla \mathit{\boldsymbol{u}}} \right) = \rho \mathit{\boldsymbol{F}}-\nabla p + \nabla \cdot \left( {2\mu \mathit{\boldsymbol{D}}} \right) $

(1)

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{u}}} \right) = 0 $

(2)

式中：u为速度向量；p为压力；D表示应变率对称张量，其各分量为D_ij=(∂_iu_j+∂_ju_i)/2；ρF代表体积力。晃荡问题坐标系如图 1所示，式(1)、(2)中各变量均在非惯性坐标系xoy上，故在该坐标系上的F各分量表示为^[13-14]:

$ {F_i} =-g\sin \theta + {x_j}\ddot \theta + {x_i}{\dot \theta ^2} + 2\dot \theta {u_j}-\ddot A $

(3)

$ {F_j} =-g\cos \theta-{x_i}\ddot \theta + {x_j}{\dot \theta ^2}-2\dot \theta {u_i} $

(4)

式中：g为重力加速度；x_i为x轴坐标；x_j为y轴坐标；u_i为x方向速度标量；u_j为y方向速度标量；θ和A(x′o′y′坐标系下)分别表示水槽外激励随时间变化的纵摇幅值和纵荡幅值(x′轴方向)，且θ=0时为图中实线位置，它们的表达式为:

$ A = {A_0}\sin \left( {{\omega _s}t + {\varepsilon _s}} \right) $

(5)

$ \theta = {\theta _0}\sin \left( {{\omega _p}t + {\varepsilon _p}} \right) $

(6)

式中A₀、θ₀、ω_s、ω_p、ε_s、ε_p分别表示纵荡和纵摇的幅值、频率与相位差。

考虑到两相流问题VOF方法中体积分数T的引入，于是，密度ρ和动力粘性系数μ可表示为：

$ \rho \left( T \right) = {\rho _1}T + {\rho _2}\left( {1-T} \right) $

(7)

$ \mu \left( T \right) = {\mu _1}T + {\mu _2}\left( {1-T} \right) $

(8)

式中下标1、2分别代表 2种不同的流体。

1.2 数值方法 1.2.1 空间离散

整个计算域的离散采用的是若干正方形有限体积单元，它们之间具有如同数据结构中四叉树(Quadtree)(八叉树Octree)的继承关系，图 2给出了一个二维情况下计算域离散后的单元对应关系，定义每一个有限体积单元(cell)的边(Edge)的长度为h，每个单元都有可能是其4个子单元(Children)(三维情况下是8个)的父单元(Parent)，需要注意的是根单元(Root cell)是结构树的根节点而叶节点处的叶单元(Leaf cell)则没有任何子单元，单元的水平(Level)是从根单元由零级开始每遇一级子单元则依次加1来确定的，每个单元C在同一水平时每个方向d(二维情况下为4个)上都有一个直接邻单元(Neighbour) N_d，它们之间通过单元的交界面(Face) C_d相连接，混合单元(Mixed cell)被用来处理嵌有固壁的计算域，控制方程中最初的变量都定义在各单元的中心位置。

为了方便在单元边界处的计算，以下约束需要注意，如图 2所示，首先，每个单元的直接邻单元的水平不能和其本身的水平相差超过1，其次，每个单元的对角线邻单元的水平不能和其本身的水平相差超过1。最后，所有和混合单元直接相邻的单元必须与混合单元具有同一水平。

目前所使用的满足数据结构算法需求的，是由Khokhlov等^[16]提出的全线程树方法，它保证了对于给定的任意单元，得到它的临单元和其自身水平及空间坐标的时间复杂度为O(1)。对于遍历所有叶单元，所有给定水平的单元集以及所有混合单元使用的是标准的基于指针的树方法，其时间复杂度为O(NlogN)。

由于树型数据结构在空间离散中的应用，对于整个自适应的算法而言就变得相对易于操作了，它实际上分为2个部分，首先，所有满足给定自适应条件的叶单元需要被细分，其次，所有满足给定自适应条件的叶单元的父单元需要被重组为叶单元，即第一步的逆。

很多变量都可以作为自适应条件的依据，以涡量为例：

$ \frac{{h\left\| {\nabla \times \mathit{\boldsymbol{U}}} \right\|}}{{\max \left\| \mathit{\boldsymbol{U}} \right\|}} > \tau $

(9)

式中：网格的大小用长度h来表示；自适应的阈值则用τ来表示(0 < τ < 1)。这样的话，一个基于局部涡量范数的自适应条件就定义了。

自适应算法较普通结构化均匀网格可以有效减少计算时所需的时间和存储空间，并且可以产生更加精确的计算结果，这是因为自适应算法可以在不需要高分辨率的网格计算区域做到网格自动加粗，而对于需要细化网格的地方进行加密，比如速度梯度较大，或者自由表面变形较大时，都需要较高的计算精度，从而保证计算的高效与准确。基于四叉树的自适应算法又继承了方形网格的优点，这样的有限体积单元形状方便对各点的变量进行插值计算和自由界面的重构，且计算精度较高。

1.2.2 时域离散

对时域的离散，本文主要使用了分数步映射法^[12]，离散后的控制方程为：

$ \begin{gathered} {\rho _{n + \frac{1}{2}}}\left[{\frac{{{u_*}-{u_n}}}{{\Delta t}} + {u_{n + \frac{1}{2}}} \cdot \nabla {u_{n + \frac{1}{2}}}} \right] = \hfill \\ \nabla \cdot \left[{{\mu _{n + \frac{1}{2}}}\left( {{D_n} + {D_{n + \frac{1}{2}}}} \right)} \right] + \rho \mathit{\boldsymbol{F}} \hfill \\ \end{gathered} $

(10)

$ {u_{n + 1}} = {u_*}-\frac{{\Delta t}}{{{\rho _{n + \frac{1}{2}}}}}\nabla {p_{n + \frac{1}{2}}} $

(11)

$ \nabla \cdot {u_{n + 1}} = 0 $

(12)

因为VOF方法中体积分数的引入，密度的对流方程(2)可等效地表示为体积分数T的对流方程，离散后的形式为：

$ {T_{n + \frac{1}{2}}}-{T_{n-\frac{1}{2}}} + \nabla \cdot \left( {{T_n}{u_n}} \right) = 0 $

(13)

式中u_*表示的是一个临时速度场，它实际体现的是以下的算法实质。

为和方程(11) ~(13)中变量区分，使用其他符号，读者可以近似地和控制方程的形式对比：

$ \frac{{{U_*}-{U_n}}}{{\Delta t}} =-{A_{n + \frac{1}{2}}} $

(14)

$ {U_*} = {U_{n + 1}} + \nabla \phi $

(15)

通过上面方程(14)和方程(15)将原来的速度场进行了分解得到了一个临时速度场U_*，实质就是对U_*运用了一个基于Hodge分解的映射法后得到的一个新的速度场U_n+1和分数步压力场p_n+1/2，注意到：

$ \nabla \cdot \mathit{\boldsymbol{U}} = 0\;\;\;\;\;{\text{in}}\;\varOmega $

(16)

$ \mathit{\boldsymbol{U}} \cdot n = 0\;\;\;\;\;\;\;\;{\text{on}}\;\;\partial \varOmega $

(17)

于是对方程(16)两边取散，就得到了压力Poisson方程：

$ {\nabla ^2}\phi = \nabla \cdot {U_*} $

(18)

满足边界条件：

$ \frac{{\partial \phi }}{{\partial n}} = {U_*} \cdot n\;\;\;\;\;\;\;\;\;在\partial \varOmega 上 $

(19)

通过求解压力Poisson方程得到无散速度场:

$ \mathit{\boldsymbol{U}} = {U_*}-\nabla \phi $

(20)

需要指出的是方程(19)的求解对于不同水平的有限体积单元，特殊的插值格式是必要的，文献[12]中有详细的讨论。

1.2.3 自由界面处理

自由界面的处理主要基于一种VOF和不可压缩四叉树自适应网格相结合的方法^[17]，其中经典的VOF几何格式被运用在树形自适应网格上，这是为了保证在自由界面上不同分辨率的空间网格可以被处理。

VOF是一种运用非常广泛且高效的处理自由表面的数值方法，最早由Debar等^[18]提出，然后由Hirt等^[19]进一步完成，VOF方法稳定且实用，其基本思想为：定义函数T和1-T，它们分别表示在计算域中水和空气的体积分数，对于每一个独立的网格单元，水和空气的体积分数之和必为1，于是存在以下3种情况：

1) T=1，网格单元内全部是水；

2) T=0，网格单元内全部是空气；

3) 0 < T < 1，网格单元内一部分是水，一部分为空气，于是水和空气之间就存在一个交界面，即自由界面或者自由表面。

这里使用的是拓展到树形空间网格的分段线性的几何VOF格式，传统的VOF主要分为2个步骤：1)界面重构；2)几何通量的计算及界面流动的处理。其主要思想^[20]是运用体积分数的梯度来确定自由表面的法向，通过方程：

$ \mathit{\boldsymbol{n}} \cdot \mathit{\boldsymbol{x}} = \alpha $

(21)

来计算出自由表面的位置的近似值。式中：n为局部自由界面的法向；x为位置向量；α和T则满足一定的函数关系：T=V(n, α)。

VOF能够用来方便地描述自由界面的变化，所以在处理有自由界面的计算流体动力学问题时，VOF是一个理想的工具。

2 计算结果分析 2.1 算例1

图 3给出了在方形容器内，宽为1 m，水深为0.7 m，纵摇幅值θ₀=4°，频率ω₀=0.5π rad/s时参考点为原点，力矩的计算结果，其中虚线为数值计算得到的数据，实线则为按照文献[2]中所给的解析解表达式得到的计算结果，从图中可以看出，计算值与解析解周期吻合，且波形也较为一致，考虑到势流理论得到的解析解的线性及刚性假设，幅值本身偏小，验证了该数值模拟对于晃荡问题力矩的计算是可靠的。

2.2 算例2

如图 4所示的容器，对其施加一个相对于水槽中心的外部纵摇简谐激励，其中θ₀分别为4°和8°，ω_p=2.0 rad/s，图 5给出了压力的计算值和Akayildiz等^[21]得到的实验值，通过对比发现，数值模拟的结果和实验测得的结果周期一致，波形大体一致，且激励幅值较小时计算值与实验值越接近，从图中数据可以看出，该数值模拟对于压力的计算也是较为准确的。

2.3 算例3

如图 6所示一矩形水槽，距离左侧0.05 m的地方有探测器来输出液面高度，外部激励则为纵荡，通过数值计算得到图 7所示。根据Abramson^[1]给出的晃荡问题一阶自然频率计算公式：

$ \omega _0^2 = \frac{{g\pi }}{l}\tanh \left( {\pi \frac{d}{l}} \right) $

(23)

通过图 7数据，有ω₀=3.767 4 rad/s。图 7 (a)和(b)为A₀=0.032 m，ω_s=1.11ω₀时液面升高数值计算结果与Faltinsen^[22]给出的解析解和实验值的比对，图 7(c)和(d)则为A₀=0.27 m，ω_s=1.28ω₀时数值结果与Faltinsen^[22]给出的解析解和实验值的对比。图 8(a)为A₀=0.032 m，ω_s=1.11ω₀时水平方向合力的数值计算结果，图 8(b)为A₀=0.27 m，ω_s=1.28ω₀时该工况下水平方向合力的数值计算结果。

通过图 8的对比，数值结果明显较解析解更符合实验所得到的结果，首先解析解的周期要和实验值有一定的误差，图 8(c)、(d)在频率偏高，振幅较小时，Faltinsen^[22]给出的理论计算方法的周期较实验值要偏小，而图 8(a)、(b)频率稍小，振幅稍微增大时，理论值的周期则又大于实验值，与此同时，无论何种情况，数值计算总能得到较为满意的结果，周期和实验完全吻合，波形上，数值方法的结果也和实验结果基本一致，外包络线基本相符，小波的峰值和周期也完全吻合，这充分说明了基于四叉树自适应网格的VOF方法具有处理二维强非线性晃荡问题的优势。

通过图 7和图 8的对比，发现水平激励下容器受到水平方向的合力与液面升高具有相同的周期，并且波形的时间历程相似，都属于拍振，是外在激励频率与自然频率相结合的结果，外包络线周期满足T=2π/|ω₀-ω_s|。

为考察自适应网格效果，图 9给出了此例计算网格模型，算例中采用了对速度、VOF函数梯度的自适应，网格此外对适应于该问题的边界层厚度进行了加密，通过图 9可以清楚的看到网格自适应的分布和变化情况。

从网格模型可见，四叉树方法能够在自由液面附近自动进行加密，晃荡液面得到了精确捕捉，提高了计算精度，流场的其他区域采用大网格计算，有效提高了计算效率。

3 结论

1) 基于四叉树自适应网格的VOF方法在处理晃荡问题时，通过对部分需要细化的网格加密自适应可以更加高效精确地得到相应的结果，从而有效地计算液面升高及力与力矩；

2) 纵摇时，对于压力的计算，数值计算结果同实验得到的数据基本一致，与此同时，激励幅值稍小时，得到的结果越符合实验结果；

3) 纵荡问题中，容器受到外部激励，当自然频率与激励频率相差不大时，水平方向的合力与液面升高具有相同的周期，且时间历程曲线都出现拍振，波形相似。