基于fluent的阻力计算

1、基于fluent的兴波阻力计算本文主要研究内容本文的工作主要涉及小型航行器在近水面航行时的绕流场及兴波模拟和阻力的数值模拟两个方面。在阅读大量文献资料的基础上，通过分析、比较上述领域所采用的理论和方法，针对目前需要解决的问题，选择合理的方法加以有机地综合运用。具体工作体现在以下几个方面： 1本人利用FLUENT软件的前处理软件GAMBIT自主建立简单回转体潜器模型，利用FLUENT求解器进行计算，得出在不同潜深下潜器直线航行的绕流场、自由面形状及阻力系数的变化情况。2通过对比潜器在不同潜深情况下的阻力系数，论证了增加近水面小型航行器的深度可以有效降低阻力。通过对模型型线的改动，为近水面小型航行

2、器的型线设计提供了一定的参考。通过改变附体形状和位置计算了附体对阻力的影响程度，为附体的优化设计提供了一定的依据。计算模型航行器粘性流场的数值计算理论水动力计算数学模型的建立根据流体运动时所遵循的物理定律，基于合理假设（连续介质假设）用定量的数学关系式表达其运动规律，这些表达式成为流体运动的数学模型，它们是对流体运动的一种定量模型化，称为流体运动控制方程组。根据控制方程组，结合预先给定的初始条件和边界条件，就可以求解反映流体运动的变量值，从而实现对流体运动的数值模拟预报，形成分析报告。基于连续介质假设的流体力学中流体运动必须满足要遵循的物理定律: 1) 质量守恒定律 2）动量守恒定律 3）能量

3、守恒定律 4）组分质量守恒方程 针对具体研究的问题，有选择的满足上述四个定律。船体的粘性不可压缩绕流运动，如果不考虑水温对水物理性质的影响，水的密度和分子粘性系数都是常数，同时没有能量的转换，就仅仅需要满足质量守恒定律、动量守恒定律。在满足这些定律下所建立的数学模型称为 Navier-Stokes方程。另外，自由液面的存在也需要建立合适的数学模型。本文是利用 FLUENT 进行数值模拟，而软件里面关于自由液面模拟是用界面追踪方法的一种流体体积法（VOF）,基于该方法所建立的数学模型称为流体体积分数方程。另外，高雷诺数下的水动力问题还需要考虑粘性不可压缩流体的湍流运动。对于湍流运动的数值模拟一直

4、是流体力学数值计算的一个难点。直接数值模拟（DNS）目前还仅仅在院校中研究，而且也仅限于二维流体问题。大涡模拟（LES）向工程应用的过渡似乎还没有完成， 并且就高雷诺数问题而言， 对计算机硬件要求很苛刻。 目前，从算法的可行性、硬件要求的可实现性、完成任务所消耗时间和人力等方面看，基于湍流模型的数值计算更为工程实际所接受。本章将会对各种湍流模型加以介绍。粘性不可压缩流体流动数学模型连续方程任何流动问题都必须满足质量守恒方程即：连续方程。根据连续介质假设，单位时间内流体微团的质量变化等于同时间间隔内进入微团的总净质量。按照这一定律，连续方程数学表达式写为：（2.1）以上是在笛卡尔直角坐标系下表示

5、，上面给出的是瞬态可压流体连续方程。由于对于潜艇粘性流场介质的不可压缩，密度 为常数，引入散度算子，则方程（2.1）变成为： （2.2）式中：速度矢量V= u ,v, w 。上式为粘性不可压缩流体运动的连续方程。动量方程动量守恒方程也是任何流动系统都必须满足的定律。根据牛顿第二定律，流体微团中流体的动量对时间的变化等于微团所受外力之和，即： （2.3） （2.4） （2.5）式中，p代表流体微团所受的压力;xx 、xy、xz等是因分子粘性作用而产生的作用在流体微团表面上的粘性应力的分量；Fx、Fy、Fz表示直角坐标系下三个方向上流体微团的体积力分量，如果体积力只有重力，且Z竖直向上，则Fx=F

6、y=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）组成了控制粘性不可压缩流体运动的基本数学模型。 对于低雷诺数的层流运动，上述方程组已经可以确切描述流体运动。但湍流流动以脉动的速度场为基本特征，各速度在时间和空间上变化很快，给流场的数值模拟带来很大困难。再则，湍流是一

