基于fluent的阻力计算

基于fluent的兴波阻力计算

本文主要研究内容

本文的工作主要涉及小型航行器在近水面航行时的绕流场及兴波模拟和阻力的数值模拟两个方面。在阅读大量文献资料的基础上，通过分析、比较上述领域所采用的理论和方法，针对目前需要解决的问题，选择合理的方法加以有机地综合运用。具体工作体现在以下几个方面：

1．本人利用FLUENT软件的前处理软件GAMBIT自主建立简单回转体潜器模型，利用FLUENT求解器进行计算，得出在不同潜深下潜器直线航行的绕流场、自由面形状及阻力系数的变化情况。

2．通过对比潜器在不同潜深情况下的阻力系数，论证了增加近水面小型航行器的深度可以有效降低阻力。通过对模型型线的改动，为近水面小型航行器的型线设计提供了一定的参考。通过改变附体形状和位置计算了附体对阻力的影响程度，为附体的优化设计提供了一定的依据。

计算模型

1

航行器粘性流场的数值计算理论

水动力计算数学模型的建立

根据流体运动时所遵循的物理定律，基于合理假设（连续介质假设）用定量的数学关系式表达其运动规律，这些表达式成为流体运动的数学模型，它们是对流体运动的一种定量模型化，称为流体运动控制方程组。根据控制方程组，结合预先给定的初始条件和边界条件，就可以求解反映流体运动的变量值，从而实现对流体运动的数值模拟预报，形成分析报告。

基于连续介质假设的流体力学中流体运动必须满足要遵循的物理定律: 1) 质量守恒定律 2）动量守恒定律 3）能量守恒定律 4）组分质量守恒方程

针对具体研究的问题，有选择的满足上述四个定律。船体的粘性不可压缩绕流运动，如果不考虑水温对水物理性质的影响，水的密度和分子粘性系数都是常数，同时没有能量的转换，就仅仅需要满足质量守恒定律、动量守恒定律。在满足这些定律下所建立的数学模型称为 Navier-Stokes方程。

另外，自由液面的存在也需要建立合适的数学模型。本文是利用 FLUENT 进行数值模拟，而软件里面关于自由液面模拟是用界面追踪方法的一种－流体体积法（VOF）,基于该方法所建立的数学模型称为流体体积分数方程。另外，高雷诺数下的水动力问题还需要考虑粘性不可压缩流体的湍流运动。对于湍流运动的数值模拟一直是流体力学数值计算的一个难点。直接数值模拟（DNS）目前还仅仅在院校中研究，而且也仅限于二维流体问题。大涡模拟（LES）向工程应用的过渡似乎还没有完成， 并且就高雷诺数问题而言， 对计算机硬件要求很苛刻。 目前，从算法的可行性、硬件要求的可实现性、完成任务所消耗时间和人力等方面看，基于湍流模型的数值计算更为工程实际所接受。本章将会对各种湍流模型加以介绍。

粘性不可压缩流体流动数学模型

连续方程

任何流动问题都必须满足质量守恒方程即：连续方程。根据连续介质假设，单位时间内流体微团的质量变化等于同时间间隔内进入微团的总净质量。按照这一定律，连续方程数学表达式写为：

（2.1）

以上是在笛卡尔直角坐标系下表示，上面给出的是瞬态可压流体连续方程。由于对于潜艇粘性流场介质的不可压缩，密度ρ 为常数，引入散度算子，则方程（2.1）变成为：

（2.2）

式中：速度矢量V= { u ,v, w }。上式为粘性不可压缩流体运动的连续方程。

动量方程

2

动量守恒方程也是任何流动系统都必须满足的定律。根据牛顿第二定律，流体微团中流体的动量对时间的变化等于微团所受外力之和，即：

（2.3） （2.4）

（2.5）

式中，p代表流体微团所受的压力;τxx 、τxy、τxz等是因分子粘性作用而产生的作用在流体微团表面上的粘性应力τ的分量；Fx、Fy、Fz表示直角坐标系下三个方向上流体微团的体积力分量，如果体积力只有重力，且Z竖直向上，则Fx=Fy=0,Fz=-ρg。

式（2.3）～（2.5）是对任何类型的流体（包括非牛顿流体）均成立的动量守恒方程。本文研究的范围属于牛顿流体，故粘性应力τ与流体的变形率成比例，有：

（2.6）

式中，μ是动力粘度系数，λ是第二粘度，一般可取λ = ?2/3，将（2.6）代入式（2.3）～（2.5）得到张量形式的动量守恒方程：

（2.7）

式（2.7）就是动量守恒方程。

方程（2.1）和（2.7）组成了控制粘性不可压缩流体运动的基本数学模型。

对于低雷诺数的层流运动，上述方程组已经可以确切描述流体运动。但湍流流动以脉动的速度场为基本特征，各速度在时间和空间上变化很快，给流场的数值模拟带来很大困难。再则，湍流是一种极度复杂的物理现象，包含无规律性，扩散性，三维涡旋波动及耗散。在实际工程计算中要对湍流进行数值模拟代价十分高昂。然而研究表明，大尺度涡在流体运动中起主要作用。由此可见，若采用时间平均、集合平均或者其他人工处理方法略去小尺度运动，将小尺度运动模型化后代入大尺度中，从而替代求解原有瞬时控制方程，就会花费较小的计算代价获得较高精度的数值解。以此为出发点，提出了将速度分解成平均值和脉动值，则瞬时速度分量u可以表达为：

（2.8）

将式（2.8）代入（2.1）和（2.7）再对时间积分就会得到下面的平均流方程。

（2.9）

3

（2.10）

（2.11）

