OpenFOAM的fvVectorMatrix——ldu矩阵学习（一）

因为研究需要，学习OpenFOAM里的fvVectorMatrix中的ldu矩阵。fvVectorMatrix是OpenFOAM中的一个数据类型，存放关于矢量的线性方程组信息。在simpleFoam/UEqn中可以看到关于速度场U的UEqn：

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + fvm::laplacian(nu, U)
      ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

UEqn存储了关于U方程的所有线性方程组的信息，线性方程组有如下形式：

那么，UEqn.diag()就是A矩阵对角线上的值，UEqn.source()则是b。在fvOptions的功能如显示多孔介质会经常用到这两个函数。

除了diag()和source()之外，fvVectorMatrix还包含upper()和lower()，分别代表A矩阵中对角线之外的上半部分和下半部分。在一整个n*n的矩阵中，因为与单个cell相邻的其他体元不会超过6个（三维，二维则是4个），所以每一行会有至少n-6个空位，也就是0值。0代表着对应位置的体元和当前行对应的体元不相邻，求解的时候没有数量上的关系。如果将A矩阵中所有的元素（n*n个）都进行存储，那么其中绝大部分都是0值，造成空间上的浪费。因此，OpenFOAM的做法是只记录和当前体元有邻里关系的若干个体元的位置，以及对应的值，并存放在upper()或lower()中。

本文所用的测试算例就是一个30x10的plug算例，在中间有一层多孔介质。边界条件是速度入口和压强出口，上下壁面设置为无滑移。如下图所示，蓝色区域为多孔介质区。为了方便看数据，计算域实际长度为3，宽度为1。横向，纵向分别划分了30和10格，每一格为边长为0.1的正方形。z方向长度为1，对后面涉及几何的计算基本无影响。入口水平速度为1。

UEqn中的每一个项都是fvVectorMatrix类型的，下面将分别学习各部分的结构。

1. 对流项fvm::div(phi,U)

1.1 对流项的diag()

在UEqn之前加这一段代码：

Info<<"start correcting ldu matrix! "<<endl;
// - save fvm::div(phi, U)'s ldu
fvVectorMatrix fvmDiv(fvm::div(phi, U));

//diag() - owner's contribution
scalarField outerDivDiag(outerCells.size(), 0.0);
scalarField innerDivDiag(outerCells.size(), 0.0);
for(label celli = 0; celli != outerCells.size(); celli++)
{
    outerDivDiag[celli] = fvmDiv.diag()[outerCells[celli]];
    Info<<"outerDivDiag: "<<outerDivDiag[celli]<<endl;
    Info<<"cell's location"<<mesh.C()[outerCells[celli]]<<endl;
    for(label i : mesh.cells()[outerCells[celli]])
    {

        if(mesh.isInternalFace(i))
        {
            Info<<"the "<<i<<"th faces is an internal face, and its flux is: "
                <<phi[i]<<endl;
        }
        else if(isA<wallFvPatch>(mesh.boundary()[mesh.boundaryMesh().whichPatch(i)]))
        {
            Info<<"the "<<i<<"th faces is a wall face, and its flux is: "
                <<phi[i]<<endl;            
        }
        else
        {
            Info<<"the "<<i<<"th faces is an empty face, and its flux is: "
                <<phi[i]<<endl;            
        }
    }
    Info<<"face's Area: "<<mesh.magSf()[interfaceList[celli]]<<endl;
    Info<<"cell's volume: "<<mesh.V()[outerCells[celli]]<<endl;
}

其中，outerCells是介质之外的体元，以此来学习ldu的过程。上段代码展示了对流项的diag()的提取，并输入相关的数据。在算例中输入求解器的命令，输出如下信息（截取）：

outerDivDiag: 0.1
cell's location(0.95 -0.45 0)
the 18th faces is an internal face, and its flux is: 0.1
the 19th faces is an internal face, and its flux is: 0
the 599th faces is a wall face, and its flux is: 0.0005
the 730th faces is an empty face, and its flux is: 0.0005
the 1030th faces is an empty face, and its flux is: 0.05
the 16th faces is an internal face, and its flux is: 0.1
face's Area: 0.1
cell's volume: 0.01

因为是均匀网格，所有体元的体积和每个面的面积都是一致的。面积为0.1*1=0.1，体积为0.1*0.1*1=0.01，输出结果与预期符合。因为一开始的速度都是1，所以面通量也就是面上速度乘面积得0.1*1=0.1。体元的对流项的diag()是0.1。看东岳网上的求解器介绍，其与当前体元UP有关的系数是：

 因为该体元有一个面是壁面的face，如果只算internal的话，那么diag=(0.1+0.1)/2=0.1，与输出的值相符合。另外empty上的通量应该也是不计入公式中的。

看另一个不与壁面接壤的体元：

cell's location(2.05 -0.35 0)
the 99th faces is an internal face, and its flux is: 0.1
the 100th faces is an internal face, and its flux is: 0
the 841th faces is an empty face, and its flux is: 0.0005
the 1141th faces is an empty face, and its flux is: 1.88979e-316
the 41th faces is an internal face, and its flux is: 0
the 97th faces is an internal face, and its flux is: 0.1
face's Area: 0.1
cell's volume: 0.01

internal faces的通量和除以二就是体元的diag()。empty上的通量数值会比较奇怪，计算的时候应该不会算上它们的通量。但是如果多跑几步，结果就会是这样：

outerDivDiag: 0.100109
cell's location(0.95 0.15 0)
the 372th faces is an internal face, and its flux is: 0.099766
the 373th faces is an internal face, and its flux is: -0.000787717
the 736th faces is an empty face, and its flux is: 1.01711
the 1036th faces is an empty face, and its flux is: 0.998378
the 314th faces is an internal face, and its flux is: -0.000343124
the 370th faces is an internal face, and its flux is: 0.0995554
flux's sum is: 0.198191
face's Area: 0.1
cell's volume: 0.01

通量之和除以二并不是等于diag，而是稍微小一些。为什么会这样？哪里认识不充分？可能需要到fvm::div的源代码中找答案。

参考：

在OpenFOAM中获取网格详细信息_姜蜉蝣的博客-CSDN博客

【翻译】OpenFOAM中的空间离散化与矩阵系数 - 知乎