目前，船舶气层减阻技术已在俄罗斯、乌克兰、日本、荷兰等国家获得了广泛应用，取得了显著的军事经济效益^[1]。我国自20世纪80年代开始，在该领域开展了相关理论及试验研究，取得了静水中模型绝对减阻率达20%~25%的效果^[2-4]，但至今尚未实现实艇应用，主要原因即是在考虑风浪作用下气泡高速艇的减阻效果及运动性能尚未明确。气泡高速艇在波浪中是否仍然具备良好的减阻效果？波浪中的艇底气层是否能保持稳定状态？气层作用下高速艇在波浪中的运动性能如何？气层与摇荡运动艇体之间的相互作用力学规律如何？这一系列问题尚需深入研究和解答。

针对以上问题，国内外公开发表的相关理论及试验研究不多。在试验研究方面，董文才等^[5-6]通过拖曳水池试验探索了波浪中的气层减阻效果以及气层对艇体运动的影响规律，但鉴于试验条件限制，未能给出艇体运动、气层形态、阻力的同步变化过程。在理论研究方面，美国的Jin-Keun Choi等^[7]联合采用面元法与求解URANS的方法研究了喷气下瘦长型船的船体兴波、气层形态及船体运动对气层分布的影响，但该方法在计算船体兴波时采用线性势流理论，不适用于计算高速艇的喷溅等非线性流动问题。

为此，本文考虑了高速艇流场的强非线性及气-液两相流的复杂性，采用粘性CFD的方法，从波浪中阻力、纵向运动及气层形态等角度进行研究，试图通过理论计算揭示波浪中气层-艇体的相互作用过程，掌握气泡高速艇在波浪中的航行性能。

1 数值波浪水池 1.1 研究对象

计算对象为一条底部设置凹槽的B.H.型气泡高速艇模型，主尺度参数如表 1所示，图 1给出了艇体及底部凹槽的三维效果图。其中，Δ为排水量，L为艇体总长，L_P为折角线长，B_PX为折角线最大宽度，β_M为中部艇底斜升角，β_T为尾部艇底斜升角；L_P/B_PX表示折角线长与折角线最大宽度之比，H_B/B_PX表示艇体凹槽深度与折角线最大宽度之比，S_BH/S_P表示水平面上的艇底凹槽投影面积与折角线投影面积之比，R_yy/L表示纵向惯性半径与艇长之比。

1.2 两相流模型

气-液两相流采用VOF模型(volume of fluid model)，界面追踪采用高分辨率HRIC方法(high-resolution interface capturing scheme)。考虑界面压缩修正的VOF方法的控制方程如式(1)~(5)所示，详细可参见文献[8]。

质量守恒方程：

$ \frac{{\partial {\alpha _q}}}{{\partial t}} + {\mathit{\boldsymbol{u}}_q} \cdot \nabla {\alpha _q} + \nabla \cdot \left( {{\mathit{\boldsymbol{u}}_c}{\alpha _q}\left( {1 - {\alpha _q}} \right)} \right) = 0 $

(1)

动量方程：

$ \begin{array}{*{20}{c}} {\frac{{\partial \left( {{\rho _q}{\alpha _q}{\mathit{\boldsymbol{u}}_q}} \right)}}{{\partial t}} + \left( {{\rho _q}{\alpha _q}{\mathit{\boldsymbol{u}}_q} \cdot \nabla } \right){u_q} = }\\ { - {\alpha _q}\nabla p + \nabla \cdot \left( {{\mu _q}{\alpha _q}\nabla {\mathit{\boldsymbol{u}}_q}} \right) + }\\ {{\rho _q}{\alpha _q}\mathit{\boldsymbol{g}} + {\mathit{\boldsymbol{F}}_{D,q}} + {\mathit{\boldsymbol{F}}_{s,q}}} \end{array} $

(2)