方程（2.9）是时均形式的连续方程，方程（2.10）是时均形式的 Navier-Stokes方程。方程（2.11）为 Reynolds 应力。由于式（2.7）采用的是 Reynolds 平均法，因此方程（2.10）被成为 Reynolds 平均 Navier-Stokes 方程（Reynolds-Averaged Navier-Stokes,简称 RANS 方程）。

有式（2.9）和（2.10）组成的方程组共有五个方程（RANS方程实际是3个）现在新增了 6 个 Reynolds 应力，再加上原来 4 个时均未知量，总共9 个未知量，因此，方程组不封闭，必须引入新的湍流模型（方程）才能使方程组（2.9）和（2.10）封闭。

湍流模型

为了使雷诺平均 N-S 方程（RANS 方程）封闭可解，要根据湍流的运动规律来寻求附加的条件和关系式，这就形成了不同的湍流模型。在 FLUENT 计算软件中可以供选择的湍流模型有：一方程模型Spalart-Allmaras(S-A)、两方程模型k?ε（包括S k?ε、RNG k?ε）和k-ω（包括S k-ω和SST k-ω）以及雷诺应力模型（RSM） ，下面将本文所用到的四种湍流模型加以介绍。

标准k ? ε模型 (S k ? ε)

标准k?ε模型是典型的两方程模型，该模型是目前应用最广泛的湍流模型。k和ε是两个基本未知量，与之相对应的输运方程为：

湍流动能k方程为：

（2.12）

湍流耗散率ε方程为：

（2.13）

式中的湍流涡粘度μ t可表示为：

（2.14）

（其中=0.09，为一常数）

式中：G k式由于平均速度梯度引起的湍动能k 的产生项，G b是由于浮力引起湍动能k的产生项，Y M代表可压缩湍流中脉动扩张的贡献，C 1ε、C 2 ε和C 3 ε为经验常数，σ k和σ ε分别是与湍动能k和耗散率ε 对应的湍流普朗特数，S k 和S ε是用户定义的源项。

模型常数C 1ε、C 2 ε、C μ、σ k、σ ε的取值为：C 1ε＝1.44, C 2 ε=1.92, Cμ =0.09, σ k =1.0, σ ε=1.3

RNG k ? ε模型

RNG k?ε湍流模型是由Yakhot 及 Orzag 提出的，该模型中的 RNG 是英文

4

“renormalization group”的缩写。在RNG k?ε湍流模型中，通过在大尺度运动和修正后的粘度项体现小尺度的影响，而使这些尺度运动有系统地从控制方程中去除。所得到的k方程和ε 方程，与标准k ?ε模型非常相似：

湍流动能k方程：

（2.15）

湍流耗散率ε方程：

（2.16）

与标准k ?ε模型相比较发现，RNG k?ε模型的主要变化： 通过修正湍动粘度，考虑了平均流动中的旋转流动情况；

在ε方程中增加了一项，从而反映了主流的时均应变率E ij，这样，RNG k?ε模型中产生项不仅与流动情况有关，而且在同一问题中也还是空间坐标的函数。

模型常数C 1ε，C 2 ε由 RNG 理论： C 1ε＝1.42，C 2 ε=1.68,

其他常数：Cμ =0.0845, σ k =1.0, σ ε=1.3

Sk ? ω模型

本文采用的Sk ?ω湍流模型是基于湍流动能k 和特殊湍流动能耗散率ω 的输运方程建立起来的经验公式。是由 Wilcox在 1998 年提出的对原k ? ω模型的改进模型。

湍流动能k方程为：

（2.17）

特殊耗散率ω方程为：

（2.18）

Γk、Γ

ω表示

k、ω的有效扩散率，表示为：

（2.19） （2.20）

σ k、σω分别为湍流动能k和湍流耗散率ω的普朗特数，湍流涡粘度μ t可表示为：

（2.21）

α?为低湍流雷诺数修正系数：

（2.22）

上式中α?=βt /3，Re t为雷诺数：

（2.23）

5

以上各式中的常数取值为：

剪切应力输运k ?ω模型(SST k ?ω)

SST k ?ω湍流模型由 Menter 提出，该模式的湍流动能方程和湍流耗散率方程与标准Sk?ω模型的形式相似：

湍流动能k方程为：

（2.24）

特殊耗散率ω 方程为：

（2.25）

Γk、Γ

ω和

μ t见式(2.18)～(2.20)，常数σ k、σω表示为：

（2.26）

（2.27）

式中F 1是混合函数：

（2.28）

（2.29）

（2.30）

式中：G k式由于平均速度梯度引起的湍动能k的产生项，Gω 是由于特殊湍流动能耗散率ω的产生，Dω为横向扩散项，Y k、Yω表示湍流k、ω的消耗，S k和Sε是用户定义项。

边界条件

边界条件类型简介

流体在运动的过程中会受到边界的限制，反映到物理模型上，就是要给控制方程加一些关于变量Ui、P、k、ε相应的边界条件。最常见的线性边界条件有两大类：第一类边界条件(Dirichlet条件)和第二类边界条件(Neumann条件)。前者描述的是计算区域的边界或部分边界上变量的值，后者则描述边界

上变量梯度的法向分量值，即：

Dirichlet条件： φ=φb 在边界上

?Neumann条件： n?ф=φn 在边界上

?式中φ为任意的物理量，n表示物体表面的单位外法线矢量，φb为给定的边界上的数

6

值，φn为给定的?ф在边界上的法向分量。

对于潜器粘性绕流，入流边界是一种人工边界，它不由物体的性质决定，因而不是固定不变的，它需要取得离潜器表面足够远才能尽量地反映真实情况。入口处边界条件属于Dirichlet条件：其速度是预先给定的，一般是均匀来流条件，湍动能k和耗散率s也是预先给定的。出流边界条件则是虚拟的，出流边界到艇尾的距离也要合理确定以消除对流场计算的影响。

