版本对应:
不管是OpenFOAM-11,还是老版本的OpenFOAM,势流求解器都是potentialFoam。
CFD: 势流算法
大量的流体力学现象,如尾涡和边界层现象,都和有旋流动密不可分。但无旋流同样大量的存在于某些细分的流体力学领域(如航空、地下水)。下图为一个圆柱势流场。很明显其是无旋的。对于无旋流的模拟,并不需要求解复杂的NS方程,求解势流方程即可,且势流方程是一个非常简单的拉普拉斯方程,求解起来非常快速。因此求解势流,可以用于快速的预估一个场来加速NS方程的迭代收敛。在势流求解算法中的主要问题是,因为速度势并不是原始变量,因此直接给速度势的边界条件在一般条件下是很困难的。因此势流求解算法,主要是研究如何通常通过速度边界条件,来隐性的施加到速度势边界条件。
普适性势流方程与求解
首先有稳态不可压缩流体连续性方程:
对于无旋流(\(\nabla\times\bfU=0\)),则必定存在一个势函数\(\Phi\),使得
为了使得速度势与压力同向,在求解的时候同样也可以定义:
本文选用方程(3)作为速度势的定义。
Warning
一维情况下方程(3)可以写为\(u=-\frac{\rd \Phi}{\rd x}\)。可见,\(x\)方向\(\Phi\)越大,\(-\frac{\rd \Phi}{\rd x}=u\)为一个负数,表示向\(x\)左侧移动,与压力作用相同。
方程(4)即为控制势流的速度势方程。可以看出在不可压缩的情况下,速度势为一个拉普拉斯方程。给定\(\Phi\)的边界条件,就可以进行求解。考虑最简单的情况,仅仅有\(x\)方向,进出口速度必然是已知的(因为要满足通量守恒)。因此进出口边界都可以给速度固定值边界条件。假定进口速度是1,那么出口速度也是1,也即\(\frac{\p\Phi}{\p x}\)在进出口都是固定值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)来进行计算。
上图左侧对于圆柱扰流。如果出口边界防止足够的远,对于无旋流动在充分发展后,其出口速度必然是与壁面垂直的(否则会发生壁面穿透)。如果出口壁面距离圆柱比较近,其也会存在一个速度但方向不太好确定。这与上图右侧的几何情况类似,即使出口边界防止足够的远,对于无旋流动在充分发展后,其出口速度也是固定的(不能发生壁面穿透)。但这个速度方向虽然是固定的,但不太好求出。那么给定速度势的出口固定梯度边界条件就比较麻烦。
可见问题主要出现在势流求解的出口边界给定上。回到方程(5),是否可以给定\(\Phi\)一个固定值边界条件,来使其迂回的实现方程(5)的效果?答案是可以的。
回到\(\frac{\p}{\p x}\left(\frac{\p\Phi}{\p x}\right)=0\),其可以想象成为一个导热方程。对于方程(5)中给定进出口都固定梯度的条件,其可以通过下述条件来施加:
出口给定\(\Phi\)的固定值边界条件,进口给定\(\Phi\)的零法向梯度边界条件(可以理解成进口去充分发展);
在计算域给定\(\frac{\p\Phi}{\p x}\)的值,使其恒定;
上述两个物理解释可以通过方程来进行解释:
方程(6)中的\(\frac{\p\Phi}{\p x}_{cells}= u^{inlet}\)的作用类似保证\(\Phi\)随着\(x\)的分布是一条直线,其斜率为\( u^{inlet}\)。这个\( u^{inlet}\)的值就是进口速度的值。方程(6)与方程(5)可以实现相同的效果。因此,方程(6)构成一种更简单的构建势流求解的方法。
OpenFOAM势流求解方法
OpenFOAM使用了另外一种“天才”的做法,将方程(6)做的更加普适性,其原则就是势流方程简单改写,然后可以将\(\Phi\)的固定梯度通过零法向梯度来施加。
首先,给定上述的网格系统。现对方程(6)进行离散,在网格点1处:
上述方程是采用方程(6)进行离散的(\(\Phi\)在进口固定梯度边界条件)。
深入来看方程(8)与(9),其也可以通过零梯度边界条件来进行。在这种情况下\(\frac{\p\Phi}{\p x}_A=0\),方程(9)可以理解为\(\Phi\)在A处 1)施加零梯度边界条件,2)附加一个离散源项。二者结合可以实现方程(9)同样的效果。
Warning
同时要注意,这个源项仅仅在零法向梯度边界存在值。否则其与方程(6)不一致。
截止至此,可以将求解方程写为:
对比方程(10)与(6),很明显,\(u^{inlet}\)要进入到源项里面才能使得二者相等。源项定义在每个网格点上。这个源项在整个网格系统,仅仅在\(\Phi\)的零法向梯度边界相邻的网格单元存在值。在一维情况下,源项可以表示为\(u_{inlet},0,0,0,0,...,0\)这样的列向量。
方程(11)是在1D情况下,OpenFOAM求解的势流方程。
现在将其进行3维拓展。将方程(11)改写为:
最后一步是,如何把\(\bfU_{inlet}\)与\(S\)建立关系。首先,将\(-\nabla\cdot(\nabla\Phi)=0\)写成离散形式:
其也可以写成:
在给定\((\nabla\Phi)_{进口}\)为固定值边界条件的时候,其等同于
移项有:
如果不对内部面以及边界面做区分,其等同于:
注意,方程(17)中\(\bfU\)只有在进口边界处存在值,其内部场的值为0。方程(17)写成连续形式即为:
方程(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;