4.1 计算流体动力学简介

        计算流体动力学（Computatonal Fluid Dynamics，简称CFD）是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所作的分析。CFD的基本思想可以归结为：把原来在时间域及空间域上连续的物理量的场，如速度场和压力场，用一系列有限个离散点上的变量值的集合来代替，通过一定的原则和方式建立起关于这些离散点上场变量之间关系的代数方程组，然后求解代数方程组获得场变量的近似值。（周雪漪 计算水力学 清华大学出版社 1995  陶文铨  数值传热学（第二版）西安交通大学出版社 2001 郭鸿志 传输过程数值模拟 冶金工业出版社 1998）

        采用CFD的方法对流体流动进行数值模拟，通常包括以下步骤：

        （1） 建立反映工程问题或物理问题本质的数学模型。具体说来就是要建立反映问题各个量之间关系的微分方程及相应的定解条件，这是数值模拟的出发点。没有正确完善的数学模型，数值模拟就毫无意义。流体的基本控制方程通常包括质量守恒方程、动量守恒方程、能量守恒方程，以及这些方程相应的定解条件。

        （2）寻求高效率、高准确度的计算方法，即建立针对控制方程的数值离散化方法，如有限差分法、有限元法、有限体积法等。这里的计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。这些内容，可以说是CFD的核心。

        （3）编制程序和进行计算。这部分工作包括计算网格划分、初始化和边界条件的输入、控制参数的设定等。这是整个工作中花时间最多的部分。由于求解的问题比较复杂，比如Navier-Stokes方程就是一个十分复杂的非线性方程，数值求解方法在理论上不是绝对完善的，所以需要通过实验加以验证。正式从这个意义上讲，数值模拟又叫数值试验。

        （4）显示计算结果。计算结果一般通过图表等方式显示，这对检查和判断分析质量和结果有重要参考意义。

        以下将以以上4点为纲领，详细叙述模型的建立和计算过程。

        本文使用商业CFD软件Fluent及其前置软件Gambit来进行模拟，使用压力基解算器求解。压力基求法原称为分离求解方法，即分别求解各个控制方程的方法。由于控制方程是非线性的，因此求解必须经过多次迭代才能获得收敛解。图4-1为求解流程图，其过程可概述如下：

        （1）流场变量更新。在第一次计算时，变量由初始化过程更新。在随后的计算中，没迭代一部既得到一个更新的解。

        （2）用当前压强和质量通量的值求解动量方程，以得到新的速度场。

        （3）因为（2）中得到的速度场的数值解无法完全满足连续方程，于是再求解压强修正方程。压强修正方程是由连续方程导出的泊松型方程，求解这个方程可以得到对压强场、速度场和质量通量的修正，进而使连续方程得到满足。

        （4）利用前面求出的解，求解湍流方程、能量方程和组元方程。

        （5）在多项流计算中那个如果考虑相间干扰，则需要通过求解弥散相轨迹计算得到连续相方程中的源项解。

        （6）检验收敛条件是否满足。如果收敛条件满足，则停止计算。如果计算没有收敛，则继续迭代过程。

图 4-1

4.2 流体的控制方程及其离散化

        设[tex]\phi[/tex]为流体系统在[tex] t [/tex]时刻所具有的某种物理量，[tex]\phi_{p}[/tex]表示单位质量流体所具有的这种物理量。

由此推导出输运方程为：

 

[tex]\frac{d\phi}{dt}=\frac{\partial}{\partial t} \int \limits_{CV}\phi_{p} \rho d V+\int \limits_{CS}\phi_{p} \rho u_{n} dA[/tex]                (4-1)

 

其中[tex]u_{n}[/tex]为速度在微元面法线方向的投影。该式表明，流体某一物理量时间的全变化率等于控制体内这种物理量的时间变化率与经过控制面的这种物理量的净通量之和。

        单位质量的流体的质量[tex]\phi_{p}=1[/tex]，质量为[tex]\phi=\int \limits_{V}\rho d V=m[/tex]因为质量守恒，所以有：

 

[tex]\frac{d \phi}{d t}=\frac{d m}{d t}=0[/tex]

 

由(4-1)式可得：

 

