第三章 控制方程及求解方法

第三章 控制方程及求解方法

本章介绍FDS 的控制方程并概述一个通用的求解方法。在后面的章节中将详细描述个别方程。控制方程是一组经过适当的简化和近似的偏微分方程。该数值计算方法基本上是由一个有限差分近似方程和一个及时更新这些方程的程序组成的。

3.1 控制方程

本节介绍了牛顿流体的基本的质量、动量和能量守恒方程。这些方程式可以在几乎任何关于流体动力学或CFD 的书籍中找到。Anderson 等人[13]的著作对方程的描述是一个特别有用的参考，方程所用的符号、各种近似都是从他们那引用的。请注意，这是一组由6个方程组成的偏微分方程组，共有6个未知数，三维空间的所有函数和时间：密度ρ，三个方向的速度分量u =[u ，v ，w] T ，温度T ，压力p 。 3.1.1 质量输运

用密度表示的的质量守恒定律：

或以单独的气态物质αY 表示的质量守恒定律：

这里∑'''='''αα，b b m

m

表示的是物种通过蒸发液滴或颗粒的生产速率。由于1Y =∑α和

0m ='''∑α

和 b

b m m '''='''∑ α，

，并且根据定义它是假定0Y D =∇∑ααρ的，因此要对所有物种产生的原始的质量守恒方程进行总结。一般来说，这最后的说法是不正确的。然而，

输运方程是求解总质量而不是单个种类的质量，这意味着未知物体种类的扩散系数的选取，使所有的扩散通量的总和为零。

3.1.2 动量输运

动量守恒方程的形式为：

uu 项表示一个双值张量。在矩阵表示法中，u =[u ，v ，w]T ，双值由张量乘以向量u 和u T 得到的。∇·ρuu 项是通过应用矢量运算子∇=⎪

⎪⎭

⎫

⎝⎛∂∂∂∂∂∂z y x ，，而形成的矢量。动量方程中的f b 项表示外力，如液滴所受到的阻力。应力张量τ

ij 定义为：

S ij 项为用传统的张量表示法表示的对称张量应变率，符号μ表示流体的动力粘度。

整体计算可以被视为直接数值模拟（DNS ）或作为一个大涡模拟（LES ），直接数值模

拟耗散项直接计算，大涡模拟的大尺度涡直接计算并对亚格子尺度耗散过程进行模拟。数值算法的设计，使LES 网格细化后变为DNS 。大多数FDS 软件的应用都是LES 。例如，在通过一个大的、多室的密闭环境中模拟的烟流，它是不可能直接解决燃烧和运输过程。然而，对于小规模的燃烧实验，可以直接计算燃烧和运输过程。

第5章详细描述了动量方程和压力方程的数值解方法。为了概述下面的求解过程，经过充分的考虑把动量方程写为：

压力方程为：

其中压力方程是通过考虑动量方程的散度得到的。 3.1.3 能量输运

以显焓h s 表示的能量守恒方程为：

显焓是温度的函数

注意物质导数的使用， D （）/ DT=∂（）/∂T +u ·∇（）。q ''' 项表示化学反应中单位体积内的热释放速率。b q

''' 表示转移到蒸发液滴中的能量。q '' 表示导电性和热辐射通量：

其中k 是热导率。

3.1.4 状态方程

在模型中应用了适合低马赫数的Navier-Stokes 方程的一种近似形式。这种近似形式包含了筛选出的声波，同时允许温度和密度发生大的变化[9]。这里给出得方程组的椭圆几何特征与低速流、热对流过程一致。在实践中，这意味着该空间分辨的压力()z y x p ，，被替换为 “平均”或“背景”的压力()t z p m ，时，这是一个在地面之上仅和时间和高度有关的函数。

考虑物质导数的背压，并把结果代入能量守恒方程中，产生了速度散度u ⋅∇，这是一个在数值算法中重要的表达式，因为它有效地消除了通过求解输运方程来得到比焓的需要。能量守恒方程的源项被纳入到散度中，这一情况也出现在质量输运方程中。温度是通过把从密度和背压代入状态方程得到的。

3.2 求解过程

FDS 对一系列连接的整流非线性网格控制方程采用了二阶精度的有限差分近似。使用一个明确的二阶Runge-Kutta 计划及时更新流动的变化。本节介绍如何使用这种算法来及时

推进密度、组分质量分数、速度分量、背景和扰动压力。令n ρ 、n

Y α、u n 、n m p 、

表示

第n 个时间步长的这些变量。

1、计算“平均斑块”的速度场（参见5.6.3节）；

2、使用明确的欧拉步长来估计下一时间步长的ρ、αY 。例如，密度的估算为：

3、在网格边界交换ρ*和*αY 的值；

4、应用ρ*和*αY 的边界条件；

5、使用估计的热力学量计算散度*

⋅∇u 。请注意，在这个阶段，没有估计在下一时刻的速度场，只计算它的散度。

6、直接求解每个网格压力波动的泊松方程：

需要注意的是矢量n

F = ()n

n

u F

，ρ是使用补丁平均速度计算的，补丁平均场的散度

已经明确计算过。

7、估计下一时间步长的速度：

请注意，估算的速度场的散度和利用热力学量估算的散度∇·u *是恒等的。 8、检查在这一点的时间步长，确保

如果时间步长太大，则降低，使其满足约束条件并将该程序返回到的开始时间步长。如果时间步长满足稳定性条件时，程序将继续校正。更详细的稳定性信息，请参见5.5节。

这就结束了时间步长的“预测”阶段。此时，

值和u n+1的组分在通过MPI 调用的

网格边界上进行交换。

1、计算“平均斑块”的速度场（参见5.6.3节）；

2、应用Runge-Kutta 的第二部分更新质量变量。例如，密度校正为：

3、在网格边界交换ρ*和Y α*的值；

4、应用ρ*和Y α*的边界条件；

5、使用校正的热力学量计算散度*

⋅∇u 。请再次注意，这些点的速度场并没有校正； 6、使用估算的量计算压力波动：

7、应用Runge-Kutta 的第二部分更新速度：

再次注意，修正的速度场的散度和前面计算的散度是恒等的。

8、时间步长结束时，

值和u n+1的组分在通过MPI 调用的网格边界上进行交换。

3.3 空间离散

控制方程中的空间导数写成直线网格上的二阶精度有限差分。整个域是一个被分为很多矩形网格的矩形框。每个单元格的指数i 、j 和k 分别表示单元格在x 、y 和z 方向上的位置。标量被分配在每个网格单元的中心，因此，n

ijk ρ表示第n 个时间步长时网格中心的密度，它的位置指标为i 、j 和k 。像速度这样的矢量被分配在它们相应的网格表面上。例如，n

ijk u 表示位置为i 、j 、k 的单元格的正面的x 向速度分量；n

jk 1-i u ，表示相同单元格的反面的x 向速度分量。