流体计算理论基础

流体计算理论基础

1 三大基本方程

1.1 连续性方程

连续性方程也称质量守恒方程，任何流动问题都必须满足质量守恒定律，该定律可表示为：单位时间内流体微元中质量的增加等于同一时间间隔内流入该微元体的净质量，其形式如下：

∂ρ∂(ρu ) ∂(ρv ) ∂(ρw ) +++=0 ∂t ∂x ∂y ∂z

可以写成：

∂ρ +div (ρu ) =0 ∂t

其中ρ密度，t 为时间，u 为速度矢量，u ，v 和w 为速度矢量在x ，y 和z 方向上的分量。 若流体不可压缩，密度为常数，于是：

∂u ∂v ∂w ++=0 ∂x ∂y ∂z

若流体处于稳态，则密度不随时间变化，可得出：

∂(ρu ) ∂(ρv ) ∂(ρw ) ++=0 ∂x ∂y ∂z

1.2 动量守恒定律

该定律可以表述为：微元体中流体的动量对时间的变化率等于外界作用在该微元体上的各种力之和，该定律实际是牛顿第二定律，按照这一定律，可导出x ，y 和z 三个方向上的动量守恒方程：

⎧∂(ρu ) ∂p ∂τxx ∂τyx ∂τzx +div (ρuu ) =-++++F x ⎪∂x ∂x ∂y ∂z ⎪∂t

⎪∂p ∂τxy ∂τyy ∂τzy ⎪∂(ρu ) +div (ρuv ) =-++++F y ⎨∂t ∂y ∂x ∂y ∂z ⎪⎪∂(ρu ) ∂p ∂τxz ∂τyz ∂τzz +div (ρuw ) =-++++F z ⎪∂t ∂z ∂x ∂y ∂z ⎪⎩

式中，p 为微元体上的压力，τxx ，τxy 和τxz 等是因分子粘性作用而产生的作用在微元体表

面上的粘性应力τ的分量。F x ，F y 和F z 是微元体上的体力，若体力只有重力，且z 轴竖直向上，则：F x =0, F y =0，F z =-ρg 。

对于牛顿流体，粘性应力τ与流体的变形率成比率，有：

∂u ∂u ∂v ⎧τ=2μ+λdiv (u ); τ=τ=μ(+) xy y x ⎪xx ∂x ∂y ∂x ⎪∂v ∂u ∂w ⎪⎨τyy =2μ+λdiv (u ); τxz =τz x =μ(+) ∂x ∂z ∂x ⎪∂w ∂v ∂w ⎪τ=2μ+λdiv (u ); τ=τ=μ(+) yz zy ⎪zz ∂x ∂z ∂y ⎩

其中，μ为动力粘度，λ为第二粘度，一般可取λ=-2，将上式代入前式中为： 3

⎧∂(ρu ) ∂p +div (ρuu ) =div (μgradu ) -+S u ⎪∂t ∂x ⎪∂p ⎪∂(ρv ) +div (ρuv ) =div (μgradv ) -+S v ⎨∂y ⎪∂t

⎪∂(ρw ) ∂p +div (ρuw ) =div (μgradw ) -+S w ⎪∂z ⎩∂t

其中：

grad ()=∂()/∂x +∂()/∂y +∂()/∂z

2μ为动力粘度(dynamic viscosity) ，λ为第二粘度(second viscosity), 一般可取：λ=-(参考3

文献：H.Schlichting,Boundary Layer Theory,8th ed,McGraw Hill, New York,1979)。S u , S v 和S w 为动量守恒方程中的广义源项，S u =F x +S x ，S v =F y +S y ，S w =F z +S z ，而其中S x ，S y 和S z 表达式为：

⎧∂∂u ∂∂v ∂∂w ∂ S =(μ) +(μ) +(μ) +(λdiv (u )) ⎪x ∂x ∂x ∂y ∂x ∂x ∂x ∂x ⎪⎪∂∂u ∂∂v ∂∂w ∂ S =(μ) +(μ) +(μ) +(λdiv (u )) ⎨y ∂x ∂x ∂y ∂y ∂x ∂y ∂y ⎪⎪∂∂u ∂∂v ∂∂w ∂ S =(μ) +(μ) +(μ) +(λdiv (u )) ⎪z ∂x ∂z ∂y ∂z ∂x ∂z ∂z ⎩

一般来讲，F x ，F y 和F z 是体积力在x ，y ，z 方向上的分量。S x ，S y 和S z 是小量，对于粘性为常数的不可压缩流体，S x =S y =S z =0，动量守恒，简称动量方程，也称N-S 方程。 关于牛顿体与非牛顿体的定义如下：

流体的内摩擦剪切力τ由牛顿内摩擦定律决定：

τ=μlim ∆u ∂u =μ ∆u →0∆n ∂n

∆u 为法向距离上∆n 其中，∆n 为沿法线方向的距离增量，∆u 对应于∆n 的流体速度的增量，

的速度变化率，所以，牛顿内摩擦定律表示：流体的内摩擦应力和单位距离上的两层流体间的相对速度成比例，比例系数μ称为流体的动力粘度，常称为粘度，单位为：N ⋅s /m 2 若μ为常数，则该类流体为牛顿流体，否则为非牛顿体，空气，水等均为牛顿体；聚合物溶液，含有悬浮粒杂质或纤维的流体为非牛顿体。

对于牛顿流体，通常用μ和[质量]密度ρ的比值ν代替动力粘度μ

ν=μ ρ

ν称为运动粘度，单位m 2/s 。

1.3 能量守恒方程

该方程可以描述为：微元体中能量的增加率等于进入微元体的净热流量加上体力与面力对微元体所做的功，实际为热力学第一定律。

∂(ρT ) k +div (ρuT ) =div (gradT ) +s T ∂t c p

k 为流体传热系数，c p 为比热容，T 为温度，S T 为流体内热源及由于粘性作用流体机械能转换为热能的部分，有时简称S T 为粘性耗散项。

以上三大基本方程参考：

《计算流体动力学分析：CFD 软件原理与应用》_王福军

2 通用控制方程

上面的基本方程可以写成下面的通用形式：

∂(ρφ) +div (ρu φ) =div (Γgrad φ) +S ∂t

展开为：

∂(ρφ) ∂(ρu φ) ∂(ρv φ) ∂(ρw φ) +++∂t ∂x ∂y ∂z ∂∂φ∂∂φ∂∂φ=(Γ) +(Γ) +(Γ) +S ∂x ∂x ∂y ∂y ∂z ∂z