对于粘性流动，在固壁边界(如艇体表面)须满足对速度和湍动能k的无 滑移边界条件，即：

u=v=w=0， k=0

然而在靠近壁面的区域，由于湍动能被强烈地耗减，耗散率达到最大值，在固壁上不易给出s的边界条件，因为它在壁面上不等于零。在与壁面相邻的粘性子层中，由于粘性的影响，局部雷诺数变得很小。由于前述k-ε模型是一种高雷诺数模型，因而对粘性子层不再适用，一般采取Launder和Spalding提出的壁函数方法来处理。

使用边界条件的注意事项

1．边界条件的组合

在CFD计算域内的流动是由边界条件驱动的。从某种意义上说，求解实际问题的过程，就是将边界线或边界面上的数据，外推扩展到计算域内部的过程。因此，提供符合物理实际且适合的边界条件是极其重要的，否则，求解过程将很难进行。CFD模拟过程中迅速发散的一个最常见的原因就是边界条件选择的不合理。例如，只给定进口边界和壁面边界，而没有给定出口边界，那么，将不可能得到计算域的稳定解，CFD将越计算越发散。这样的边界条件组合显然是不合理的。

在使用出口边界时需要特别注意，该边界只在进入计算域的流动是以进口边界条件给定(如在进口给定速度和标量)时才使用，而且仅推荐在只有一个出口的计算域中使用。物理上，出口压力控制着流体在多出口间的分流情况，因此，在出口给定压力值要比给定出口条件(梯度为零)合理。将出口条件和一个或多个恒压边界结合使用是不允许的，因为零梯度的出口条件不能指定出口的流量，也不能指定出口的压力，这样将使问题不可解。

2．流动出口边界的位置选取

如果流动出口边界太靠近固体障碍物，流动可能尚未达到充分发展的状态(在流动方向上梯度为零)，这将导致相当大的误差。一般来讲，为了得到准确的结果，出口边界必须位于最后一个障碍物后10倍于障碍高度或更远的位置。对于更高的精度要求，还要研究模拟结果对出口位于不同距离时的影响的敏感程度，以保证内部模拟不受出口位置选取的影响。

3．近壁面网格

在CFD模拟时，为了获得较高的精度，常需要加密计算网格，而另一方面，在近壁面处为快速得到解，就必须将k-ε模型与结合了准确经验数据的壁面函数法一起使用。要保证壁面函数法有效，就需要使离壁面最近的以内节点位于湍流的对数律层之中，即Y+必须大于11.63(最好是在30～500之间)。这就相当于给最靠近壁面的网格到壁面的距离设定了一个下限。但是，在流动的任意位置都使上述要求得到保证常常不太可能，典型的例子就是包含回流的流动。

4．随时间变化的边界条件

这类边界条件是针对非稳定问题而言的。就是说边界上的有关流动变量并不是一成不变，而是随着时间变化的。对于这类边界条件，需要将边界条件离散成与时间步长相应的离散结果，然后存储起来，供计算到相应的时间步时调用。这类边界条件一般是与初始条件一

7

同给定的。

k和ε的计算公式

在入口、出口或远场边界流入流域的流动，FLUENT需要指定输运标量的值。

在某些情况下流动流入开始时，将边界处的所有湍流量指定为统一值是适当的。在大多数湍流流动中，湍流的更高层次产生于边界层而不是流动边界进入流域的地方，因此这就导致了计算结果对流入边界值相对来说不敏感。然而必须注意的是要保证边界值不是非物理边界。非物理边界会导致你的解不准确或者不收敛。对于外部流来说这一特点尤其突出，如果自由流的有效粘性系数具有非物理性的大值，边界层就会找不到了。

湍流强度I定义为相对于平均速度u_avg的脉动速度的均方根。

小于或等于1%的湍流强度通常被认为低强度湍流，大于10%被认为是高强度湍流。从外界测量数据的入口边界，你可以很好的估计湍流强度。例如：如果你模拟风洞试验，自由流湍流强度通常可以从风洞指标中得到。在现代低湍流风洞中自由流湍流强度通常低到0.05。

对于内部流动，入口的湍流强度完全依赖于上游流动的历史，如果上游流动没有完全发展者没有被扰动，你就可以使用低湍流强度。如果流动完全发展，湍流强度可能就达到了百之几。完全发展的管流的核心的湍流强度可以用下面的经验公式计算：

I=0.16×(雷诺数Re) -1/8

雷诺数Re=速度×当量直径×（密度÷粘度）

湍流尺度l是和携带湍流能量的大涡的尺度有关的物理量。例如在完全发展的管流中，l被管道的尺寸所限制，因为大涡不能大于管道的尺寸。L和管的物理尺寸之间的计算关系如下： L=0.07×l

其中L为管道的相关尺寸。因子0.07是基于完全发展湍流流动混合长度的最大值的。 在速度入流，压力边界条件中，需要输入如下两个物理量，其计算公式如下：

k?(Iwwall)2

??0.093/4Lk3/2

llaw为最大速度，I为湍流强其中0.09是湍流模型中指定的经验常数（近似为0.09）。w度，L为相当尺寸。

数值计算方法

网格生成

用CFD方法进行流场计算时，首先要将计算区域离散化，即划分网格。网格是CFD模型的几何表达形式，也是模拟与分析的载体。计算网格的好坏直接影响到数值计算的可行性、收敛性以及计算精度。对于复杂的CFD问题，网格生成极为耗时，且极易出错，生成网格所需时间常常大于实际CFD计算的时间。因此，有必要对网格生成方式给以足够的关注。

网格(grid)分为结构网格和非结构网格两大类。把节点看成是控制体积的代表。在离散过程中，将一个控制体积上的物理量定义并存储在该节点处。若节点排列有序，即当给出了一

8

