CFD有限体积2-离散方程

参考文献： 

F. Moukalled, L. Mangani, M. Darwish. The finite volume method in computational fluid dynamics. An advanced introduction with Openfoam and Matlab[M]. Chapter 5.

网页资料：

数值分析8-高斯型求积公式.ppt (book118.com)

FVM: vertex-centered, cell-centered, or cell-vertex? -- CFD Online Discussion Forums (cfd-online.com)

FVM的优势在于可以直接在物理域上进行数值计算，而不需要向计算域映射，包括两步：

1. 在每个网格单元上对微分控制方程进行积分，生成半离散化方程

2. 假设参数分布型线，生成离散化方程

        1. 生成半离散化方程

       对于标量phi，其通用守恒方程可写作

 从左至右依次为瞬态项、对流项、扩散项和源项。定常形式为

对单元C

积分得

根据散度定理，得半离散化方程

        2. 生成离散化方程         

        对对流通量和扩散通量进行积分

记

总通量为

 根据散度定理利用面积分代替体积积分

上述面通量在单元界面上预估而不是在单元体内预估。根据高斯积分定理对界面f有

该式表示界面f上共有ip(f)个积分点，wip为权函数，积分点和权函数的选择决定了数值积分精度。 

 以最简单的梯形积分规则为例，做梯形积分时，控制单元每个界面上只有一个积分点，该点处于界面中心，权函数为1.该方法具有二阶精度，适用于2D和3D问题。对2D问题，有两个积分点（ip(f)=2)，分别位于

L为积分线总长度，对应权函数均为1/2，该方法具有三阶精度。积分点越多，计算成本越高。

        最终得

         对源项进行积分，同理

        对FVM，习惯上采用一个积分点构成二阶计算精度。当采用一个积分点时，对下式

有 

将上式化为代数方程组。对通量进行线性化

C和F为界面两侧控制单元。右式前两项为两个控制单元对界面通量的线性贡献，最后一项为非线性项。

对所有界面通量

 对控制单元源项，也可将其划分为线性项和非线性项

若源项为恒定源项，则有 

将以上线性化结果代入半离散积分方程中，整理为

有

        对于非定常问题，在体积分外还需进行时间积分

         3. 边界条件

        两个最常见的边界条件为Dirichlet条件和Neumann条件，前者给定边界上参数值，后者给定边界上参数导数值。

        3.1 Dirichlet边条

 此时对边界上控制单元，边界界面上对流通量为

 因此

         3.2 Neumann边条

    给定    

        因此

         4. 离散精度

        界面上通量和源项值由界面/单元中心点参数值决定，参数的分布为哪类型线，这两者决定了数值积分精度。例如，基于网格单元中心xc参数值和线性参数分布时，x处的phi值为

         4.1 参数空间变化假设

        x处phi值对xc处phi值的泰勒展开式为

 与

比较，可得二者误差与 abs(x-xc)的平方成比例，因此空间精度为二阶（second order spatial accuracy）。

        4.2 网格单元参数平均值假设

当C位于网格单元中心点时，第二项为0，因此有上式

代表利用网格单元中心点处phi值代替网格单元平均值的精度为二阶。

对于二阶方法，界面积分和参数值分布都应为二阶精度。当需要提高方法精度时，可以从界面积分方法和参数值分布形式两个方面进行改善。

        5. 离散方程性质

        5.1 守恒性

        界面两侧通量值相等但大小相反。守恒性为FVM固有属性。

        5.2 精度

        精度严格来讲应指数值解与精确解之间误差，但通常精确解未知，因此以离散格式的截断误差（离散精度上文第4小节）来表征精度。截断误差并不代表数值解与精确解之间的误差，但其代表了当网格加密时，离散误差的下降速度。如二阶精度代表当网格加密一倍时，离散误差会下降4倍。

        5.3 收敛性

       当两个迭代步之间数值解相差值小于某一规定值时，称计算收敛。

        5.4 一致性

        当时间步和网格大小均趋于无限小时，数值解趋于微分方程精确解。因此，当离散误差应是delta(t)和delta(x)的指数型函数。当离散误差是delta(x)/delta(t)的函数时，delta(x)应比delta(t)更快趋于0。

        5.5 稳定性

        衡量代数方程组求解器的性质，当求解器可以求解足够多初始和边界条件组合下的代数方程组时，称该求解器是稳定的。例如，显式求解器和隐式求解器的稳定性不同。

        5.6 经济性

        求解一个问题所需时间。

        5.7 输运性（transportiveness）

        流动传热问题的输运性可由Pe描述

当Pe=0时，为扩散问题，某点参数Phi受边界及边界内所有点同样权重影响；当Pe大于0时，某点参数Phi受上游点影响的权重大于下游点影响的权重。因此输运性与迎风格式构建有关。选择数值格式应考虑物理问题的输运性，否则会出现不稳定解。

        5.8 有界性（boundedness）

        当求解域内无源项时，内部节点参数值应在边界节点参数值之间：

Fi为C相邻节点。 

        6. Vertex-centered FVM / Cell-centered FVM

         两种参数值布局方法：网格中心or网格顶点

         6.1 基于网格顶点的FVM

        流动参数值被存储于网格顶点（网格点），此时需要用到双网格系统在网格顶点周围构造控制单元。两种双网格生成方法：连接原始网格中心点生成第二套网格；将原始网格中心点与其界面中心点相连生成第二套网格，如下图所示。

        基于网格顶点的FVM精度较低，通量守恒型较差（见下文描述）。

FVM: vertex-centered, cell-centered, or cell-vertex? -- CFD Online Discussion Forums (cfd-online.com)        6.2 基于网格中心点的FVM

        最常用的FVM为基于网格中心点的FVM，即参数值存储于网格单元中心点，通常该方法具有二阶精度。主要缺点：在处理non-conjunctional单元（？具体含义）精度较低，在非正交网格上计算扩散项时鲁棒性较差。

例如，对下图C和F的平均值为f '处值，而非界面中心点f处值，此时离散精度达不到二阶。

对该方法，离散误差与网格光滑程度有很大关系。在生成网格的过程中，应保证网格斜率连续。

总结：只有在扭曲网格的情况下，vertex-centered格式优于cell-centered格式。在其他情况下，应优先选择cell-centered格式。

        7. 隐式格式/显式格式

        在CFD中，更倾向于优先选择隐式格式。

        8. 网格

        网格生成模块需要向求解器模块提供的信息包括：

        网格单元相关，网格单元的：索引、中心点、界面、相邻单元、界面法向量定义

        界面相关，界面的：索引，中心点，面向量、相邻单元

        计算域的边界