[tex]\frac{\partial}{\partial t}\int \limits_{CV}\rho d V+\int\limits_{CS}\rho u_{n} d A=0[/tex]                (4-2)

 

此即是连续方程的积分形式，其微分形式为：

 

[tex]\frac{\partial \rho}{\partial t}+\nabla \cdot(\mathbf{\rho u})=0[/tex]                (4-3)

 

在直角坐标系中表示为：

 

[tex]\frac{\partial \rho}{\partial t}+\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=0[/tex]        (4-4)

 

        单位质量的流体动量为[tex]\phi_{p}=u[/tex]，流体系统动量为[tex]p=\int \limits_{V}\rho \mathbf{u} d V[/tex]，由(4-1)式得：

 

[tex]\frac{d}{d t}\int \limits_{V}\mathbf{u}\rho d V=\frac{\partial}{\partial t}\int \limits_{CV}\mathbf{u}\rho d V+\int \limits_{CS}\mathbf{u}\rho u_{n}d A[/tex]                (4-5)

 

        根据质点系动量定理，流体系统动量的时间变化率等于作用在系统上的外力之和。作用在流体流体微元上的外力有质量力和表面力之分。定义[tex]f[/tex]为作用在单位质量流体上的质量力分布函数，定义[tex]p_{n}[/tex]为作用在单位质量流体上的表面力分布函数。那么运动方程的完整积分表述方式为：

 

 [tex]\frac{\partial}{\partial t}\int \limits_{CV}\mathbf{u}\rho d V+\int \limits_{CS}\mathbf{u}\rho u_{n} d A=\int \limits_{CV}f \rho d V+\int \limits_{CS}p_{n} d A[/tex]         (4-6)

 

 [tex]\because \quad \frac{d}{d t}\int \limits_{V}\mathbf{u}\rho d V=\int \limits_{V}\rho \frac{d \mathbf{u}}{d t}d V[/tex]

 

又根据奥高定理：

 

[tex]\int \limits_{CS}p_{n} d A=\int \limits_{CS}\mathbf{n} \cdot P d A=\int \limits_{CV}\nabla \mathbf{P} d V[/tex]

 

于是式(4-6)变为：

 

[tex]\int \limits_{CV}(\rho \frac{d \mathbf{u}}{d t}-\rho \mathbf{f}-\nabla \mathbf{P}) d V=0[/tex]

 

[tex]\therefore \quad \rho\frac{d \mathbf{u}}{d t}=\rho\mathbf{f}+\nabla\mathbf{P}[/tex]               (4-7)

 

此即是运动方程的微分形式

在直角坐标系中，变为：

 

[tex]\rho(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z})=\rho f_{x}+\frac{\partial p_{xx}}{\partial x}+\frac{\partial p_{xy}}{\partial y}+\frac{\partial p_{xz}}{\partial z}[/tex]                (4-8a)

[tex]\rho(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z})=\rho f_{y}+\frac{\partial p_{yx}}{\partial x}+\frac{\partial p_{yy}}{\partial y}+\frac{\partial p_{yz}}{\partial z}[/tex]                (4-8b)

[tex]\rho(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z})=\rho f_{z}+\frac{\partial p_{zx}}{\partial x}+\frac{\partial p_{zy}}{\partial y}+\frac{\partial p_{zz}}{\partial z}[/tex]                (4-8c)

 

由本构方程：

 

[tex]p_{xx}=-p-\frac{2}{3}\mu \nabla \mathbf{u}+2\mu\frac{\partial u}{\partial x}[/tex]

[tex]p_{yy}=-p-\frac{2}{3}\mu \nabla \mathbf{u}+2\mu\frac{\partial v}{\partial y}[/tex]

[tex]p_{zz}=-p-\frac{2}{3}\mu \nabla \mathbf{u}+2\mu\frac{\partial w}{\partial z}[/tex]

[tex]p_{xy}=p_{yx}=\mu(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})[/tex]

[tex]p_{zx}=p_{xz}=\mu(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})[/tex]

[tex]p_{yz}=p_{zy}=\mu(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})[/tex]

 

代入式（4-8）得到N-S方程：

 

