0 引　言

随着人口增长和能源消耗，海洋资源的开发和利用日益重要，水下航行器在海洋探测和研究、资源开发以及军事应用等领域发挥着关键作用，其设计与优化也成为重中之重。计算流体力学方法（CFD）基于许多优秀软件，如Fluent、CFX、SHIPFLOW、STAR-CCM+等不仅能够获得与试验数据相符的计算结果，同时计算周期短，成本低廉^[1]。粘性流域中，固体壁面附近的边界层在垂直壁面方向存在很大的速度梯度，影响CFD计算结果的准确性。因此针对边界层不同第一次高度、最少网格层数、最终层高度分别对计算结果产生的影响有助于提高计算精度。

目前，对近壁区域的处理方法主要有2种：一种是壁面函数法，可跳过粘性底层和过渡层直接在对数律层上进行求解，其本质是在湍流核心区配合${k}-{\varepsilon }$湍流模型，在近壁区域采用壁面函数法进行求解^[2]。另一种是近壁模型法，使用该方法时需要在壁面附近划分足够细的网格，使之可以求解粘性影响的区域^[3]。

在${k}-{\varepsilon }$湍流模型与壁面函数相结合的模式下，资丹等^[4]利用3种不同边界层网格尺度，对比分析泵站流场模拟计算精度，证明了边界层网格尺度对模拟计算结果产生较大影响。Griffin等^[5]采用局部重构方法计算了大范围雷诺数非平衡流中边界层的厚度，为非平衡流条件下边界层的选取提供了理论依据。Hasanuzzaman等^[6]研究了在中、高雷诺数条件下，在湍流动能的传输和耗散过程中，壁面湍流边界层及其尺度对相关摩擦阻力的作用。Medjnoun等^[7]使用多尺度粗糙表面研究了粗糙度等级对湍流边界层的影响，实验表明改变粗糙度边界层会发生实质性变化。

众多研究探索了边界层的选取以及边界层整体对模拟计算产生的影响，早在2016年Posa等^[8]使用30亿网格对SUBOFF模型绕流场进行计算，以捕捉湍流边界层内的流动，但并未对边界层第一次高度、最少网格层数、最终层高度对计算精度的影响展开深入研究，并且在实际模拟计算中，受到计算资源、效率、时间成本等因素的制约往往不能完全采用计算结果的第一层高度与层数。

本文基于STAR-CCM+软件，采用Standard ${k}-{\varepsilon }$模型对SUBOFF航行阻力进行预测计算，并将不同边界层因素的计算结果与实验结果进行对比分析。

1 数值计算理论 1.1 控制方程 1.1.1 质量守恒方程

质量守恒是流体质点在流场中运动的必然定律。质量守恒定律是指流体质点在流场运动过程中质量始终保持不变，即流过流体微团表面的质量通量等于流体微团质量的变化率。

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0 。$

(1)

式中：$\rho$为流体质点密度；$t$为时间变量；$u、v、w$为流体质点在笛卡尔坐标系下$i、j、k$方向上的速度。此次假设流体为定常流动的理想流体，因此式(1)可表示为：

$\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0。$

(2)

1.1.2 动量守恒方程

流体质点在流动过程中需要满足牛顿第二定律即动量守恒方程，在流体力学中常用纳维-斯托克斯方程描述动量守恒定律，如下式：

$\frac{{{\mathrm{d}}}V}{{{\mathrm{d}}}t}=F-\frac{1}{{\rho }}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{P}+\frac{1}{{\rho }}\mathrm{d}\mathrm{i}\mathrm{v}\left(2{\mu }S\right)-\frac{2}{3}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\left({\mu }\mathrm{d}\mathrm{i}\mathrm{v}V\right)。$

(3)

式中：$V$为流体质点的速度；$t$为时间变量；$F$为流体质点所受体力；$\mathrm{\mu }$为流体质点的粘度系数；$P$为流体质点压力；$S$为应变率张量。

1.2 边界层计算

