FVM in CFD 学习笔记_第8章_空间离散之扩散项

最新推荐文章于 2023-02-15 09:36:13 发布

梅冠华

最新推荐文章于 2023-02-15 09:36:13 发布

阅读量3.3k

点赞数 11

分类专栏： FVM in CFD 学习笔记

本文链接：https://blog.csdn.net/meiguanhua/article/details/104789091

版权

FVM in CFD 学习笔记专栏收录该内容

14 篇文章 135 订阅

订阅专栏

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 8 Spatial Discretization: The Diffusion Term

本章将详细讲解由Laplace算子所表示的扩散项的空间离散方法。流体控制方程中，扩散项和对流项反映了两种截然不同的物理现象，所以两者的离散方法也是不同的，需要分开来讲解。本章首先讲解二维矩形计算域在笛卡尔网格上的含源项的扩散方程的离散方法，然后讲解Dirichlet、Von Neumann、mixed和symmetry边界条件的添加方法。接下来，介绍下在非正交的笛卡尔网格上的离散方法，并详细讲讲非正交的结构和非结构网格上的离散方法、非正交交叉扩散，还会提及非各向同性扩散的处理方法、对高度非线性系统的欠松弛处理方法。最后，讲讲uFVM中对扩散项离散过程的代码实现。

1 在矩形域上的二维扩散