[tex]\rho \frac{Du}{Dt}\!=\!\rho f_{x}\!-\!\frac{\partial p}{\partial x}\!+\!\frac{\partial}{\partial x}[\mu(2\frac{\partial u}{\partial x}\!-\!\frac{2}{3}\nabla \mathbf{u})]\!+\!\frac{\partial}{\partial y}[\mu(\frac{\partial v}{\partial x}\!+\!\frac{\partial u}{\partial y})]\!+\!\frac{\partial}{\partial z}[\mu(\frac{\partial u}{\partial z}\!+\!\frac{\partial w}{\partial x})][/tex]        (4-9a)

 

[tex]\rho \frac{Dv}{Dt}=\rho f_{y}-\frac{\partial p}{\partial y}+\frac{\partial}{\partial x}[\mu(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})]+\frac{\partial}{\partial y}[\mu(2\frac{\partial v}{\partial y}-\frac{2}{3}\nabla \mathbf{u})]+\frac{\partial}{\partial z}[\mu(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z})][/tex]        (4-9b)

 

[tex]\rho \frac{Dw}{Dt}=\rho f_{z}-\frac{\partial p}{\partial z}+\frac{\partial}{\partial x}[\mu(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})]+\frac{\partial}{\partial z}[\mu(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z})]+\frac{\partial}{\partial z}[\mu(2\frac{\partial w}{\partial z}-\frac{2}{3}\nabla \mathbf{u})][/tex]        (4-9c)

 

当流体为均质不可压缩常粘度时，N-S方程简化为：

 

[tex]\rho(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z})=\rho f_{x}-\frac{\partial p}{\partial x}+\mu(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}})[/tex]        (4-10a)

[tex]\rho(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z})=\rho f_{y}-\frac{\partial p}{\partial y}+\mu(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}+\frac{\partial^{2}v}{\partial z^{2}})[/tex]        (4-10b)

[tex]\rho(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z})=\rho f_{z}-\frac{\partial p}{\partial z}+\mu(\frac{\partial^{2}w}{\partial x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}+\frac{\partial^{2}w}{\partial z^{2}})[/tex]        (4-10c)

 

        单位质量的流体的能量为[tex] \phi_{p}=e+u^{2}/2[/tex]，其中[tex]e[/tex]为单位质量流体的热力学能。则[tex]\phi=\int \limits_{V}(e+u^{2}/2)\rho d V[/tex]，由（4-1）式得：

 

[tex]\frac{d}{dt}\int \limits_{V}(e+u^{2}/2)\rho d V=\frac{d}{dt}\int \limits_{CV}(e+u^{2}/2)\rho d V+\int \limits_{CS}(e+u^{2}/2)\rho v_{n}d V[/tex]        (4-11)

 

根据能量守恒和转换定律，流体系统中能量的时间全变化率等于作用在系统上的质量力和表面力所作的功率和与外界换热率之和。则得：

 

[tex]\frac{d}{dt}\int \limits_{V}(e+u^{2}/2)\rho d V=\int \limits_{V}f \rho \mathbf{u} d V+\int \limits_{CA}p_{n} \mathbf{u} d A+Q[/tex]                (4-12)

 

则能量守恒方程积分形式为：

 

[tex]\frac{d}{dt}\int \limits_{CV}(e+u^{2}/2)\rho d V+\int \limits_{CS}(e+u^{2}/2)\rho v_{n}d V=\int \limits_{V}f \rho \mathbf{u} d V+\int \limits_{CA}p_{n} \mathbf{u} d A+Q[/tex]                (4-13)      

 

其微分形式为：

 

[tex]\frac{\partial \rho T}{\partial t}+\frac{\partial \rho u T}{\partial x}+\frac{\partial \rho v T}{\partial y}+\frac{\partial \rho w T}{\partial z}=\frac{\partial}{\partial x}(\frac{\lambda}{c_{p}}\frac{\partial T}{\partial x})+\frac{\partial}{\partial y}(\frac{\lambda}{c_{p}}\frac{\partial T}{\partial y})+\frac{\partial}{\partial z}(\frac{\lambda}{c_{p}}\frac{\partial T}{\partial z})+Q[/tex]                (4-14)

 

对于不可压缩，常物性无内热源可简化为：

 

