CFD动态网格技术在水下航行体非定常操纵运动预报中的应用研究

李迎华，吴宝山，张 华

（中国船舶科学研究中心，江苏 无锡214082）

1 引 言

对于水下滑翔机[1-2]这一类以剩余浮力为主要驱动方式的水下低速航行体而言，其在进行纵倾与浮力调节时，所作的运动是非定常的，因此对其非定常操纵运动的预报研究，是指导其设计开发的基础。

由于水下低速航行体非定常操纵运动中其雷诺数低于临界雷诺数，水动力系数会随时间发生变化，因而采用传统的拘束模型试验[3-5]，以及基于准定常假设的数值方法[6]对其进行运动预报时，则需要对各种雷诺数下的水动力分别进行测试和计算，这显然需要大量的工作，而且对于拘束模型试验来说，其代价也是极其昂贵的。近年来，随着商业软件Fluent的不断发展，基于RANS求解器的CFD数值预报方法已成为一种经济、有效的手段。因此，研究非定常操纵运动的CFD数值模拟技术对水下低速航行体设计研究具有重要的意义。

对于数值模拟时水动力及运动参数需要实时求解的这类操纵运动，Fluent专门提供了一种研究手段—动网格技术。在模拟时将运动方程通过用户自定义函数（UDF）导入至Fluent中，即可实现水动力和运动参数实时自主求解。因此，动网格技术适合于水下低速航行体非定常操纵运动的直接模拟。

目前动网格技术大多应用于周期性操纵运动的数值计算，如垂荡、升沉等，所采用的动网格划分形式一般是全非结构网格[7]，由于这种网格划分形式在近壁面粘性边界层内难以生成大展弦比的非结构网格，且网格数量较结构网格多出许多，因而将其推广到非周期性操纵运动预报上有一定的局限性，因此，近年来，结合结构网格和非结构网格优势的混合网格技术得到了重视和发展，张来平等[8]提出了三棱柱/四面体/矩形三维动态混合网格划分形式，并用无粘流的圆球超声速绕流算例进行了验证，得到了与定常计算相一致的结果。

本文在此基础上对三维动态混合网格模型进行了发展，提出了六面体/四面体三维动态混合网格模型，并以此建立了基于动网格技术的CFD数值预报方法。本文首先以带攻角均匀直航的滑行状态为研究对象，应用动网格技术对一典型水下低速航行体缩比模型主体及带附体的水动力进行了计算以验证其水动力的计算精度，为便于比较，同时给出了静态网格计算结果，进而将其应用于一水下滑翔机自航模非定常滑翔运动的预报，以验证动网格技术的应用能力。

2 基于动网格技术的数值计算方法

2.1 控制方程

时均的不可压缩连续性方程和N-S方程（RANS方程）如（1）、（2）式所示：

2.2 湍流模型

对于本文研究的水下低速航行体而言，其运动过程中纵向雷诺数较低。RNG k-ε湍流模型提供了针对低雷诺数的有效粘性的微分解析式，且具备数值稳定性好和处理流线弯曲程度较大的流动更优等优点[9]，本文数值计算即采用RNG k-ε湍流模型。

2.3 几何守恒律

动网格需满足几何守恒律，控制体积的时间导数由（3）式计算：

式中：nf是控制面的数目，是j面的面积矢量。

式中：δVj是控制面j在一个时间步长内的体积更新量。

2.4 边界条件

（1）远场边界条件：采用零剪应力壁面边界条件处理。

（2）物面条件：满足壁面黏附条件，壁面处流体速度与运动边界速度相同。

2.5 动态混合网格生成方法

本文的六面体/四面体动态混合网格模型如下建立：在近壁面处依次生成贴体六面体网格和四面体网格共同构成运动区域随运动边界一起运动 （四面体网格与六面体网格之间采用交界面过渡），外场生成四面体网格并时刻保持静止，在运动区域和静止区域之间采用四面体网格过渡构成变形区域。当物体位移较小时，采用“类弹簧原理”对变形区域网格节点进行松弛，当物体位移较大导致变形区域网格质量较差时，则在局部重新生成网格。

本文采用上述改进的动态混合网格模型对一典型水下低速航行体缩比模型表面及其周围流动区域进行了网格划分，此模型主体为细长回转体，长2 150mm，直径300mm，在主体中后部布置有一对后掠的水平翼，翼两端间距1 200mm。网格划分的具体形式为：主体C-O型（纵向-周向），附体C-O型（展向—弦向）；计算域大小为：沿航行体首前端、航行体径向分别延伸3L距离，沿航行体尾后端延伸4L距离（L为航行体纵向长度）。单独主体网格数65万，带附体网格数120万。图1、2为带附体表面及近壁面流场区域网格划分示意图。为便于对比分析，也给出了静态网格划分示意图，如图3所示。单独主体的静态网格数量32万，带附体的静态网格数量57万。

2.6 动网格更新算法

