版本对应:

不管是OpenFOAM-11,还是老版本的OpenFOAM,势流求解器都是potentialFoam。

CFD: 势流算法

大量的流体力学现象,如尾涡和边界层现象,都和有旋流动密不可分。但无旋流同样大量的存在于某些细分的流体力学领域(如航空、地下水)。下图为一个圆柱势流场。很明显其是无旋的。对于无旋流的模拟,并不需要求解复杂的NS方程,求解势流方程即可,且势流方程是一个非常简单的拉普拉斯方程,求解起来非常快速。因此求解势流,可以用于快速的预估一个场来加速NS方程的迭代收敛。在势流求解算法中的主要问题是,因为速度势并不是原始变量,因此直接给速度势的边界条件在一般条件下是很困难的。因此势流求解算法,主要是研究如何通常通过速度边界条件,来隐性的施加到速度势边界条件。

_images/potential.gif

普适性势流方程与求解

首先有稳态不可压缩流体连续性方程:

(1)\[ \nabla \cdot \bfU =0 \]

对于无旋流(\(\nabla\times\bfU=0\)),则必定存在一个势函数\(\Phi\),使得

(2)\[ \nabla \Phi=\bfU \]

为了使得速度势与压力同向,在求解的时候同样也可以定义:

(3)\[ -\nabla \Phi=\bfU \]

本文选用方程(3)作为速度势的定义。

Warning

一维情况下方程(3)可以写为\(u=-\frac{\rd \Phi}{\rd x}\)。可见,\(x\)方向\(\Phi\)越大,\(-\frac{\rd \Phi}{\rd x}=u\)为一个负数,表示向\(x\)左侧移动,与压力作用相同。

对于不可压缩流动,将方程(3)代入到方程(1),有:

(4)\[ -\nabla\cdot(\nabla \Phi)=0 \]

方程(4)即为控制势流的速度势方程。可以看出在不可压缩的情况下,速度势为一个拉普拉斯方程。给定\(\Phi\)的边界条件,就可以进行求解。考虑最简单的情况,仅仅有\(x\)方向,进出口速度必然是已知的(因为要满足通量守恒)。因此进出口边界都可以给速度固定值边界条件。假定进口速度是1,那么出口速度也是1,也即\(\frac{\p\Phi}{\p x}\)在进出口都是固定值1。因此其可以用方程描述:

(5)\[ -\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=0, \frac{\p\Phi}{\p x}_{inlet}=1, \frac{\p\Phi}{\p x}_{outlet}=1 \]

方程(5)由于相容性问题存在无穷多个解。在实际操作过程中,可以任意的固定某个网格点的\(\Phi\)值,也即给定参考\(\Phi\)值。方程(5)可以求出。

Warning

在不存在固定值边界条件的时候,压力方程也存在相容性问题。因此在OpenFOAM里的pRefValue以及pRefCell就是为了解决相容性问题需要设置的参数。其给定一个参考压力,以及给定一个参考网格。更详细的信息请参考无痛苦ns方程笔记

还是考虑这个1D问题。一般情况下在CFD中给定的是进口速度固定值,出口速度让其自由的去发展,因此出口速度的边界条件为\(\frac{\p u}{\p x}=0\)。在这种情况下,出口的\(\Phi\)不再是固定值。其边界条件变成了\(\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)_{outlet}=0\)。然而其就退化成了方程(5),会导致求解困难。但因为势流是求解无旋的流动,因此为了避免这个问题,可以将出口设置成足够的远,在这种情况下,已知进口速度,已知足够远的出口,为了满足连续性方程,出口速度也必然是已知的。因此还可以通过方程(5)来进行计算。

_images/poten.JPG

上图左侧对于圆柱扰流。如果出口边界防止足够的远,对于无旋流动在充分发展后,其出口速度必然是与壁面垂直的(否则会发生壁面穿透)。如果出口壁面距离圆柱比较近,其也会存在一个速度但方向不太好确定。这与上图右侧的几何情况类似,即使出口边界防止足够的远,对于无旋流动在充分发展后,其出口速度也是固定的(不能发生壁面穿透)。但这个速度方向虽然是固定的,但不太好求出。那么给定速度势的出口固定梯度边界条件就比较麻烦。

可见问题主要出现在势流求解的出口边界给定上。回到方程(5),是否可以给定\(\Phi\)一个固定值边界条件,来使其迂回的实现方程(5)的效果?答案是可以的。