个节点的编号后，立即可以得出其相邻节点的编号。这种网格称之为结构网格(structured grid)。结构网格是一种传统的网格形式，网格自身利用了几何体的规则形状。近几年来，还出现了非结构网格(unstructured grid)。非结构网格的节点以一种不规则的方式布置在流场中。这种网格虽然生成过程比较复杂，但却有着极大的适应性，尤其对具有复杂边界的流场计算问题特别有效。

无论是结构网格还是非结构网格，都需要按下列过程生成网格：

1．建立几何模型。几何模型是网格和边界的载体。对于二维问题，几何模型是二维面；对于三维问题，几何模型是三维实体。

2．划分网格。在所生成的几何模型上应用特定的网格类型、网格单元和网格密度对面或体进行划分，获得网格。

3．指定边界区域。为模型的每个区域指定名称和类型，为后续给定模型的物理属性、边界条件和初始条件做好准备。

网格生成是一个“漫长而枯燥”的工作过程，经常需要进行大量的试验才能取得成功。因此，出现了许多商品化的专业网格生成软件。如GAMBIT、T Grid、Geo Mesh、pre BFC和ICEMCFD等。此外，一些CFD或有限元结构分析软件，如ANSYS、I—DEAS、NASTRAN、PATRAN和ARIES等，也提供了专业化的网格生成工具。

这些软件或工具的使用方法大同小异，且各软件之间往往能够共享所生成的网格文件，例如FLUENT就可读取上述各软件所生成的网格。有一点需要说明，由于网格生成涉及几何造型，特别是3D实体造型，因此，许多网格生成软件除自己提供几何建模功能外，还允许用户利用CAD软件(AutoCAD、Pro／ENGINEER)先生成几何模型，然后再导入到网格软件中进行网格划分。因此，使用前处理软件，往往需要涉及CAD软件的造型功能。

方程离散

数值计算是将描述物理现象的偏微分方程在一定的网格系统内离散，用网格节点处的场变量值近似描述微分方程中各项所表示的数学关系，按一定的物理定律或数学原理构造与微分方程相关的离散代数方程组。引入边界条件后求解离散代数方程组，得到各网格节点的场变量分布，用这一离散的场变量分布近似代替原微分方程的解析解。当前求解流体流动和传热方程的数值计算方法比较多，如有限差分法(Finite Difference Method)，有限元法(Finite Element Method)、有限体积法(Fini te Volume Method)、边界元法、特征线法、谱方法、有限分析法和格子类方法等。每种数值计算方法有各自的特点和各自的适用范围，其中通用性比较好、应用比较广泛的是前3种。

1 有限差分法

有限差分法(Finite Difference Method,简称 FDM)是数值解法中最经典的方法。它是将求解域划分成差分格式，用有限个网格节点代替连续的求解域，然后将偏微分方程（控制方程）的导数用差商代替，推导出含有离散点上有限个未知数的差分方程组。求差分方程组（代数方程组）的解，就是微分方程定解问题的数值近似解，这是一种直接将微分问题变成代数问题的近似数值解法。

这种方法发展较早，比较成熟，较多的用于求解双曲型和抛物型问题。用它求解边界条件复杂、尤其是椭圆问题不如有限元或有限体积法方便。

2 有限元法

有限元法（Finite Element Method,简称 FEM）与有限差分法都是广泛应用的流体动力学数值计算方法。 有限元是将一个连续的求解域任意分成适当形状的许多微小单元，并与各小单元分片构造插值函数，然后根据极值原理（变分或加权余量法），将问题的控制方程转

9

化为所有单元上的有限元方程， 把总体的极值作为各单元极值之和，即将局部单元总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程就得到各节点上待求的函数值。

有限元法的基础是极值原理和划分插值，它吸收了有限差分法中离散处理的内核，又采用了变分计算中选择逼近函数并对区域进行积分的合理方法，是这两类方法互相结合、取长补短发展的结果。它具有很广泛的适应性，特别适用于几何及物理条件比较复杂的问题，而且便于程序的标准化。对椭圆型方程问题又更好的适应性。

有限元法因求解速度较有限差分法和有限体积法慢，因此，在商用CFD 软件中应用不普遍。

3 有限体积法

有限体积法（Finite Volume Method）又称控制体积法（Control Volume Method，CVM）。其基本思想是：将计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体积；将待解微分方程（控制方程）对每个控制体积积分，从而得出一组离散方程。其中的未知量是网格点上因变量φ。为了求出控制体积的积分，必须假定因变量φ在网格点之间的变化规律。从积分区域的选取方式来看，有限体积法属于加权余量法中的子域法，从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简而之，子域法加离散，就是有限体积法的基本方法。

就离散方法而言，有限体积法可视作有限元法和有限差分法的中间物。有限元法必须假定φ值在网格节点之间的变化规律（即插值函数），并将其作为近似解。有限差分法只考虑网格点上φ的数值而不考虑φ值在网格点之间如何变化。有限体积法只寻求φ值节点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定φ值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。

综上所述，有限体积法是目前在流体流动和传热问题求解中最有效的数值计算方法。有限体积法也称为控制容积积分法，是20世纪六七十年代逐步发展起来的一种主要用于求解流体流动和传热问题的数值计算方法。有限体积法是在有限差分法的基础上发展起来的，同时它又吸收了有限元法的一些优点。有限体积法与有限元法和有限差分法一样，要对求解域进行离散，将其分割成有限大小的离散格式。在有限体积法中每一网格点按一定的方式形成一个包围该节点的控制容积矿，有限体积法的关键步骤是将控制微分方程式在控制容积内进行积分。有限体积法获得的离散方程，物理上表示的是控制容积的通量平衡，方程中各项有明确的物理意义。有限体积法区域离散的节点网格与进行积分的控制容积分立。由于有限体积法的诸多优点，当今大多数CFD软件都采用它来离散求解，如PHOENICS、FLUENT、STAR-CD、NUMECA等。

