我国陆上易采油气资源的日渐枯竭，海上油气的开采工作迫在眉睫，对于其中诸多深水气田，浮式液化天然气平台(floating liqueied natural gas facility，FLNG)具有良好的应用前景。为保证FLNG大型液舱的安全，为其液舱设计有效的减晃方案，对避免液体晃荡破坏平台结构和保证作业安全有重要意义。目前常见的液舱减晃技术主要针对常规液货舱设计，大致分为主动式和被动式两类。

主动式减晃技术依赖外部电力控制，需设置传感器实时反馈液货晃荡状态，进而采用气泡注入^[1]或动态隔板^[2]等技术，实时调整减晃机构的运动或阻尼特性等参数，实现减弱晃荡的效果。被动式减晃技术则无需依赖电力控制，主要利用舱内固定结构增加流体运动时的能量耗散，来达到减晃目的。Wang等^[3]采用半解析方法分析了圆柱液舱内环形隔板的减晃效果；Xue等^[4]对方形液舱中不同隔板形式的减晃特性开展了试验研究；Jung等^[5]研究了舱底垂直隔板的减晃效果；Wei等^[6]考虑了舱底中线上垂直栅形减晃结构；Yu等^[7]对自由面上双垂直隔板的减晃方案进行了试验研究。另有部分研究者提出了结合舱内运动结构的“半被动式”减晃技术。如采用浸没毯式结构^[8]、浸没的膜状结构^[9]、毯式覆盖结构^[10]、在自由面上方填充泡沫^[11]。值得注意的是，Koh等^[12]提出了一种基于浮式结构的减晃概念。Hosseini等^[13]采用悬浮式环形板结构减弱浮顶运动。浮式结构可沿舱内滑轨或铰接于舱壁，随着液面晃荡自由浮动。舱内液体共振主要发生在液体的自然晃荡频率与外部激励频率接近的条件下，而该方案中舱内液体晃荡的自然频率随浮式结构在液体中的位置变化而动态调整，理论上液体晃荡的自然频率始终偏离外部激励频率，这可以很大程度地避免液体共振条件的满足，进而减少液舱共振的发生几率。但对于液舱内浮式结构的运动对液体自然晃荡频率的改变范围，目前尚缺乏量化研究。因此，本文将三维方形液舱为例，研究浮动柱形浮子式减晃结构位于液舱内不同位置时对舱内液体晃荡的自然频率和模态的影响。

1 数学模型和数值方法

针对上述浮子式减晃荡结构，本文给出该结构数值模型的建立方法。建立如图 1所示的直角坐标系O-xyz，其中O位于舱内平均自由面中心，且Ox、Oy和Oz分别与方形液舱的长、宽、高方向平行。

基于线性势流理论假设，速度势ϕ(x, y, z, t)在平均流域V内满足拉普拉斯方程:

$ {\nabla ^2}\phi = 0 $

(1)

对于固定液舱，速度势在舱内平均湿表面${{\bar{S}}_{\text{B}}}$和平均自由面${{\bar{S}}_{\text{F}}}$上分别满足边界条件：

$ \begin{array}{*{20}{c}} {\partial \phi /\partial n,}&{在\;{{\bar S}_{\rm{B}}}\;上} \end{array} $

(2)

$ \begin{array}{*{20}{c}} {\partial \phi /\partial z = \partial \eta /\partial t,}&{在\;{{\bar S}_{\rm{F}}}\;上} \end{array} $

(3)

$ \begin{array}{*{20}{c}} {\partial \phi /\partial t = - g\eta ,}&{在\;{{\bar S}_{\rm{F}}}\;上} \end{array} $

(4)

式中:∂ϕ/∂n表示速度势在流域边界一点处的法向偏导数；η(x, y, t)为自由表面升高。

对于该线性问题，速度势ϕ可构造为各自然晃荡模态的线性叠加:

$ \phi \left( {x,y,z,t} \right) = \sum\limits_{i = 1}^\infty {{\xi _i}\left( t \right){{\bar \phi }_i}\left( {x,y,z} \right)} $

