[转载]OpenFOAM中的神奇方程定义方式的背后
没有用过或者刚开始用OpenFOAM的人会感觉到OpenFOAM中的方程定义紧凑性所惊奇，我开始也是被OpenFOAM微分方程定义的精炼性所吸引的，从而放弃了mfix而转向OpenFOAM的。闲言少叙。OpenFOAM中的方程类为fvMatrix，该类为一个基于有限容积的稀疏矩阵类，更确切的说她是一个类模板。对于定义标量方程 fvMatrix<scalar>，也可以写成fvScalarMatrix，两者一样，因为openfoam运用了 typedef进行了类型别名定义
typedef fvMatrix<scalar> fvScalarMatrix ; //标量
typedef fvMatrix<vector> fvVectorMatrix ; //矢量
搞cfd的人都知道，微分方程是通过离散从而转化为代数方程组进行求解。对于一个代数方程组主要有两部分：系数矩阵A，右边值b。构造成类似Ax＝b的代数方程形式。
方程的离散有显式和隐式,隐式离散得到方程稀疏矩阵A的，显式离散得到右边的值b。一旦通过某种离散形式得到A和b就可以求解线性方程组了。
fvMatrix中有两个基本变量A 和 b初始化为0，用以存储系数矩阵A和右边值b。
OpenFOAM的所有的显式离散都在fvc空间中，所以以fvc开头的都为显式，比如fvc::div(phi).隐式离散在fvm空间中，所有的隐式离散都是以fvm开头，比如fvm::ddt, fvm::div(phi,U)等。所以，以fvm开头的项声称A，以fvc开头的项生成b。下面以的动量方程为例对他们的实现过程进行详述
fvVectorMatrix UEqn //定义一个fvmatrix对象，向量的方程
(
fvm::ddt(U)
+ fvm::div(phi, U)
-fvm::laplacian(nu, U)
==
-fvc::grad(p)
);
fvm::ddt(U)隐式离散生成矩阵类，该项会被加到UEqn中的A上
fvm::div(phi, U)隐式离散矩阵类，该项加到UEqn中的A上
-fvm::laplacian(nu, U)隐式离散，该项从UEqn中的A上减下来，因为前面是负号。
==是c＋＋优先级最小的符号，在OpenFOAM中他的功能和“－”相同。
-fvc::grad(p)，显式离散，该项加入到UEqn中的b上，前面“－”又有一个“＝＝”负负得正
这样稀疏矩阵A和右边的b在UEqn中已全部构建好，也就是得到方程了，求解即可。
solve(UEqn);
or
UEqn.solve();
功能一样，前面是全局函数，后面这个是UEqn的成员函数罢了。
也许你会感觉到纳闷，OpenFOAM怎么知道他该加到A上还是b上。这点很简单，fvm返回的是fvmatrix, 而fvc返回的是GeometryField;类型不一样，重载一下就行了。呵呵。
还有一个问题,如何避免下面的加和，U1和U2属于不同的场
fvm::ddt(U1)+fvm::div(phi,U2)
所以在每次相加之前，他都会check一下，看看两个相加或者相减的两者是否是一个对象的操作（比较他们地址是否一样即可），因此fvMatrix里面还有一个重要变量psi_，来记录当前求解变量的地址（实际上里面存有求解场的一个引用）.
现在明朗点了吗？
转自OpenFOAM研究：http://blog.sina.com.cn/openfoamresearch
页:
[1]