7、种极度复杂的物理现象，包含无规律性，扩散性，三维涡旋波动及耗散。在实际工程计算中要对湍流进行数值模拟代价十分高昂。然而研究表明，大尺度涡在流体运动中起主要作用。由此可见，若采用时间平均、集合平均或者其他人工处理方法略去小尺度运动，将小尺度运动模型化后代入大尺度中，从而替代求解原有瞬时控制方程，就会花费较小的计算代价获得较高精度的数值解。以此为出发点，提出了将速度分解成平均值和脉动值，则瞬时速度分量u可以表达为： （2.8）将式（2.8）代入（2.1）和（2.7）再对时间积分就会得到下面的平均流方程。 （2.9） （2.10） （2.11）方程（2.9）是时均形式的连续方程，方程（2.10）是时

8、均形式的 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

9、-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可表

10、示为： （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 11.44, C 2 =1.92, Cµ =0.09, k =1.0, =1.3RNG k 模型RNG k湍流模型是由Yakhot 及 Orzag 提出的，该模型中的 RNG 是英文“renormalization gr

11、oup”的缩写。在RNG k湍流模型中，通过在大尺度运动和修正后的粘度项体现小尺度的影响，而使这些尺度运动有系统地从控制方程中去除。所得到的k方程和 方程，与标准k 模型非常相似：湍流动能k方程： （2.15）湍流耗散率方程： （2.16）与标准k 模型相比较发现，RNG k模型的主要变化： 通过修正湍动粘度，考虑了平均流动中的旋转流动情况； 在方程中增加了一项，从而反映了主流的时均应变率E ij，这样，RNG k模型中产生项不仅与流动情况有关，而且在同一问题中也还是空间坐标的函数。 模型常数C 1，C 2 由 RNG 理论： C 11.42，C 2 =1.68, 其他常数：Cµ =

12、0.0845, k =1.0, =1.3Sk 模型本文采用的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）以上各式中的常数取值为：剪切应力输运k 模型(SST k ) SST k 湍流模

13、型由 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是用户定义项。边界条件边界条件类型简介流体在运动的过程中会受到边界的限制，反映到物理模型上，就是要给控制方程加一些关于变量U i、P、k、

14、相应的边界条件。最常见的线性边界条件有两大类：第一类边界条件(Dirichlet条件)和第二类边界条件(Neumann条件)。前者描述的是计算区域的边界或部分边界上变量的值，后者则描述边界上变量梯度的法向分量值，即：Dirichlet条件： =b 在边界上Neumann条件： =n 在边界上式中为任意的物理量，表示物体表面的单位外法线矢量，b为给定的边界上的数值，n为给定的在边界上的法向分量。对于潜器粘性绕流，入流边界是一种人工边界，它不由物体的性质决定，因而不是固定不变的，它需要取得离潜器表面足够远才能尽量地反映真实情况。入口处边界条件属于Dirichlet条件：其速度是预先给定的，一般是均

15、匀来流条件，湍动能k和耗散率s也是预先给定的。出流边界条件则是虚拟的，出流边界到艇尾的距离也要合理确定以消除对流场计算的影响。对于粘性流动，在固壁边界(如艇体表面)须满足对速度和湍动能k的无滑移边界条件，即：u=v=w=0， k=0 然而在靠近壁面的区域，由于湍动能被强烈地耗减，耗散率达到最大值，在固壁上不易给出s的边界条件，因为它在壁面上不等于零。在与壁面相邻的粘性子层中，由于粘性的影响，局部雷诺数变得很小。由于前述k-模型是一种高雷诺数模型，因而对粘性子层不再适用，一般采取Launder和Spalding提出的壁函数方法来处理。使用边界条件的注意事项1边界条件的组合在CFD计算域内的流动是

16、由边界条件驱动的。从某种意义上说，求解实际问题的过程，就是将边界线或边界面上的数据，外推扩展到计算域内部的过程。因此，提供符合物理实际且适合的边界条件是极其重要的，否则，求解过程将很难进行。CFD模拟过程中迅速发散的一个最常见的原因就是边界条件选择的不合理。例如，只给定进口边界和壁面边界，而没有给定出口边界，那么，将不可能得到计算域的稳定解，CFD将越计算越发散。这样的边界条件组合显然是不合理的。在使用出口边界时需要特别注意，该边界只在进入计算域的流动是以进口边界条件给定(如在进口给定速度和标量)时才使用，而且仅推荐在只有一个出口的计算域中使用。物理上，出口压力控制着流体在多出口间的分流情况，