(5)

式中：ξ_i(t)为时间相关系数；${{{\bar{\phi }}}_{i}}$为液体晃荡的自然模态。令ξ_i(t)=cos(ω_it)，可得的关于${{{\bar{\phi }}}_{i}}$的稳态问题，将式(5)代入式(1)~(4)可得到特征值问题：

$ \begin{array}{*{20}{c}} {{\nabla ^2}\overline {{\phi _i}} = 0,}&{在\;\bar V\;上} \end{array} $

(6)

$ \begin{array}{*{20}{c}} {\partial {{\bar \phi }_i}/\partial n = 0,}&{在\;{{\bar S}_B}\;上} \end{array} $

(7)

$ \begin{array}{*{20}{c}} {\frac{{\partial {{\bar \phi }_i}}}{{\partial z}} = \frac{{\omega _i^2}}{g}{{\bar \phi }_i},}&{在\;{{\bar S}_F}\;上} \end{array} $

(8)

式中特征值ω_i即为液体晃荡的自然频率，对应的特征向量${{{\bar{\phi }}}_{i}}$即为晃荡的自然模态。

为了简化表述式，下文将省略下标i。利用格林第三公式将关于ϕ的拉普拉斯方程转换为边界积分方程：

$ \begin{gathered} c\left( {{x_{\text{p}}}} \right)\bar \phi \left( {{x_{\text{p}}}} \right) = \iint\limits_S {\left[ {G\left( {{x_{\text{p}}},x} \right)\frac{{\partial \bar \phi \left( x \right)}}{{\partial n}} - } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\bar \phi \left( x \right)\frac{{\partial G\left( {{x_{\text{p}}},x} \right)}}{{\partial n}}} \right]{\text{d}}S \hfill \\ \end{gathered} $

(9)

式中：S为流体边界；下标p代表边界上的任意配置点；c(x_p)为配置点处的固角系数；G(x_p, x)为格林函数：

$ \begin{array}{l} G\left( {{x_{{p}}},x} \right) = \left[ {{{\left( {x - {x_p}} \right)}^2} + {{\left( {y - {y_p}} \right)}^2} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left. {{{\left( {z - {z_p}} \right)}^2}} \right]^{ - 1/2}} \end{array} $

(10)

采用三角或者四边形单元离散流域边界，并假定各单元上速度势和速度势法向导数为常数，从而将边界积分方程离散为数值积分的形式：

$ 2{\rm{ \mathsf{ π} }}{{\bar \phi }_i} + \sum\limits_{e = 1}^{{N_e}} {{{\bar \phi }_e}{H_{ie}}} = \sum\limits_{e = 1}^{{N_e}} {{{\left( {\frac{{\partial \bar {\mathit{\boldsymbol{\phi}}} }}{{\partial n}}} \right)}_e}{G_{ie}}} ,i = 1,2, \cdots ,{N_e} $

(11)

$ \left\{ \begin{gathered} {H_{ie}} = \iint\limits_{Se} {\left[ {\frac{{\partial G\left( {{x_i},x} \right)}}{{\partial n}}} \right]{\text{d}}S} \hfill \\ {G_{ie}} = \iint\limits_{Se} {\left[ {G\left( {{x_i},x} \right)} \right]{\text{d}}S\left[ {{G_{ba}}} \right]{{\bar \phi }_a}} \hfill \\ \end{gathered} \right. $

(12)

式中：N_e为单元的总数；S_e为第e个单元的面积；x_i为i个单元的几何中心。在单元序列中，假定前N_a个单元为自由液面单元，后N_b=N_e-N_a个单元为液舱湿表面单元，可将方程(11)可写为：

$ \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_{aa}}}&{{\mathit{\boldsymbol{H}}_{ab}}}\\ {{\mathit{\boldsymbol{H}}_{ba}}}&{{\mathit{\boldsymbol{H}}_{bb}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\left\{ {\bar \phi } \right\}}_a}}\\ {{{\left\{ {\bar \phi } \right\}}_b}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{G}}_{aa}}}&{{\mathit{\boldsymbol{G}}_{ab}}}\\ {{\mathit{\boldsymbol{G}}_{ba}}}&{{\mathit{\boldsymbol{G}}_{bb}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\left\{ {{{\bar \phi }_n}} \right\}}_a}}\\ {{{\left\{ {{{\bar \phi }_n}} \right\}}_b}} \end{array}} \right) $

