请教个openfoam的关于fvmatrix问题

BeatNBBZ

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
               fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)
                );
                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();

                if (nonOrth == nNonOrthCorr)
                {
                    phi -= pEqn.flux();
                }
            }
看不懂这段代码，我怎么觉得这个for循环没有用呢， fvScalarMatrix pEqn( fvm::laplacian(rUA, p) == fvc::div(phi) );这语句跟p有什么关系？fvm::laplacian(rUA, p)返回的不就是个矩阵类嘛？矩阵的系数难道还和p的当前值有关？
我觉得每次执行pEqn.solve();的结果都是一样的，不解。请指教

su_junwei

这里的for循环是进行非正交修正的。 当相邻的网格不正交时，这是将拉普拉斯离散过程中的压力梯度项的表面流律计算分为两部分：正交＋非正交，而非正交部分是作显示处理的，因此当前的计算结果不会对正交部分产生影响，但会影响到非正交部分，上面的for循环就是为了让这部分在迭代过程中更趋近真值。  
fvScalarMatrix pEqn( fvm::laplacian(rUA, p) == fvc::div(phi) );
这句话就是定义压力波松方程的，求解新的压力。
非正交修正完之后，进行界面流律更新。

BeatNBBZ

多谢指教。但是我的疑惑之处是：除非fvm::laplacian(rUA, p)生成的矩阵的系数跟p有关，否则反复执行fvScalarMatrix pEqn( fvm::laplacian(rUA, p) == fvc::div(phi) );和pEqn.solve();之后求解的p的值是不变的。因为rUA在循环中没有变化，phi只在nonOrth == nNonOrthCorr之后才变化，也就是对于nonOrth < nNonOrthCorr时每次执行pEqn.solve();我觉得p的值都不变。难道fvm::laplacian(rUA, p)生成的矩阵的系数跟p有关吗？按理说laplacian离散后，系数只跟网格几何结构有关啊

su_junwei

呵呵，fvm::laplacian(rUA,p)  生成两部分，一部分是系数矩阵，另外一部分是显示项b（非正交部分） fvc::div(phi)生成系数矩阵b的另外一部分。也就是
Ax+b1==b2
修正过程中对A没有影响，但会对非正交部分b1有影响，迭代过程中会使得b1不断接近真值。

liuhuafeifei

非常好的一个问题！！！

rUa是扩散相系数，在循环中为变化，fvc::div(phi)也为变化，系数矩阵Ap,AL,Au也没有变化，但是非正交源项，即 fvc::linearInterpolation(rUa)*CorrectionVector()*fvc::linearInterpolation(grad(p))，扩散系数乘以非正交向量*压力梯度，是随着p的迭代变化的，所以需要迭代，但是迭代仅仅需要改变的是源项。

这里可能触及了openfoam需要改进的地方，它的非正交源项必须用fvm::laplacian(rUA, p) 才能得到。但是每次迭代都要重新计算，AP AL AU SU，只有SU是必需的，其它都是无用功。再往深处探究，线形方程求解器，多数是要预条件的，预条件需要的D-1的计算也是无用功，因为它只与AP，AL，AU相关。

如果是其它方程，浪费一点时间也无所谓，但是压力方程是分离式算法求解的主要耗时步骤，如果openfoam真如上面所分析，还真值得改进。

仅供参考
liuhuafei