式中：ρ_q、α_q和u_q分别表示第q相的密度、体积分数和流动速度；F_{D, q}和F_{s, q}分别为相间曳力与表面张力，本文计算过程中忽略了这两项的影响；u_c为界面压缩速度，其表达式为u_c=C_α|u|(∇α/∇α)，式中|u|为流动速度，C_α为锐化因子，取值范围为0~1。

界面追踪采用高分辨率HRIC格式，其基础是NVD方法(normalized variable diagram)，原理示意图如图 2所示。图中，α_D、α_A、α_U分别表示施主单元、受主单元及迎风单元的体积分数，α_f表示控制体边界的体积分数，其计算表达式为：

$ {\alpha _f} = \left( {1 - {{\tilde \beta }_f}} \right){\alpha _D} + {{\tilde \beta }_f}{\alpha _A} $

其中${{\tilde \beta }_f} = ({{\tilde \alpha }_f} - {{\tilde \alpha }_D})/(1 - {{\tilde \alpha }_D})、{{\tilde \alpha }_f} = ({\alpha _f} - {\alpha _U})$/$({\alpha _A} - {\alpha _U})、{{\tilde \alpha }_D} = ({\alpha _D} - {\alpha _U})/({\alpha _A} - {\alpha _U})$。

HRIC方法在求解过程中考虑了界面与网格夹角、当地库朗数等因素的影响，其求解的主要方程和步骤如下^[9-10]：

1) 采用NVD方法计算归一化体积分数${{\tilde \alpha }_f}$。

$ {{\tilde \alpha }_f} = \left\{ \begin{array}{l} {{\tilde \alpha }_D},\;\;\;\;\;\;{{\tilde \alpha }_D} < 0,{{\tilde \alpha }_D} \ge 1\\ 2{{\tilde \alpha }_D},\;\;\;\;0 \le {{\tilde \alpha }_D} < 0.5\\ 1,\;\;\;\;\;\;\;\;\;0.5 \le {{\tilde \alpha }_D} < 1 \end{array} \right. $

(3)

2) 考虑角度θ_f的影响，通过权重因子γ_f对界面伪变形进行修正。

$ \tilde \alpha _f^ * = {\gamma _f}{{\tilde \alpha }_f} + {{\tilde \alpha }_D}\left( {1 - {\gamma _f}} \right) $

(4)

式中：θ_f定义为施主单元D指向受主单元A的单元向量d与界面法向n之间的夹角(如图 3所示)，γ_f=(cos θ)^C_θ，C₀称为角度因子。

3) 考虑当地库朗数的影响，对计算的不稳定性进行修正。

$ \tilde \alpha _f^{ * * } = \left\{ \begin{array}{l} \tilde \alpha _f^ * \;\;\;{{\tilde \alpha }_D} < 0,{{\tilde \alpha }_D} > 1\;\;\;\;\;\;\;{C_f} < {C_{UL}}\\ {{\tilde \alpha }_D}\;\;\;0 < {{\tilde \alpha }_D} < 0.5\;\;\;\;\;\;\;\;\;{C_f} > {C_{UU}}\\ {{\tilde \alpha }_D} + \left( {\tilde \alpha _f^ * - {{\tilde \alpha }_D}} \right)\frac{{{C_{UU}} - {C_f}}}{{{C_{UU}} - {C_{UL}}}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{C_{UL}} \le {C_f} \le {C_{UU}} \end{array} \right. $

(5)

式中：C_UL为库朗数下界值，C_UU为库朗数上界值，C_f为当地库朗数。

1.3 湍流模型

湍流模型采用在船舶流体性能计算中被广泛应用和验证的Realizable k-ε模型^[11-12]，其基本数学表达形式如下。

湍动能方程(k方程)：

$ \begin{array}{*{20}{c}} {\frac{{\partial k}}{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho k{u_i}} \right) = }\\ {\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + {G_k} - \rho \varepsilon - {Y_M}} \end{array} $

(6)

耗散率方程(ε方程)：

$ \begin{array}{*{20}{c}} {\frac{{\partial k}}{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \varepsilon {u_i}} \right) = }\\ {\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + \rho {C_1}S\varepsilon - \rho {C_2}\frac{{{\varepsilon ^2}}}{{k + \sqrt {\upsilon \varepsilon } }}} \end{array} $