其中φ为通用变量，可以代表u ，v 和w 以及T 等求解变量，Γ为广义扩散系数，S 为广义源项。

2 几种数值求解方法

2.1 有限差分法

主要的思路是用差商代替微商, 来近似的表示微分方程. 其形式简单, 对任意复杂的偏微分方程都可以写成其对应的差分方程, 但是微分方程中各项的物理意义和微分方程所反映的物理定律(如守恒定律) 在差分方程中所表现的特点, 在差分方程中没有得到体现. 只是微分方程的数学近似, 没有反映物理特性, 计算结果可能表现出某些不合理现象.

2.2 有限元法

20世纪60年代出现, 离散方程获得的方法主要有:直接刚度法, 虚功原理推导, 泛函原理推导或加权余量法推导.

有限元法的优点是解题能力强, 可以较精确的模拟各种复杂的曲线或曲面边界, 网格划分比较随意, 可以统一处理多种边界条件, 离散方程形式规范, 便于编写通用程序, 但在应用流体流动和传热原理中却遇到了一些困难, 其原因可归结为按加权余量法推导出的有限元离散方程也只是对原微分方程的数学近似, 当处理流动和传热问题的守恒性, 强对流, 不可压缩条件等方面的要求时, 有限元离散方程中的各项还无法给出合理的物理解释, 对计算中出现的一些误差也难以进行改进, 因此有限元法在流体力学和传热学中的应用还存在一些问题.

2.3 有限体积法FVM(或控制体积法CVM)

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

有限体积法得到的离散方程，要求因变量的积分守恒对于任意的一组控制体积都满足，对于整个计算区域，自然满足，这个是有限体积法吸引人的优点。有些离散方法，例如有限差分，仅当网格极其细密时，离散方程才满足积分守恒，而有限体积法在粗网格的情况下，也

显示出准确的积分守恒。

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

插值方式称之为离散格式，有：中心差分，一阶迎风格式，混合格式，指数格式，乘法格式。

二阶迎风格式，QUICK 格式，

有限体积法有四项基本原则：

（1）控制体积界面上的连续性原则

（2）正系数原则

（3）源项的负斜率线性化原则

（4）系数αp 等于相邻节点系数之和原则

2.4 谱方法

基本思想是考虑热传导方程的初值问题, 当微分方程的解足够光滑时, 谱方法给出的近似解将以很高的精度逼近微分方程的精确解, 而且该方法得到的近似解应用于整个区域而不是局部区域. 这个是区别有限元法的重要特征.

2.5 边界元法

20世纪70年代发展起来, 针对有限差分和有限元占计算机内存资源过多而发展起来的求解偏微分方程的数值方法. 最大特点是降维, 只在求解区域的边界进行离散就能求解整个流场的解. 这样三维问题降维二维, 二维问题降为一维, 可用小机器计算大问题. 基本思想是用边界积分方程将求解域的边界条件与域内任意一点的待求变量值联系起来, 然后求解边界积分方程即可. 但是若流体描述方程本身比较复杂时, 如粘性的N-S 方程, 则对应的权函数算子基本解

不一定能找到, 因此应用受到很大限制.

3 离散控制方程求解概述

建立了离散方程之后, 所生成的离散方程不能直接用来求解, 还必须对离散方程进行某种调整, 并且对各未知量(速度，压力，温度等) 的求解顺序及方式进行特殊处理，对于这些计算方法有：

耦合式解法的基本过程：

（1）假定初始压力和速度等变量，确定离散方程的系数及常数项。

（2）联立求解连续方程，动量方程，能量方程

（3）求解湍流方程及其它标量方程

（4）判断当前时间步上的计算是否收敛，若不收敛，返回到第（2）步，迭代计算，若收敛，重复上述步骤，计算下一时间步的物理量。

若所有变量整场联立求解，称为隐式解法，部分变量整场联立求解称为显隐式解法，在局部地区（如一个单元上）对所有变量联立求解称为显式解法。对于显式解法，是在一个单元上求解所有变量后，逐一的在其它单元上求解所有的未知量，这种方法在求解某个单元时，要求相邻单元的变量解是已知的。

分离式解法的基本思路是：不直接联立方程组，而是顺序地，逐个的求解各变量代数方程组。依据是否直接求解原始变量u,v,w 和p ，分离式解法可分为原始变量法和非原始变量法。 涡量—流函数法不直接求解原始变量u,v,w 和p ，而是求解旋度ω和流函数ψ。

涡量—速度法不直接求解原始变量p ，而是求解旋度ω和u,v,w ，此两种方法共同特点是不求解压力项，避免压力项带来的问题。

缺点是：不易扩展到三维情况，当存在压力时，需要单独求解压力，对于固定壁面，其上旋度极难确定，往往使得涡量方程的数值解发散或者不合理。

原始变量法中，解压力泊松方程法需要采用对方程取散度等方法将动量方程转变为泊松方程，然后对泊松方程进行求解，与这种方法对应的是著名的MAC 方法和分布法。

人为压缩法主要是受可压的气体可以通过联立求解速度分量与密度的方法来求解的启发，引入人为压缩性和人为状态方程，以此对不可压流体的连续方程施加干扰，将连续方程写为包含有人为密度的项，而人为密度前有一个极小的系数，这样，方程可以转化为求解人为密度的基本方程，但是这种方法要求的时间步长必须很小。因此限制了它的广泛应用。

目前工程上最为广泛的流场数值解法是压力修正法，实质是迭代法，在每一时间步长中，

先给出压力场的初始猜测值，据此求出猜测的速度场，再求解根据连续方程导出的压力修正方程，对猜测的压力场和速度场进行修正，如此循环往复，可得出压力场和速度场的收敛解，其基本思路是：

（1）假定初始压力场

（2）利用压力场求解动量方程，得到速度场

（3）利用速度场求解连续方程，使压力场得到修正

（4）如果需要，求解湍流方程及其它标量方程

（5）判断当前时间步上的计算是否收敛，若不收敛，返回到第（2）步，迭代计算，若收敛，重复上述步骤，计算下一时间步的物理量。

