于向阳 姚凌虹 孟庆昌 刘巨斌 张志宏

（1海军航空大学青岛校区 青岛 266000）（2海军工程大学 武汉 430033）

1 引言

潜艇等水下航行体粘性流场和水动力对其操纵性预报至关重要。而潜艇等水下航行体作操纵运动时，其水动力主要集中在主体部分，同时复杂的三维流动（潜艇非线性水动力的主要来源）也主要发生在主体上，所以主体水动力预报的准确程度备受关注。6：1椭球体类似于潜艇等水下航行体的主体，其外形虽然简单，但其带攻角的绕流流动却是一个复杂、典型的三维分离和涡流流动，而对其流动预报的结果可以间接检验数值计算模型对潜艇等水下航行体操纵性预报的能力［1～4］。

本文应用CFD软件对6：1椭球体粘性流场和水动力进行数值模拟。采用基于非结构化网格的有限体积法离散控制方程和湍流模式，速度-压力耦合采用SIMPLE（Semi-implicitMethodforPressure-LinkedEquations）方法处理，对流项采用二阶迎风格式，采用不同湍流模型，包括k-ε湍流模型，SA一方程湍流模型，k-ω湍流模型，SST两方程湍流模型，得到5°、10°、15°、20°、25°、30°这6个攻角的水动力数值计算结果，同时分析了不同计算域范围，不同网格密度等对数值计算结果的影响，并与文献［5］中的实验结果进行了比较。

2 数值计算方法

2.1 数值离散方法

采用基于非结构化网格的有限体积法，数值求解定常的RANS方程，在SA模式中，生成项所采用的形式为

其中

2.2 计算对象与边界条件设置

计算对象为长细比6：1的椭球体，采用直角坐标系，其原点位于椭球体的几何中心，x轴与椭球体的对称轴重合，方向与无攻角时来流方向一致，z轴竖直向上，y轴水平。考虑到对称性，仅对半椭球体流动进行数值计算；计算域为半圆柱体区域，范围：-a≤x/L≤b（a，b取值在各节中不同），-2≤z/L≤2，0≤ y/L≤2，其中L为椭球体长度，取值L=1.37，计算域如图1所示。

图1 模型的边界条件设置

计算域的边界条件：半圆柱面的下半部分AEFD（z≤0）和半圆柱体的左半圆面AEB（x/L=a处）定义为速度入口，给定速度Uo，方向沿x轴正向；半圆柱面的上半部分EBCF（z＞0）和半圆柱体的右半圆面CFD（x/L=b处）定义为压力出口，考虑到在迭代过程中，有可能存在数值上的反向流动，设置出口回流方向的方式为垂直于出口面；物面为无滑移壁面，速度分量和湍动能为零；对称面上，法向速度分量为零，平行于对称面的速度分量的法向导数和所有标量的法向导数为零。

在所定义的坐标系下，规定来流的z方向速度分量为正时为正攻角，计算雷诺数Re=4.2×106，收敛准则为所有求解变量的无量纲残差下降4个量级。

3 边界层处理

3.1 边界层理论

边界层问题表现为近壁湍流问题，在工程实际中流体具有粘性，因而流体在物面上的速度相对于物面应为零，即无滑移条件。而在紧贴物面的薄层区域内，分子输运现象起主导作用，致使一般湍流模型的很多基本假设不再成立；同时近壁处的速度、压力等变量沿物面法向方向梯度很大。在数值计算中，为了获得可靠的数值计算结果，不得不考虑上述问题。一种解决方法是在物面附近，使用壁面函数法来处理，但针对具有高雷诺数的流动问题；另一种解决方法是在物面附近，增加网格密度，以取得较密的节点，但这样会增加计算量。

在数值计算中，引入无量纲值u+和 y+，进行近壁处网格的布置。

其中，τω是壁面剪切力，y是沿物面法向的距离。

如采用近壁模型法，要求第一层网格y+＜5，而且在法向雷诺数Re＜200的范围内，至少布置10层网格，但这势必要大大增加计算量。

图2 近壁处理方法

3.2 网格划分

计算域网格由边界层网格和体网格两部分组成。

CFD对网格选取有一定要求，一方面考虑到粘性效应近壁面采用较密的贴体网格，另一方面网格的疏密程度与流场参数（如压力、速度等）的变化梯度趋于一致［6～11］。为捕捉壁面附近的强分离流动，使椭球体表面附近的压力，速度等物理参数有较好的模拟效果，采用三角形网格划分底层物面，沿垂直于壁面方向，边界层进行加密处理，以求得高质量的贴体网格，如图3所示；体网格采用基于四面体的非结构化网格形式，以适应椭球体表面的曲率变化［12～15］。估算 y+并确定第一层网格控制点离壁面的距离△y/L，控制在1×10-5。

