OpenFOAM 中的有限体积离散
Contents
介绍
这篇文章中你将会学到:
- 基于任意多面体网格的有限体积离散方法
- 时间项的处理方式
- 对流项的处理方式
- 源项的处理方式
- 边界条件的处理方式
- OpenFOAM 中对有限体积离散的代码实现
基本概念
OpenFOAM 采用基于同位网格(collocated grid)的有限体积法进行离散,所有的变量都存储在控制体网格内部的同一个点上。
(图片来自 Hrvoje Jasak1)
考虑有如图所示的任意多面体网格单元 $V_P$ ,$P$ 为其体心,定义如下
\[ \begin{equation} \int_{V_P} (\mathbf x - \mathbf x _P) dV = 0 \end{equation} \]
基于以上定义,假设物理量在控制体内线性变化, 位于单元内 $\mathbf x$ 处的任意物理量可表示为
\[ \begin{equation} \phi (\mathbf x) = \phi_P + (\mathbf x - \mathbf x_P) \cdot (\nabla \phi)_P \end{equation} \]
将上式对单元进行体积分,可得
\[ \begin{equation} \begin{aligned} \int_{V_P} \phi (\mathbf x) dV &= \int_{V_P} [\phi_P + (\mathbf x - \mathbf x_P) \cdot (\nabla \phi)_P ] dV \\ &= \phi_P V_P + (\nabla \phi)_P \cdot \int_{V_P} (\mathbf x - \mathbf x_P) dV \\ & = \phi_P V_P \end{aligned} \end{equation} \]
对流扩散方程
考虑一般形式的对流扩散方程,其形式可以表示如下
\[ \begin{equation} \underbrace{\frac{\partial \phi}{\partial t}}_{时间项} + \underbrace{\nabla\cdot (\phi \mathbf u)}_{对流项} - \underbrace{\nabla \cdot (\gamma \nabla \phi)}_{扩散项或粘性项} = \underbrace{S_\phi(\phi)}_{源项} \end{equation} \]
将其对单元做体积分,可得
\[ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV + \int_{V_P} \nabla\cdot (\phi \mathbf u) dV - \int_{V_P} \nabla \cdot (\gamma \nabla \phi) dV = \int_{V_P} S_\phi(\phi) dV \end{equation} \]
可通过高斯定理将对流项和扩散项的体积分转为面积分,上式变为
\[ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV + \sum_{f} [\mathbf S_f \cdot (\phi \mathbf u)_f] - \sum_{f} [\mathbf S_f \cdot (\gamma \nabla \phi)_f] = \int_{V_P} S_\phi(\phi) dV \label{eq:int_diff_conv} \end{equation} \]
时间项的处理
为了表述方便,我们定义 $\phi^o$ 和 $\phi^n$ 为旧时刻和新时刻的值(通常省略 $n$)。以Euler格式为例, $\eqref{eq:int_diff_conv}$ 中的时间项可以写作
\[ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV =\frac{\phi^n - \phi^o}{\Delta t} V_P \end{equation} \]
对流项的处理
式 $\eqref{eq:int_diff_conv}$ 中的下标 $f$ 代表面单元上的物理量,这些值通过与相邻单元插值得到。考虑下图所示的两个相邻控制体单元 $V_P$ 和 $V_N$ ,它们共享面 $f$, $\mathbf S_f$ 为面法向向量,其大小为 $f$ 的面积。
(图片来自 OpenFOAM Foundation2)
面心物理量可由体心物理量插值得到,以二阶精度的中心差分格式为例
\[ \begin{equation} \phi_f = \lambda \phi_P + (1 - \lambda) \phi_N \end{equation} \]
这里 $\lambda = \overline{fN} / \overline{PN}$ 。
于是 $\eqref{eq:int_diff_conv}$ 中的对流项可以表示为
\[ \begin{equation} \sum_{f} [\mathbf S_f \cdot (\phi \mathbf u)_f] = \underbrace{\sum_{f} [\lambda (\mathbf S_f \cdot \mathbf u_P)] \phi_P}_{P 单元} + \underbrace{\sum_{f} [(1 - \lambda) (\mathbf S_f \cdot \mathbf u_N) \phi_N]}_{与 P 相邻的所有单元} \label{eq:conv_discre} \end{equation} \]
扩散项的处理
扩散项由于包含面法向 $\mathbf S_f$ 与梯度项 $(\nabla \phi)_f$ 的点乘,处理起来比对流项更为复杂。
先看简单的情况,考虑正交网格, $\mathbf S_f$ 与 $\overline{PN}$ 平行,扩散项可以写作
\[ \begin{equation} -\sum_{f} [\mathbf S_f \cdot (\gamma \nabla \phi)_f] = -\sum_{f} \gamma_f | \mathbf S_f | \frac{\phi_N - \phi_P}{|\overline{PN}|} \end{equation} \]
对于非正交网格,则需对上式进行修正,考虑如下图所示的非正交网格,面法向量表示为 $\mathbf S_f = \Delta + \mathbf k$ ,其中 $\Delta$ 与 $\overline{PN}$ 平行,$\mathbf k$ 与面 $f$ 平行
(图片来自 Hrvoje Jasak)
于是有
\[ \begin{equation} \mathbf S_f \cdot (\nabla \phi)_f = \underbrace{|\Delta| \frac{\phi_N - \phi_P}{|\overline{PN}|}}_{正交部分} + \underbrace{\mathbf k \cdot (\nabla \phi)_f}_{非正交修正部分} \end{equation} \]
非正交修正将产生两部分:正交部分和非正交修正部分。OpenFOAM 将正交部分进行隐式处理,将其放入系数矩阵中,而将非正交部分显式处理,放入源项中。这样一来便涉及到两个问题
- 显式处理的非正交部分采用的是旧迭代步(或时间步)的值,存在信息滞后
- 显式处理的非正交部分放入源项中,将产生一个负的源项,若太大将导致求解崩溃
对于第一个问题,通常采用迭代方式来减小旧值和新值的差距。如在 system/fvSchemes
里通过指定关键字 nNonOrthogonalCorrectors
来指定求解压力泊松方程时需要进行非正交修正的次数。
对于第二个问题,则通过对放入源项的非正交部分 $\mathbf k \cdot (\nabla \phi)_f$ 施加限制器 $\lambda$ 来保证数值稳定,定义如下
\[ \begin{equation} \mathbf k \cdot (\nabla \phi)_f < \lambda|\mathbf S_f| \frac{\phi_N - \phi_P}{|\overline{PN}|} \end{equation} \]
源项的处理
源项通常随时间和空间变化,较为复杂。为了确保数值求解的稳定性和有界性,通常对其做线性化处理
\[ \begin{equation} S_\phi(\phi) = Su (\phi) + Sp (\phi) \phi \end{equation} \]
$\eqref{eq:int_diff_conv}$ 中的源项可以写作
\[ \begin{equation} \int_{V_P}S_\phi (\phi) dV = Su V_P + Sp V_P \phi_P \label{eq:linearized_source_term} \end{equation} \]
注意到这里为止,对流项和扩散项离散后的物理量($\phi$ 和 $\mathbf u$)均为新时刻的值,即省略了上角标 $n$ 。
将 $\phi_N$ 和 $\phi_P$ 的系数都相加,同时将时间项中 $\phi^o$ 移至方程右侧,记 $r_P = Su V_P + Sp V_P \phi_P + \frac{\phi^o V_P}{\Delta t}$ ,于是方程 $\eqref{eq:int_diff_conv}$ 变成了
\[ \begin{equation} a_P \phi_P + \sum_f a_N \phi_N = r_P \end{equation} \]
联立所有单元的离散方程得到一个大型稀疏线性方程组,写作
\[ \begin{equation} \mathbf A \pmb \phi = \mathbf R \label{eq:discre_coeffs_matrix} \end{equation} \]
$\mathbf A$ 被称为系数矩阵 ,通常被分为三部分,分别是下三角元素 $\mathbf L$ 、对角元素 $\mathbf D$ 和上三角元素 $\mathbf U$ 。
边界条件的处理
常见的边界条件有两种:一种是给定值的 Dirichlet 边界条件,另一种是给定梯度的 Neumann 边界条件。
对于 Dirichlet 边界,边界上的面单元的值等于给定值: $\phi_f = \phi_b$ 。
- 对流项
\[ \begin{equation} \mathbf S_f \cdot (\phi \mathbf u)_f = \mathbf S_b \cdot (\phi \mathbf u)_b \end{equation} \]
这项只影响 $\mathbf R$ 。
- 扩散项
\[ \begin{equation} \mathbf S_f \cdot (\gamma \nabla \phi)_f = \mathbf S_b \cdot (\gamma \nabla \phi)_b = \gamma_b |\mathbf S_b |\frac{\phi_b - \phi_P}{\mathbf d_b} \end{equation} \]
$\phi_P$ 影响 $\mathbf D$,而 $\phi_b$ 影响 $\mathbf R$ 。
对于 Neumann 边界,边界上的面单元的梯度等于给定值:$\frac{\mathbf S_b}{|\mathbf S_b|}(\nabla \phi)_b = g_b$ 。
- 对流项
边界面上的值 $\phi_b = \phi_P + \mathbf d_b \cdot (\nabla \phi)_b = \phi_P + | \mathbf d_b | g_b$ 。
\[ \begin{equation} \mathbf S_f \cdot (\phi \mathbf u)_f = \mathbf S_b \cdot \mathbf u_b (\phi_P + |\mathbf d_b| g_b) \end{equation} \]
对流项将影响 $\mathbf D$ 和 $\mathbf R$ 。
- 扩散项
\[ \begin{equation} \mathbf S_f \cdot (\gamma \nabla \phi)_f = \gamma_b |\mathbf S_f| g_b \end{equation} \]
扩散项只影响 $\mathbf R$ 。
代码实现
系数矩阵
描述系数矩阵的类为 lduMatrix
及其子类 fvMatrix
,涉及到的主要函数有:
Foam::lduMatrix<Type>::diag()
Foam::lduMatrix<Type>::lower()
Foam::lduMatrix<Type>::upper()
Foam::fvMatrix<Type>::D()
Foam::fvMatrix<Type>::A()
Foam::fvMatrix<Type>::source()
lower()
返回下三角元素 $\mathbf L$diag()
返回不是对角元素 $\mathbf D$,而是未考虑边界条件影响的对角元素upper()
返回上三角元素 $\mathbf U$D()
返回的是对角元素 $\mathbf D$,通过addCmptAvBoundaryDiag
施加边界条件影响
|
|
A()
返回的不是 $\eqref{eq:discre_coeffs_matrix}$ 中的 $\mathbf A$ ,是 $\mathbf D$ 除以网格体积之后的值。
|
|
时间项
时间项产生的系数矩阵由 fvm::ddt
系列函数计算得到。以 icoFoam
中使用的 fvm::ddt(U)
为例,对应的实现如下
|
|
这里的 New
函数利用 RTS 机制构造字典文件中 ddt(U)
指定的时间离散格式,获得系数矩阵。
对流项
对流项产生的系数矩阵由 fvm::div
系列函数计算得到。同样以 icoFoam
中使用的 fvm::div(phi, U)
为例,对应的实现如下
|
|
实际上调用的是 fvm::div(phi, U, "div(phi, U)")
,对应的实现如下
|
|
和时间项一样,这里也利用 RTS机制,读取字典中 "div(phi, U)"
对应的关键字,构造相应的对流格式,获得系数矩阵。
扩散项
扩散项产生的系数矩阵有 fvm::laplacian
系列函数得到。以 icoFoam
中使用的 fvm::laplacian(nu, U)
为例,对应的实现如下
|
|
这里做了一次转换,用类型为 dimensionedScalar
的 nu
构造了一个 surfaceScalarField
类型的场量,最终调用了以下实现
|
|
同样是利用 New
函数,读取字典中 "laplacian(nu, U)"
对应的关键字,构造相应的 laplacian 格式,获得系数矩阵。
源项
式 $\eqref{eq:linearized_source_term}$ 中线性化后的源项右侧有包括 $\phi_P$ 的导数项,需要做隐式处理。OpenFOAM 采用 fvm::su
, fvm::sp
和 fvm::susp
将源项的影响施加到系数矩阵上。其中 fvm::su
用来处理纯显式项 $Su V_P$,fvm::sp
用来处理纯隐式项 $Sp V_P \phi_P$,fvm::susp
则用来处理既包含显式项又包含隐式项的情况。
fvm::su
将 su
的值乘以网格体积后施加到 source()
上:
|
|
fvm::sp
将 sp
的值乘以网格体积后施加到对角元素 diag()
上:
|
|
fvm::susp
较复杂,它为了数值稳定性考虑,把源项里的正负成分(对应源和汇)分开处理,将大于 0 的部分乘以网格体积后施加到对角元素 diag()
上,小于 0 的部分乘以网格体积和 $\phi_P^o$ 后施加到 source()
上:
|
|