压力修正法有很多方式，其中压力耦合方程组的半隐式方法(simple算法) 应用最为广泛，其首先使用一个猜测的压力场来求解动量方程，得到速度场，接着求解通过连续方程所建立的压力修正方程，得到压力场的修正值，然后利用压力修正值更新速度场和压力场，最后检查结果是否收敛，若不收敛，以得到的压力场作为新的猜测的压力场，重复该过程，为了启动该迭代过程，需要提供初始的，带有猜测性的压力场和速度场，随着迭代的进行，这些猜测的压力场和速度场不断改善，得到的压力与速度分量值逐渐逼近真解。SIMPLE 算法及其改进是算法一般都依赖于交错网格。

具体有：

SIMPLE,SIMPLER,SIMPLEC,PISO 等算法。

3.1 SIMPLE

全称：Semi-Implict Method for Pressure-Linked Equations的缩写，意为“求解压力耦合方程组的半隐式方法”于1972年提出，核心是采用“猜测—修正过程”，在交错网格的基础上计算压力场，从而达到求解动量方程的目的。

基本思路如下：

3.2 SIMPLER

是SIMPLE Revised的缩写，是SIMPLE 算法的改进。

基本过程：

3.3 SIMPLEC

是英文SIMPLE Consistent的缩写，意为协调一致的SIMPLE 算法。

3.4 PISO

全称：Pressue Implicit with Splitting of Operators的缩写，意为压力的隐式算子分割算法，1986年提出。

5 各种方程的在非结构网格离散形式

由于fluent 在5.5版本后，就采用非结构网格，故介绍一下非结构网格。 通用方程如下：

∂(ρφ)

+div (ρu φ) =div (Γgrad φ) +S ∂t

对任意控制体积的积分有：

∂(ρφ)

+div (ρu φ) dV =⎰⎰⎰div (Γgrad φ) dV +⎰⎰⎰SdV ⎰⎰⎰∆v ∂t ⎰⎰⎰∆v ∆v ∆v

依据散度定理

∂(ρφ) ∂φ+ρu v dS =Γ⎰⎰⎰∆v ∂t ⎰⎰∆s i i ⎰⎰∆s ∂x i v i dS +⎰⎰⎰∆v SdV

其中，∆v 为控制体的体积，∆s 为控制体的表面积（二维中为多边形的边长），x i 表坐标方向，x 1=x , x 2=y ，而v i 表示控制体积各边的单位法向量，v 1=v x , v 2=v y , u i 表示速度

矢量，u 1=u , u 2=v 。 具体各项形式如下 瞬态项：

(ρφ) p -(ρφ) 0p ∂(ρφ)

∆V ⎰⎰⎰∆v ∂t =∆t

φp 是变量φ在控制体积中心点P 的值，上标表示前个时间步的值，∆t 是时间步长。

源项

⎰⎰⎰

∆v

SdV =S ∆V =(S C +S P φp ) ∆V =S C ∆V +S P φp ∆V

其中，

S C 是常数，S P 是随时间和物理量φ变化的项。

扩散项

，

N S ⎧⎫∂φφ-φ⎪

Γv dS =[Γ(v ∆y -v ∆x )]⎬+C diff x y ⎰⎰∆s ∂x i i ∑E =0⎪E ⎭其中N S 为控制体积P 的总面数，也是相邻控制体积的数量，E 表示与控制体积有公共界面的控制体积，v x 和v y 表示控制体积各界面的单位法向矢量的分量，∆x 和∆y 表示界面的外法线矢量的分量，δx 和δy 是两个控制体积之间节点P 到节点E 的矢量分量，C diff 为公共界面上的变叉扩散项。当矢量N 和界面e 垂直时, 通过该界面的变叉扩散量等于0, 对于一般的准正交网格, C diff 是小量, 可按0处理, 若网格高度奇异, 则C diff 不可忽略, 但是目前还没有办法准确计算C diff 这一项, 因此, 为了避免计算, 构建网格的时候尽量选择正交网格. 对流项

⎰⎰

∆s

ρu i v i dS =∑[ρφ(u ∆y -v ∆x )]]

E =1

E

N i

注意, 上式中界面处的φ值要通过插值公式(空间离散格式) 计算.

将上式综合在一起, 再在世界域上引入全隐式时间积分方案, 得到其在非结构网格上的离散

格式.

a p φp =∑a E φE +b p

E

N s

其中,

a p =∑a E +

E

N s

(ρp ∆V ) 0

∆t

-S p ∆V

b p =

(ρp φp ∆V )

∆t

+S c ∆V

系数a E 取决于对流项所使用的离散格式, 例如, 若对流项使用一阶迎风格式

a E =D e +max(0,-F e )

对于中心差分格式有:

a E =D e -

F e

2

符号e 表示控制体积P 与E 相邻的界面, F e 和D e 分别是界面e 上的对流质量流量与扩散传导性，计算公式如下：

F e =ρu S

S N

D e =Γφ 2

N

各个量见下图：

6 各种离散方程组的解法

无论采用何种离散格式，无论采用什么算法，最终都要生成离散方程组，除非对瞬态问题采用显式解法，都需要求解离散方程组。 解法分为直接解法和迭代法两类 6.1 直接解法

6.1.1 Cramer矩阵求逆法

其只适用于方程组规模非常小的情况。 6.1.2 Gauss消去法

Gauss 消去法先把系数矩阵通过消元而化为上三角阵，然后逐一回代，从而得到方程组的解。比Cramer 矩阵求逆法能够适应较大规模方程组，但是不如迭代法效率高。

6.2 迭代法

Jacobi 迭代法和Gauss-Seidel 迭代法。 6.3 TDMA解法

Tomas 在较早以前开发了一种能快速求解三对角方程组的解法TDMA(Tri-Diagonal Matrix Algorithm) ，目前得到了广泛的应用，对于一维的CFD 问题，其实际是一种直接解法，但是它可以迭代使用，从而用于二维和三维问题中的非三对角方程组。

7 湍流

7.1 湍流的定义

当流速很小时，流体分层流动，互不混合，相邻的流体层彼此有序地流动，称为层流(laminar FLOW) ；逐渐增加流速，流体的流线开始出现波浪状的摆动，摆动的频率及振幅随流速的增加而增加，此种流况称为过渡流；当流速增加到很大时，流线不再清楚可辨，流场中有许多小漩涡，层流被破坏，相邻流层间不但有滑动，还有混合。这时的流体作不规则运动，有垂直于流管轴线方向的分速度产生，这种运动称为湍流(turbulent flow)，又称为乱流、扰流或紊流。