17、因此，在出口给定压力值要比给定出口条件(梯度为零)合理。将出口条件和一个或多个恒压边界结合使用是不允许的，因为零梯度的出口条件不能指定出口的流量，也不能指定出口的压力，这样将使问题不可解。2流动出口边界的位置选取如果流动出口边界太靠近固体障碍物，流动可能尚未达到充分发展的状态(在流动方向上梯度为零)，这将导致相当大的误差。一般来讲，为了得到准确的结果，出口边界必须位于最后一个障碍物后10倍于障碍高度或更远的位置。对于更高的精度要求，还要研究模拟结果对出口位于不同距离时的影响的敏感程度，以保证内部模拟不受出口位置选取的影响。3近壁面网格在CFD模拟时，为了获得较高的精度，常需要加密计算网格，而另

18、一方面，在近壁面处为快速得到解，就必须将k-模型与结合了准确经验数据的壁面函数法一起使用。要保证壁面函数法有效，就需要使离壁面最近的以内节点位于湍流的对数律层之中，即Y+必须大于11.63(最好是在30500之间)。这就相当于给最靠近壁面的网格到壁面的距离设定了一个下限。但是，在流动的任意位置都使上述要求得到保证常常不太可能，典型的例子就是包含回流的流动。4随时间变化的边界条件这类边界条件是针对非稳定问题而言的。就是说边界上的有关流动变量并不是一成不变，而是随着时间变化的。对于这类边界条件，需要将边界条件离散成与时间步长相应的离散结果，然后存储起来，供计算到相应的时间步时调用。这类边界条件一般

19、是与初始条件一同给定的。k和的计算公式在入口、出口或远场边界流入流域的流动，FLUENT需要指定输运标量的值。在某些情况下流动流入开始时，将边界处的所有湍流量指定为统一值是适当的。在大多数湍流流动中，湍流的更高层次产生于边界层而不是流动边界进入流域的地方，因此这就导致了计算结果对流入边界值相对来说不敏感。然而必须注意的是要保证边界值不是非物理边界。非物理边界会导致你的解不准确或者不收敛。对于外部流来说这一特点尤其突出，如果自由流的有效粘性系数具有非物理性的大值，边界层就会找不到了。湍流强度I定义为相对于平均速度u_ avg的脉动速度的均方根。小于或等于1%的湍流强度通常被认为低强度湍流，大于1

20、0%被认为是高强度湍流。从外界测量数据的入口边界，你可以很好的估计湍流强度。例如：如果你模拟风洞试验，自由流湍流强度通常可以从风洞指标中得到。在现代低湍流风洞中自由流湍流强度通常低到0.05。对于内部流动，入口的湍流强度完全依赖于上游流动的历史，如果上游流动没有完全发展者没有被扰动，你就可以使用低湍流强度。如果流动完全发展，湍流强度可能就达到了百之几。完全发展的管流的核心的湍流强度可以用下面的经验公式计算：I=0.16×(雷诺数Re) -1/8雷诺数Re=速度×当量直径×（密度÷粘度）湍流尺度l是和携带湍流能量的大涡的尺度有关的物理量。例如在完全发展的管

21、流中，l被管道的尺寸所限制，因为大涡不能大于管道的尺寸。L和管的物理尺寸之间的计算关系如下： L=0.07×l其中L为管道的相关尺寸。因子0.07是基于完全发展湍流流动混合长度的最大值的。在速度入流，压力边界条件中，需要输入如下两个物理量，其计算公式如下：其中0.09是湍流模型中指定的经验常数（近似为0.09）。为最大速度，I为湍流强度，L为相当尺寸。数值计算方法网格生成用CFD方法进行流场计算时，首先要将计算区域离散化，即划分网格。网格是CFD模型的几何表达形式，也是模拟与分析的载体。计算网格的好坏直接影响到数值计算的可行性、收敛性以及计算精度。对于复杂的CFD问题，网格生成极为耗

22、时，且极易出错，生成网格所需时间常常大于实际CFD计算的时间。因此，有必要对网格生成方式给以足够的关注。网格(grid)分为结构网格和非结构网格两大类。把节点看成是控制体积的代表。在离散过程中，将一个控制体积上的物理量定义并存储在该节点处。若节点排列有序，即当给出了一个节点的编号后，立即可以得出其相邻节点的编号。这种网格称之为结构网格(structured grid)。结构网格是一种传统的网格形式，网格自身利用了几何体的规则形状。近几年来，还出现了非结构网格(unstructured grid)。非结构网格的节点以一种不规则的方式布置在流场中。这种网格虽然生成过程比较复杂，但却有着极大的适应性