在这里插入图片描述
如上图所示，将简单的矩形区域用笛卡尔网格剖分，并在该网格上对下面的稳态扩散方程进行离散
$-\nabla\cdot(\Gamma^\phi\nabla\phi)=Q^\phi$
其中 $\phi$ 代表一个标量（比如温度、湍动能等）， $Q^\phi$ 为单位体积 $\phi$ 的生成项（源项）， $\Gamma^\phi$ 为扩散系数。该方程可写成更加通用的形式，定义diffusion flux扩散通量 $\bold J^{\phi,D}$ 为
$\bold J^{\phi,D}=-\Gamma^\phi\nabla\phi$
则原方程变为
$\nabla\cdot\bold J^{\phi,D}=Q^\phi$
根据前面章节提到的离散方式和Gauss梯度定理，可得单元 $C$ 的离散形式
$\sum_{f \sim nb(C)}(-\Gamma^\phi\nabla\phi)_f \cdot \bold S_f = Q_C^\phi V_C$
上式的展开形式为
$(-\Gamma^\phi\nabla\phi)_e \cdot \bold S_e + (-\Gamma^\phi\nabla\phi)_w \cdot \bold S_w + (-\Gamma^\phi\nabla\phi)_n \cdot \bold S_n + (-\Gamma^\phi\nabla\phi)_s \cdot \bold S_s = Q_C^\phi V_C$
对于图中的均匀笛卡尔网格，单元的构成面的面积矢量为
$\begin{aligned} \bold S_e &= +(\Delta y)_e \bold i=||\bold S_e|| \bold i=S_e \bold i, \quad\quad \bold S_w =-(\Delta y)_w \bold i=-||\bold S_w|| \bold i=-S_w \bold i \\ \bold S_n &=+(\Delta x)_n \bold j=||\bold S_n|| \bold j=S_n \bold j, \quad\quad \bold S_s =-(\Delta x)_s \bold j=-||\bold S_s|| \bold j=-S_s \bold j \end{aligned}$
那么东边面的扩散通量为
$\begin{aligned} J^{\phi,D}_e &= -(\Gamma^\phi\nabla\phi)_e \cdot \bold S_e \\ &= -\Gamma^\phi_e S_e \left( \frac{\partial\phi}{\partial x} \bold i+ \frac{\partial\phi}{\partial y} \bold j \right)_e \cdot \bold i \\ & = -\Gamma^\phi_e (\Delta y)_e \left( \frac{\partial\phi}{\partial x} \right)_e \end{aligned}$
该扩散通量的离散形式可写成
$J^{\phi,D}_e = FluxT_e=FluxC_e\phi_C+FluxF_e\phi_E+FluxV_e$
为了确定系数 $FluxC_e,~FluxF_e,~FluxV_e$ ，需要给出变量\phi在该面两侧单元形心之间的变化规律，以便计算面上的梯度值。假设 $\phi$ 在两单元形心之间是线性变化的，如图
在这里插入图片描述
那么在面 $e$ 上沿着 $\bold i$ 方向的梯度可以写成
$\left( \frac{\partial\phi}{\partial x} \right)_e=\frac{\phi_E-\phi_C}{(\delta x)_e}$
那么面 $e$ 处扩散通量的最终离散形式为
$\begin{aligned} FluxT_e &=-\Gamma^\phi_e (\Delta y)_e\frac{\phi_E-\phi_C}{(\delta x)_e} \\ &=-\Gamma^\phi_e \frac{(\Delta y)_e}{(\delta x)_e} (\phi_E-\phi_C) \\ &=FluxC_e\phi_C+FluxF_e\phi_E+FluxV_e \end{aligned}$
令
$gDiff_e=\frac{(\Delta y)_e}{(\delta x)_e} = \frac{||\bold S_e||}{||\bold d_{CE}||} =\frac{S_e}{d_{CE}}$
其中 $\bold d_{CE}$ 为单元 $C$ 和 $E$ 形心之间的距离矢量，系数为
$\begin{aligned} FluxC_e &= \Gamma^\phi_e gDiff_e \\ FluxF_e &= - \Gamma^\phi_e gDiff_e \\ FluxV_e &= 0 \end{aligned}$
采用同样的方法，也可以推得w界面的扩散通量离散形式
$FluxT_w =FluxC_w\phi_C+FluxF_w\phi_W+FluxV_w$
其中
$\begin{aligned} FluxC_w &= \Gamma^\phi_w gDiff_w \\ FluxF_w &= - \Gamma^\phi_w gDiff_w \\ FluxV_w &= 0 \end{aligned}$
而
$gDiff_w=\frac{(\Delta y)_w}{(\delta x)_w} = \frac{||\bold S_w||}{||\bold d_{CW}||} =\frac{S_w}{d_{CW}}$
同样也可得到 $n$ 和 $s$ 面上的扩散通量离散形式
$\begin{aligned} FluxT_n &= FluxC_n\phi_C+FluxF_n\phi_N+FluxV_n \\ FluxT_s &= FluxC_s\phi_C+FluxF_s\phi_S+FluxV_s \end{aligned}$
其中
$\begin{aligned} FluxC_n &= \Gamma^\phi_n gDiff_n, \quad FluxF_n = - \Gamma^\phi_n gDiff_n, \quad & FluxV_n = 0 \\ FluxC_s &= \Gamma^\phi_s gDiff_s, \quad FluxF_s = - \Gamma^\phi_s gDiff_s, \quad & FluxV_s = 0 \end{aligned}$
而
$gDiff_n=\frac{(\Delta x)_n}{(\delta y)_n} = \frac{||\bold S_n||}{||\bold d_{CN}||} =\frac{S_n}{d_{CN}} \\ \quad \\ gDiff_s=\frac{(\Delta x)_s}{(\delta y)_s} = \frac{||\bold S_s||}{||\bold d_{CS}||} =\frac{S_s}{d_{CS}}$
将各个面的对流通量离散形式汇总代入到原方程的半离散格式，可得
$a_C\phi_C + a_E\phi_E + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=FluxF_e = - \Gamma^\phi_e gDiff_e \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_e+FluxC_w+FluxC_n+FluxC_s \\ &=-(a_E+a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-(FluxV_e+FluxV_w+FluxV_n+FluxV_s) \end{aligned}$
或者，写成紧致格式
$a_C\phi_C+\sum_{F\sim NB(C)} a_F \phi_F=b_C$
其中
$\begin{aligned} a_F &=FluxF_f = - \Gamma^\phi_f gDiff_f \\ a_C &=\sum_{f\sim nb(C)} FluxC_f \\ b_C &= Q_C^\phi V_C-\sum_{f\sim nb(C)} FluxV_f \end{aligned}$
其中的 $F$ 代表单元 $C$ 的邻近单元 $E, W, N, S$ ，小写 $f$ 代表单元 $C$ 的构成面（邻近面） $e, w, n, s$ 。

2 离散方程的解释

合适的离散方法得到的离散代数方程应该能反映原守恒方程的特性。所以推出的离散方程的系数需要满足如下两个准则。

2.1 零加和准则

即无源项的时候，各系数加和应为0。
$a_C+\sum_{F\sim NB(C)} a_F =0$
或
$a_C=-\sum_{F\sim NB(C)} a_F$
或
$\sum_{F\sim NB(C)} \frac{a_F}{a_C} =-1$

2.2 反号准则

即 $a_C$ 和 $a_F$ 是反号的，这样才能保证当 $\phi_F$ 增加/减小的时候， $\phi_C$ 是增加/减小的，从而保证了有界性boundedness。

3 边界条件

边界条件的添加至关重要，尽管控制方程相同，然而不同的边界条件所得到的解是不同的，对应的物理意义也不一样，添加错误的边界条件会导致无法得到解，或者得到的解与原物理问题毫不相干。

扩散方程，即，拉普拉斯方程的边界条件通常有Dirichlet、Von Neumann、mixed和symmetry，边界条件是施加在边界单元上的，它们有1个或多个面位于边界上，变量 $\phi$ 既存储在单元形心上，也存储在面形心上。

在这里插入图片描述
如图， $C$ 代表边界单元的形心，该单元有1个边界面， $b$ 为该边界面的形心， $\bold S_b$ 为该边界面的面积矢量。跟前面一样地，也可以推出单元 $C$ 的离散格式
$\sum_{f \sim nb(C)} (\bold J^{\phi,D} \cdot \bold S)_f = Q_C^\phi V_C$
内部面的通量还跟之前一样离散，而边界通量的离散则是要构建跟 $\phi_C$ 相关的线性格式
$\begin{aligned} \bold J^{\phi,D}_b \cdot \bold S_b &=FluxT_b \\ &= -\Gamma_b^\phi (\nabla \phi)_b \cdot \bold S_b \\ &= FluxC_b\phi_b+FluxV_b \end{aligned}$
边界条件的添加，要么给定边界值 $\phi_b$ ，要么是给出边界通量 $\bold J^{\phi,D}_b$ 。不同边界条件的添加方法如下。

3.1 Dirichlet边界条件

Dirichlet边界条件是直接给定边界上的变量值 $\phi$ ，即
$\phi_b=\phi_{specified}$
那么
$\begin{aligned} FluxT_b &= -\Gamma_b^\phi (\nabla \phi)_b \cdot \bold S_b \\ &= -\Gamma_b^\phi \frac{||\bold S_b||}{||\bold d_{Cb}||} (\phi_b - \phi_C) \\ &= FluxC_b\phi_C+FluxV_b \end{aligned}$
推得
$\begin{aligned} FluxC_b &= \Gamma_b^\phi gDiff_b=a_b \\ FluxV_b &= -\Gamma_b^\phi gDiff_b \phi_b = -a_b \phi_b \end{aligned}$
其中
$gDiff_b=\frac{S_b}{d_{Cb}}$
对于图中所示的 $C$ 单元， $a_E$ 系数变成0了（这个边界单元也没有 $E$ 单元，所以系数肯定是0了），其离散方程变为
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_b+\sum_{f\sim nb(C)}FluxC_f\\ &=FluxC_b+(FluxC_w+FluxC_n+FluxC_s) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$
注意：

系数 $a_b$ 比其它邻居系数都要大，因为 $b$ 距离 $C$ 更近，所以其对 $\phi_C$ 的影响更大。
系数 $a_C$ 仍旧是所有邻居系数的加和，包括 $a_b$ 。这也就是说，边界单元的 $\sum_{f\sim nb(C)}|a_F|/|a_C|<1$ ，不能保证离散方程满足零系数加和准则，有 $|a_C|>\sum_{f\sim nb(C)}|a_F|$ （严格对角占优了），从而在使用迭代方法求解线性方程的时候，是能获得收敛解的。
$a_b\phi_b$ 项产生了 $FluxV_b$ ，现在位于方程的右端项，是 $b_C$ 的一部分，因为其不含未知量。

3.2 Von Neumann边界条件

在这里插入图片描述
如上图，如果边界上给定的是 $\phi$ 的通量（或垂直于面的梯度值），那么这种边界条件就是Neumann边界条件。此时，这个给定的通量是
$-(\Gamma^\phi\nabla\phi)_b\cdot\bold i=q_b$
这样的话，通量 $\bold J_b^{\phi,D}$ 可得
$\bold J_b^{\phi,D}\cdot\bold S_b=-(\Gamma^\phi\nabla\phi)_b \cdot ||\bold S_b|| \bold i=q_b||\bold S_b||=FluxC_b\phi_C+FluxV_b$
其中
$\begin{aligned} FluxC_b &=0 \\ FluxV_b &=q_bS_b=q_b(\Delta y)_C \end{aligned}$
当与所使用的坐标系统方向相同的时候，通量分量为正值。

这样，单元 $C$ 的离散方程为：
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=\sum_{f\sim nb(C)}FluxC_f=-(a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$
注意：

Von Neumann边界条件生成的系数 $a_C$ 并不是占优的（也就意味着迭代解法是不靠谱的了）。
如果 $q_b$ 和 $S_C^\phi$ 为0，即边界不流入通量，且单元内部没有源项，则 $\phi_C$ 的值是完全由其周围邻近单元所限定的。而在其它情况下（ $q_b$ 或 $S_C^\phi$ 不0，即有边界通量流入，或者是单元内部有源项生成的时候）， $\phi_C$ 则可以超过（或是低于）周围邻近值 $\phi$ ，这是合情合理的。举个例子，如果 $\phi$ 是温度，则 $q_b$ 代表的就是施加在边界上的热流通量，那么如果边界上有热量流入的话，则靠近边界区域的温度是会比内部区域的温度要高的。
关于边界上 $\phi$ 值的计算，一旦单元 $\phi_C$ 算出来后，边界的 $\phi_b$ 可以用下面的式子来计算
$\phi_b=\frac{\Gamma_b^{\phi}gDiff_b\phi_C-q_b}{\Gamma_b^{\phi}gDiff_b}$
Von Neumann边界条件可以认为是FVM的自然边界条件，因为当某个边界面所给定的通量为0时，离散后跟该面相关的啥项也没有了。可是，Dirichlet边界条件中，即便边界面的给定 $\phi$ 为0，依然在离散方程中出现了相关项。

3.3 混合边界条件

在这里插入图片描述
混合边界条件如上图所示，边界条件是通过对流传热系数（ $h_\infin$ ）和 $\phi$ 的周围环境值 $\phi_\infin$ 来给出的
$\bold J_b^{\phi,D}\cdot\bold S_b=-(\Gamma^\phi\nabla\phi)_b \cdot ||\bold S_b|| \bold i = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C$
重新写成如下形式
$-\Gamma_b^\phi S_b \left( \frac{\phi_b - \phi_C}{\delta x_b} \right) = - h_\infin(\phi_\infin-\phi_b)S_b$
从这个方程可以解出 $\phi_b$
$\phi_b=\frac{h_\infin\phi_\infin + (\Gamma_b^\phi/\delta x_b)\phi_C}{h_\infin+(\Gamma_b^\phi/\delta x_b)}$
将 $\phi_b$ 代入到 $\bold J_b^{\phi,D}\cdot\bold S_b = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C$ 中，可得
$\begin{aligned} \bold J_b^{\phi,D}\cdot\bold S_b & = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C \\ & = - \underbrace{\left[ \frac{h_\infin (\Gamma_b^\phi/\delta x_b)}{h_\infin+(\Gamma_b^\phi/\delta x_b)} S_b\right]}_{R_{eq}}(\phi_\infin-\phi_C) \\ & = FluxC_b\phi_C + FluxV_b \end{aligned}$
其中
$\begin{aligned} FluxC_b &=R_{eq} \\ FluxV_b &=-R_{eq}\phi_{\infin} \end{aligned}$
那么，对于单元 $C$ 的离散方程，其最终形式为
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_b+\sum_{f\sim nb(C)}FluxC_f=FluxC_b+(FluxC_w+FluxC_n+FluxC_s) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$

3.4 对称边界条件

沿着对称边界条件，标量 $\phi$ 的法向通量为0，因此，对称边界条件等效于Neumann边界条件（将通量值设为0），即， $FluxC_b=FluxV_b=0$ ，那么直接把之前推出来的Neumann边界条件离散格式中的 $q_b$ 设成0就行了，即
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=\sum_{f\sim nb(C)}FluxC_f=-(a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-\sum_{f\sim nb(C)}FluxV_f \end{aligned}$

4 界面扩散系数

前面的离散方程中，还需要用到扩散系数 $\Gamma_e^\phi,~\Gamma_w^\phi,~\Gamma_n^\phi,~\Gamma_s^\phi$ ，它们分别代表着扩散系数 $\Gamma^\phi$ 在 $e, w, n, s$ 面处的值。当扩散系数 $\Gamma^\phi$ 随空间位置变化时（比如在可压缩流动问题中，粘性系数是温度的函数，随温度不同粘性也不相同），扩散系数同样也是存储在单元形心 $E, W, N, S . . .$ 上的，那么就需要去计算界面上的扩散系数值了。实际上，湍流粘性、热传导系数等都是非均匀分部的。当然，如果扩散系数整场是定值的话，就不用这么折腾了。

计算界面扩散系数最简单的方法莫过于直接用两侧单元值做线性插值了，即
$\Gamma_e^\phi=(1-g_e)\Gamma_C^\phi+g_e\Gamma_E^\phi$
其中插值因子 $g_e$ 是距离的比值，即
$g_e=\frac{d_{Ce}}{d_{Ce}+d_{eE}}$
如果是笛卡尔网格，且界面刚好位于两侧单元中间的话， $g_e=0.5$ ，那么 $\Gamma_e^\phi=(\Gamma_C^\phi+\Gamma_E^\phi)/2$ 。同样也可以定义其它插值系数
$g_w=\frac{d_{Cw}}{d_{Ce}+d_{wW}} \\ \quad \\ g_n=\frac{d_{Cn}}{d_{Cn}+d_{nN}} \\ \quad \\ g_s=\frac{d_{Cs}}{d_{Cs}+d_{sS}}$
然而，在某些情况下，这种处理方法会导致错误，比如当复合材料中的传导突变的情况下。幸运的是，有个更好的且相对简单易于实现的方法。在该方法中，界面的局部传导系数不是重点，它主要关注的是如何获得界面上的扩散通量 $\bold J^{\phi,D}$ 。
在这里插入图片描述
如图所示的1维问题，假设单元 $C$ 是由导热系数 $\Gamma_C^\phi$ 的材料构成，而单元 $E$ 是由导热系数为 $\Gamma_E^\phi$ 的材料构成，则对于C和E之间的非均匀分片，稳态1维分析（不含源项）有（界面e的任何一边的通量假设是相同的）
$\begin{aligned} \bold J_e^{\phi,D}\cdot \bold S_e &=\frac{\phi_C-\phi_e}{\displaystyle \frac{(\delta x)_{Ce}}{\Gamma_C^\phi}} = \frac{\phi_e-\phi_E}{\displaystyle \frac{(\delta x)_{eE}}{\Gamma_E^\phi}} \\ &=\frac{\phi_C-\phi_E}{\displaystyle \frac{(\delta x)_{Ce}}{\Gamma_C^\phi} + \frac{(\delta x)_{eE}}{\Gamma_E^\phi}} = \frac{\phi_C-\phi_E}{\displaystyle \frac{(\delta x)_{CE}}{\Gamma_e^\phi}} \end{aligned}$
于是推得了该片的导热系数
$\frac{(\delta x)_{CE}}{\Gamma_e^\phi} = \frac{(\delta x)_{Ce}}{\Gamma_C^\phi} + \frac{(\delta x)_{eE}}{\Gamma_E^\phi} \Rightarrow \frac{1}{\Gamma_e^\phi} = \frac{1-g_e}{\Gamma_C^\phi} + \frac{g_e}{\Gamma_E^\phi}$
当界面刚好位于 $C$ 和 $E$ 的中间时， $g_e=0.5$ ，上式变为
$\Gamma_e^\phi=\frac{2\Gamma_C^\phi\Gamma_E^\phi}{\Gamma_C^\phi+\Gamma_E^\phi}$
这是 $\Gamma_C^\phi$ 和 $\Gamma_E^\phi$ 的调和平均（harmonic mean），而非算术平均。

需要注意的是，对于不连续扩散系数的调和平均插值方法，仅对1维扩散问题是精确的。然而，其对高维问题的应用也很有优势，比如流体和固体区域可以不加区分，当成一个区域来计算传热系数。

例1
在这里插入图片描述

如图热传导问题，2维矩形区域由两种材料构成，如下偏微分方程为其控制方程
$\nabla \cdot (k \nabla T) = 0$
其中 $T$ 代表温度，热传导系数 $k$ 和边界条件如图所示
a.推导出所有单元的代数方程
b. 使用Gauss-Seidel迭代方法求解系统方程获得单元 $T$ 值
c. 计算底部、顶部、右侧的 $T$ 值
d. 计算通过底部、顶部、左侧边界的净热流通量，并检查能量守恒是否满足。

解
求解过程较为繁琐，不再详述。

5 非笛卡尔正交网格

在这里插入图片描述
如图，网格还是正交的，但是它和 $x$ 与 $y$ 轴不是对齐的，这样的网格相当于是把笛卡尔网格转了个角度。

这种网格上的离散方程和笛卡尔网格上的离散方程式完全一样的，边界条件也是一样的，解法也是一样的。

比如
$\nabla\cdot\bold J^{\phi,D}=Q^\phi$
半离散形式
$\sum_{f \sim nb(C)}(-\Gamma^\phi\nabla\phi)_f \cdot \bold S_f = Q_C^\phi V_C$
考虑到界面 $e$ ，有
$J^{\phi,D}_e = -\Gamma^\phi_e (\nabla\phi \cdot \bold n)_e S_e = -\Gamma^\phi_e \left( \frac{\partial\phi}{\partial n} \right)_e S_e$
其中
$(\nabla\phi \cdot \bold n)_e = \left( \frac{\partial\phi}{\partial n} \right)_e$
为变量 $\phi$ 在界面 $e$ 上沿着 $\bold n$ 方向的梯度，再次对 $\phi$ 使用线性假设，该梯度可得
$\left( \frac{\partial\phi}{\partial n} \right)_e=\frac{\phi_E-\phi_C}{d_{CE}}$
其它项也将得出和笛卡尔网格同样的离散方程。

6 非正交非结构网格

6.1 非正交

在这里插入图片描述

前面所讲的情况下，通量是与面垂直的。然而，在一般情况下，结构化曲线网格和非结构网格是非正交的，即，面矢量 $\bold S_f$ 和连接面两侧单元形心的矢量 $\bold {CF}$ 是不共线的，如上图。这样，垂直于面的梯度就无法写成 $\phi_F$ 和 $\phi_C$ 的函数，而它实际上有了个垂直于 $\bold {CF}$ 的分量。

这样，在正交网格上梯度沿着界面法向为
$(\nabla\phi\cdot\bold n)_f=\left(\frac{\partial\phi}{\partial n}\right)_f=\frac{\phi_F-\phi_C}{||\bold r_F - \bold r_C||}=\frac{\phi_F-\phi_C}{d_{CF}}$
因为 $\bold {CF}$ 和 $\bold n$ （面单位法矢量）是共线的。而在非结构网格上，上面这个式子所表示的梯度方向是包含 $\phi_F$ 和 $\phi_C$ 且沿着连接点 $C$ 和 $F$ 的连线方向的（跟面法向的梯度是不同的，它仅仅是真正的面法向梯度的一个分量而已）。

如果 $\bold e$ 代表沿着连接节点 $C$ 和 $F$ 单位矢量，那么
$\bold e=\frac{\bold r_F - \bold r_C}{||\bold r_F - \bold r_C||}=\frac{\bold d_{CF}}{d_{CF}}$
这样，在 $\bold e$ 方向的梯度可以写成
$(\nabla\phi\cdot\bold e)_f = \left(\frac{\partial\phi}{\partial e}\right)_f=\frac{\phi_F-\phi_C}{||\bold r_F - \bold r_C||}=\frac{\phi_F-\phi_C}{d_{CF}}$
这样为了获得非正交网格上的线性化通量，面矢量 $\bold S_f$ 应该写成两个矢量 $\bold E_f$ 和 $\bold T_f$ 的加和形式，即
$\bold S_f = \bold E_f + \bold T_f$
其中 $\bold E_f$ 是沿着 $\bold {CF}$ 方向的，以保证扩散通量中的一部分可以写成节点值 $\phi_C$ 和 $\phi_F$ 的函数，即
$\begin{aligned} (\nabla\phi)_f\cdot\bold S_f &= \underbrace{(\nabla\phi)_f\cdot\bold E_f}_{orthogonal-like~contributation} + \underbrace{(\nabla\phi)_f\cdot\bold T}_{non-orthogonal~like~contributation} \\ &=E_f\left(\frac{\partial\phi}{\partial e}\right)_f + (\nabla\phi)_f\cdot\bold T_f \\ &=E_f\frac{\phi_F-\phi_C}{d_{CF}} + (\nabla\phi)_f\cdot\bold T_f \end{aligned}$
上式中右侧第一项与正交网格上的式子一样，其代表了正交贡献量。而第2项则是由于网格的非正交性所产生的交叉扩散或非正交扩散。对 $\bold S_f$ 的分解可以有不同的选择，如下所述。

6.2 最小修正方法Minimum Correction Approach

在这里插入图片描述
如图所示， $\bold S_f$ 的分解方法是要保证非正交修正尽可能的小，即让 $\bold E_f$ 和 $\bold T_f$ 正交。伴随着非正交修正的增加，从 $\phi_F$ 到 $\phi_C$ 的扩散贡献将减小。矢量 $\bold E_f$ 是这样计算的
$\bold E_f=(\bold e \cdot \bold S_f) \bold e=(S_f \cos \theta)\bold e$

6.3 正交修正方法Orthogonal Correction Approach

在这里插入图片描述这种方法让包含 $\phi_F$ 到 $\phi_C$ 的项和它们在正交网格上的贡献是同样的，即， $\bold S_f$ 和 $\bold E_f$ 的幅值是一样的，那么矢量 $\bold E_f$ 选择为
$\bold E_f=S_f\bold e$

6.4 过松弛方法 Over-Relaxed Approach

在这里插入图片描述
该方法中，包含 $\phi_F$ 到 $\phi_C$ 的项将随着非正交性的增加而增加，如上图所示。通过将 $\bold T_f$ 选择为垂直于 $\bold S_f$ 来实现，这样 $\bold E_f$ 是这样计算的
$\bold E_f=\left( \frac{S_f}{\cos \theta} \right)\bold e=\left( \frac{S_f^2}{S_f \cos \theta} \right) \bold e = \frac{\bold S_f \cdot \bold S_f}{\bold e \cdot \bold S_f}\bold e$
总结一下，在非正交网格上的单元面的扩散通量不能简单地写成该面两侧单元值的形式，必须添加上一个非正交修正项。该修正项在文献中常称为“交叉扩散”，其算法如下
$\begin{aligned} (\nabla\phi)_f \cdot \bold T_f &= (\nabla\phi)_f \cdot (\bold S_f - \bold E_f) \\ & = \begin{cases} (\nabla\phi)_f \cdot (\bold n - \cos \theta \bold e)S_f & minimum~correction \\ (\nabla\phi)_f \cdot (\bold n - \bold e)S_f & normal~correction \\ (\nabla\phi)_f \cdot (\bold n - \displaystyle \frac{1}{\cos \theta} \bold e)S_f & over-relaxed \end{cases} \end{aligned}$
对于正交网格而言， $\bold n$ 和 $\bold e$ 是共线的，角度 $\theta$ 是0，因此交叉扩散项是0。当交叉扩散项非0时，因为其不能写成 $\phi_F$ 和 $\phi_C$ 的函数，所以它将作为源项添加在单元的代数方程中。
上面的方法都是正确的，它们的差别是在非正交网格上的准确性和稳定性（accuracy和stability）。其中over-relaxed方法已经被发现是最稳定的，即便是网格高度非正交的时候。然而，扩散项离散后的最终通有形式，对于三种方法来说都是一样的。

6.5 交叉扩散项的处理方法

交叉扩散项无法展开成节点值的形式，因为其含有 $(\nabla\phi)_f$ ，所以，它只好作为延迟修正方式来处理，通过使用当前梯度场来算出交叉扩散项的值，并将其作为源项添加到代数方程的右端项中去。如下所述，梯度是在主网格（单元形心）上计算的，即 $(\nabla\phi)_C$ ，并且插值到界面上去的，即 $(\nabla\phi)_f$ 。（注意，这里的交叉扩散项指的是 $(\nabla\phi)_f\cdot\bold T_f$ ，即该项该如何计算，而在第9章第4节中所讲的是 $(\nabla\phi)_f$ 在非正交网格上式如何计算的，两者是不一样的！）

6.6 梯度计算

在正交网格上，单元形心的梯度值可以直接由单元两侧单元的形心值给出，面形心的梯度值也可由其两侧单元的形心值给出。然而，在非正交网格和非结构网格上，梯度的计算较为复杂，常用的方法是Green-Gauss或梯度定理来计算，即对于任何封闭体积 $V$ ，其由面 $\partial V$ 所包围而成，则有
$\int_V\nabla\phi~\text dV=\oint_{\partial V}\phi \text d \bold S$
其中 $\text d\bold S$ 为指向外侧的面矢量。为获取上式的离散形式，使用中值定理，得
$\int_V\nabla\phi~\text dV=\overline{\nabla\phi}V$
于是，在单元 $C$ 上的平均梯度为
$\overline{\nabla\phi}_C = \frac{1}{V_C}\oint_{\partial V_C}\phi_f \bold S_f$
接下来，对于单元面的积分可近似为面形心值与面面积的乘积，这样 $\overline{\nabla\phi}_C$ ，或简写成 $\nabla\phi_C$ ，可由下式计算
${\nabla\phi}_C = \frac{1}{V_C}\sum_{f\sim nb(C)}\phi_f \bold S_f$
算得单元形心的梯度后，面形心的梯度可由面两侧单元形心梯度的加权平均来获取，即
$\nabla\phi_f=g_C\nabla\phi_C+g_F\nabla\phi_E$
其中的 $g_C$ 和 $g_F$ 为几何插值系数，其与单元 $C$ 和单元 $F$ 以及面 $f$ 的位置相关。

6.7 非正交网格的代数方程

将面矢量 $\bold S_f$ 分解为 $\bold E_f$ 和 $\bold T_f$ 矢量，并将其等效表达式代入到半离散扩散通量方程中去，可得
$\begin{aligned} \sum_{f\sim nb(C)}\left( \bold J_f^{\phi,D} \cdot \bold S_f\right) &= \sum_{f\sim nb(C)} \left( -(\Gamma^\phi\nabla\phi)_f \cdot (\bold E_f+\bold T_f) \right) \\ &= \underbrace{\sum_{f\sim nb(C)} \left( -(\Gamma^\phi\nabla\phi)_f \cdot \bold E_f \right)}_{Orthogonal~Linearizeable~Part} + \underbrace{\sum_{f\sim nb(C)} \left( -(\Gamma^\phi\nabla\phi)_f \cdot \bold T_f \right)}_{Non-orthogonal~Non-Linearizeable~Part} \\ &= \sum_{f\sim nb(C)} \left( -\Gamma^\phi_f E_f \frac{\phi_F-\phi_C}{d_{CF}} \right) + \sum_{f\sim nb(C)} \left( -(\Gamma^\phi\nabla\phi)_f \cdot \bold T_f \right) \\ &= \sum_{f\sim nb(C)} \Gamma^\phi_f gDiff_f (\phi_C-\phi_F) + \sum_{f\sim nb(C)} \left( -(\Gamma^\phi\nabla\phi)_f \cdot \bold T_f \right) \\ &=\left( \sum_{f\sim nb(C)} FluxC_f \right)\phi_C+ \sum_{f\sim nb(C)} (FluxF_f\phi_F) + \sum_{f\sim nb(C)} FluxV_f \end{aligned}$
其中 $gDiff_f$ 为几何扩散系数，即
$gDiff_f=\frac{E_f}{d_{CF}}$
使用上述的扩散通量形式，在非结构/结构的非正交网格上的扩散方程的最终离散形式为
$a_C\phi_C+\sum_{f\sim nb(C)}a_F\phi_F=b_C$
其中
$\begin{aligned} & a_F=FluxFf=-\Gamma_f^\phi gDiff_f \\ & a_C=\sum_{f\sim nb(C)}FluxC_f=-\sum_{f\sim nb(C)}FluxF_f=\sum_{f\sim nb(C)}\Gamma_f^\phi gDiff_f \\ & b_C=Q_C^\phi V_C - \sum_{f\sim nb(C)}FluxV_f=Q_C^\phi V_C+ \sum_{f\sim nb(C)}\left( (\Gamma^\phi\nabla\phi)_f \cdot \bold T_f \right) \end{aligned}$
注意非正交项的符号有所变化，因为其作为源项整到了方程的右端项中去。

例2
对于如图所示的多边形单元 $C$ 和其邻居单元 $F$ ，在任何点处的解满足 $\phi=x^2+y^2+x^2y^2$ ，如果单元 $C$ 的体积是 $V_C=8.625$ ，那么计算：

单元 $C$ 的梯度值，即 $\nabla\phi_C$ ，使用数值方法求解，并和解析解对比；
解析解 $\nabla\phi_F$ ；
用数值解 $\nabla\phi_C$ 和解析解 $\nabla\phi_F$ 插值计算 $\nabla\phi_{f_1}$ 的近似解，并与解析解 $\nabla\phi_{f_1}$ 对比；
将 $\nabla\phi_{f_1}\cdot\bold S_{f_1}$ 展开为 $\phi_C$ 和 $\phi_F$ 的形式，使用
a. 最小修正方法Minimum Correction Approach
b. 正交修正方法Orthogonal Correction Approach
c. 过松弛方法 Over-Relaxed Approach

在这里插入图片描述

解
较为繁琐，省去不表，感兴趣的话可自行计算。

6.8 非正交网格的边界条件

非正交网格的边界条件处理与正交网格是非常类似的，只是在非正交扩散的处理上略有不同。

6.8.1 Dirchlet 边界条件

在这里插入图片描述
在Direchlet边界条件上，变量值 $\phi$ 是已知的，与在结构网格上的离散过程类似，只是需要考虑交叉扩散，只不过这时是把边界面当成内部面来处理就好了。当面矢量与连接边界面紧邻边界单元的形心与该边界面的形心的矢量不共线的时候，需要考虑交叉扩散。扩撒通量沿着边界面可被离散为
$\begin{aligned} \bold J_b^{\phi,D} \cdot \bold S_b &= -(\Gamma^\phi\nabla\phi)_b \cdot \bold S_b = -(\Gamma^\phi\nabla\phi)_b \cdot (\bold E_b+\bold T_b) \\ &= -\Gamma^\phi_b \left(\frac{\phi_b-\phi_C}{d_{Cb}}\right) E_b - \Gamma^\phi_b(\nabla\phi)_b \cdot \bold T_b \\ & = FluxC_b\phi_C + FluxV_b \end{aligned}$
其中
$\begin{aligned} FluxC_b &= \Gamma_b^\phi gDiff_b \\ FluxV_b &= -\Gamma_b^\phi gDiff_b \phi_b - \Gamma^\phi_b(\nabla\phi)_b \cdot \bold T_b \\ gDiff_b &= \frac{E_b}{d_{Cb}} \end{aligned}$
把上面式子代入到单元的离散方程中去，可得各项修正的系数为
$\begin{aligned} a_F&=FluxF_f \quad\quad a_C=\sum_{f\sim nb(C)}FluxC_f+FluxC_b \\ b_C&=Q_C^\phi V_C - FluxV_b - \sum_{f\sim nb(C)}FluxV_f \end{aligned}$

6.8.2 Neumann 边界条件

Neumann边界条件对于非结构网格的处理与结构网格是一样的，边界上给定的通量将作为源项添加。边界的代数方程是一样的，修正系数也是一样的，并没有任何差别和改变。

6.8.3 Mixed边界条件

对于混合边界条件，其由对流导热系数 $h_\infin$ 和周围变量值 $\phi_\infin$ 来定义，扩散通量可写成
$\begin{aligned} \bold J_b^{\phi,D}\cdot\bold S_b &=-\Gamma^\phi_b\left( \frac{\phi_b - \phi_C}{d_{Cb}} \right) E_b - \Gamma_b^\phi (\nabla \phi)_b \cdot \bold T_b\\ &= -h_\infin(\phi_\infin-\phi_b)S_b \end{aligned}$
可得到 $\phi_b$ 的表达式如下
$\phi_b=\frac {h_\infin S_b \phi_\infin + \displaystyle\frac{\Gamma_b^\phi E_b}{d_{Cb}}\phi_C - \Gamma_b^\phi (\nabla \phi)_b \cdot \bold T_b} {h_\infin S_b + \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}}$
将上式待会到上上式中，可得
$\begin{aligned} \bold J_b^{\phi,D}\cdot\bold S_b &= \underbrace{-\left[ \frac{h_\infin S_b \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}}{h_\infin S_b + \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}} \right] }_{a_b} (\phi_\infin - \phi_C) - \underbrace{\frac{h_\infin S_b \Gamma_b^\phi(\nabla\phi)_b\cdot\bold T_b}{h_\infin S_b + \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}}}_{S_b^{\phi,CD}} \\ &= FluxC_b \phi_C + FluxV_b \end{aligned}$
其中
$\begin{aligned} FluxC_b &= \frac{h_\infin S_b \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}}{h_\infin S_b + \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}} \\ FluxV_b &= -FluxC_b\phi_\infin - \frac{h_\infin S_b \Gamma_b^\phi(\nabla\phi)_b\cdot\bold T_b}{h_\infin S_b + \displaystyle \frac{\Gamma_b^\phi E_b}{d_{Cb}}} \end{aligned}$
对于边界单元的离散方程的修正系数为
$\begin{aligned} a_F&=FluxF_f \quad\quad a_C=\sum_{f\sim nb(C)}FluxC_f+FluxC_b \\ b_C&=Q_C^\phi V_C - FluxV_b - \sum_{f\sim nb(C)}FluxV_f \end{aligned}$