动网格更新算法有弹簧近似光滑算法（Spring Based Smoothing）、动态分层算法（Dynamic Layering）和局部重构算法（Local Remeshing）。对于本文研究的四面体/六面体混合网格，选用弹簧近似光滑算法和局部重构算法作为动网格更新算法。

2.7 数值离散和求解

（1）时间项的离散：采用直接一阶隐式离散。

（2）空间项的离散：其中扩散项以中心差分格式差分，对流项采用二阶迎风格式。

应用SIMPLE法处理压力速度耦合问题，离散方程以Gauss-Seidel迭代法求解。

2.8 时间步长的选取

由于时间项采用一阶隐式离散，是无条件稳定的，因此时间步长的选取可以不用考虑离散格式的稳定性，理论上可以取较大的时间步长，但过大的时间步长将导致网格更新过快，进而出现负体积。本文从网格更新角度推导出时间步长计算公式如下：

不考虑控制面上的面积更新量则由（4）式可得对于单元j上的当地时间步长：

式中：αj=hj/Δhj，Δhj为单元j上的高度更新量，hj为单元j上更新前的高度。

全场时间步长应该取全动区域最小的即

式中：h1表示动区域最小的第一层网格高度；vel为边界运动的合速度；α为一系数。

2.9 双振荡圆板验证算例

为对动网格技术的可行性进行验证，本文选取文献[10]中的双强迫振荡圆板算例进行了计算，并与文献给出的结果进行了对比。采用全非结构网格划分形式对其进行网格划分，网格数量12万，如图4所示。

采用弹簧近似光滑算法和局部重构算法作为动网格更新算法，将双板振荡速度以用户自定义函数（UDF）形式导入到Fluent中，时间步长 Δt=0.1s。

本文计算得到的阻力系数变化周期T=25s，阻力系数幅值A=13.7；而文献给出的阻力系数变化周期T=25s，阻力系数幅值A=14.2，可见本文采用的动网格技术可以得到相对比较精确的计算结果，因此可以应用于水下滑翔机非定常运动的预报。

3 水下滑翔机非定常运动的CFD计算研究

本文以中国船舶科学研究中心自主开发的“前哨”号水下滑翔机自航模[2]为对象，应用动网格技术对其非定常操纵运动进行了预报，以进一步验证动网格技术的应用。

操纵运动预报的基础在于水动力的预报，因此本文首先以带攻角均匀直航的滑行状态为研究对象，应用动网格技术对前述典型水下低速航行体缩比模型的主体及带附体水动力的预报精度进行了验证，为便于比较分析，同时给出了静态网格计算结果。为分析水动力计算精度对运动预报的影响，本文还基于带附体动网格计算的水动力对定常滑翔运动参数进行了分析。

3.1 基于动网格技术的水动力预报验证

通过用户自定义函数 （UDF）将水下低速航行体运动速度实时导入到Fluent中，时间步长Δt=0.002s，应用基于动态混合网格模型的动网格技术依次计算了雷诺数Re=2.4×106，攻角分别为6°、8°、10°、12°四种工况下主体及带附体模型的水动力，得到了无因次水动力系数和水动力中心作用点位置lα′（lα′=-M′/Z′），并与风洞模型试验结果[12]进行了对比。

（1）主体模型水动力计算结果

主体模型静态网格和动网格计算结果如表1所示，水动力随攻角变化曲线如图5-8所示。

表1 水下航行体主体模型水动力计算结果（Re=2.4×106）Tab.1 The calculated results of the hydrodynamic forces for the main body of the underwater vehicle model（Re=2.4×106）

由计算结果可知，除个别工况外，主体模型静态网格计算的轴向力系数相差16%以内，垂向力系数相差8%以内，俯仰力矩系数相差5%以内，水动力作用中心位置相差5%以内；动网格计算的轴向力系数相差5%以内，垂向力系数、俯仰力矩系数均相差10%以内，水动力作用中心位置相差13%以内。

（2）带附体模型水动力计算结果及对运动预报的影响

带附体模型操纵性水动力计算结果如表2所示，水动力计算结果随攻角变化曲线如图9-12所示（由于静态网格在12°时失速，因此只给出前三个攻角对应的水动力值）。

由计算结果可知，除个别工况外，带附体模型静态网格计算的轴向力系数相差16%以内，垂向力系数和俯仰力矩系数相差14%以内，水动力作用中心位置相差5%以内；动网格计算的轴向力系数相差10%以内，垂向力系数相差9%以内，俯仰力矩系数相差13%以内，水动力作用中心位置lα′相差10%以内。

表2 水下航行体带附体模型水动力计算结果（Re=2.4×106）Tab.2 The calculated results of the hydrodynamic forces for the main body with appendages of the underwater vehicle model（Re=2.4×106）

由主体及带附体模型水动力计算结果可知，基于动态混合网格模型的动网格技术水动力预报结果与试验吻合较好，且可以达到与静态网格相类似的精度。