观测表明湍流带有旋转流动结构这就是湍流涡turbulent eddies简称涡eddy 。从物理结构上看，可以把湍流看成是由各种不同尺寸的涡叠合而成的流动，这些涡的大小和旋转轴的方向分布是随机的。大尺度的涡主要是由流动的边界条件所决定，其尺寸可以与流场的大小相比拟，它主要受惯性影响而存在是引起低频脉动的原因。小尺度的涡主要是由粘性力所决定的其尺寸可能只是流场尺度的千分之一量级，是引起高频脉动的原因。大尺寸的涡不断地从主流中获得能量通过涡间相互作用能量逐渐向小尺寸的涡传递。最后由于流体粘性的作用，小尺度的涡就不断消失，机械能就耗散为流体的热能。同时由于边界的作用扰动及速度梯度的作用，新的涡又不断产生，构成了湍流运动。对某些简单的均匀时均流场，如果湍流脉动是均匀的、各向同性的，可以用经典的统计理论进行分析。但实际上湍流是不均匀的。由于湍流的存在速度脉动量，在流线方向的分量和垂直于流线方向的分量之间建立了关联量，它代表着一种横向交换通量，也可以认为是由于湍流流动引起的一种附加剪切应力——影响动量的输运过程。湍流的存在使传热和传质通量提高。由于湍流会促进这些基本过程，因此对某些物理现象就会产生强烈的影响，如脉动过程的消衰、均相化学反应率的增加以及液滴蒸发的强化。某些因素会影响湍流的形成。如当湍流定性尺度和脉动强度非常小时流体的粘度会直接影响当地的湍流度。当马赫Mach 数达到5以上时密度的脉动量与当地的湍流有密切的关系。强烈的化学反应、气流的旋转流动、颗粒的存在以及浮力或电磁场的作用都会影响当地的湍流结构。

一般认为，无论湍流多么复杂，非稳态的连续方程和N-S 方程，对于湍流的瞬时运动仍然是适用的。

7.2 湍流的数值解法

目前湍流的数值模拟方法分为：直接数值模拟方法和非直接数值模拟方法。所谓直接数值模拟方法就是直接求解瞬时湍流流动方程，而非直接数值模拟方法就是不直接计算湍流的脉动特性，而是设法对湍流作某种程度的近似和简化处理。

7.2.1 直接数值模拟(DNS)

直接数值模拟(Direct Numerical Simulation) 方法就是直接用瞬时的N-S 方程对湍流进行计算。DNS 最大的好处就是无需对湍流流动作任何简化或近似，理论上可以得到相对准确的计算结果。

但是，实验测试表明，在一个0.1×0.1m 2大小的流动区域内，在高Reynolds 数的湍流中包含尺度为10um~100um的涡，要描述所有尺度的涡，则计算的网格节点数将高达109~1012. 同时，湍流脉动的频率约为10KHZ, 因此，必须将时间的离散步长取为100us 以下。对计算机能力提出了很大的挑战。因此，目前还无法用于真正意义上的工程计算。 7.2.2 大涡模拟(LES)

放弃对全尺度范围上涡的运动模拟，而只将比网格尺度大的湍流运动通过N-S 方程直接计算出来，对于小尺度的涡对大尺度运动的影响则通过建立模型来模拟，从而形成了目前的大涡模拟法(Large eddy simulation,简称LES) 。

LES 方法基本思想可概括为：用瞬时的N-S 方程直接模拟湍流中的大尺度涡，不直接模拟小尺度涡，而小涡对大涡的影响通过近似的模型来考虑。

总体而言，LES 方法对计算机内存及CPU 要求仍然比较高，但是低于DNS 。 7.2.3 Reynolds平均法(RANS)

Reynolds 平均法的核心是不直接求解瞬时的N-S 方程，而是想办法求解时均化的Reynolds 方程，这样，不仅可以避免DNS 方法的计算量大的问题，而且对工程实际应用可以取得很

好的效果。根据对Reynolds 应力作出的假定或处理方式不同，目前常用的湍流模型有两大类：Reynolds 应力模型和涡粘模型。 Reynolds 时均格式的N-S 方程（RANS ）：

标量 的时均输送方程：

7.2.3.1 Reynolds应力模型

直接构建表示Reynold 应力的方程，然后联立求解，通常情况下，Reynolds 应力方程是微分形式，称为Reynolds 应力方程模型，若将Reynolds 应力方程的微分形式简化为代数方程的形式，则称这种模型为代数应力方程模型。 7.2.3.1.1Reynolds 应力方程模型 Reynold 应力输送方程如下：

第一项为瞬态项，其它为：

然后对各项进行计算，综合后，可采用SIMPLE 算法进行计算。

Reynold 应力方程模型是高Re 数的湍流模型, 可是, 在近壁区内的流动,Re 数较低, 湍流发展并不充分, 湍流的脉动影响不如分子粘性的影响大, 湍流应力几乎不起作用, 这样在这个区域就不能使用k -ε相关的模型, 必须进行特殊处理. 如壁面函数或低Re 数模型。

同时，计算实践表明，RSM 虽能考虑一些各向异性效应，但是并比一定比其它模型效果更好，在计算突扩流动分离区和计算湍流输送各向异性较强的流动时，RSM 优于双方程模型，但是对于一般的回流流动，其不一定比k -ε模型要好。同时，就三维计算而言，采用RSM 意味要多求解6个Reynolds 应力的微分方程，计算量大，对计算机要求高。因此，RSM 不如k -ε应用更广，但是RSM 是一种更加有潜力的湍流模型。

7.2.3.1.2代数应力方程模型

由于RSM 过于复杂，计算量大，有许多学者从RSM 出发，建立Reynolds 应力的代数方程模型，即将RSM 中包含Reynolds 应力的微商的项用不包含微商的表达式去代替，就形成了代数应力方程模型（Algebraic Stress equation Model简称ASM ），简化时，重点集中在对流项和扩散项的处理上。

一种简化方案是采用局部平衡假定，即Reynolds 应力的对流项和扩散项之差为0；另外一种方案是假定Reynolds 应力的对流项和扩散项之差正比与湍动能k 的对流项和扩散项之差。 ASM 是将各向异性的影响合并到Reynolds 应力中进行计算的一种经济算法，但是，因要解6个代数方程组，计算量还是远大于k -ε模型。

ASM 虽然不像k -ε广泛，但是可用于k -ε不能满足的场合以及不同的传输假定对计算精度影响不是十分明显的场合。例如：对于像方形管道和三角形管道内的扭曲和二次流的模拟，由于流动特征是由Reynolds 正应力的各向异性造成的，因此使用标准的k -ε模型得不到理想的结果，而使用ASM 就非常有效。