(13)

式中{ }_a和{ }_b分别表示自由液面和液舱内壁单元上的变量。将式(7)和(8)代入方程(13)可以得到：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{H}}_{aa}}{\left\{ {\bar \phi } \right\}_a} + {\mathit{\boldsymbol{H}}_{ab}}{\left\{ {\bar \phi } \right\}_b} = \frac{{{\omega ^2}}}{g}{\mathit{\boldsymbol{G}}_{aa}}{\left\{ {\bar \phi } \right\}_a}\\ {\mathit{\boldsymbol{H}}_{ba}}{\left\{ {\bar \phi } \right\}_a} + {\mathit{\boldsymbol{H}}_{bb}}{\left\{ {\bar \phi } \right\}_b} = \frac{{{\omega ^2}}}{g}{\mathit{\boldsymbol{G}}_{ba}}{\left\{ {\bar \phi } \right\}_a} \end{array} \right. $

(14)

进一步消除{ϕ}_b项可以得到一个标准的特征方程：

$ \begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{G}}_{aa}} - {\mathit{\boldsymbol{H}}_{ab}}{\mathit{\boldsymbol{H}}_{bb}}^{ - 1}{\mathit{\boldsymbol{G}}_{ba}}} \right)}^{ - 1}} \cdot }\\ {\left( {{\mathit{\boldsymbol{H}}_{aa}} - {\mathit{\boldsymbol{H}}_{ab}}{\mathit{\boldsymbol{H}}_{bb}}^{ - 1}{\mathit{\boldsymbol{H}}_{ba}}} \right){{\left\{ {\bar \phi } \right\}}_a} = \frac{{{\omega ^2}}}{g}{{\left\{ {\bar \phi } \right\}}_a}} \end{array} $

(15)

该方程的特征值为ω_i²/g(i=1, 2, …, N_a)，特征向量为ϕ_i，即得到液舱晃荡的自然频率ω_i和自然模态ϕ_i。应该指出的是，液体晃荡有无穷多个自然频率和自然模态，而本数值求解可至多得到前N_a个自然频率，即0 < ω₁≤ω₂≤…≤ω_Na(等号表示不同自然模态对应的自然频率可能具有相同的值)。

应该指出的是，液舱晃荡属于强非线性问题，本文所采用的线性理论通常仅适用于准确预报舱内液体的自然频率和模态^[14]。

2 收敛性数值验证

在本节中，液舱的长、宽分别用L、B表示，舱内液体的平均深度设为H。本文以方形液舱(L=B=1, H=0.2)为例，对数值收敛性进行验证。以L和g为基底对所有变量进行无量纲化处理，无量纲后的长度和频率分别表示为(x, y, z, B, H)→(x, y, z, B, H)L和ω→ω$\sqrt{g/L}$。考虑3种边界网格：Mesh-1:10×10×4，Mesh-2:20×20×8和Mesh-3:30×30×12，且网格单元从Mesh-1到Mesh-3逐渐加密。对于Mesh-2, 数组20、20和8分别表示沿方形流域长、宽、深方向的单元数，该网格的单元分布如图 2所示。对于复杂形状的自由面，基于三角形单元的网格相对规则，有利于提高计算的准确性，因此本文采用三角形单元离散自由面边界。基于各个网格，可通过式(15)计算得到从小到大排列的前N_a个自然晃荡频率。图 3对比了基于3网格模型得到的前20个自然频率，并与方形液舱自然晃荡频率的解析解进行比较。三维方形液舱自然晃荡频率的无量纲解析解为

$ {\omega _{mn}} = \sqrt {\frac{{\rm{ \mathsf{ π} }}}{L}\sqrt {{m^2} + {n^2}} \tanh \left( {\frac{{H{\rm{ \mathsf{ π} }}}}{L}\sqrt {{m^2} + {n^2}} } \right)} $