同时从计算结果也可看出：数值计算结果与试验仍存在一定差距，这在很大程度上是因为实际流动处于层流和湍流的混合区，而本文采用适用于充分发展湍流（Re>1×107）的湍流模型去解算，因而会带来一定的误差。

为分析水动力计算精度对运动预报的影响，本文基于上述动网格计算的水动力和试验测得的水动力，分别计算了给定雷诺数Re=2.4×106下带附体模型稳定滑翔时不同剩余浮力系数k（k=ΔB/B，ΔB为剩余浮力，B为水下静均衡后的浮力） 对应的潜浮角 χ （χ=θ-α，θ为纵倾角，α 为攻角）、水平速度分量uξ和垂直速度分量uζ，计算结果曲线如图13-15所示。

由计算结果来看：动网格预报的结果与试验预报结果相比，潜浮角相差0.4°以内，水平速度和垂直速度相差0.007m/s以内，表明了上述动网格水动力计算误差对运动预报的影响比较小。

3.2 基于动网格技术的非定常运动的CFD预报

对“前哨”号水下滑翔机自航模垂直面内滑翔运动的预报可采用如图16所示的两类坐标系：地面坐标系E-ξζ和运动坐标系o-xz，其中原点o位于水下航行体浮心处，ox轴指向首部为正，oz轴垂直ox轴指向下方为正。图中X，Z，M分别为此水下航行体所受的轴向、垂向水动力和纵倾力矩。

此水下滑翔机垂直面内非定常滑翔运动方程可简化为如（7）～（12）式所示的形式：

式中：W（t）为水下滑翔机重量；xg（t）为水下滑翔机重心纵向坐标；zg（t）为水下滑翔机重心垂向坐标；B为水下滑翔机所受浮力。

将滑翔运动方程（7）～（12）进行时间一阶向前差分离散，并把离散后的方程以UDF形式导入到Fluent中，即可实现Fluent中每一时间步的水动力和运动参数实时自主求解。其求解流程如图17所示。

本文采用基于动态混合网格模型的动网格技术对“前哨”号水下滑翔机非定常下潜运动中的一段进行了数值模拟，结果如图18、19所示。

从数值计算结果与试验的比较来看：基于动网格技术数值方法的预报结果与试验符合较好，且变化趋势一致。

4 结论和展望

本文基于传统的动态混合网格模型，发展了一种改进的三维动态混合网格模型，以此建立了非定常操纵运动数值预报方法，并以带攻角均匀直航的滑行状态为研究对象，应用动网格技术对一典型水下低速航行体缩比模型主体及带附体的水动力进行了计算，结果与风洞模型试验结果吻合较好，从而验证了基于动态混合网格模型的动网格技术对操纵运动水动力的预报精度；本文进而将其应用于一水下滑翔机自航模非定常运动的预报，预报结果与试验符合较好，初步验证了动网格技术的应用能力和本文建立的方法的有效性。同时动网格技术具有更广阔的应用前景，主要表现在：

（1）动网格技术可以应用于水下多体相对运动，潜艇的坐底、靠岸等的数值模拟。

（2）动网格技术可以应用于非定常操纵运动过程中流体动力的精细研究。

[1]马冬梅,马 铮.水下滑翔机运动及水动力特性研究[R].无锡:中国船舶科学研究中心科技报告,2007.

[2]李 龙,张 华,陈季军,马 铮,吴宝山.水下滑翔机自航模演示试验研究[C].中国造船工程学会学术论文集,2008,No.200804.

[3]Ma Ling,Cui Weicheng.Simulation of dive motion of a deep manned submersible[J].Journal of Ship Mechanics,2004,8(3):31-38.

[4]张 华,吴宝山.扁平型潜水器空间大机动数学模型研究[C].中国造船工程学会学术论文集,2004,No.200405.

[5]谢俊元,马 岭,胡 震.载人深潜器运动操纵模拟[J].中国造船,2006,47(2):62-69.

[6]王长涛,俞建成,吴利红,封锡盛.水下滑翔机器人运动机理仿真与实验[J].海洋工程,2007,25(1):64-69.

[7]黄昆仑,庞永杰,苏玉民,朱 军.潜器线性水动力系数计算方法研究[J].船舶力学,2008,12(5):697-703.

[8]张来平,张涵信.三维动态混合网格技术及非定常计算方法[C].全国流体力学青年研讨会论文集,2005.

[9]Fluent 6.3.26 User’s Guide.Fluent Inc.[CP].

[10]顾 正.二维单圆柱、双圆柱绕流问题及三维振荡板运动的数值模拟[D].上海:上海交通大学,2007.

[11]施生达.潜艇操纵性[M].北京:国防工业出版社,1995.

[12]倪歆韵,潘子英.低速水下航行体风洞模型试验[R].无锡:中国船舶科学研究中心科技报告,2008.