23、，尤其对具有复杂边界的流场计算问题特别有效。无论是结构网格还是非结构网格，都需要按下列过程生成网格：1建立几何模型。几何模型是网格和边界的载体。对于二维问题，几何模型是二维面；对于三维问题，几何模型是三维实体。2划分网格。在所生成的几何模型上应用特定的网格类型、网格单元和网格密度对面或体进行划分，获得网格。3指定边界区域。为模型的每个区域指定名称和类型，为后续给定模型的物理属性、边界条件和初始条件做好准备。网格生成是一个“漫长而枯燥”的工作过程，经常需要进行大量的试验才能取得成功。因此，出现了许多商品化的专业网格生成软件。如GAMBIT、T Grid、Geo Mesh、pre BFC和ICEM

24、CFD等。此外，一些CFD或有限元结构分析软件，如ANSYS、IDEAS、NASTRAN、PATRAN和ARIES等，也提供了专业化的网格生成工具。这些软件或工具的使用方法大同小异，且各软件之间往往能够共享所生成的网格文件，例如FLUENT就可读取上述各软件所生成的网格。有一点需要说明，由于网格生成涉及几何造型，特别是3D实体造型，因此，许多网格生成软件除自己提供几何建模功能外，还允许用户利用CAD软件(AutoCAD、ProENGINEER)先生成几何模型，然后再导入到网格软件中进行网格划分。因此，使用前处理软件，往往需要涉及CAD软件的造型功能。方程离散数值计算是将描述物理现象的偏微分方程

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

26、较好、应用比较广泛的是前3种。1 有限差分法 有限差分法(Finite Difference Method,简称 FDM)是数值解法中最经典的方法。它是将求解域划分成差分格式，用有限个网格节点代替连续的求解域，然后将偏微分方程（控制方程）的导数用差商代替，推导出含有离散点上有限个未知数的差分方程组。求差分方程组（代数方程组）的解，就是微分方程定解问题的数值近似解，这是一种直接将微分问题变成代数问题的近似数值解法。 这种方法发展较早，比较成熟，较多的用于求解双曲型和抛物型问题。用它求解边界条件复杂、尤其是椭圆问题不如有限元或有限体积法方便。2 有限元法 有限元法（Finite Element M

27、ethod,简称 FEM）与有限差分法都是广泛应用的流体动力学数值计算方法。 有限元是将一个连续的求解域任意分成适当形状的许多微小单元，并与各小单元分片构造插值函数，然后根据极值原理（变分或加权余量法），将问题的控制方程转化为所有单元上的有限元方程， 把总体的极值作为各单元极值之和，即将局部单元总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程就得到各节点上待求的函数值。 有限元法的基础是极值原理和划分插值，它吸收了有限差分法中离散处理的内核，又采用了变分计算中选择逼近函数并对区域进行积分的合理方法，是这两类方法互相结合、取长补短发展的结果。它具有很广泛的适应性，特别适用于几何及物理条件

28、比较复杂的问题，而且便于程序的标准化。对椭圆型方程问题又更好的适应性。 有限元法因求解速度较有限差分法和有限体积法慢，因此，在商用CFD 软件中应用不普遍。3 有限体积法 有限体积法（Finite Volume Method）又称控制体积法（Control Volume Method，CVM）。其基本思想是：将计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体积；将待解微分方程（控制方程）对每个控制体积积分，从而得出一组离散方程。其中的未知量是网格点上因变量。为了求出控制体积的积分，必须假定因变量在网格点之间的变化规律。从积分区域的选取方式来看，有限体积法属于加权余量法中的子域法，从

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

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

31、分的控制容积分立。由于有限体积法的诸多优点，当今大多数CFD软件都采用它来离散求解，如PHOENICS、FLUENT、STAR-CD、NUMECA等。SIMPLE算法流场计算方法的本质是对离散后的控制方程组的求解。目前各种商用CFD软件普遍采纳的算法是压力耦合方程组的半隐式方法(SIMPLE算法)，它属于压力修正法的一种。SIMPLE是英文SemiImplicit Method for Pressure-Linked Equations的缩写。SIMPLE方法由Patankar与Spalding于1972年提出，是种主要用于求解不可压流场的数值方法，也可用于求解可压流动。它的核心是采用“猜测一