[tex]\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}=\frac{\lambda}{\rho c_{p}}(\frac{\partial^{2}T}{\partial x^{2}})+\frac{\partial^{2}T}{\partial y^{2}}+\frac{\partial^{2}T}{\partial z^{2}}[/tex]        (4-15)

      

        比较基本的控制方程，可以看出，尽管这些方程中因变量各不相同，但它们均反映了单位时间单位体积内物理量的守恒性质。用 [tex]\phi[/tex]表示通用变量，则上述控制方程可以表示成以下通用形式：

 

[tex]\int \limits_{V}\frac{\partial \rho \phi}{\partial t}dV+\oint \rho \phi \mathbf{u}\cdot d\mathbf{A}=\oint \Gamma_{\phi}\nabla \phi \cdot d\mathbf{A}+\int \limits_{V}S_{\phi}dV[/tex]                (4-16)

 

其中[tex]\Gamma_{\phi}[/tex]是变量[tex]\phi[/tex]的扩散系数

        [tex]S_{\phi}[/tex]是[tex]\phi[/tex]的单位体积源项

转化为微分形式为：

 

[tex]\frac{\partial(\rho \phi)}{\partial t}+\nabla \cdot (\rho \mathbf{u}\phi)=\nabla \cdot (\Gamma_{\phi}\nabla \phi)+S_{\phi}[/tex]        (4-17)

 

式中从左到右依次代表瞬态项、对流项、扩散项和源项。

 

        对反映物理问题的控制方程离散化是数值模拟里核心的 步骤，离散化方式的选择决定着计算的稳定性和精度。目前成熟和应用较多的离散化方式有：有限差分法（Finite Difference Method）、有限元法（Finite Element Method）和有限体积法（Finite Volume Method）。  

        有限体积法的基本思路是：将计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体积；将待解的微分方程对每一个控制体积积分，从而得到一组离散方程。其中的未知数是网格点上的因变量[tex]\phi[/tex]。为了求出控制体积的积分，必须假定[tex]\phi[/tex]值在网格点之间的变化规律。在有限体积法的网格中，计算节点（node）被看成是控制体积（cell）的代表。在离散化过程中，将一个控制体积上的物理量定义并储存在该节点处。在有限体积法中，计算节点的值即因变量[tex]\phi[/tex]处于控制体积的中心，有异于网格交点。其数值的确定方式和在节点间的数值传递（变化）方式即是离散化方式。 

        Fluent对于时间上的偏导数，有一阶和二阶的显式或隐式格式可选。由于本文讨论的问题为定常流动，因此控制方程中没有瞬态项。对于空间上的离散方式即对对流项和扩散项的[tex]\phi[/tex]值的离散，通常有一阶迎风、二阶迎风、中心差分和QUICK等格式。一阶迎风格式将控制体积界面上的未知量[tex]\phi_{f}[/tex]恒取上游的计算节点的值，即：

 

[tex]\phi_{f}=\phi[/tex]

 

而二阶迎风格式在一阶迎风的基础上考虑了物理量在节点间分布曲线的曲率影响，其表示为：

 

[tex]\phi_{f}=\phi+\nabla \phi \cdot \mathbf{r}[/tex]

 

其中[tex]\phi \quad \nabla \phi[/tex]分别为计算节点的值和其在上游节点的梯度

        [tex]\mathbf{r}[/tex]为上游节点质心道界面质心的方向向量

 

中心差分格式计算控制体积界面上的值[tex]\phi_{f}[/tex]为：

 

[tex]\phi_{f}=\frac{1}{2}(\phi_{0}+\phi_{1})+\frac{1}{2}(\nabla \phi_{0}\cdot \mathbf{r_{0}}+\nabla \phi_{1}\cdot \mathbf{r_{1}})[/tex]

 

其中[tex]\phi_{0} \quad \phi_{1}[/tex]指共面于[tex]f[/tex]的控制体的计算节点值

 

        本文所建立的模型流动方向与网格方向基本一致，且为结构化网格，因此选用一阶迎风格式。一阶迎风格式具有稳定性高，计算速度快的优点。对于[tex]\nabla \phi[/tex]的离散化，在Fluent中使用Green-Gauss理论，分别有Cell-Based和Node-Based两种方法。对于控制体中心[tex]c0[/tex]的值[tex]\phi[/tex]，其梯度的离散化形式为：

 