运用壁面函数是以对数律为基础计算边界层规律，其完全忽略了粘性底层和过渡层。因此在CFD模拟过程中需要保证边界层第一层网格高度位于对数律的范围内，从而保证计算结果的准确性。边界层第一层高度可由下列公式计算：

$\mathit{Re}=\frac{\rho ud}{\mu } ，$

(4)

${C}_{f}=0.0576{\left(Re\right)}^{\frac{-1}{5}} ，$

(5)

${\tau }_{w}=\frac{1}{2}{C}_{f}\rho {U}_{\infty }^{2}，$

(6)

${v}^{*}=\sqrt{\frac{{\tau }_{w}}{\rho }}，$

(7)

$y=\frac{{y}^+v}{{v}^{*}} 。$

(8)

式中：$Re$为雷诺数；$\rho$为流体质点密度；$d$为所计算物体特征长度；$u$为流体质点的速度；$\mu$为流体质点的粘度系数；${C}_{f}$为壁面摩擦系数；${\tau }_{w}$为壁面剪切应力；${U}_{\infty }$为流体质点的动力粘度；${v}^{*}$为边界层内流体速度；$y$为边界层第一层高度；${y}^{+}$为到壁面处的无量纲距离；$v$为流体质点运动粘度。

2 水力性能计算 2.1 模型选择

对航行器阻力进行分析时^[9]，多数学者以SUBOFF模型为例^[10]，且不同航速下的实验阻力值已经给出，其全附体模型主要包括1个水滴形的主体部分、1个围壳和4个对称的尾翼，其长度为L=4 356 m，最大直径为D=508 m。本次同样选择SUBOFF模型，根据实验数据确定模拟计算的准确性。其外形及具体尺寸数据如表1及图1所示。

2.2 计算域的确定及网格划分

本次模拟SUBOFF直航状态下的航行阻力无其他复杂运动，故确定计算域为圆柱形的静域。由于SUBOFF模型尺寸较大，雷诺数较高且本次讨论边界层对计算结果的影响，需要足够大的静域以避免因计算域的大小产生误差。确定计算域的长度为11L，半径为2.5L，入口距航行器首部3L，出口距航行器尾部7L。其中L为SUBOFF模型的长度。静域由速度入口、压力出口、壁面、SUBOFF模型外壁面组成，如图2所示。

将该模型导入STAR-CCM+中采用非结构网格进行网格划分，选择收敛速度更快的棱柱层加切割体网格^[11]。为保证计算精度，在确定模型表面合适尺寸的同时，采用非常慢体积增长率细化模型周围区域网格尺寸，以更高分辨率捕捉模型周围流体流动状态。取SUBOFF模型表面$Y+$为30，根据边界层计算章节所述方法计算得边界层第一层高度为0.00035 m，层数为20层。整体及边界层网格划分情况如图3和图4所示。

2.3 边界条件及湍流模型的选择

静域边界条件设置如下：入口为速度入口并设置来流方向及速度，出口为压力出口，静域壁面为滑移壁面与来流相对静止。SUBOFF模型表面设置为壁面。因Standard ${k}-{\varepsilon }$湍流模型对高雷诺数和复杂几何周围外部流动求解效果较好^[12]，且具有很好的收敛速率和较低的内存需求，故选择Standard ${k}-{\varepsilon }$湍流模型结合壁面函数进行定常水力计算及分析。

2.4 网格无关性验证

在仿真模拟计算之前需进行网格无关性验证，以排除网格数量对计算结果的影响，确保仿真计算结果的可靠性。因此对SUBOFF模型以表2所示4种不同数量的网格进行网格无关性验证，并且前3种网格数量之差在2.5倍以上^[13]，计算SUBOFF在10 kn下航行的阻力。

通过对比4种网格的计算结果可知，除粗网格外其余网格计算结果变化不大，变化率在3%以内，且对比细和细+网格变化率更不足1%。如图5所示，将模拟结果和实验结果进行对比分析，中、细和细+3等种网格均与实验结果相差不大，中网格的计算误差最低为2.64%。根据模拟结果的准确性以及计算成分等方面考虑，最终选择中密度网格进行计算。