32、修正”的过程，在交错网格的基础上来计算压力场，从而达到求解动量方程(NavierStokes方程)的目的。SIMPLE方法的基本思想是对于给定的压力场(它可以是假定的值或者是上一次迭代计算所得到的结果)，求解离散形式的动量方程，得出速度场。因为压力场是假定的或不精确的，由此得到的速度场一般不满足连续方程，所以，必须对给定的压力场加以修正。修正的原则是与修正后的压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把由动量方程的离散形式所规定的压力与速度的关系代入连续方程的离散形式，从而得到压力修正项，由压力修正方程得出压力修正值。接着，根据修正后的压力场，求得新的速度场。然后检查速度场

33、是否收敛，若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算。如此反复，直到获得收敛的解。VOF方法由于本文研究的是近水面潜体周围的绕流场，所以必须考虑自由面问题。20世纪50、60年代，自由面的数值模拟方法有了实质的进展，其中美国Los Alamos的科学家们提出和发展的格子类方法，应用效果较好。著名的格子类方法主要有PIC方法、FLIC方法和MAC方法等。从20世纪70年代开始，自由面追踪的数值方法有了进一步发展，比较著名的有标高法、线段法、VOF方法等。VOF(Volume of Fluid)方法是使用固定网格系统捕捉两种或两种以上不相掺混流体交界面的一种数值方法。在运动界面

34、追踪问题的数值模拟方法中，VOF方法是最为重要的方法之一，它的特点是将运动界面在空间网格内定义成一种流体体积函数，并构造这种流体体积函数的发展方程，从而界面追踪问题的目的就是如何随着主场的模拟过程，通过流体输运，精细地确定该运动界面的位置、形状和变形方向，达到追踪的目的。VOF方法是美国学者Hirt和Nichols等人在MAC方法基础上提出的，它是一种可以处理任意自由面的方法。其基本原理是利用计算网格单元中流体体积量的变化和网格单元本身体积的比值函数，来确定自由面的位置和形状。VOF方法追踪的是网格单元中流体体积的变化，而非追踪自由液面流体质点的运动，这与Harlow和疆welch提出的MAC

35、方法不同，后者则是从流体质点入手。相对于MAC方法，VOF法可以处理自由面重入等强非线性现象，所需计算时间更短、存储量更少，但在处理网格单元中体积比函数，的变化时，稍显繁琐，而且有一定的人为因素。VOF法同MAC法一样，以压力P和速度“，V作为独立原始变量，边界条件易处理，为计算程序的编制提供了很大方便，对于研究多相流体交界面的运动变化有着非常大的吸引力。在VOF方法中，所有流体满足同一组动量方程，在整个计算区域上跟踪每一种流体在每个计算单元中的体积分数，根据各个时刻流体在网格单元中所占体积分数来构造和追踪自由面。假设第q种流体在单元中的体积分数为F q，第q种流体标量函数定义为f q，存在第

36、q种流体空间点的f q值等于1，其它不被第q种流体占据点的f q值为0。在各网格单元中对f q值积分，并把此积分值除以单元的体积，得到单元的f q平均值，即网格单元中第q种流体所占据的体积分数F q。若在某时刻网格单元中F q=1，则说明单元中充满第q种流体；若F q=0，则单元中不含第q种流体；当0<F q<1时，则该单元为存在不同流体的交界面。VOF法将流体体积分数设定在单元中心，根据相邻网格的流体体积分数和网格单元四边上的流体速度来计算流过指定单元网格的流体体积，借此来确定单元内下时刻的流体体积分数，并根据相邻网格单元的流体体积分数来确定自由面单元内自由面的位置和形状。本章小

37、结本章首先较详细地介绍了流体流动所遵守的流体动力学控制方程，包括粘性流动的基本方程、本构方程，NS方程和RANS方程等。在控制方程的基础上加入湍流模型以构成封闭方程组用以求解各速度分量和压力。对于湍流模型详细介绍了应用广泛的k-湍流模型，这也是本文所用的湍流模型。其后介绍了潜体绕流中应用到的几种边界条件，包括：入口边界条件、出口边界条件、固壁边界条件压力边界条件等，对使用边界条件时的注意事项作了说明。对于数值计算方法介绍了CFD的网格生成，方程离散所用到的有限体积法的基本思想以及控制方程组离散后的求解所用到的SIMPLE算法。最后对于自由面问题，介绍了运动界面追踪问题的VOF法。不同潜深的近水