7.2.3.2 涡粘模型

不直接处理Reynolds 应力项，而是引入湍动粘度（turbulent viscosity），或称涡粘系数（eddy viscosity ），然后把湍流应力表示成湍动粘度的函数，整个计算的关键在于确定这种湍动粘度。 湍动粘度的提出来源于Boussinesq 提出的涡粘假定，把因湍流引起的，由脉动速度相关联的剪切应力τ，模仿层流由以时间平均速度的梯度来表达，即：

μt 为湍流粘度，u i 为时均速度，δij 是”Kronecker delta”符号（当i=j时，δij =1；当i ≠j 时，

，k 为湍流动能(turbulent kinetic energy ) δij =0）

1

k =(u '2+v '2+w '2)

2

湍动粘度μt 是空间坐标的函数，取决于流动状态，而不是物性参数。引入Boussinesq 假定后，计算湍流流动的关键在于如何确定μt 。

所谓的涡粘模型，就是把μt 与湍流时均参数联系起来的关系式，依据确定μt 的微分方程的数目多少，涡粘模型包括：零方程模型，一方程模型，两方程模型。目前两方程模型应用最广，最基本的两方程模型是标准的k -ε模型。还有各种改进的k -ε模型，如：RNG k -ε模型和Realizable k -ε模型。

7.2.3.2.1 零方程模型

指不使用微分方程，而是用代数关系式，把湍动粘度与时均值联系起来的模型。 零方程模型方案有很多种，最著名的是：Prandtl 提出的混合长度模型（mixing length model）。Prandtl 假定湍动粘度正比于时均速度的梯度和混合长度的乘积。但是只有在简单的流动中才比较容易给定混合长度，对复杂流动很难确定。

只能解释某些简单的流动过程，其在实际工程中很少使用。 7.2.3.2.2 一方程模型

为了弥补混合长度假定的局限性，人们在时均连续方程和Reynolds 方程基础上，再建立一个湍动能的输送方程。

由Kolmogorov-Prandtl 表达式，有：

μt =ρC μ

一方程模型中如何确定长度比尺仍为不易解决的问题，因此很难得到推广。 7.2.3.2.3 标准k -ε两方程模型

在一方程基础上，新引入一个关于湍流耗散率ε的方程形成。 湍流耗散率定义为：

在标准的k -ε模型中，k 和ε是两个基本未知量，与之对应的输送方程为：

其中，G k 是由于平均速度梯度引起的湍动能k 的产生项，G b 是由于浮力引起的湍动能产生项，Y M 代表可压湍流中脉动扩张的贡献，C 1ε, C 2εC 3ε为经验常数，σk k 和耗散率ε对应的Prandtl 数，S k 和S ε是用户定义的源项。 各个参数计算如下：

σε分别是与湍动能

（对于不可压缩流体，G b =0）

Pr t 是湍动Prandtl 数，可取Pr t =0.85，g i 是重力加速度在第i 方向的分量，β是热膨胀系

数，可由压流体的状态方程求出：

1∂ρβ=-*

ρ∂T

对于不可压缩体，Y M =0，对于可压缩体，

，其中，M t 是湍流Mach

数，

M t =

是声速，a =

若采用标准的k -ε模型求解流动及换热问题时，控制方程包括连续性方程，动量方程，能量方程，k 方程，ε方程与k -ε模型中对应的输送方程。 如果考虑传质或者化学变化情况，则应该加上组分方程，这些方程可以表示为下面的通用形式：

即：

∂(ρφ)

+div (ρu φ) =div (Γgrad φ) +S ∂t

关于k -ε模型，有如下的说明：

（1）模型中的有关系数，主要根据一些特殊条件下的实验结果而确定，在不同的文献讨论不同的问题时，这些值可能有出入。在数值计算的过程中，针对特定的问题，参考相关文献，寻求更合理的取值。

（2）上述k -ε模型，是针对湍流发展非常充分的湍流流动来建立的，是一种针对高Re 数的湍流计算模型，而当Re 数较低时，例如，在近壁区内的流动，湍流发展并不充分，湍流的脉动影响可能不如分子粘性的影响大，在更贴近壁面的底层内，流动可能处于层流。因此，对于Re 数较低的流动使用上面建立的k -ε模型，就会出现问题。这时，必须采用特殊的处理方式，以解决近壁区内的流动计算及低Re 数时的流动问题。使用上面的k -ε模型可能就会出现问题。常用解决方法有壁面函数法和低Re 数的k -ε模型。

（3）标准的k -ε比零方程模型和一方程模型有了很大的改进，但是对于强旋流，弯曲壁面流动或弯曲曲线流动时，会产生一定失真，原因在标准k -ε中，对于Reynolds 应力的各个分量，假定粘度系数μt 是相同的，即假定μt 是各向同性的标量，而在弯曲流线的情况下，湍流是各向异性的，μt 应该是各向异性的张量。

7.2.3.2.4 RNG k -ε模型

RNG k -ε模型由Yakhot 及Orzag 提出，RNG 全称：renormalization group缩写。 其基本思想是：

在谱空间内对N-S 方程引入了所谓的“对应原理”，利用Gauss 统计法在平衡态展开，经过一序列移去小尺度部分对余下部分重新标度的运算，得到一针对大尺度运动的方程。其中小尺度对大尺度的影响在方程中以涡粘性的方式体现。若移去的仅是那些最小的尺度就得到大涡模拟中的亚格子模型，若移去的尺度继续增大，最终得到就是涡粘性模型，如代数模式，两方程模式，非线性模式。在高雷洛数极限情况下，所得模式（称RNG 模式）与标准模式形式上完全一样，仅在系数上有所差别。

所得到的k -ε方程为：

其中：

RNG c 模型能够更好的处理高应变率及流线弯曲程度较大的流动。仍然是对充分发展的湍流是有效的。

此外，Fluent 将RNG 模型所引入的反映主流的时均应变率E ij 一项，ε归入了方程的系数C 2b 中，且表达式多了一个系数C μ，而不是归入C 1E 。

7.2.3.2.5Realizble k -ε模型

标准的k -ε模型对于时均应变率特别大的情况, 可能出现负的正压力. 为使流动符合湍流的物理规律, 需要对正应力进行某种数学约束, 为了保证这种约束的实现, 认为湍流粘度计算式中的系数不应是常数, 而应与应变率联系起来, 从而出现了Realizble k -ε模型.

其模型形式如下

:

其中

:

Realizble k -ε模型已被有效的用于各种不同类型的流动模拟, 包括旋转均匀剪切流, 包含有射流和混合流的自由流动, 管道内流动, 边界层流动以及带有分离的流动等.

7.2.3.3 关于近壁区的流动

以上的k -ε模型都是高Re 数的湍流模型, 可是, 在近壁区内的流动,Re 数较低, 湍流发展并不充分, 湍流的脉动影响不如分子粘性的影响大, 湍流应力几乎不起作用, 这样在这个区域就不能使用k -ε相关的模型, 必须进行特殊处理.

解决这个问题有两种途径:一个是采用壁面函数法, 一个是采用低Re 数的k -ε模型.

7.2.3.3.1 壁面函数法

基本思想是:对于湍流核心区的流动使用k -ε模型求解, 而在壁面区不进行求解, 直接使用半经验公式将壁面上的物理量与湍流核心区内的求解变量联系起来. 这样, 不需要对壁面区内的流动进行求解, 就可以直接得到与壁面相邻控制体积的节点变量值.

可以将壁面区分为三个子层:粘性底层, 过渡层, 对数律层.

壁面函数法是FLUENT 选用的默认方法，它对各种壁面流动都非常有效。相对低Re 数的k -ε模型，壁面函数法计算效率高，工程实用性强。当然壁面函数法无法像低Re 数k -ε模型那样得到粘性底层和过渡层内的”真实”速度分布。

壁面函数法还有一定的局限，当流动分离过大或近壁面的流动处于高压之下时，此方法不是很理想，为此，fluent 还提供了非平衡的壁面函数法及增强的壁面函数法。

7.2.3.3.2 低Re 数的k -ε模型

低Re 数的流动主要体现在粘性底层，流体的分子粘性起着绝对支配的地位，因此必须对高Re 数的c 模型进行三个方面的修改：

（1）为体现分子粘性影响，扩散方程的扩散系数项必须同时包含湍流扩散系数与分子扩散系数两部分。

（2）控制方程的有关系数必须考虑不同流态的影响，即系数计算公式中引入湍流雷洛数Re t ，Re t =ρk 2/(ηε) 。

（3）在k 方程中应该考虑壁面附近湍动能的耗散不是各向同性这一因素。

因此得到其输送方程为：

7.3 fluent中的粘度模型

7.3.1 Inviscid

进行无粘计算。

7.3.2Laminar 模型

用层流的模型进行流动模拟，不需要用户输入任何与计算模型有关的参数。

7.3.3Spalart-Allmaras(1 eqn)模型

这个是用于求解动力涡粘输送方程相对简单的一种模型，Spalart-Allmaras 模型是专门用于求解航空领域的壁面限制流动，对于受逆压力梯度作用的边界层流动，已取得很好的效果，在透平机械中的应用也越来越普遍。

原始的Spalart-Allmaras 模型实际是一种低雷洛数模型，要求在近壁区的网格划分很细，但是在fluent 中，由于引入了壁面函数法，这样，Spalart-Allmaras 模型用在较粗的壁面网格也可以。当精确的湍流计算并不是十分需要时，这种模型是最好的选择。

7.3.4 k-epsilon(2 eqn)模型

使用k -ε双方程模型进行湍流计算，就像上面所述，又分为标准的k -ε，RNG k -ε模型和Realizable k -ε模型三种。初次使用FLUENT 时，可暂时使用默认值。

这三种模型均是针对充分发展的湍流才有效，均是高Re 数的湍流模型。只能用于求解湍流核心区域的流动。

对于壁面的流动，fluent 采用了壁面函数法来弥补不足。

7.3.5 k-omega(2 eqn)模型

标 使用k-ω双方程模型进行湍流计算。k-ω双方程模型分为标准的k-ω和SST k-ω模型。

准的k-ω模拟基于Wilcox k-ω模型，在考虑低雷洛数，可压缩性和剪切流特性的基础上修改而成。Wilcox k-ω模型在预测自由剪切流传播速率时，取得了很好的效果，成功应用于尾迹流，混合层流动，平板绕流，圆柱绕流和放射状喷射。因而可以说该模型能够应用于壁面约束流动和自由剪切流动。SST k-ω模型全称是剪切应力输运(shear-stress transport) k-ω模型，是为了使标准k-ω模型在近壁面有更好的精度和算法稳定性而发展起来的，可以说是将k -ε模型转换到k-ω模型的结果，因此，sst k-ω模型在很多时候比标准的k-ω模型更加有效。

7.3.6 Reynolds Stress 模型

使用RSM 进行湍流计算，在fluent 中，Reynolds 应力模型是最精细制作的湍流模型，它放弃了各向同性的涡粘假定，直接求解Reynolds 应力方程。由于它比单方程和双方程模拟更加严格地考虑了流线弯曲，旋涡，旋转和张力快速变化，它对于复杂的流动总体上有更高的预测精度。但是，为了事Reynolds 方程封闭而引入了附加模型，也会使得这种方法的预测结果的真实性收到挑战。

总体来讲，Reynolds 应力模型的计算量很大，当要考虑Reynolds 应力的各向异性时，例如飓风流动，燃烧室高速旋转流，管道中二次流，必须采用Reynolds 应力模型。

7.3.7 Large Eddy Simulation 模型

使用大涡(LES)模拟，该模型只对三维问题有效，是目前比较有潜力的湍流模型。

在上述的几种湍流模型中，从计算的角度看，Spaimn-Allmaras 模型在fluent 中是最经济的湍流模型，标准的k -ε模型比它耗费更多的计算资源。Realizable k -ε模型比标准的k -ε模型需要多一点的计算资源，RSM 因为增加了Reynolds 应力方程而需要更多的内存和SPU 时间，而FLUENT 的有效程序设计，使得RSM 算法并没有在CPU 时间方面增加很多，在每个迭代法中，RSM 比k -ε模型要多耗50~60%的CPU 时间，以及多需要15~20%的内存。

8 边界条件

9 Fluent软件思路

9.1 fluent软件简介

本质上,fluent 只是一个求解器. 其提供英制(British),国际单位(SI)和厘米-克-秒(CGS)等单位制，这些单位之间可以相互转换，但是fluent 规定，对于边界特征，源项，自定义流场函数，外部创建的x-y 图散点图的数据文件数据，必须采用国际单位，对于网格文件，不管在创建时用的什么单位制，读入时一律采用国际单位制。因此需要进行缩放处理。

fluent4及以前版本，都是采用结构网格，而fluent5之后使用非结构网格。

但是兼容传统的结构网格和块结构网格。

9.1.1 Scheme表达式

Scheme 是Lisp 的一个分支，有着非常统一而又简单的命令格式：

（commandname argument1 argument2…）

每一个命令调用都是一个函数调用，因此也就会输出一个结果，命令名和变量名不区分大小写，但是只能以字母开头。

注释采用；；开头。

9.1.1.1 RP变量

获得变量的值：比如或者模拟时间

(rpgetvar ‘flow-time)

设置变量的值：

(rpsetvar ‘flow-time 0)

9.1.1.2 CX-变量

读取变量，比如说颜色历程表

(cxgetvar ‘cmap-list)

设置变量

(cxsetvar ‘def-cmap “rgb ”)

9.2 fluent软件数值计算方法

fluent 提供了三种数值计算方法

非耦合隐式算法、耦合显式算法、耦合隐式算法，分别适用于不可压、亚音速、跨音速、超音速乃至高超音速流动。

9.2.1 求解器的选择

fluent 有两类求解器：分离式和耦合式。默认使用分离式求解器，但是对于高速可压流动，由强体积力（如浮力或者旋转力）导致的强耦合流动，或者在非常精细的网格上求解的流动，需要考虑耦合式求解器。耦合隐式求解器所需内存大约是分离式求解器的1.5~2倍。如果计算机内存不够，可采用分离式或耦合显式，但是收敛性稍微差点。

需要注意的是：在分离式求解器中提供的几个物理模型，在耦合式求解器中是没有的，这些物理模型包括：VOF 模型，多相混合模型，欧拉模型，PDF 燃烧模型，预混合燃烧模型，部分预混合燃烧模型，烟灰和NOx 模型，Rosseland 辐射模型，熔化和凝固等相变模型，指定质量流量的周期流动模型，周期性热传导模型和壳传导模型等。而下列模型只在耦合式求解器中有效，在分离式求解器中无效：理想气体模型，用户定义的理想气体模型，NIST 理

想气体模型，非反射边界条件和用于层流火焰的化学模型。

9.2.1基于压力的非耦合算法

9.3 fluent湍流参数计算公式

采用不同的湍流模型，需要的参数如下：

9.3.1 湍流强度Turbulence Intensity

湍流强度小于1％时，可以认为湍流强度是比较低的，而在湍流强度大于10％时，则 可以认为湍流强度是比较高的。在来流为层流时，湍流强度可以用绕流物体的几何特征粗 略地估算出来。

(9-3-1)

式中，D H 为水力直径，Re D H 雷诺数是以水力直径为特征长度求出的。

9.3.2湍流的长度尺度l 与水力直径Hydraulic Diameter

湍流能量主要集中在大涡结构中，而湍流长度尺度l 则是与大涡结构相关的物理量。 在充分发展的管流中，因为漩涡尺度不可能大于管道直径，所以l 是受到管道尺寸制约的 几何量。湍流长度尺度l 与管道物理尺寸L 关系可以表示为：

l = 0.07L (9-3-2)

式中的比例因子0.07 是充分发展管流中混合长的最大值，而L 则是管道直径。在管道 截面不是圆形时， L 可以取为管道的水力直径。

湍流的特征长取决于对湍流发展具有决定性影响的几何尺度。在上面的讨论中，管道直径是决定湍流发展过程的唯一长度量。如果在流动中还存在其他对流动影响更大的物体，比如在管道中存在一个障碍物，而障碍物对湍流的发生和发展过程起着重要的干扰作用。在这种情况下，湍流特征长就应该取为障碍物的特征长度。

从上面的分析可知，虽然上式对于大多数管道流动是适用的，但并不是普遍适

用的，在某些情况下可以进行调整。

在FLUENT 中选择特征长L 或湍流长度尺度l 的方法如下：

1）对于充分发展的内流，可以用Intensity and Hydraulic Diameter（湍流强度与水力直 径）方法定义湍流，其中湍流特征长度就是Hydraulic Diameter（水力直径）D H 。

2）对于导向叶片或分流板下游的流场，可以用Intensity and Hydraulic Diameter（湍流 强度与水力直径）定义湍流，并在Hydrauli Diameter（水力直径）中将导向叶片或分流板 的开口部分的长度L 定义为特征长度。

3）如果进口处的流动为受到壁面限制且带有湍流边界层的流动，可以在Intensity and Length Scale 面板中用边界层厚度δ99通过公式l=0.4δ99计算得到湍流长度尺度l 。最后在 Turbulence Length Scale（湍流长度尺度）中输入l 的值。

9.3.3 湍流粘度比μt /μ

湍流粘度比μt /μ与湍流雷诺数Re t 成正比。湍流雷诺数的定义为：

Re t =k 2

ευ (9-3-3)

剪切层和充分发展的管道流动中的数值较大，其量级大约在100 到Re t 在高雷诺数边界层、

1000 之间。而在大多数外部流动的自由流边界上，μt /μ的值很小。在典型情况下，μt /μ的值在1 到10 之间。

用湍流粘度比定义流动时，可以使用Turbulent viscosity Ratio（湍流粘度比）或Intensity and Viscosity Ratio（湍流强度和粘度比）进行定义。前者适用于Spalart-Allmaras 模型，后 者适用于k −ε模型、k −ω模型和RSM 模型。

9.3.4 修正的湍流粘度

在使用Spalart-Allmaras 模型时，可以用湍流强度I 和长度尺度l 求出修正的湍流粘度，具

体公式如下：

=υ(9-3-4) avg Il

在使用FLUENT 时，如果在Spalart-Allmaras 模型中选择Intensity and Hydraulic

Diameter （湍流强度与水力直径）选项，则修正的湍流粘度就用这个公式求出。

9.3.5 湍流动能k

3k =(u avg I ) 2(9-3-5) 2

如果在使用FLUENT 时没有直接输入湍流动能k 和湍流耗散率ε 的值，则可以使用Intensity and Hydraulic Diameter（湍流强度与水力直径）、Intensity and Length Scale（湍流强度与长度尺度）或Intensity and Viscosity Ratio（湍流强度与粘度比）等方法确定湍流动能。

9.3.6 湍流耗散率

ε=C 3/4

μk 3/2(9-3-6) l

C μ为湍流模型中的一个经验常数，约为0.09.

在没有直接输入湍流动能k 和湍流耗散率ε 的情况下，可以用Intensity and Viscosity

Ratio （湍流强度与粘度比）定义湍流变量，实际上就是利用上述公式算出湍流耗散率ε 。 湍流耗散率ε 与湍流粘度比μ / μ t 和湍流动能k 的关系如下：

k 2μt -1ε=ρC μ() μμ(9-3-7)

式中C μ为湍流模型中的一个经验常数，其值约等于0.09。

如果计算风洞阻尼网下游试验段中的流场，可以用下式求出湍流耗散率ε ： ε≈∆kU ∞(9-3-8) L ∞

式中Δk 是湍流动能k 的衰减量，比如可以设为入口处k 值的10％，U ∞是自由流速度， L ∞是自由流区域的长度。该式是对高雷诺数各向同性湍流衰减指数律的线性近似， 其理论基础是衰减湍流中湍流动能k 的方程：

U ∂k =-ε(9-3-9) ∂x

如果用这种方法计算ε，还需要用（9-3-7）式检验计算结果，以保证湍流粘度比μ / μ t 不过大。虽然这种方法在FLUENT 中没有使用，但是可以用这种方法估算出自由流中的湍 流耗散率ε ，然后再用（9-3-5）式确定k ，最后在Turbulence Specification Method（湍流定义方法）下拉列表中选择K and Epsilon（ k 和ε ）并k 和ε 的计算结果输入到相应的栏目中。

9.3.7 比耗散率

如果知道湍流长度尺度l ，可以用下式确定ω ：

k 0.5

ω=0.25C μl (9-3-9)

式中C μ和长度尺度l 的取法与前面段落中所述相同。在使用Intensity and Hydraulic

Diameter （湍流强度与水力直径）或Intensity and Length Scale（湍流强度与长度尺度）定 义湍流时,fluent 就是采用这种方法.

ω 的值还可以用μ / μ t 和k 通过下式计算得出

ω=ρμμk μt -1()

在使用Intensity and Viscosity Ratio（湍流强度与粘度比）方法定义湍流时，FLUENT 就是使用上述关系式对湍流进行定义的。

9.3.8 雷洛应力分量

在使用RSM （雷诺应力模型）时，如果用户没有在Reynolds-Stess Specification Method （雷诺应力定义方法）的Reynolds-Stress Components（雷诺应力分量）选项中直接定义雷 诺应力的值，则雷诺应力的值将由给定的k 值计算得出。假定湍流是各向同性的，即： u i ' u ' j =0

且：

' ' u αu α=2k 3

如果用户在Reynolds-Stress Specification Method（雷诺应力定义方法）下拉列表中选择 K or Turbulence Intensity（ k 或湍流强度）时，FLUENT 就用这种方法定义湍流。

9.3.9在大涡模拟LES 中定义进口湍流

在使用速度进口条件时，可以将湍流强度作为对LES 进口速度场的扰动定义在边界条 件中。在实际计算中，根据湍流强度求出的随机扰动速度分量与速度场叠加后形成LES 算

法边界上的、随机变化的速度场。

9.4 Fluent中凝固熔化模型的理论

enthalpy-porosity 技术被应用于fluent 中的凝固和熔化处理过程中。在这项技术中，深化界面没有被很明显的跟踪。取而代之，大量的液体分数（液体分数显示那些流体组成的单元体积）被联合到每个单元在整个区域内。这个液体所占份额在热平衡的基础上通过反复计算被估计。

糊状的地方是流体分数在0-1之间的区域。这些糊状的地方在模拟过程中假冒为在多孔性从1-0递减的凝固材料中的多孔介质。当这些材料完全凝固成一个单元时，多孔性变为0因此速度也降为0。

9.4.1 凝固熔化模型中的能量方程

材料的焓能用显焓值h 和潜热ΔH 来计算：

H =h +∆H

其中,h

h =h ref +⎰c p dT T ref T

其中

h ref 参考焓

T ref 参考温度

c p 定压比热

流体分数β, 被定义为:

⎧0 (T

⎪T liquidus -T solidus

⎪1 (T >T liquidus )⎩

潜热内容可以查阅相关的材料的潜热性能表，L ：

∆H =βL

潜热的值从0（对固体）到L （对液体）之间变化。

在关于组分分离的多组分的凝固实例中; 例如，在组分传输的凝固或熔化

过程中，固相线和液相线被用来代替组分计算

T solidus =T melt +

T liquidus =T melt +solutes ∑K mY i i i

solutes ∑mY i i

Ki 是溶质i 的分离系数，是固体和与液体界面的浓度比率，Yi 是溶持i 的质量分数，mi 是液相线表面考虑Yi 之后的梯度。它被假设为混合物的最后一种组分材料是溶剂并且其它的组分是溶质。 对于凝固/熔化的问题，能量方程可以写作：

∂ (ρH ) +∇⋅(ρuH ) =∇⋅(k ∇T ) +S ∂t

其中: u 为流体速度, S 为源项.

9.4.2 凝固熔化模型中的动量方程

在enthalpy-porosity 技术把糊状区域（部分凝固的区域）看作为多孔介质。每个单元的多孔性在单元中设置相等的流体阻力。对于全凝固的区域，多孔性为0，这些区域的速度也为0。动量的损失是由于在糊状区域的多孔性的减少造成的如下式：

(1-β) 2 S =3A mush (u -u p ) β+ε

β是流体体积分数，ε是一个小于0.0001数，防止被0除，Amush 是一个糊状区域的连续数， 是固体速度由于拉凝固材料出范围（也指牵连速度）。

糊状区域的常数是测试阻尼的振幅的尺度；这个值越大，当凝固时的速度到0越是梯度大。非常大的值可能造成凝固振荡。在连续的铸造处理过程中，牵连速度包括来说明凝固材料从主体区的连续性运动。在上式出现一个项，允许一个新的凝固材料在牵连速度下运动。

如果凝固材料没有被从空上主体区拉出，即

9.4.3 凝固熔化模型中的湍流方程

Sink 更能说明在固体区域出现糊状和凝固区。Sink 项和动量方程的sink 项非常相似。 u p =0

(1-β) 2

S =3A mush φ β+ε

Φ代表被（k, ε, ω, 等）求解的湍流数量，在糊状区是一个常数