SIMPLE算法

流场计算方法的本质是对离散后的控制方程组的求解。目前各种商用CFD软件普遍采纳的算法是压力耦合方程组的半隐式方法(SIMPLE算法)，它属于压力修正法的一种。SIMPLE是英文Semi—Implicit Method for Pressure-Linked Equations的缩写。SIMPLE方法由Patankar与Spalding于1972年提出，是～种主要用于求解不可压流场的数值方法，也可用于求解可压流动。它的核心是采用“猜测一修正”的过程，在交错网格的基础上来计算压力场，从而达到求解动量方程(Navier—Stokes方程)的目的。SIMPLE方法的基本思想是对于给定的压力场(它可以是假定的值或者是上一次迭代计算所得到的结果)，求解离散形式的动量方程，得出速度场。因为压力场是假定的或不精确的，由此得到的速度场一般不满足连

10

续方程，所以，必须对给定的压力场加以修正。修正的原则是与修正后的压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把由动量方程的离散形式所规定的压力与速度的关系代入连续方程的离散形式，从而得到压力修正项，由压力修正方程得出压力修正值。接着，根据修正后的压力场，求得新的速度场。然后检查速度场是否收敛，若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算。如此反复，直到获得收敛的解。

VOF方法

由于本文研究的是近水面潜体周围的绕流场，所以必须考虑自由面问题。20世纪50、60年代，自由面的数值模拟方法有了实质的进展，其中美国Los Alamos的科学家们提出和发展的格子类方法，应用效果较好。著名的格子类方法主要有PIC方法、FLIC方法和MAC方法等。从20世纪70年代开始，自由面追踪的数值方法有了进一步发展，比较著名的有标高法、线段法、VOF方法等。

VOF(Volume of Fluid)方法是使用固定网格系统捕捉两种或两种以上不相掺混流体交界面的一种数值方法。在运动界面追踪问题的数值模拟方法中，VOF方法是最为重要的方法之一，它的特点是将运动界面在空间网格内定义成一种流体体积函数，并构造这种流体体积函数的发展方程，从而界面追踪问题的目的就是如何随着主场的模拟过程，通过流体输运，精细地确定该运动界面的位置、形状和变形方向，达到追踪的目的。

VOF方法是美国学者Hirt和Nichols等人在MAC方法基础上提出的，它是一种可以处理任意自由面的方法。其基本原理是利用计算网格单元中流体体积量的变化和网格单元本身体积的比值函数，来确定自由面的位置和形状。VOF方法追踪的是网格单元中流体体积的变化，而非追踪自由液面流体质点的运动，这与Harlow和疆welch提出的MAC方法不同，后者则是从流体质点入手。相对于MAC方法，VOF法可以处理自由面重入等强非线性现象，所需计算时间更短、存储量更少，但在处理网格单元中体积比函数，的变化时，稍显繁琐，而且有一定的人为因素。VOF法同MAC法一样，以压力P和速度“，V作为独立原始变量，边界条件易处理，为计算程序的编制提供了很大方便，对于研究多相流体交界面的运动变化有着非常大的吸引力。

在VOF方法中，所有流体满足同一组动量方程，在整个计算区域上跟踪每一种流体在每个计算单元中的体积分数，根据各个时刻流体在网格单元中所占体积分数来构造和追踪自由面。假设第q种流体在单元中的体积分数为F q，第q种流体标量函数定义为f q，存在第q种流体空间点的f q值等于1，其它不被第q种流体占据点的f q值为0。在各网格单元中对f q值积分，并把此积分值除以单元的体积，得到单元的f q平均值，即网格单元中第q种流体所占据的体积分数F q。若在某时刻网格单元中F q=1，则说明单元中充满第q种流体；若F q=0，则单元中不含第q种流体；当0

本章小结

本章首先较详细地介绍了流体流动所遵守的流体动力学控制方程，包括粘性流动的基本方程、本构方程，NS方程和RANS方程等。在控制方程的基础上加入湍流模型以构成封闭方程组用以求解各速度分量和压力。对于湍流模型详细介绍了应用广泛的k-ε湍流模型，这也是本文所用的湍流模型。其后介绍了潜体绕流中应用到的几种边界条件，包括：入口边

11

界条件、出口边界条件、固壁边界条件压力边界条件等，对使用边界条件时的注意事项作了说明。对于数值计算方法介绍了CFD的网格生成，方程离散所用到的有限体积法的基本思想以及控制方程组离散后的求解所用到的SIMPLE算法。最后对于自由面问题，介绍了运动界面追踪问题的VOF法。

12

不同潜深的近水面无附体小型航行器的数值模拟

引言

当潜体在近水面航行时，潜体位于水面下航行，必然遭受到水的反作用力，产生一种与潜体运动方向相反的流体作用力，也就是潜体阻力。同时，由于潜体在近水面航行，潜体与自由面之间相互作用，自由面将产生较大变形，甚至在自由面上有波破碎等现象出现。传统的方法是不考虑流体的粘性，以势流理论为基础数值求解拉普拉斯方程。这种方法很难揭示近水面潜体直航的运动规律，因为当潜体很接近自由面时，粘性起了很重要的作用，直接数值求解Ns方程是现行的方法。Yoon&Jung利用标高法模拟了二维正方形物体在自由面附近运动的流场，同时考察了阻力系数与潜深的关系。洪方文使用VOF法对方柱在自由面附近运动的三维流场进行了数值模拟，考察了自由面与方柱在不同距离时的流场、不同长度方柱的流场以及它的阻力系数。在三维流场数值模拟中，回转体是最基本的潜体形状。本章以商用CFD软件FLUENT作为工具，利用VOF法结合k-ε湍流模型数值模拟近水面两端光滑过度的圆柱体直航的运动情形。