38、面无附体小型航行器的数值模拟引言当潜体在近水面航行时，潜体位于水面下航行，必然遭受到水的反作用力，产生一种与潜体运动方向相反的流体作用力，也就是潜体阻力。同时，由于潜体在近水面航行，潜体与自由面之间相互作用，自由面将产生较大变形，甚至在自由面上有波破碎等现象出现。传统的方法是不考虑流体的粘性，以势流理论为基础数值求解拉普拉斯方程。这种方法很难揭示近水面潜体直航的运动规律，因为当潜体很接近自由面时，粘性起了很重要的作用，直接数值求解Ns方程是现行的方法。Yoon&Jung利用标高法模拟了二维正方形物体在自由面附近运动的流场，同时考察了阻力系数与潜深的关系。洪方文使用VOF法对方柱在自由面

39、附近运动的三维流场进行了数值模拟，考察了自由面与方柱在不同距离时的流场、不同长度方柱的流场以及它的阻力系数。在三维流场数值模拟中，回转体是最基本的潜体形状。本章以商用CFD软件FLUENT作为工具，利用VOF法结合k-湍流模型数值模拟近水面两端光滑过度的圆柱体直航的运动情形。计算模型为了建模简单，选择一个中间为圆柱体的简单回转体作为计算模型，圆柱体前端加上一个半球，末端采用光滑过度的椎体。圆柱体的直径为0.4m，模型的总长为3m。本人用SolidWorks建立了模型如图3.1所示。尺寸标注见图3.2。图3.1 本人用SolidWorks 建立的计算模型图3.2 SolidWorks智能尺寸标注

40、潜器在近水面运动时，潜器的运动与自由面是互相耦合的。为了使问题简化，研究方便，去除潜器的所有附体并假设它作一个自由度的运动即直航，且在运动过程中与静水面之间的距离保持不变，来流方向平行于潜器（回转体）的轴线。对于简单回转体而言，由于其几何形状规则，可以直接用FLUNT前处理软件GABIT来完成几何建模。GAMBIT是目前被广泛应用的计算流体力学(CFD)的前处理器。其建立几何模型的过程遵循点、线、面、体的顺序。首先设定回转体的各个控制点，然后将点连成线，画出圆弧，最后只需保留一个封闭的半轮廓线，将此半轮廓线形成面再绕一个坐标轴（本人是x轴）旋转360度即形成一个回转体，为了研究潜器模型与自由面

41、的相互作用，选择潜器潜深H分别为H/D=l，H/D=2，H/D=3三种情况对潜器周围绕流场进行数值模拟。其中潜深H为静水面到潜器上沿的垂向距离。潜深需在定义入流边界条件时控制边界长度，由自由面函数定义。对于深水潜器，可以假定为无限水深。在划分网格时，首先划分旋转面的面网格，采用Tri这个类型里的默认类型pave。然后划分网格疏密控制线上的网格，划分成一定比例，靠近潜器的部分密一些。最后划分体网格，采用默认的TGRID网格。根据计算要求的精确程度来定义网格大小。网格划分情况见图3.3和3.4。设置控制区域及划分网格完成几何建模之后，就要建立合适的控制区域。控制区域大小的选择将会影响到数值模拟的真

42、实性和计算结果的准确性，因此必须要根据模型的尺度及计算的要求建立恰当的控制区域。控制区域的选取原则是控制区域为其边界对潜器周围绕流场的影响可以被忽略的最小区域。由于潜器近水面直航，考虑到自由面的存在，控制区域被分为上、下两部分，上部分的流体为空气，下部分的流体为水。根据无界流场的概念和计算的需要，参考CFD粘流计算的文献，将控制区域选为长方体，其边界范围是-15m<x<15m,-5m<y<5m,-5m<z<5m。水面位置为y=2m。模型距离左边（x负方向）入流边界7m。这个控制区域是足够大的。由于潜器直航时流场左右对称，因此求解区域只需为控制区域的一半，这里

43、取Z0的部分。这样计算所需的网格数减少了，大大缩短了计算所需的时间，提高了计算效率。控制区域建立以后，与模型进行布尔并运算即可形成计算区域，对计算区域进行网格划分后需要定义边界条件。下面是网格的划分情况。图3.3 中间密的TGrid网格（四条蓝色的线为网格疏密控制线）图3.4 放大观察对称面附近的网格（可以看到许多漂亮的六面体网格）边界条件Fluent的计算是从边界条件开始的，边界条件的设置首先在GAMBIT软件中初步完成，再在FLUENT软件中对边界条件的参数进行具体设定。由于自由面的存在，在GAMBIT软件中需要将控制区域的入口和出口分别划分成空气入口、水入口和空气出口、水出口。根据计算模

