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方程的所有线性方程组的信息,线性方程组有如下形式:

A\\bold{x}=b

那么,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()中。

本文所用的测试算例就是一个30×10的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中的空间离散化与矩阵系数 – 知乎

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片