(16)

式中：m, n=0, 1, 2, …；mn≠0。由图 3可见，随着网格单元密度的增加，数值计算得到的自然频率值逐渐收敛于解析解，本文所采用数值方法的收敛性得到验证。同时，由于Mesh-2和Mesh-3两种网格产生的数值结果几乎重合，为保证计算精度和减少计算时间，后文将选取Mesh-2开展数值计算。

3 数值结果与分析

本节重点研究舱内柱形浮子式减晃结构的位置、吃水和尺寸对液体自然晃荡频率和模态的影响。考虑3个典型的浮子半径r为0.05、0.1和0.2，以及3个典型浮子吃水h为0.05、0.1和0.15。浮子在液舱中的测试位置如图 4中黑点所示，图 4中黑点的覆盖区域实际上就是浮子的水平运动范围。为保证柱形浮子与舱壁之间有足够大的空隙，对于浮子直径较大的情况(r=0.2)，浮子的运动范围从图 4(a)调整为图 4(b)。液舱中液体深度分别设H为0.2、0.3和0.4。

以r=0.1，h=0.15和H=0.2这一情况为例。当浮子处于图 4(a)中A、B、C3个不同位置时，舱内液体晃荡的前3个基本模态分别如图 5、6和7所示。实际海况中，通常最小的自然频率w₁最接近真实的海浪频率w，即最容易引起共振^[14]，故后文仅关注w₁。为方便表达，使用“自然晃荡频率w”来代指“最小自然晃荡频率w₁”。

3.1 浮子位置对液舱自然晃荡频率的影响

本节考虑水深H为0.2的液舱。图 8给出了半径r为0.05、吃水分别为h为0.05、0.1和0.15的浮子处于不同位置时，舱内液体的自然晃荡频率分布。图中各坐标点处的颜色表示浮子位于该位置时，舱内液体晃荡的自然频率值。由图 8可见:1)浮子位于液舱角落位置时，舱内液体的自然晃荡频率达到最大值；2)当浮子靠近且沿着液舱边界运动时，舱内液体自然晃荡频率的变化范围最大；3)当位于液舱中心的浮子平行于液舱边界运动时，由浮子位置变化所引起的自然晃荡频率变化的幅值最小。图 9和10考虑了浮子半径r增大至0.1和0.2的情况，从图中可以得到与图 8相似的规律。另外，可以比图 8更明显地看出，当浮子位于靠近液舱边界中点位置时，舱内液体的自然晃荡频率达到最小值。

3.2 浮子吃水和半径对液舱自然晃荡频率变化范围的影响

本节将研究不同吃水和半径的浮子靠近且沿液舱边界移动时，舱内液体自然晃荡频率的变化范围。考虑水深H为0.2的液舱。

图 11(a)给出了半径r为0.05且吃水分别为h=0.05、0.1和0.15的浮子沿液舱壁移动时，舱内自然晃荡频率的变化情况。图 11中，横坐标λ表示柱形浮子中轴线到相邻舱壁中点的水平距离(如图 4所示)，纵坐标表示有浮子液舱的自然晃荡频率与无浮子液舱的自然晃荡频率的偏差(即(w-w′)/w′×100%，其中w₁为式(15)给出的有浮子液舱的自然晃荡频率，w′为式(16)给出的无浮子液舱的自然晃荡频率)。由图中可以看出，浮子从液舱角落舱壁中点移动时，液舱的自然晃荡频率逐渐减小。然而对于半径r为0.05和吃水h为0.05的浮子，由于其排水体积极小，导致其水平运动引起的液舱自然晃荡频率的变化范围很小(小于1%)。将吃水h增加至0.1和0.15后，液舱自然晃荡频率的变化范围增大，不过依然小于2%。图 11(b)和(c)进一步将浮子半径r增加为0.1和0.2，液舱的自然晃荡频率受浮子运动的影响明显增大。对于吃水h为0.15的浮子，当其半径r从0.05增加至0.1和0.2时，液舱自然晃荡频率的变化范围从2%增加至4%和13%。同时，比起小直径r为0.05浮子，大直径浮子r=0.2的吃水变化对液舱自然晃荡频率的影响更为明显。然而，对于某一位置固定半径的浮子，在其吃水h从0.05增加至0.1和吃水h从0.1增加至0.15两种情况时，液舱自然晃荡频率的变化幅值相当，这显示液舱自然晃荡频率变化与吃水变化近似成线性关系。