44、型的坐标可知，来流方向为x轴正方向，将来流的入口设为速度入口(VELOCITY-INLET)，该边界条件专门用于不可压流动。一般假定在出口处来流未受到潜体的扰动影响，从而在出口处也可给定来流的速度分布。但实际上，在控制区域出口处的边界条件是未知的情形，故应该将控制区域的出口设定为自由出流(OUTFLOW)。但在大量计算练习时发现定义OUTFLOW边界条件时流体回流问题较为明显，即使网格细化后也难以解决此类问题，所以将出流边界重新设定为速度入流，速度设为负值，这样就相当于速度出流。总长仅仅3m的模型在30m长的距离内，对外面的流体影响已经可以近似为零，所以将出流边界设置为速度出流并非完全不可。控

45、制区域的对称面设定为对称边界(SYMMETRY)。鉴于出流没有设置为OUTFLOW，所以可以将充满空气的那部分体积设置一个101325pa的压强（压力边界和自由出流边界不可以同时存在）。由于计算控制区域的范围取得特别大，所以可以将控制区域的其它外边界包括潜器表面设定为无滑移的壁面(WALL)。计算结果表明这样设置是完全可以的。计算条件计算中模型潜深分别取HD=1、2、3，直航航速分别为v=1.5，2，2.5，2.7，3时的情况下进行计算。由于按照经验当Fr=0.5时阻力系数会出现峰值，故加算一个速度为2.7时的情况。以便观察此时的阻力变化情况。对于具有自由面的流动问题，计算是非稳态的，采用FL

46、UENT求解器求解流动控制RANS方程。选择分离式求解器Segregated Solver，采用隐式方案Implicit Formulation进行控制方程的线性化，选择非定常Unsteady Time状况。用VOF模型，自由面重构格式采用几何重构Geo-Reconstruct。设置空气为第一相，水为第二相。其中空气为FLUENT求解器中默认的介质，无需再定义介质的材料参数；对于水需选择FLUENT求解器材料数据库中的液态水(water-liquid)，即温度为20度时的水，密度为9982kgm3，动力粘性系数为0001003kgm·S，也可以根据需要改变有关物理参数。考虑到网格数量

47、和质量以及计算时间等多方面原因，对于湍流模型选用工程上常用的标准k-模型Standard k-epsilon Model，对于近壁区域的处理采用壁面函数法Standard Wall Functions。在运行环境中设置参考压力值为101325Pa，并计及重力影响，设置工作流体密度为1225kgm3，这样设置工作流体密度为较轻相的密度可以排除在较轻相中建立静水压力的计算，从而改善动量平衡计算的精度。FLUENT求解器中对边界条件的设定进行修改和进一步设置参数时，主要就是设置入口的不同来流速度和相关的湍流参数，即湍动能k和耗散率的初始值。对于k-模型，其计算与湍流参数即湍动能k和耗散率的初始值无关

48、，但是湍动能k和耗散率对计算时间和计算结果的收敛性是有影响的，因此k和初始值的合理设定将使计算更好的收敛。根据相关参考书的规定，采用如下计算公式：其中0.09是湍流模型中指定的经验常数（近似为0.09）。为最大速度，I为湍流强度，L为相当尺寸。本文中壁面条件将壁面的粗糙度全部设置为默认值0.5。在solution controls中分别采用了PISO和SIMPLE两种算法，结果相差不大。其中的Pressure相选用Body Force Weighted。Report相中的area填写潜器的湿面积。在求解计算的过程中，可以通过检查变量的残差、统计值、力的收敛趋势等，随时动态地监视计算的收敛性和当

49、前计算结果。计算中时间步长取0.002s，每个时间步内迭代计算次数设定为20次，运算800次。结果显示阻力系数是接近收敛的。由于本人所掌握的知识有限，计算机的配置偏低，不同计算条件下Gambit网格的划分情况并不完全相同等原因，对于这样计算的正确性和精确度不能予以绝对保证。但与其他参考文献中的结果相近。三维运算结果十分符合常理而且与参考文献相当。可见FLUENT软件的计算结果非常好。计算结果本计算主要考虑到近水面小型航行器的兴波问题和阻力系数问题。当潜器距离自由面较近时，自由面对潜器流场的影响是不能忽略的，并且这时自由面本身也发生变形，不再为静水平面，这样它对流场的影响就会更大。当自由面上发生