计算模型

为了建模简单，选择一个中间为圆柱体的简单回转体作为计算模型，圆柱体前端加上一个半球，末端采用光滑过度的椎体。圆柱体的直径为0.4m，模型的总长为3m。本人用SolidWorks建立了模型如图3.1所示。尺寸标注见图3.2。

图3.1 本人用SolidWorks 建立的计算模型

13

图3.2 SolidWorks智能尺寸标注

潜器在近水面运动时，潜器的运动与自由面是互相耦合的。为了使问题简化，研究方便，去除潜器的所有附体并假设它作一个自由度的运动即直航，且在运动过程中与静水面之间的距离保持不变，来流方向平行于潜器（回转体）的轴线。

对于简单回转体而言，由于其几何形状规则，可以直接用FLUNT前处理软件GABIT来完成几何建模。GAMBIT是目前被广泛应用的计算流体力学(CFD)的前处理器。其建立几何模型的过程遵循点、线、面、体的顺序。首先设定回转体的各个控制点，然后将点连成线，画出圆弧，最后只需保留一个封闭的半轮廓线，将此半轮廓线形成面再绕一个坐标轴（本人是x轴）旋转360度即形成一个回转体，为了研究潜器模型与自由面的相互作用，选择潜器潜深H分别为H/D=l，H/D=2，H/D=3三种情况对潜器周围绕流场进行数值模拟。其中潜深H为静水面到潜器上沿的垂向距离。潜深需在定义入流边界条件时控制边界长度，由自由面函数定义。对于深水潜器，可以假定为无限水深。在划分网格时，首先划分旋转面的面网格，采用Tri这个类型里的默认类型pave。然后划分网格疏密控制线上的网格，划分成一定比例，靠近潜器的部分密一些。最后划分体网格，采用默认的TGRID网格。根据计算要求的精确程度来定义网格大小。网格划分情况见图3.3和3.4。

设置控制区域及划分网格

完成几何建模之后，就要建立合适的控制区域。控制区域大小的选择将会影响到数值模拟的真实性和计算结果的准确性，因此必须要根据模型的尺度及计算的要求建立恰当的控制区域。控制区域的选取原则是控制区域为其边界对潜器周围绕流场的影响可以被忽略的最小区域。由于潜器近水面直航，考虑到自由面的存在，控制区域被分为上、下两部分，上部分的流体为空气，下部分的流体为水。

根据无界流场的概念和计算的需要，参考CFD粘流计算的文献，将控制区域选为长方体，其边界范围是-15m

由于潜器直航时流场左右对称，因此求解区域只需为控制区域的一半，这里取Z≥0的部分。这样计算所需的网格数减少了，大大缩短了计算所需的时间，提高了计算效率。

控制区域建立以后，与模型进行布尔并运算即可形成计算区域，对计算区域进行网格划分后需要定义边界条件。下面是网格的划分情况。

14

图3.3 中间密的TGrid网格（四条蓝色的线为网格疏密控制线）

图3.4 放大观察对称面附近的网格（可以看到许多漂亮的六面体网格）

边界条件

Fluent的计算是从边界条件开始的，边界条件的设置首先在GAMBIT软件中初步完成，再在FLUENT软件中对边界条件的参数进行具体设定。

由于自由面的存在，在GAMBIT软件中需要将控制区域的入口和出口分别划分成空气入口、水入口和空气出口、水出口。根据计算模型的坐标可知，来流方向为x轴正方向，将来流的入口设为速度入口(VELOCITY-INLET)，该边界条件专门用于不可压流动。一般假定在出口处来流未受到潜体的扰动影响，从而在出口处也可给定来流的速度分布。但实际上，在控制区域出口处的边界条件是未知的情形，故应该将控制区域的出口设定为自由出流(OUTFLOW)。但在大量计算练习时发现定义OUTFLOW边界条件时流体回流问题较为明显，即使网格细化后也难以解决此类问题，所以将出流边界重新设定为速度入流，速度设为负值，这样就相当于速度出流。总长仅仅3m的模型在30m长的距离内，对外面的流体影响已经可以近似为零，所以将出流边界设置为速度出流并非完全不可。控制区域的对称面设定为对称边界(SYMMETRY)。鉴于出流没有设置为OUTFLOW，所以可以将充满空气的那部分体积设置一个101325pa的压强（压力边界和自由出流边界不可以同时存在）。由于计算控制区域的范围取得特别大，所以可以将控制区域的其它外边界包括潜器表面设定为无滑移

15

用量纲分析法分析二维圆柱体绕流阻力FD与相关物理量ρ、V、d、μ的关系，可得

上式表明圆柱绕流阻力系数由流动Re数（ρVd/μ）唯一确定。图3.27为二维光滑圆柱体绕流的CD-Re关系曲线。根据阻力与速度的关系及阻力系数变化特点，可将曲线分为6个区域。

图3.27 二维光滑圆柱体绕流的CD-Re关系曲线

（1）Re ＜＜1，称为低雷诺数流动或蠕动流。几乎无流动分离，阻力以摩擦阻力为主，且与速度一次方成比例。

（2）1≤Re≤500，有流动分离。当Re=100，圆柱后部有一对驻涡。当Re >100时从圆柱后部交替释放出旋涡，组成卡门涡街。阻力由摩擦阻力和压差阻力两部分组成，且大致与速度的1.5次方成比例。

45

（3）500≤Re<2×10，流动分离严重，大约从Re=10起，边界层甚至从圆柱的前部就开始分离，涡街破裂成为湍流,形成很宽的分离区。阻力以压差阻力为主，且与速度的二次方成比例，即CD几乎不随Re数变化。

55

