[转载]OpenFOAM中的神奇方程定义方式的背后

su_junwei

没有用过或者刚开始用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