[tex](\nabla \phi)_{c0}=\frac{1}{V}\sum \limits_{f}\overline{\phi_{f}}\mathbf{A_{f}}[/tex]        (4-18)

 

Cell-Based方法：

 

[tex]\overline{\phi_{f}}=\frac{\phi_{c0}+\phi_{c1}}{2}[/tex]        (4-19)

 

Node-Based方法：

 

[tex]\overline{\phi_{f}}=\frac{1}{N_{f}}\sum \limits_{n}^{N_{f}}\overline{\phi_{n}}[/tex]        (4-20)

 

 其中[tex]N_{f}[/tex]——界面上的节点数

         [tex]\overline{\phi_{n}}[/tex]——包围节点的单元的值的加权平均

Node-Based方法比默认的Cell-Based方法精度更高，因此本模型选用Node-Based方法计算[tex]\nabla \phi[/tex]。

 

 4.3 几何模型的建立及网格化

        本文用Gambit建立两个平直型翅片的模型，其尺寸参数分别为：

                                        模型一                                        模型二

水力直径mm                    2.45                                             2.79

翅片间距mm                    2.0                                                2.1

翅片厚度mm                    0.3                                                0.3

翅片高度mm                     4.7                                               6.5

隔板厚度mm                     0.5                                               0.5

流动长度mm                     50                                                50

 

        由于翅片很薄，其与板片接触部分以钎焊连接，换热量少。为了简化模型，减少网格数量节约计算资源将此部分去除。此外，参考1947年美国海军研究署的实验装置，取板翅式换热器一个完整流道和其相邻的两个流道的一半，建立模型结构如图4-2。

图 4-2

        图中1-2-3部分为流体区域，上下板片及翅片被分割为12部分。注意在Gambit中对体的切割有连接（connected）和非连接（unconnected）之分，非连接的体在其交接面出产生连个重叠的面，在载入Fluent时默认为两个壁面（wall）。如果两个壁面间有参数的传递，比如流体/固体内部或流固耦合处。则这种设置将阻断参数的传递，使其彼此成为孤立的体。若要使其之间有参数的传递，即使其网格结构不同或属于流固两种区域，则需要使其连接或分别定义为interface边界条件。interface边界条件为不同网格结构的区域提供了一种插值的参数传递方式，然而其设置繁杂，会增加计算量且存在一定误差。此处体的分割均为连接。即两个分割的体公用一面，在载入Fluent时会自动生成wall和wall.shadow两个壁面。其中一个为虚拟壁面，Fluent自动耦合两壁面间参数。

        由于平直型翅片几何结构较为规整，因此使用六边形的网格单元以map的方式分别网格化图 4-2中的各区域。考虑到网格密度对计算精度的影响，先行对各待网格化的区域在流动方向上的面和流体入口面的边网格化，设置合适的网格数量。得到整体网格如下图 4-3

         图 4-3

        为了减少入口效应和出口效应，在流体的进口处和出口处分别建立延长段。进口延长段长度为进口水力直径的1.5倍，出口延长段长度为出口水力直径的5倍，如上图。

        

4.4 解的设置

4.4.1 边界条件

        此模型考虑到计算资源有限，采用压力基隐式求解方式。因此将流体近似认为不可压缩流体，密度取为常数。因为速度入口条件应用于不可压缩流体，且比质量流量入口条件稳定，因此流体入口定义为速度入口边界条件。入口流体的温度保持为300K，以考察各模型对流体的加热情况。速度方向设置为垂直于边界（Normal to Boundary），所以忽略了导流片对流动分配的影响，近似认为入口流动是均匀的。

        模型一预期取得十个实验数据点，雷诺数从500～9500以1000递增。模型二取雷诺数500和9500两个数据点与模型一对比。模型一速度入口的速度值与雷诺数对应为：

速度m/s                雷诺数

1.4                          500

4.2                          1500

6.95                        2500

9.7                          3500

12.5                        4500

15.3                        5500

18                           6500

20.87                      7500

23.65                      8500

26.4                        9500

模型二速度与雷诺数对应为：