3 边界层影响分析

设置速度入口进速为10 kn，以计算所得边界层为基础改变第1层高度、层数、边界层总厚度和表面第一层Y+值，将模拟计算所得结果与实验结果对比，进而分析边界层不同因素对阻力预报的影响。相关计算结果见表3、表4、图6～图8。表3为仅改变边界层网格第1层高度的计算结果；表4为仅改变边界层最少网格层数的计算结果。A₀表示计算所得基础边界层。

计算结果表明，改变边界层确实对仿真准确度产生影响。其中，A₀边界层与实验数据的误差为2.64%。在仅改变第1层高度的情况下，如图6(a)所示，相比降低第1层高度，增加第1层高度与误差曲线的斜率明显更低，对计算的影响更小。当第1层高度增加到0.006约为A₀边界层第1层高度的17倍时，计算结果误差对比A₀边界层仅增加了2.01%，且误差仍在5%以内；当边界层第1层高度降低到0.0001约为A₀边界层第1层高度的1/3时，误差上升至10.56%。如图6(b)所示，在仅改变层数的情况下，适当增加或减少最少网格层数对计算误差的影响不大，26层到14层范围内误差均在5%以下，但最少网格层数过低到10层时，误差上升到8.11%。

当边界层第1层高度和层数发生改变时，边界层总厚度随之改变。如图7所示，将上述2种方案的边界层总厚度进行对比分析，总体来说边界层总厚度过低对于2种方案计算误差都会大幅上升，并且第一层高度降低所导致的总厚度减少对计算的影响更加剧烈，误差更高。当第1层高度相比A₀边界层降低了0.00005 m，总厚度相比A₀边界层降低了0.0003 m时，误差由2.64%到5.56%，误差上升了110.6%；当最少网格层数由A₀的20层降低到14层，总厚度相比A₀边界层降低了0.0356 m，误差由2.64%到3.63%，误差上升了37.5%。然而，增加边界层总厚度对于2种方案而言皆对计算结果影响较低，当边界层总厚度为1.12 m时，误差仅上升2%。2种方案近似厚度下的误差并不完全相同，边界层总厚度分别为0.0934、0.0949时，误差为0.52%、2.5%。

对于表面第1层Y+而言，如图8所示，配合${k}-{\varepsilon }$模型的Y+分布位于30～300之间误差相对较低，误差在0.52%～8.11%之间。结合2种方案分析，Y+与误差之间无线性关系，但Y+在30～300区间外误差上升明显，Y+由31.4到20.1误差上升幅度为66.7%。

由以上计算结果表明，增加第一层高度以及降低最少网格层数对计算结果的影响相对较低，为探究增加第一层高度同时降低层数产生的影响，制定如表5方案。

此方案均增加了第1层高度以及减少了最少网格层数，并且保证边界层总厚度大于A₀边界层总厚度或小幅度低于A₀边界层总厚度。最低边界层总厚度为0.0296相比A₀边界层降低了47.4%。如图9所示，各类型边界层计算误差均低于2.64%，可见以此方向改变边界层较为理想。

4 结　语

本文基于STAR-CCM+软件，运用${k}-{\varepsilon }$湍流模型结合壁面函数探究了边界层网格对水下航行器阻力预报的影响，为航行器阻力预报提供了一定参考。

1）仅改变第1层高度的情况下，增加第1层高度对计算的影响更小。当第1层高度增加到0.006约为A₀边界层第1层高度的17倍时，计算误差仍在5%以内。

2）仅改变最少网格层数的情况下，最少网格层数在20%以内浮动对计算的影响较低，计算误差在5%以内。

3）计算误差由第1层高度与最少层数共同影响，边界层总厚度不起决定性作用，且总体上表面第1层Y+值在30～300区间，计算误差相对较低在10%以内。

4）增加第1层高度的同时降低最少网格层数，并且保证边界层总厚度与计算总厚度相似或高度计算总厚度，计算误差较为理想皆在3%以下。