图3 边界层网格

4 不同湍流模型下的粘性流场和水动力数值计算结果比较

湍流模型的实质是为RANS方程中出现的雷诺应力寻求计算方法。考虑到求解问题的复杂性及计算条件的局限性，必须选择合适的湍流模型进行数值计算。

4.1 湍流模型的选择

由于湍流流动问题的差异性，没有一种湍流模型能够对所有工程实际问题具有通用性，即使比较常用的湍流模型也只是对某一类或某一些特定问题适用，因此需要对所研究的特定问题选取不同湍流模型，以验证其对计算对象的适用性。对湍流模型的选择通常注意以下几个方面：1）所计算问题的特殊性，物理问题是否可以简化；2）所拥有的计算资源，对计算时间的限制；3）对计算结果的要求，如计算精度的要求等［16～18］。

本文采用的湍流模型如表1所示，SA模型适合应用于存在翼形、壁面边界层等场合的湍流流动问题，不适合射流类等具有自由剪切流动的湍流问题，本文中的椭球体绕流问题属于其适用范围；标准k-ε模型应用较为广泛，具有较好的稳定性，同时也具有较高的计算精度，适合高雷诺数湍流流动，但对旋流等各向异性较强的流动具有局限性；标准k-ω模型，包含低雷诺数影响和可压缩性影响，适用于受壁面限制的流动附着边界层湍流流动计算；剪切应力输运（SST），综合了k-ω模型在近壁区计算和k-ε模型在远场计算的优越性，适用于带逆压梯度的流动计算。湍流模型的计算量与求解的方程数量有关，进行数值模拟时，需要考虑计算资源的合理利用；各湍流模型的计算速度依次为SA一方程模型，k-ε和k-ω双方程模型。

表1 采用的湍流模型

4.2 计算结果分析

本节适度选取网格间距（网格参数如表2），对不同湍流模型的水动力数值计算结果进行比较，结果包括力矩系数、升力系数CL和阻力系数随攻角的变化情况，并与国外实验结果进行了比较，结果表明表1中列举的几种湍流模型都能够反映流场的水动力和粘性分离特性。

表2 网格参数

图4～图6（图中exp为实验值，下文不再另行说明）分别给出了椭球体升力系数CL和阻力系数CD和力矩系数CM随攻角的变化情况，其中力矩的参考点取在椭球体的中心。与实验值比较，升力系数均小于实验值，但SA（Cυ=20）和k-ω模型的计算值与实验值较为接近，且较好地模拟了实验值的趋势。

图4 升力系数与攻角的关系

图5 阻力系数与攻角的关系

依据实验本身曲线的总体趋势来分析，随着攻角的增加，非线性作用明显；在力矩系数的比较中，SST和k-ε模型与实验值符合均较好，SA（Cυ=20）模型存在偏差；从阻力系数的计算结果比较中，可以看到除k-εstand模型外，其余四种模型在低攻角时，变化较为平直，大攻角时曲线陡直，非线性作用增强。综合比较5种湍流模型，对椭球体的水动力数值计算效果，SA和k-ω模型具有较好的计算效果。

图6 力矩系数与攻角的关系

图7 ～图11分别给出了椭球体20°攻角时x/L=0.272截面速度向量图，在背风面可以清晰地看到有分离涡出现，反映了不同湍流模型下流场的粘性分离特性。

图7 SST湍流模型20°攻角时x/L=0.272截面速度向量图

图8 k-ωstand湍流模型20°攻角时x/L=0.272截面速度向量图

图9 k-εstand湍流模型20°攻角时x/L=0.272截面速度向量图

图10 SA（Cυ=20）湍流模型20°攻角时x/L=0.272截面速度向量图

图11 SA（vorticity）湍流模型20°攻角时x/L=0.272截面速度向量图

5 不同计算域范围对水动力数值计算结果的影响

5.1 计算域范围的选取

采用数值方法求解控制方程时，需要在空间区域上进行离散，以得到离散方程组。空间上离散控制方程，需要进行网格划分，在划分网格时，需要选取适当的计算域，防止边界条件产生不利影响。经过查找大量的相关领域的文献，发现对计算域范围并没有定量的讨论，本文在参考相关文献的基础上选择不同的计算域范围进行了数值计算，并对不同计算域的数值计算结果进行了分析和比较。

表3 不同计算域参数

5.2 计算结果分析

图12～图14分别给出了椭球体升力系数CL和阻力系数CD和力矩系数CM随攻角的变化情况，其中力矩的参考点取在椭球体的中心。与实验值比较，升力系数的数值计算结果均小于实验值，但较好地模拟了实验值的趋势，不同的计算域对应的数值计算结果基本一致，但随着攻角的增加，特别是大攻角处，存在差异性，计算域选取-2.5≤＜x/L≤4.5时，已收敛；在力矩系数的比较中，不同计算域对应的计算结果与实验值符合均较好，计算域选取-2.5≤＜x/L≤4.5时，在大攻角处，显示出较好的模拟效果；对阻力系数的数值计算中，大攻角时曲线陡直，不同计算域都捕捉到了其非线性作用。