7 Skewness（偏斜）

构成通用离散方程的很多项在计算的时候都需要，单元面上的 $\phi$ 值。可以用面上的平均值来代替面上的值，那么如果面上各处的 $\phi$ 值不同的话，变量 $\phi$ 在面上的平均值就是位于面形心处的 $\phi$ 值了。

然而一般是用面两侧单元形心的值来做线性插值获取面上的值，但是这样子获得的是面两侧单元形心连线与面交点处的值，如下图所示，这个交点是 $f^{'}$ ，这个点跟面形心 $f$ 点往往是不重合的（除非是在正交网格上面两点才重合），也就是说获取的面上点的值，并非是面形心处的值！！

在这里插入图片描述
为了保证离散方法的2阶精度，所有面积分的积分点应在面形心 $f$ 处，这就需要采取修正措施将 $f^{'}$ 处的值修正到 $f$ 处，一般用Taylor展开来实现，即
$\phi_f=\phi_{f'}+(\nabla\phi)_{f'}\cdot\bold d_{f'f}$
其中 $\bold d_{f'f}$ 为交点 $f^{'}$ 到面形心点 $f$ 的距离矢量。

这里只是简要提及一下而已，关于该部分修正的详细内容可参考第9章。

8 各向异性（非各向同性）扩散

目前的扩散方程中都是假设材料对于 $\phi$ 的传导是没有方向性的，即对于所有的方向扩散系数都是一样的，即，材料是各向同性的。但是，对于某些材料而言，扩散系数是跟方向相关的，扩散是各向异性的！某些固体是各向异性的，其半离散扩散方程为
$\sum_{f\sim nb(C)}(-\boldsymbol \kappa^\phi \cdot \nabla \phi )_f\cdot\bold S_f = Q_C^\phi V_C$
其中 $\boldsymbol \kappa^\phi$ 为2阶对称张量，对于最通常的三维情况，左端项经过一些数学处理，可写成
$(-\boldsymbol \kappa^\phi \cdot \nabla \phi )_f\cdot\bold S_f = -\begin{bmatrix} \boldsymbol \kappa^\phi_{11} & \boldsymbol \kappa^\phi_{12} & \boldsymbol \kappa^\phi_{13} \\ ~\\ \boldsymbol \kappa^\phi_{21} & \boldsymbol \kappa^\phi_{22} & \boldsymbol \kappa^\phi_{23} \\ ~\\ \boldsymbol \kappa^\phi_{31} & \boldsymbol \kappa^\phi_{32} & \boldsymbol \kappa^\phi_{33} \end{bmatrix}_f \begin{bmatrix} \displaystyle \frac{\partial\phi}{\partial x} \\ ~\\ \displaystyle \frac{\partial\phi}{\partial y} \\ ~\\ \displaystyle \frac{\partial\phi}{\partial z} \end{bmatrix}_f \cdot \bold S_f$
经过矩阵相乘，上式变为
$(-\boldsymbol \kappa^\phi \cdot \nabla \phi )_f\cdot\bold S_f = -\begin{bmatrix} \boldsymbol \kappa^\phi_{11}\displaystyle \frac{\partial\phi}{\partial x} + \boldsymbol \kappa^\phi_{12} \displaystyle \frac{\partial\phi}{\partial y} + \boldsymbol \kappa^\phi_{13} \displaystyle \frac{\partial\phi}{\partial z} \\ ~\\ \boldsymbol \kappa^\phi_{21} \displaystyle \frac{\partial\phi}{\partial x} + \boldsymbol \kappa^\phi_{22}\displaystyle \frac{\partial\phi}{\partial y} + \boldsymbol \kappa^\phi_{23} \displaystyle \frac{\partial\phi}{\partial z}\\ ~\\ \boldsymbol \kappa^\phi_{31} \displaystyle \frac{\partial\phi}{\partial x} + \boldsymbol \kappa^\phi_{32}\displaystyle \frac{\partial\phi}{\partial y} + \boldsymbol \kappa^\phi_{33}\displaystyle \frac{\partial\phi}{\partial z} \end{bmatrix}_f \cdot \begin{bmatrix} S^x \\ ~\\ S^y \\ ~\\ S^z \end{bmatrix}_f$
进一步处理得到
$\begin{aligned} (-\boldsymbol \kappa^\phi \cdot \nabla \phi )_f\cdot\bold S_f &= -\begin{bmatrix} \displaystyle \frac{\partial\phi}{\partial x} & \displaystyle \frac{\partial\phi}{\partial y} & \displaystyle \frac{\partial\phi}{\partial z} \end{bmatrix}_f \begin{bmatrix} \boldsymbol \kappa^\phi_{11} & \boldsymbol \kappa^\phi_{21} & \boldsymbol \kappa^\phi_{31} \\ ~\\ \boldsymbol \kappa^\phi_{12} & \boldsymbol \kappa^\phi_{22} & \boldsymbol \kappa^\phi_{32} \\ ~\\ \boldsymbol \kappa^\phi_{13} & \boldsymbol \kappa^\phi_{23} & \boldsymbol \kappa^\phi_{33} \end{bmatrix}_f \cdot \begin{bmatrix} S^x \\ ~\\ S^y \\ ~\\ S^z \end{bmatrix}_f \\ &= -(\nabla\phi)_f\cdot[(\boldsymbol\kappa^\phi)^T\cdot S]_f=-(\nabla\phi)_f\cdot\bold S'_f \end{aligned}$
于是，扩散方程的离散格式变成
$\sum_{f\sim nb(C)}-(\nabla\phi)_f\cdot\bold S'_f = Q_C^\phi V_C$
通过以上分析发现，只需要在离散过程中把原本的扩散系数 $\Gamma$ 替换为 $1$ ，把面矢量 $\bold S_f$ 替换为 $\bold S'_f$ ，而其余的部分是完全相同的，就可以求解各向同性和各向异性扩散问题了。

9 迭代求解问题的欠松弛处理

对于更加通常的扩散问题， $\Gamma^\phi$ 可能是未知量 $\phi$ 的函数，网格可能高度的非正交，有很大的交叉扩散项，这时就需要做延迟修正处理了。因为 $\phi$ 在迭代中的较大变化会导致较大的源项和较大的系数变化，这些会导致迭代迭代求解过程的发散。该发散是由于系数和交叉扩散项的非线性效应引起的，即，源项被当前流场影响较大而无法收敛。为了获得收敛解，需要减慢 $\phi$ 在迭代过程中的变化，这种技术就是欠松弛。有许多方法来引入欠松弛，这里介绍一种，其余的可参考第14章内容。在通有离散方程上进行处理
$a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=b_C$
上述方程可以写成
$\phi_C=\frac{-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F+b_C}{a_C}$
让 $\phi_C^*$ 代表上次迭代的 $\phi_C$ 值，将其从右端项中减去并加上，即
$\phi_C=\phi_C^* + \left( \frac{-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F+b_C}{a_C} - \phi_C^* \right)$
上式中括号里的项代表了当前迭代步对 $\phi_C$ 的改变，该变化可以通过引入松弛因子 $\lambda^\phi$ 来修正，即
$\phi_C=\phi_C^* + \lambda^\phi \left( \frac{-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F+b_C}{a_C} - \phi_C^* \right)$
或
$\frac{a_C}{\lambda^\phi}\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=b_C+\frac{(1-\lambda^\phi)a_C}{\lambda^\phi}\phi_C^*$
即，收敛的时候， $\phi_C$ 和 $\phi_C^*$ 是一样的。

根据 $\lambda^\phi$ 的不同，方程可以是欠松弛（ $0<\lambda^\phi<1$ ）的，也可以是超松弛（ $\lambda^\phi>1$ ）的。在CFD中，一般常用的都是欠松弛的处理方式。 $\lambda^\phi$ 接近1意味着轻微的欠松弛，而其接近0则意味着较为强烈的欠松弛，即 $\phi_C$ 每次迭代的变化非常小。

最适宜的欠松弛因子是因问题而定的，并没有通用的规则可寻。影响松弛因子的因素包括求解问题的类型、系统方程的规模（计算域中网格点数目）、网格尺度、网格增长率、采用的迭代方法等。通常来说， $\lambda^\phi$ 是凭借之前计算中所积攒的经验来选取的。而且，没必要再整个计算域中采用同样的松弛因子，也没必要在每次迭代时采用同样的松弛因子。（这话说的，跟没说也差不多了，总之就是不收敛了你就调小点，收敛慢了你就调大点，大概就是这么个意思了）

加了松弛因子的离散方程也能归并成统一的格式，即，让 $a_C$ 和 $b_C$ 项做替换就可以了
$a_C \leftarrow \frac{a_C}{\lambda^\phi}$
$b_C \leftarrow b_C+\frac{(1-\lambda^\phi)a_C}{\lambda^\phi}\phi_C^*$
实际上，松弛处理的方法是求解非线性问题中常用的一种稳定求解的方法。

10 代码讲解

10.1 uFVM

在uFVM中，对于内部面的扩散项离散的实现是由函数cfdAssembleDiffusionTermInterior来完成的，代码和个人理解展示如下

% 将gamma从单元形心值插值到面形心值
gamma_f = cfdInterpolateFromElementsToFaces('Average',gamma);
gamma_f = [gamma_f(iFaces)];

% 提取面几何扩散系数 gDiff_f = Ef / CF
geoDiff_f = [theMesh.faces(iFaces).geoDiff]';
% 提取面面积矢量Sf
Sf = [theMesh.faces(iFaces).Sf]';
% 提取面的非正交修正矢量Tf = Sf - CF
Tf = [theMesh.faces(iFaces).T]';

% 面左单元标识
iOwners = [theMesh.faces(iFaces).iOwner]';
% 面右单元标识
iNeighbours = [theMesh.faces(iFaces).iNeighbour]';

% 与8.6.7（第8章第6.7节）中的离散代数方程形式相对应
% 面扩散通量FluxCf，加和后为phiC的系数，正值，也是该面的own单元系数
theFluxes.FLUXC1f(iFaces,1) = gamma_f .* gDiff_f;
% 面扩散通量FluxFf，phiF的系数，负值，也是该面的neighbour单元系数
theFluxes.FLUXC2f(iFaces,1) = -gamma_f .* gDiff_f;

% 非正交扩散项FluxVf，加和后为源项bC的第2项
theFluxes.FLUXVf(iFaces,1) = gamma_f .* dot(grad_f(:,:)',Tf(:,:)')';
% 经过这个面的所有通量加和，即，经过这个面的净通量
theFluxes.FLUXTf(iFaces,1) = theFluxes.FLUXC1f(iFaces) .* ...
    phi(iOwners) + theFluxes.FLUXC2f(iFaces) .* phi(iNeighbours) + ...
    theFluxes.FLUXVf(iFaces);

注意，上述代码并没有做面循环（by the way，matlab里面做循环的话效率会低到让人想撞墙），这里是整场来做运算的，把全部向量一起计算了。此外，面几何扩散系数 $gDiff_f = Ef / CF$ 的计算是在函数cfdProcessOpenFoamMesh.m函数中进行的，相关代码如下

% 面右单元形心F到左单元形心C的矢量CF
theMesh.faces(iFace).dCF = element2.centroid - element1.centroid;

% CF矢量单位化，的eCF矢量
eCF = dCF/cfdMagnitude(dCF);

% 在矢量eCF上整一个跟面面积一样的长度矢量E
E = theFace.area*eCF;

% 面的几何扩散系数gDiff=Ef/dCF
theMesh.faces(iFace).gDiff = cfdMagnitude(E)/cfdMagnitude(dCF);

% 非正交扩散矢量T = 面积矢量Sf - 正交矢量E(Sf在CF方向的分量)
theMesh.faces(iFace).T = theFace.Sf - E;

% 注意，这里的非正交修正用的是Orthogonal Correction Approach
% 即，让Ef和Sf的幅值相同，即正交项和正交网格上的值相同。

边界条件对方程的影响在边界单元上应充分考虑，下面给出Dirichlet边界条件的情况下扩散项在边界面的处理代码

% 取出网格、单元总数、内部面总数
theMesh = cfdGetMesh;
numberOfElements = theMesh.numberOfElements;
numberOfInteriorFaces = theMesh.numberOfInteriorFaces;

% 取出第iPatch个边界Patch，及其所含边界面数目
theBoundary = theMesh.boundaries(iPatch);
numberOfBFaces = theBoundary.numberOfBFaces;

% 该片边界patch的起始、结束、范围面标识信息
iFaceStart = theBoundary.startFace;
iFaceEnd = iFaceStart+numberOfBFaces-1;
iBFaces = iFaceStart:iFaceEnd;

% 该片边界对应的起始单元、结束单元、单元范围标识信息
iElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces;
iElementEnd = iElementStart+numberOfBFaces-1;
iBElements = iElementStart:iElementEnd;

% 取出面的几何扩散系数、非正交修正矢量、面的所属单元标识（边界面对应一个边界单元哈）
geodiff = [theMesh.faces(iBFaces).geoDiff]';
Tf = [theMesh.faces(iBFaces).T]';
iOwners = [theMesh.faces(iBFaces).iOwner]';

% 计算扩散方程离散形式的相关系数
% 计算边界条件在 aC系数中的添加项 FluxCb = gamma * gDiff
theFluxes.FLUXC1f(iBFaces) = gamma(iBElements).*geodiff;
% 计算边界条件在 bC系数中的添加项 FluxVb1 = - FluxCb = - gamma * gDiff
theFluxes.FLUXC2f(iBFaces) = - gamma(iBElements).*geodiff;
% 计算边界条件在 bCx系数中的添加项 FluxVb2 = - gamma * (grad(phi)) dot T
theFluxes.FLUXVf(iBFaces) = - ...
    gamma(iBElements).*dot(grad(iBElements,:)',Tf(:,:)')';

其它边界条件的添加方式可参考cfAssembleDiffusionTerm.m函数。

当这些线性系数对于内部面和边界面都算好之后，就可以把它们组装成整体（稀疏）矩阵了。这部分操作是在cfdAssembleIntoGlobalMatrixFaceFluxes函数中进行的，内部面系数组装成左端项LHS矩阵和右端项RHS矢量的代码展示如下：

%
% Assemble fluxes of interior faces
%
for iFace = 1:numberOfInteriorFaces
    theFace = theMesh.faces(iFace);
    iOwner = theFace.iOwner;
    iOwnerNeighbourCoef = theFace.iOwnerNeighbourCoef;
    iNeighbour = theFace.iNeighbour;
    iNeighbourOwnerCoef = theFace.iNeighbourOwnerCoef;
    %
    % assemble fluxes for owner cell
    ac(iOwner) = ac(iOwner) ...
        + vf_f(iFace)*theFluxes.FLUXC1f(iFace);
    anb{iOwner}(iOwnerNeighbourCoef) = anb{iOwner}(iOwnerNeighbourCoef) ...
        + vf_f(iFace)*theFluxes.FLUXC2f(iFace);
    bc(iOwner) = bc(iOwner) ...
        - vf_f(iFace)*theFluxes.FLUXTf(iFace);
    %
    % assemble fluxes for neighbour cell
    ac(iNeighbour) = ac(iNeighbour) ...
        - vf_f(iFace)*theFluxes.FLUXC2f(iFace);
    anb{iNeighbour}(iNeighbourOwnerCoef) = anb{iNeighbour}(iNeighbourOwnerCoef) ...
        - vf_f(iFace)*theFluxes.FLUXC1f(iFace);
    bc(iNeighbour) = bc(iNeighbour) ...
        + vf_f(iFace)*theFluxes.FLUXTf(iFace);
end

对于每个内部面来说，其系数同时组装到owner和neighbour两个单元的方程中去。Owner的系数是ac(iOwner)、anb(iOwner)（iOwnerNeighbourCoef)和bc(iOwner)，而Neighbour的系数则是ac(iNeighbour)、anb(iNeighbour)(iNeighbourOwnerCoef)和bc(iNeighbour)。组装俩方程时候不同的正负号是因为面矢量是指向neighbour单元，而指出owner单元的。

10.2 OpenFOAM

to be continued…

梅冠华

关注

11
点赞
踩
16

收藏

觉得还不错? 一键收藏
11
评论
FVM in CFD 学习笔记_第8章_空间离散之扩散项

本章将详细讲解由Laplace算子所表示的扩散项的空间离散方法。流体控制方程中，扩散项和对流项反映了两种截然不同的物理现象，所以两者的离散方法也是不同的，需要分开来讲解。本章首先讲解二维矩形计算域在笛卡尔网格上的含源项的扩散方程的离散方法，然后讲解Dirichlet、Von Neumann、mixed和symmetry边界条件的添加方法。接下来，介绍下在非正交的笛卡尔网格上的离散放法，并详细讲讲非正交的结构和非结构网格上的离散方法、非正交交叉扩散，还会提及非各向同性扩散的处理方法、对高度非线性系统的....
复制链接

扫一扫