回到\(\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=0\),其可以想象成为一个导热方程。对于方程(5)中给定进出口都固定梯度的条件,其可以通过下述条件来施加:

  1. 出口给定\(\Phi\)的固定值边界条件,进口给定\(\Phi\)的零法向梯度边界条件(可以理解成进口去充分发展);

  2. 在计算域给定\(\frac{\p\Phi}{\p x}\)的值,使其恒定;

上述两个物理解释可以通过方程来进行解释:

(6)\[\begin{split} &-\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=0 \\\\ &\frac{\p\Phi}{\p x}_{inlet}= u^{inlet}, \Phi_{outlet}=\mathrm{fixedValue} \end{split}\]

方程(6)中的\(\frac{\p\Phi}{\p x}_{cells}= u^{inlet}\)的作用类似保证\(\Phi\)随着\(x\)的分布是一条直线,其斜率为\( u^{inlet}\)。这个\( u^{inlet}\)的值就是进口速度的值。方程(6)与方程(5)可以实现相同的效果。因此,方程(6)构成一种更简单的构建势流求解的方法。

OpenFOAM势流求解方法

OpenFOAM使用了另外一种“天才”的做法,将方程(6)做的更加普适性,其原则就是势流方程简单改写,然后可以将\(\Phi\)的固定梯度通过零法向梯度来施加。

_images/pt2.JPG

首先,给定上述的网格系统。现对方程(6)进行离散,在网格点1处:

(7)\[ -\left( \frac{ \frac{\p\Phi}{\p x}_B-\frac{\p\Phi}{\p x}_A }{\Delta x} \right) =0 \]
(8)\[ -\left( \frac{\Phi_2-\Phi_1}{\Delta x} -\frac{\p\Phi}{\p x}_A \right) =0 \]
(9)\[ -\left( \frac{\Phi_2-\Phi_1}{\Delta x} \right) =\frac{\p\Phi}{\p x}_A=u_{inlet} \]

上述方程是采用方程(6)进行离散的(\(\Phi\)在进口固定梯度边界条件)。

深入来看方程(8)(9),其也可以通过零梯度边界条件来进行。在这种情况下\(\frac{\p\Phi}{\p x}_A=0\),方程(9)可以理解为\(\Phi\)在A处 1)施加零梯度边界条件,2)附加一个离散源项。二者结合可以实现方程(9)同样的效果。

Warning

同时要注意,这个源项仅仅在零法向梯度边界存在值。否则其与方程(6)不一致。

截止至此,可以将求解方程写为:

(10)\[\begin{split} &-\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=源项 \\\\ &\frac{\p\Phi}{\p x}_{inlet}= 0, \Phi_{outlet}=\mathrm{fixedValue} \end{split}\]

对比方程(10)(6),很明显,\(u^{inlet}\)要进入到源项里面才能使得二者相等。源项定义在每个网格点上。这个源项在整个网格系统,仅仅在\(\Phi\)的零法向梯度边界相邻的网格单元存在值。在一维情况下,源项可以表示为\(u_{inlet},0,0,0,0,...,0\)这样的列向量。

(11)\[\begin{split} &-\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=S \\\\ &\frac{\p\Phi}{\p x}_{inlet}= 0, \Phi_{outlet}=\mathrm{fixedValue} \\\\ &S=[u_{inlet},0,0,0,0,...,0]^T \end{split}\]

方程(11)是在1D情况下,OpenFOAM求解的势流方程。

现在将其进行3维拓展。将方程(11)改写为:

(12)\[\begin{split} &-\nabla\cdot(\nabla \Phi)=S \\\\ &\bfn\cdot\nabla\Phi_{inlet}=0, \Phi_{outlet}=\mathrm{fixedValue} \end{split}\]

最后一步是,如何把\(\bfU_{inlet}\)\(S\)建立关系。首先,将\(-\nabla\cdot(\nabla\Phi)=0\)写成离散形式:

(13)\[ -\sum (\nabla\Phi)_f\cdot\bfS_f=0 \]

其也可以写成:

(14)\[ -\sum (\nabla\Phi)_{内部面}\cdot\bfS_{内部面} -\sum (\nabla\Phi)_{进口}\cdot\bfS_{进口} =0 \]