(7)

式中：${\mu _t} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon },{G_k} = {\mu _t}{S^2},S = \sqrt {2{S_{ij}}{S_{ij}}} ,{Y_M} = $$2\rho \varepsilon M_t^2,{M_t} = \sqrt {k/{a^2}} ,{C_1} = {\rm{max}}\left[ {0.43,\frac{\eta }{{\eta + 5}}} \right],\eta = Sk/\varepsilon $壁面附近采用两层模型进行处理，单方程模型为切向力驱动的Wolfstein模型，其尺度函数和湍流粘度系数函数的表达式为

$ {l_\varepsilon } = {c_l}y\left[ {1 - \exp \left( { - \frac{{\mathit{R}{\mathit{e}_y}}}{{{A_\varepsilon }}}} \right)} \right] $

(8)

$ \frac{{{\mu _t}}}{\mu } = \mathit{R}{\mathit{e}_y}C_\mu ^{1/4}\kappa \left[ {1 - \exp \left( { - \frac{{\mathit{R}{\mathit{e}_y}}}{{{A_\mu }}}} \right)} \right] $

(9)

式中：$R{e_y} = \rho \sqrt {ky} /\mu $(y表示至壁面的距离)，当Re_y＜Re^*(临界值，本文Re^*=200)时采用单方程模型，Re_y＞Re^*时采用Realizable k-ε两方程模型；A_ε=2c_l；c_l=κC_μ^-3/4；C_μ=0.09；κ=0.42；A_μ=70。

1.4 数值造波

数值波浪采用线性波，无限水深下的波浪参数满足：