50、波破碎时，这种影响还会增加。自由面对流场的影响，势必会对潜器的受力造成影响。鉴于三维网格难以细化到本文所要求的水平及计算机硬件条件的限制，又为了非常清楚地体现兴波情况，本文中兴波图形均采用二维网格计算。立体云图及阻力系数等均采用三维网格精确计算。兴波图中的网格特征尺寸大小为0.004-0.02m，模型周围的网格尤为细密。图中绿色线条为波面，上面为空气，下面为水。由图中那条细细的绿色线条可以看出网格划分的非常细密。图3.53.10分别是v=1，1.5，2，2.5，2.7，3时潜深分别为H/D=1、2、3时的潜器周围绕流场所引起的自由面变形后的相图。图3.5 v=1时不同潜深的潜器引起的兴波情况图

51、3.6 v=1.5时不同潜深的潜器引起的兴波情况图3.7 v=2时不同潜深的潜器引起的兴波情况图3.8 v=2.5时不同潜深的潜器引起的兴波情况图3.9 v=2.7时不同潜深的潜器引起的兴波情况图3.10 v=3时不同潜深的潜器引起的兴波情况从这些图中可以看出不同情况下潜器周围绕流场所引起的自由面变形有很大区别。在同样的速度情况下，随着潜器潜深的增加，潜器周围绕流场对自由面的影响越来越小，自由面所受到的扰动会逐渐减小。当H/D=3时，自由面的变化进一步减小，按照经验可知此时可以认为潜器周围绕流场对自由面的影响是可以忽略的。对比上述几张图片可知，当速度较小时，自由面基本维持水平状态，随着速度的升

52、高，自由面将受到很大的影响而形成表面波，并且所激发的表面波随着速度的增加变得更加明显。当潜器周围绕流场对自由面产生干扰时，自由面不再为水平面，而是在潜器尾部附近形成一个凹坑。由于粘性的作用，潜器尾部的压力值比潜器前方相对应点的压力值低，这一方面使潜器产生形状阻力，另一方面使潜器尾部上方的自由面下陷，形成凹坑，激发自由面产生表面波。这种影响只产生于潜器附近的自由面，随着表面波向下游传播，其在自由面上扩展开来，幅度会越来越小。可见潜器周围绕流场对于潜器附近的自由面扰动最大，随着离潜器的距离越来越大，其绕流场对自由面的扰动会越来越小。图3.11-3.16列举了H/D=2时控制区域对称面的动压力云图。

53、图3.17-3.22列举了H/D=2时控制区域对称面的速度云图。从这些图中可以看出：来流速度越大，潜器周围绕流场的影响范围就越大。图3.11 v=1时的动压力云图 图3.12 v=1.5时的动压力云图图3.13 v=2时的动压力云图 图3.14 v=2.5时的动压力云图图3.15 v=2.7时的动压力云图 图3.16 v=3时的动压力云图图3.17 v=1时的速度云图 图3.18 v=1.5时的速度云图图3.19 v=2时的速度云图 图3.20 v=2.5时的速度云图图3.21 v=2.7时的速度云图 图3.22 v=3时的速度云图 图3.23 V=1.5时水平面上的动压力和总压分布云图图3.

54、23是V=1.5时水平面上的动压力和总压分布云图，由图中可以看到截面突变处压力变化比较明显，远离潜器的控制区域压力变化和缓。图3.24是立体图。可以看出立体图与水平和竖直平面图对应的非常好。 图3.24 V=1.5时潜器表面的的动压力和总压分布云图 图3.25 有时潜器尾部会产生漩涡脱落（unsteady模型计算所得）在FLUENT软件的后处理中可以得到如上的很多参数图，而且从图中可以得到每一点的参数值，此处不一一列举。下面主要分析潜器的阻力情况。图3.26为不同潜深潜器的总阻力系数曲线。从图中可看出阻力系数总体呈下降趋势且随着深度的增加阻力系数有所下降。潜深比较浅时（潜深0.4-0.8），阻力系数在速度为0.5时会增加，超过0.7时又有所下降。可能是因为有兴波阻力的存在。不过从总的变化趋势来看，兴波阻力的趋势（在Fr=0.5（本模型速度为2.7）前很小，当Fr接近0.5时急剧增加然后很快下降.）并不明显，可以由此推测虽然潜深比较浅时有兴波存在，但兴波阻力所占的成分并不是主要的。当潜深为1.2（三倍于模型直径）时，虽然有兴波存在，但兴波阻力的成分就更小了。