（4）2×10≤Re≤5×10，层流边界层变为湍流边界层，分离点向后推移，阻力减小，CD5

下跌，至Re = 5×10时，CD=0.3达最小值，此时的分离区最小。

56

（5）5×10≤Re≤3×10，分离点又向前移，CD回升。

6

（6）Re ＞3×10，CD与Re无关，称为自模区。

6

本文中模型特征尺寸采用3m，速度为1m/s时雷诺数为2.6×10，速度为3m/s时雷诺数

6

约为7.9×10。计算区域基本全部处于自模区。所以阻力系数大体不变。

2 粘压阻力

粘压阻力主要是由于界层分离，其次由于边界层排挤厚度的存在使流线受挤压造成压力降低而产生。影响粘压阻力的因素有：物体后端形状，后端收缩缓和者粘压阻力较低；前端型线变化亦不宜太陡；物面型线应光顺；分离点之前的流态对分离点位置影响很大；层流中法向流速分布曲线较瘦，易发生分离，分离点

21

偏前，导致分离区域扩大，Rpv增大；当层流转变为紊流时，分离点突然后移，将使阻力突然下降。

在同一流态下；若界层未发生分离，则随速度的增加,Cpv逐渐下降；在超过界分离的临界雷诺数之后，Cpv与速度的增加或Re的增加基本无关,即Cpv=常数。减小Rpv的重要方法是改进船型,使之流线型化,避免界层分离。

小结

本章模拟了不同潜深的近水面潜器直航时的绕流情况，得到了它的兴波图形和阻力系数并与相关资料对比，发现计算结果正确的可能性很大。本章还介绍了本人使用成功的Gambit建模、选择控制区域及网格划分，边界条件的选择和初始条件的计算。得到了许多直观的数据图像，如速度云图，总压力云图等，为分析阻力提供了直观的依据。

由本节的计算结果可以得到如下结论：

1 近水面航行器在水下航行时兴波阻力并不占绝对优势，即使是Fr=0.5（v=2.7m/s）附近时阻力系数也没有明显表现出兴波阻力的变化趋势。由此可知完全在水下航行的近水面航行器阻力的主要成分依然是粘性阻力（摩擦阻力和粘压阻力）。

2 当潜深大于3倍于潜器直径时，兴波阻力的影响将很小。增加潜深可以减小近水面航行器的阻力系数。

3 由压力云图可以看出圆柱体的前端与半圆球体接触的地方压力变化明显，这一点对压阻力的影响很大。由此可以推断更加光滑的流线型潜器的阻力情况可以有所改善。本文将继续研究型线方面的优化。

22

对潜器型线的优化和对附体形状的优化

型线优化

流线型的计算模型及边界条件和计算条件

在上一章发现圆柱形潜器首部压力变化明显的研究基础上，本文又建立了如下流线型的3D计算模型。该模型针对上述研究做出很大优化。模型及网格划分情况如图4.1。

图4.1 模型对称面附近的六面体网格

该模型型线大体为流线型，总长仍然是3米，最宽处仍然是0.4米。

模型的控制区域仍然同前一章一样，边界条件方面来流采用速度入口，出流采用自由出流（outflow），其余均采用了壁面条件，计算介质为水，用来模拟深水直航时的情况。与之形成对比的前一章所提到的模型采用同样的边界条件和初始条件进行精确计算。

计算条件大体与前一章一样，采用unsteady状况k-ε模型并考虑了重力影响，用SIMPLE算法。壁面粗糙度设置仍然是默认的0.5。其余参数由计算公式计算得到。

计算结果

本章将在这一节里对两种型线的潜器的计算结果进行对比，得到阻力系数更小的型线并分析原因，为设计更加快速的潜器提供依据。首先对比潜器的压力情况，分析两种型线的潜器表面压力有什么不同，从而得到阻力不同的原因。然后对比计算所得的阻力系数，确认结果与事实相符。两种模型的压力云图如图4.2-4.5。

23

图4.2 v=1m/s时的动压力云图对比

图4.3 v=2.5m/s时的动压力云图对比

图4.4 v=1m/s时的总压力云图对比

图4.5 v=2.5m/s时的总压力云图对比

由两个模型的动压力云图对比可知，流线型的潜器并没有像圆柱形的潜器那样在首部有很明显的压力变化，它的压力分布比较均匀，压力变化较为缓和，因此粘压阻力减小了。由总压云图的对比可知，流线型的潜器只有在首部那一小部分压力比较大，尾部的一小部分压力比较小，其它地方总压力分布比较均匀，变化不大。首尾压差也不大。圆柱形的潜器首部较大面积里压力很大，而且压力较流线型潜器更趋向于指向后方。所以压力的水平分量比较大。过了半球体以后，圆柱体前半截的压力变小，这样就形成了更大的压差，压阻力显然比较大些。所以总阻力也比流线型大一些。

24

图4.6 v=2.5m/s时垂直面上的动压力云图

图4.6是垂直面上的动压力图。由于重力的影响，云图并不是对称的。这张图反映了模型周围流体的压力变化情况。由图可知压力流畅内压力分布均匀且变化（形成压差）的区域不大。

计算结果表明，当速度为1m/s时，圆柱形潜器的阻力系数为0.0063，而流线型潜器的阻力系数仅为0.003，是圆柱形潜器阻力系数的47.6%，还不及圆柱形潜器阻力系数的一半；当速度为2m/s时，圆柱形潜器的阻力系数为0.0061，流线型潜器的阻力系数仅为0.0027，是圆柱形潜器阻力系数的44.3%；当速度为2.5m/s时，圆柱形潜器的阻力系数为0.00596，而流线型潜器的阻力系数仅为0.00262，是圆柱形潜器阻力系数的44.0%，同样不及圆柱形潜器阻力系数的一半。由此可见型线对潜器阻力系数的影响之大。这并不奇怪，因为潜器的阻力系数主要是粘压阻力系数，摩擦阻力系数所占比例不是很大。另外由上述数据可以看出潜器速度越高，流线型潜器阻力系数减小越明显。