速度m/s                雷诺数

2.2                          500

39.36                      9500

         

        对于雷诺数由2500～9500的情形，虽然实际上流体处于过度流动阶段。然而在模拟中使用湍流模型，调整湍流强度来模拟其流动状况。界面上的湍流强度定义为：

 

[tex]I=\frac{\sqrt{u^{'2}+v^{'2}+w^{'2}}}{u_{avg}}[/tex]        (4-21)

 

其中[tex]u^{'},v^{'},w^{'}[/tex]是速度脉动量，[tex]u_{avg}[/tex]是平均速度。

 

        湍流强度小于1%时，可以认为湍流强度是比较低的，而在湍流强度大于10%时，则可以认为湍流强度是比较高的。内流问题的湍流强度取决于上游的流动状态。如果上游是没有充分发展的未受扰动的流动，则进口处可以使用低湍流强度。如果上流是充分发展的湍流，则进口处湍流强度可以达到几个百分点。鉴于此，考虑到本模型的实际情况，选用较小的湍流强度。其数值在1%到10%浮动，以水力直径为特征长度。

        压力出口边界条件指定出口边界处的静压，而静压值仅在流场为亚音速时使用。如果在出口边界上流场达到超音速，则边界上的压强将从流场内部通过插值得到。其他流场变量均从流场内部插值获得。定义压力出口的静压值均为1000Pa（表压），无回流，回流温度均为零。

        定义板片的上下两个面为定热流密度（热通量）的壁面，热流密度均设置为[tex]20000W/m^{2}[/tex]。定义模型的左右各个面为对称边界条件。其余的面均定义为绝热壁面。

 

4.4.2 物理材质的物性参数定义

        对于Fluent中的物理材质，还需根据本模型的具体情况进行调整。

        空气的密度设置为常数，其数值为[tex]1.5 kg/m^{3}[/tex]。定压比热容以温度的多项式计算，设置为式（2-24）。导热系数以运动理论计算，设置为式（2-22）。运动粘度系数以Sutherland三系数公式计算，系数取默认值。

        固体材质为铝，其密度、比热容和导热系数均认为常数。其值分别设置为[tex]\rho=2770ka/m^{3} \quad c_{p}=963J/(kg\cdot K) \quad \lambda=175W/(m\cdot K)[/tex]。

         

4.4.3 湍流模型及壁面函数定义

        本文模型选用Standard[tex]k-\epsilon[/tex]湍流模型，标准的[tex]k-\epsilon[/tex]模型通过解两个独立的输运方程来求取湍动粘度[tex]\mu_{t}]/tex]和湍流尺度[tex]l[/tex]。湍动能[tex]k[/tex]和湍动耗散率[tex]\epsilon[/tex]的输运方程分别为：

[tex]\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{i}}(\rho k u_{i})=\frac{\partial}{\partial x_{j}}[(\mu+\frac{\mu_{t}}{\sigma_{k}})\frac{\partial k}{\partial x_{j}}]+G_{k}+G_{b}-\rho \epsilon-Y_{M}+S_{k}[/tex]        (4-22a)

 

[tex]\frac{\partial}{\partial t}(\rho \epsilon)+\frac{\partial}{\partial x_{i}}(\rho \epsilon u_{i})=\frac{\partial}{\partial x_{j}}[(\mu+\frac{\mu_{t}}{\sigma_{\epsilon}})\frac{\partial \epsilon}{\partial x_{j}}]+C_{1\epsilon}\frac{\epsilon}{k}(G_{k}+C_{3\epsilon}G_{b})-C_{2\epsilon}\rho \frac{\epsilon^{2}}{k}+S_{\epsilon}[/tex]        (4-22b)

 

其中：[tex]k[/tex]——湍动能

            [tex]\epsilon[/tex]——湍动耗散率

            [tex]G_{k}[/tex]——平均速度梯度对湍动能的影响系数

                        [tex]G_{k}=-\rho \overline{u_{i}^{'}u_{j}^{'}}\frac{\partial u_{i}}{\partial x_{i}}[/tex]

                        根据Boussinesq 假设：

                                [tex]G_{k}=\mu_{t}S^{2}[/tex]        (4-23)

                                其中[tex]S[/tex]为平均应力张量系数

            [tex]G_{b}[/tex]——浮力（弹性力）对湍动能的影响系数

                        [tex]G_{b}=\beta g_{i}\frac{\mu_{t}}{P_{rt}}\frac{\partial T}{\partial x_{i}}[/tex]        (4-24)

                        其中[tex]P_{rt}[/tex]——湍动能的普朗特数，默认设置为0.85

                                [tex]g_{i}[/tex]——重力加速度在[tex]i[/tex]方向的分量

                                [tex]\beta[/tex]——热力扩散系数

                                                [tex]\beta=-\frac{1}{\rho}(\frac{\partial \rho}{\partial T})_{p}[/tex]        (4-25)

             [tex]Y_{M}[/tex]——流体的可压缩性对湍动耗散率的影响系数

                        [tex]Y_{M}=2\rho \epsilon M_{t}^{2}[/tex]

                        其中[tex]M_{t}[/tex]——湍流马赫数

                                [tex]M_{t}=\sqrt{\frac{k}{s}}[/tex]        (4-26)

                                        其中[tex]s[/tex]——当地声速

              [tex]C_{1\epsilon},C_{2\epsilon},C_{3\epsilon}[/tex]——常数，默认设置值为1.44，1.92。

              [tex]\sigma_{k},\sigma_{\epsilon}[/tex]——分别是湍动能和湍动耗散率的普朗特数，默认设置值为1.0，1.3。

              [tex]S_{k},S_{\epsilon}[/tex]——用户自定义源项。

              [tex]\mu_{t}[/tex]——湍流粘度

                        [tex]\mu_{t}=\rho C_{\mu} \frac{k^{2}}{\epsilon}[/tex]        (4-27)        

                        其中[tex]C_{\mu}[/tex]——常数，默认设置值为0.09。

        标准的[tex]k-\epsilon[/tex]模型，是针对湍流发展非常充分的湍流流动来建立的，也就是说，它是一种针对高Re数的湍流计算模型，而当Re数比较低时，比如在近壁区域内的流动，湍流发展并不充分，湍流的脉动影响可能不如分子粘性的影响大，在更贴近壁面的底层内，流动可能处于层流状态。因此，对Re数较低的流动使用此模型计算就需要对壁面进行处理，以解决近壁面区域内的流动计算及低Re数时的流动计算问题。常用的解决办法是采用壁面函数法。

        流体近壁面可分为三个层：最里面一层称为”粘性子层“，流动几乎为层流，粘性是主要作用力；外面一层称为”充分湍流层“，湍流是主要作用力；它们之间有一个过渡层，粘性力和湍流的影响相当。本模型用标准壁面函数来模拟近壁面的流动。

动量：

 

[tex]U^{\ast}=\frac{1}{\kappa}\ln(E y^{\ast})[/tex]        (4-28)

 

其中[tex]\kappa[/tex]——冯·卡门系数，等于0.4187。

        [tex]E[/tex]——经验系数，等于9.793。

 

[tex]y^{\ast}=\frac{\rho C_{\mu}^{1/4} k_{p}^{1/2} y_{p}}{\mu}[/tex]        (4-29)

 

其中[tex]k_{p}[/tex]——p点的湍动能

        [tex]y_{p}[/tex]——p点离壁面的距离

当[tex]30\le y^{\ast} \le 300[/tex]时，标准壁面函数是正确的。在Fluent中，应在[tex]y^{\ast} \ge 11.225[/tex]时应用。当[tex]y^{\ast} \le11.225[/tex]时，则[tex]U^{\ast}=y^{\ast}[/tex]。

 

能量：

 

[tex]T^{\ast}=\frac{(T_{w}-T_{p})\rho c_{p} C_{\mu}^{1/4} k_{p}^{1/2}}{q}[/tex]        (4-30)

 

其中[tex]T_{w}[/tex]——壁面温度

        [tex]T_{p}[/tex]——p点温度

        至此，模型的建立和计算前的设置均已经完成。

[取消回复评论]

昵称 登录 E-mail: *
Web:
Twitter:

当有新评论通过 E-mail 通知我

(输入验证码)

or Ctrl+Enter