版本对应：
不管是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;