在给定\((\nabla\Phi)_{进口}\)为固定值边界条件的时候,其等同于

(15)\[ -\sum (\nabla\Phi)_{内部面}\cdot\bfS_{内部面} -\sum (\bfU)_{进口}\cdot\bfS_{进口} =0 \]

移项有:

(16)\[ -\sum (\nabla\Phi)_{内部面}\cdot\bfS_{内部面} =\sum (\bfU)_{进口}\cdot\bfS_{进口} \]

如果不对内部面以及边界面做区分,其等同于:

(17)\[\begin{split} &-\sum (\nabla\Phi)_f\cdot\bfS_f=\sum \bfU_f\cdot\bfS_f, \\\\ &(\nabla\Phi)_{进口}=零法向梯度, \bfU_f^{边界}=固定值, \bfU_f^{内部}=0 \end{split}\]

注意,方程(17)\(\bfU\)只有在进口边界处存在值,其内部场的值为0。方程(17)写成连续形式即为:

(18)\[ -\nabla\cdot(\nabla \Phi)=\nabla\cdot\bfU, (\nabla\Phi)_{进口}=零法向梯度, \bfU_{内部}=0,\bfU_{进口}=固定值 \]

方程(18)为OpenFOAM最终植入的方程。

速度更新

在有了\(\Phi\)场之后,可以通过方程(3),也即\(\bfU=-\nabla \Phi\)来计算速度。我们的目的是通量守恒。在这里再一次强调这个概念,有限体积法的目的是通量守恒,而不是速度守恒,速度本身并不是一个守恒量。类似的,有限体积法强调的是能量守恒,而不是温度守恒。温度也不是一个守恒量。所以目前看起来可能有两个路径:

  • 通过\(\bfU=-\nabla \Phi\)来计算速度;

  • 通过一种方法来计算守恒通量,然后从通量来构建速度;

很明显,通过\(\bfU=-\nabla \Phi\)来计算速度,然后再从速度插值到通量这是不能满足通量守恒的。这种方法不会满足任何特性。仅仅是满足\(\bfU=-\nabla \Phi\)

另一方面,如何构建守恒的通量?回看方程(17),会发现其实\(-(\nabla\Phi)_f\cdot\bfS_f\)中的\((\nabla\Phi)_f\)就是\(\Phi\)插值到面上的值,\((\nabla\Phi)_f\cdot\bfS_f\)就是面上的\(\Phi\)的通量,因为\(\sum \bfU_f\cdot\bfS_f=0\),因此方程(17)表示的就是面上的\((\nabla\Phi)_f\)的通量是守恒的。同时,\(-\nabla\Phi=\bfU\),也就是表示面上的\((\nabla\bfU)_f\)的通量是守恒的。这自然而然就是通量守恒。\((\nabla\bfU)_f\cdot\bfS_f\)就是守恒通量。在有了守恒通量后,可以通过守恒通量来重组速度场。

势流求解也可以判断是否满足连续性方程。在这种情况下,其判断\((\nabla\bfU)_f\cdot\bfS_f\)是否等于\(0\)

代码分析

针对方程(18)中的\(\bfU_{内部}=0,\bfU_{进口}=固定值\),通过下述代码进行:

volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

U = dimensionedVector(U.dimensions(), Zero);

方程(18)中的\(-\nabla\cdot(\nabla \Phi)=\nabla\cdot\bfU\),通过下述代码进行:

fvScalarMatrix PhiEqn
(
    fvm::laplacian(dimensionedScalar(dimless, 1), Phi)
  ==
    fvc::div(phi)
);

方程(18)中的\((\nabla\Phi)_{进口}=零法向梯度\),通过压力边界条件来进行(不需要显性的指定\(\Phi\)的边界条件)。例如在算例中,压力p还是通过常规的进口零法向梯度,出口固定值边界条件来设置。代码直接将压力的边界转移到\(\Phi\)边界。

下面的代码计算守恒的通量,以及重组速度场:

phi -= PhiEqn.flux();
U = fvc::reconstruct(phi);
U.correctBoundaryConditions();

下面的代码用来判断连续性方程误差:

Info<< "Continuity error = "
    << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
    << endl;

下面的代码,用来判断从守恒的通量场,重组后获得的速度场,继续插值到面上获得的通量场,与原本守恒通量场的误差:

Info<< "Interpolated velocity error = "
    << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
    << endl;