图 12为r=0.2保持不变但不同半径的浮子沿舱壁移动时，舱内自然晃荡频率的变化情况。由图 12可见，对于同样吃水的浮子，增大其浮子半径后，液舱自然晃荡频率随空间位置变化的曲线变陡。而对于同样吃水的浮子，其半径r从0.05增加至0.1后液舱自然晃荡频率变化幅值的2倍，半径从0.1增加至0.2后液舱自然晃荡频率变化幅值，显示液舱自然晃荡频率变化与浮子半径变化呈现非线性关系。相较于小吃水h为0.05的情况，大吃水h为0.15浮子的半径变化对液舱自然晃荡频率的改变更加明显。

在实际设计中，浮子的位置和吃水均在液舱中动态改变，液舱自然晃荡频率也随之动态调整。图 13给出了浮子沿舱壁移动且吃水在0.05~0.15变化时，液舱自然晃荡频率的变化范围(由图中阴影区域表示)。当浮子在该空间范围运动时，液舱自然晃荡频率一定出现在图中阴影区域内。换言之，图中阴影面积的大小意味着运动浮子对液舱自然晃荡频率改变范围的大小。图中3块阴影区域分别对应着3个不同的浮子半径，可见采用大直径运动浮子，能在更大频率范围内改变液舱的自然晃荡频率。应该指出的是，图中每块阴影区域的上下边界分别对应着h=0.05和h=0.15两种等吃水情况，即阴影上下边界处浮子的排水体积保持不变。图中将各阴影区域上下边界对应的排水体积标记出来。值得注意的是，尽管r=0.1阴影区域下边界对应的浮子排水体积(4.71×10^-3)小于r=0.2阴影区域上边界对应的浮子排水体积(6.28×10^-3)，但是液舱的自然晃荡频率在r=0.1阴影区域下边界条件下的改变范围更大。这意味着排水体积更大的浮子，并不一定会使液舱的自然晃荡频率发生更大改变。另外，与无浮子液舱的自然晃荡频率相比，采用浮子结构后液舱的自然晃荡频率普遍降低。

3.3 水深对液舱自然晃荡频率的影响

本节将进一步分析水深对液舱自然晃荡频率的影响。图 14(a)显示了半径r=0.2的浮子在不同水深条件下，液舱自然晃荡频率随浮子吃水的变化。图 14(b)显示了吃水为h=0.15的浮子在不同水深条件下，液舱自然晃荡频率随浮子半径的变化。由图可见，水深越小，液舱自然晃荡频率随浮子吃水或半径的变化越大。图 14(c)给出了不同水深条件下半径为r=0.2、吃水为h=0.15的浮子沿舱壁运动时，液舱自然晃荡频率的变化。从图中可以明显看出，浮子运动对浅水液舱的自然晃荡频率影响更大，即圆柱浮子式结构对浅水液舱的减晃效果更好。

4 结论

1) 与无浮子液舱的自然晃荡频率相比，采用浮子结构后液舱的自然晃荡频率普遍降低。

2) 浮子位于液舱角落位置时，舱内液体的自然晃荡频率达到最大值；当浮子位于靠近液舱边界中点位置时，舱内液体的自然晃荡频率达到最小值；当浮子靠近且沿着液舱边界运动时，舱内液体自然晃荡频率的变化范围最大。

3) 在其他条件不变的情况下，大半径、大吃水浮子运动时，对液舱自然晃荡频率的影响更为明显。

4) 浮子运动对浅水液舱的自然晃荡频率影响更大，即圆柱浮子式结构对浅水液舱的减晃效果更好。