$ \left\{ \begin{array}{l} {\omega ^2} = gk\\ U = a\omega \sin \left( {\left( {kx - \omega t} \right){{\rm{e}}^{kz}}} \right) \end{array} \right. $

(10)

式中：ω为圆频率，k为波数，a为波幅，U为边界入口处的流体速度，x和z分别为流体的纵向和垂向坐标位置。此外，为保持压力出口边界条件的稳定性，在距离压力出口1倍艇长区域内添加阻尼项，动量源的函数表达形式为

$ S_z^d = \rho \left( {{f_1} + {f_2}\left| w \right|} \right)\frac{{{{\rm{e}}^k} - 1}}{{{\rm{e}} - 1}}w,k = \left( {\frac{{x - {x_{sd}}}}{{{x_{ed}} - {x_{sd}}}}} \right) $

(11)

式中：S_z^d为动量源函数，x_sd为消波开始坐标，x_ed为消波终止坐标，f₁、f₂和n_d为模型参数，w为流体质点垂向速度。

1.5 计算域及网格

考虑到流动关于艇体中纵剖面对称，仅对流场的一半进行网格离散。计算流域分成随艇体一起运动的重叠网格区域和固定不动的背景区域，如图 4所示。这两个区域均采用剪切型网格进行离散，区域之间的数据传递采用重叠网格技术。

重叠网格区域前方距艇艏0.25L，后方距艇艉0.3L，上表面距甲板为0.3L，下表面距艇底0.25L，侧面距舷侧为0.3L，均设置为重叠网格边界。中纵剖面设置为对称面，艇体表面设置为无滑移物面，喷气口位于艇底凹槽的艏部，设置为速度入口。

背景区域在计算过程中保持静止，该区域总长5.5L、总宽2L、总高3.1L，上表面距离艇体甲板为1.3L，下表面距离艇底部为1.7L，侧面距离艇体舷侧为1.8L，均设置为对称面；入口距离艇艏为1.5L，设置为速度入口；出口距离艇艉为3.0L，设置为压力出口；中纵剖面设置为对称面。

图 5给出了重叠网格区域及背景区域的网格划分情况。对艇底凹槽、艏部喷溅及艉部兴波等区域进行了网格加密，艇体附近的网格尺度为5 mm，喷溅处的网格尺度为3 mm，边界层网格厚度为3 mm，设置4层网格，沿法向的增长率为1.25。对自由液面附近的网格进行了加密，垂直自由液面方向的网格尺度为4 mm，同时为保证艇体姿态变化过程中重叠区的网格尺度基本相同，自由液面加密区的高度在重叠网格区域中为0.06 m，在背景区域中为0.1 m，且在两者相互重叠的位置处进行适当增加(如图 5(b)所示)。网格总数为120万，其中重叠网格区域的网格数为90万，背景区域的网格数为30万。

2 数值方法验证

艇体纵摇和垂荡的计算采用6-DOF刚体运动模型，仅计及沿垂直方向的平动及绕横轴方向的转动，限制了其余自由度的运动。非稳态时间离散采用二阶格式，时间步长Δt=5.0×10^-4 s。空气的密度为1.184 kg/m³，动力粘性系数为1.855×10^-5 Pa·s；水的密度为997.561 kg/m³，动力粘性系数为8.887×10^-4 Pa·s。

在计算过程中为避免气-液相在艇体表面发生伪扩散现象，考虑模型尺度、航速、时间步长等因素的影响，HRIC-VOF方法的主要参数取值为C_а=0.7、C_θ=0.5、C_UL=0.5、C_UU=6.0，其中模型试验在中国特种飞行器的高速拖曳水池中进行，试验结果可参见文献[14]。

2.1 数值波浪

图 6给出了波长λ=4.0 m、波高h_w=0.06 m时充分发展后的波形图像以及距离区域入口1.5L处(艇艏所在位置)数值计算波升与理论结果的对比。从图 6可以看出：采用本文的数值造波方法可以清晰表达所需的数值波浪，计算所得波形与理论值吻合较好，验证了本文数值造波方法的有效性。

2.2 静水阻力

图 7给出了饱和喷气下的阻力、纵倾计算值与试验值的对比，图中R_t/Δ为单位排水量阻力。从图 7可以看出：采用该数值方法计算所得阻力、纵倾值与试验结果吻合较好，阻力最小相对偏差可达1.48%，最大相对偏差为4.59%；纵倾值的最小相对偏差可达3.10%，绝对偏差值为0.14°。可见，该数值计算方法计算艇体阻力及航行姿态且具备较好的精度。图 8给出了F_V=2.80时，艇体兴波及气层形态图像的计算结果。通过试验过程中的观察发现，计算所得兴波及艇底气层在形态上与试验结果基本一致。

2.3 顶浪中的纵向运动

图 9给出了喷气流量Q=12 m³/h、容积傅汝德数F_V=2.24、波长船长比λ/L_p=2.171、波高船长比h_w/L_p=0.021 7时艇体垂荡与纵摇的时历过程。图中k为入射波的波数，ζ_a为波幅，Z_G为垂荡值，θ为纵摇角。从图中可以看出，当时间t＞15 s后，艇体运动处于稳定波动状态。

图 10为艇体运动稳定状态下无因次垂荡、纵摇与试验结果的对比。从图中可以看出，垂荡与纵摇的计算结果与试验结果吻合较好，垂荡幅值与试验结果的偏差为6.4%，纵摇幅值与试验结果的偏差为2.8%。说明该数值计算方法求解艇体在波浪中的摇荡运动具备较好的可靠性。

3 气层对波浪中阻力及纵向运动性能的影响 3.1 气层对波浪中阻力的影响

图 11给出了F_V=2.24时波浪中喷气前后的总阻力变化，此时波长λ/L_P=2.171、波高h_w/L_P=0.021 7、气流量Q=12 m³/h。从图中可以看出：在波浪作用下，艇体阻力随时间呈周期性波动，气层可使波浪中的阻力大幅度降低。

表 2给出了不同波长下喷气前后艇体平均阻力及减阻率的变化，从表中可以看出：在不同波长下，气层均可使阻力大幅度降低，减阻率可达27.53%~30.41%；从表 2中还可以看出，随着波长增长，减阻率呈增大趋势。

图 12给出了不喷气及喷气下波浪中阻力相对静水阻力的增加量随波长的变化，图中ΔR/R表示波浪中的阻力增加与静水阻力之比。从图中可以看出：喷气及不喷气条件下艇体阻力增值随波长的变化规律相类似；但气层作用下艇体在波浪中的阻力增值有所增大；随着波长增加，两者阻力增值的百分比趋于相同。

3.2 气层对纵向运动的影响

图 13给出了喷气前后艇体垂荡、纵摇时历曲线的对比，此时航速F_V=2.237，波长λ/L_P=2.17，波高h_w/L_P=0.021 7，喷气时的气流量Q=12 m³/h。图 14给出的是艇体运动处于稳定波动状态下时中心化时历曲线(去均值)的对比情况。

从图 13看出：喷气使得垂荡和纵摇运动的波动中心发生变化，气层使得垂荡运动的平均值增大，而纵摇运动的平均值略有降低。去除波动运动的中心值后发现：气层对垂荡运动幅值的影响不大，使纵摇运动幅值略有降低(如图 14所示)。

图 15给出了不同波长下，喷气前后艇体垂荡与纵摇响应的变化。从图 15中可以看出：顶浪规则波中，艇体气层对垂荡运动响应的影响不大，对纵摇运动响应有所影响；喷气使得纵摇运动响应有所减小，且随着波长增长，气层导致的纵摇响应减小量呈增大趋势。

4 艇体运动对气层形态及其稳定性的影响 4.1 气层面积

图 16给出了稳定运动状态下，艇底气层面积覆盖率随时间的变化，此时的波高、喷气量、波浪等参数与上一节计算工况完全相同。其中，气层面积覆盖率是指气层面积S_A与艇底凹槽面积S_BH之比，表征了气层在凹槽中所占有的面积比率。

从图 16可以看出：在波浪作用下，气层面积覆盖率呈周期性波动，这表征了艇体运动对气层形态的干扰作用；波长不同，气层面积覆盖率的波动情况略有不同，但差别不大。

表 3给出了不同波长下，气层平均覆盖率、波动幅值的统计结果，从表中可以看出：在不同波长下，气层在凹槽中的平均覆盖率可达78.67%~83.76%，且波动幅值仅为6.92%~9.54%，说明气层在纵向规则波中未发生严重破碎现象，具备良好的稳定性，这正是B.H.型气泡高速艇在顶浪中仍然具备良好减阻效果的主要原因。

4.2 气层形态

图 17给出了波长λ/L_P=2.171下艇底气层在一个运动周期中的变化过程，航速、波长、波高、喷气量等参数与上一节的计算相一致。图中，t₁时刻表示纵摇角位于波谷位置，t₂时刻表示纵摇角位于从下往上穿过平衡位置，t₃表示纵摇位于波峰位置，t₄表示纵摇位于从上往下穿过平衡位置。从图中可以看出：在一个纵摇运动周期过程中，气层基本上均能填满整个底部凹槽，其形态与艇底凹槽形状接近；除了气层尾部发生小量变化外，其余变化不大。这说明采用设置凹槽的艇底构型，可较好地保持气层在顶浪航行状态下的稳定性。

5 结论

1) 采用HRIC-VOF方法、Overset网格技术及刚体6-DOF模型相结合，可有效计算气层作用下高速艇波浪中的阻力、垂荡、纵摇及气层形态。

2) B.H.型气泡高速艇在顶浪中仍可取得良好的气层减阻效果，波长对减阻率有所影响，但影响不大。

3) 气层对顶浪条件下B.H.型气泡高速艇垂荡的影响甚微；气层对纵摇有所影响，喷气可使纵摇响应减小。

4) 顶浪条件下，气层在凹槽中的覆盖率可达78.67%~83.76%，气层面积变化不大，这是该型艇在波浪中仍具备良好喷气减阻效果的主要原因。