不同附体对潜器阻力系数的影响

计算模型、边界条件和计算条件

为了使潜器能够在水中按照人们的意愿行驶，潜器要有完成控制方向等功能的附体。本文仅对两种不同形状的“舵”对阻力系数的影响做了简单的计算研究。一种是方形舵，另一种是较方形舵超过潜器直径较少的三角形舵。虽然两种舵形状不同，但面积相差不大。三角形舵沿着潜器轴线布置得更长一些。两种模型的网格划分情况与前面几个雷同，这里不再给出。其形状将在计算结果的压力云图中看到。

边界条件及计算条件和上一节中完全相同。

计算结果

依然按照前面所述的方法，首先对比压力云图。图4.7-4.10。

25

图4.7 v=1m/s时的动压力云图对比

图4.8 v=2.5m/s时的动压力云图对比

图4.9 v=1m/s时的总压力云图对比

图4.10 v=2.5m/s时的总压力云图对比

由压力云图可知，方形附体的潜器压力分布均匀，说明附体的压力与潜器主体的压力比较接近。而三角形附体附体部分的压力较主体小，主体部分压力图类似于无附体潜器。附体部分不难看出附体前端的压力大，后端的压力小，附体前后存在压力差。可见附体同样存在压阻力。但由于附体厚度很小，压阻力不大。

图4.11 v=2.5m/s时水平面上的总压力云图对比

图4.11是潜器外部流场的压力情况。

26

计算结果表明，当速度为1m/s时，方附体潜器的阻力系数为0.007，比无附体潜器大了0.0007。而三角附体潜器的阻力系数仅为0.0063。这并不意味着三角附体的潜器阻力下降了，而是因为计算作用的湿面积大了。所以虽然阻力系数并没有增加，但总阻力却是增加的。这可能是因为三角附体处在潜器的尾流中，并没有超出潜器的半径太多。而附体本身由于宽度很小，压力差所形成的阻力不大。主要是摩擦阻力起作用。当速度为2.5m/s时，方附体潜器的阻力系数为0.00653，比无附体潜器大了0.00057 。而三角附体潜器的阻力系数仅为0，00553，比无附体潜器低了0.00043 。而在fluent计算中，这一数字甚至可以比得上波动幅值，所以不能认为三角附体潜器的阻力系数比无附体的小。如刚刚论述的那样，虽然阻力系数没有增加，但阻力是增加的。因为计算湿面积增大了。

由以上分析可以得出结论：超出潜器半径较多的方附体的附体阻力系数比相对超出潜器半径很少的附体阻力系数大。但由于附体宽度很小，主要受摩擦阻力影响。即使是附体增加了阻力系数，增加的数值也不算大。所以考虑到“舵效”，采用方形舵可能更为合适。

小结

本章通过对两种不同型线的无附体潜器的阻力系数计算，压力云图的对比得出了潜器的型线对潜器阻力系数影响很大的结论。因为潜器在水下航行时，不仅仅会受到摩擦阻力的影响，还会受到粘压阻力的影响并且影响很大。粘压阻力很容易受到形状的影响，所以又称为形状阻力。流线型的物体阻力系数小是个常识，本文只是验证了这一点。单从快速性考虑，潜器采用流线型的型线要比采用普通圆柱形好得多。

本章的第二个任务是论证了附体的形状（相对于潜器主体的布置而言）对潜器快速性的影响。由计算结果可以看到将附体放在潜器尾部并尽量少的超出其直径是合适的。因为这样对改善潜器的快速性有效。这一点是符合常识的，也是很容易想得到的。不过就潜器整体的阻力系数而言，本文所采用的附体对整体阻力系数的影响并不大，所以当考虑到附体的功能本身时，采用上述布置和形状可能并不合适，但就快速性而言，附体尽量不要超出主体的特征尺寸太多为好。

27

结论

本文研究的是近水面航行器的快速性，使用商用软件fluent模拟了不同潜深的潜器周围绕流场的不同和兴波情况的变化，并通过计算不同潜深的航行器的阻力系数得到了近水面航行器可以通过增加潜深的方法来降低阻力。这是因为航行于水下的航行器在近水面航行时，自由面将受到扰动而兴起波浪，兴波是需要能量的，这部分能量当然要由潜器来提供。所以近水面航行的潜器的阻力会因为自由面的扰动消耗能量而增加。但是，潜深越深，兴波就越不明显，因此波浪消耗的能量也会越小。所以增加潜深可以减小航行器的阻力系数。

本文还通过对比不同型线的潜器的绕流场，压力云图和阻力系数论证了型线对于潜器快速性的重要性以及原因。这虽是常识，但依然值得研究。通过本文的论证，发现流线型的潜器阻力系数得到了大大的改善，对于要求快速的近水面航行器而言，这一点是万万不可忽略的。对于并不要求十分快速的稍大型潜器，从节能的角度考虑，采用流线型或减少过于剧烈的型线变化是完全应该的。

另外本文还考虑到潜器附体形状对潜器阻力系数的影响。通过计算发现将附体放在潜器尾部并尽量少的超出其直径是合适的。因为这样对改善潜器的快速性是有效的。当然考虑到附体的功能本身，采用上述布置和形状可能并不合适，但就快速性而言，附体尽量不要超出主体的特征尺寸太多为好。

总之，本文综合考虑了多种增加潜器快速性的方法并利用fluent软件做了计算，论证等工作，得到了若干种行之有效的降低近水面小型潜器阻力系数的方法。

28

本文来源：https://www.bwwdw.com/article/7ob.html