图12 升力系数与攻角的关系

图13 阻力系数与攻角的关系

图14 力矩系数与攻角的关系

6 不同网格密度对水动力数值计算结果的影响

随着CFD技术应用的不断深入，许多复杂的工程问题得到解决，特别是几何外形复杂的实际问题更多的提到日程上来，在这里我们不得不提到CFD前处理软件，为生成高效率的满足计算精度的网格所做的贡献。

6.1 网格无关解

用有限体积法离散得到的代数方程与原偏微分方程比较，存在截断误差。截断误差的大小与网格的间距（或网格数目）有关，网格间距越小（网格数目越多）截断误差就越小。因此就减小截断误差而言，网格的数目越多越好。

事实上，由于计算机内存和计算时间的限制，我们不可能无止境地减小网格间距，增加网格数目，同时，由于数值计算中总会有舍入误差的出现，没有必要片面地对截断误差作太苛刻的要求。

工程实际中，为了获得网格无关解，划分网格时，需要对同一个计算问题划分多套网格，通过逐渐增加网格数目，来考察网格数目对计算结果的影响。

如果在两个不同数目的网格上获得的解基本相同，说明此时进一步加大网格数目，增加网格密度，已无必要，或者说我们已经得到了与网格无关的解。

6.2 网格密度的选择

为了数值模拟得到与网格无关的收敛解，应该选取密度不同的网格进行数值计算，作比较以检验网格无关性。网格密度的选取一般要考虑所使用计算机的运行速度和内存容量，网格的生成往往占用整个CFD任务70%以上的工作量，所以在不影响计算精度的前提下，尽量节省计算资源，选取适度的网格密度。通过参考国内外文献，确定以下三种网格进行数值计算（如表4所示）。

表4 不同网格密度

6.3 计算结果分析

图15～图17分别给出了椭球体升力系数CL和阻力系数CD和力矩系数CM随攻角的变化情况，其中力矩的参考点取在椭球体的中心，即坐标系原点处。数值计算方法如前所述，湍流模式均采用SST模式。与实验值比较，升力系数与力矩系数的计算值均小于实验值，在小攻角下与实验值符合较好，只是在30°时，不同的网格密度在数值模拟结果存在偏差，不同网格密度计算结果

图15 升力系数与攻角的关系

总体趋势保持一致，网格密度较小（Δx/L=0.05）的计算值在本例中并没有表现出众的模拟效果，网格密度取Δx/L=0.07时，已经取得较好的计算值，与网格密度为Δx/L=0.05时相比，相对偏差在1.5%左右，表明此时计算结果已收敛，网格密度选取已合理；阻力系数的比较中，不同网格密度的结果总体趋势保持一致。

图16 阻力系数与攻角的关系

图17 力矩系数与攻角的关系

综合考虑计算资源与数值计算精度的要求，网格密度取Δx/L=0.07是可行的，以相对偏差1.5%的结果换取百万网格的代价，解决现有问题是可以考虑的。

7 结语

6：1椭球体类似于潜艇等水下航行体的主体，对其流动预报的结果可以间接检验数值计算模型对潜艇等水下航行体操纵性预报的能力。本文应用CFD软件对长细比为6：1椭球体水动力进行数值模拟。采用基于非结构化网格的有限体积法离散控制方程和湍流模式，对流项采用二阶迎风格式，速度-压力耦合采用SIMPLE方法处理，得到5°、10°、15°、20°、25°、30°这6个攻角的水动力数值计算结果，并与已知的实验结果进行了比较。

分析了不同湍流模式对数值计算结果的影响，综合比较5种湍流模型对椭球体的水动力数值计算结果，认为SA和SST模型具有较好的计算效果。

分析了不同计算域，不同网格密度等对数值计算结果的影响。不同的计算域对应的数值计算结果基本一致，但随着攻角的增加，特别是大攻角处，存在差异性，计算域选取-2.5≤＜x/L≤4.5时，略显优势；不同的网格密度在数值模拟上存在微小偏差，较小网格密度（Δx/L=0.05）的计算结果在本例中并没有表现特别出众的模拟效果，网格密度取Δx/L=0.07时，已经取得较好的计算结果，与网格密度为Δx/L=0.05的计算结果相比，相对偏差在1.5%左右，本文认为以相对偏差1.5%的结果换取百万网格的代价，解决现有问题是可以考虑的。本文的计算结果和结论为后续湍流模式、计算域范围和网格密度的研究提供了参考依据。