首发于有所思

OpenFOAM中的数据结构-icoFoam为例【未完成】

以前我不是很喜欢写这种博客,一方面觉得费时间,一方面担心被大佬批评。后来发现这种写东西整理的方法确实对提高工作效率很有用,并且也没有被大佬痛批。所以今天再开一篇,讲讲OpenFOAM中的数据结构。


今年暑假我就想做这个题目了,但是一直没时间,实习呀回国呀开会呀报告呀,总是像赶火车一样,头痛医头脚痛医脚,现在正好大雪封校,又是周五晚上,我今晚就要把这个问题搞清楚。


事先声明,我写的一切文章,动笔的时候都是还没搞清楚的。我是通过写文章的方式整理思路,所以看到不对劲的地方可以告诉我。我随时改。我写的时候也会不断完善。


那么言归正传,这篇文章要说什么呢?

这篇文章希望从数据结构,也就是矩阵,数列的角度,详细分析一下,OpenFOAM是如何对一个case进行运算的。我将从一个最简单的例子开始(1*2*3的网格),number by number的演示openFOAM的运算方式。我将举两个例子,一个是不带turbulence model的icoFoam例子,一个是带turbulence model的pisoFoam例子。这篇文章先讲icoFoam,以后有空了再讲pisoFoam。

这篇文章将:

  1. 【数学公式】先讲讲icoFoam的数学公式;
  2. 观察代码】然后讲讲代码结构;
  3. 【数据流动】最后结合一个例子讲讲其中数据的流动过程。

  • 【数学公式】

在讨论code之前,必须先了解其数学形式,才好一一对应。

关于icoFoam的数学形式我主要参考李东岳博士的这篇文章:

icoFoamdyfluid.com

基本上icoFoam是一个特别简单的NS方程,没有湍流项,所以也就算不了湍流模型。

一切的开始是动量方程:

我们假设忽略压力梯度项(在后文中即对应动量预测步骤):

因为OpenFOAM是FVM(有限体积法),所以你需要对每一项进行体积积分来离散化,最后这个方程就会变成这样:

Vp表示网格单元体积,Ff为通量,SS表示网格单元的各个面的面矢量,ν为动力粘度,d表示网格单元体心之间的距离。上标n表示为当前的时间步(已知),上标r表示预测时间步(待求)。下标N表示相邻网格单元的速度,下标P表示当前网格单元的速度。

现在再把之前忽略的压力梯度项加回来:

这里的p就是p/rho,单位是 m^2/s^2 (这也是OpenFoam的一个特点,OpenFoam里p文件定义的压力默认就是p/rho。但是也不一定所有solver都是,有的solver里p就是p,所以最好用的时候检查一下量纲,否则怎么错的都不知道。)

可以看出来,这是一个线性方程,在某个时间n上,除了U^r不知道,其他变量都是知道的。那么这就是一个ax=b的问题,列一个矩阵就可以把这些U^r都算出来。

第二步是压力泊松方程:

现在我们有下一时刻的速度U^r了,我们还需要什么?想一想我们还缺什么?就光有下一时刻的速度U就能完成循环了吗?还缺什么?

对,我们还缺下一时刻的压力。所以循环还不能完成。

同时,我们也应当注意到,这样算出来的下一时刻的速度U^r应当是不准确的,因为只用到了动量方程,没有用到连续性方程,那么新预测出来的U^r就不会满足连续性方程,这肯定不行,所以我们要对U^r进行修正,使其满足连续性方程。

所以,piso循环的任务,就是通过循环不断修正U^r和p^r,并加入动量方程和连续性方程两大物理constraint,最终得到满足物理constraint的U^r和p^r,带入下一个timestep进行计算。

那么piso循环内是如何更新速度和压力的呢?

将动量方程带进连续性方程,我们可以得到一个浓缩的2合1方程式,即压力泊松方程(推导过程略):

其中HbyA是一个在openfoam中会经常出现的变量,所以我列出公式:

上标n代表第几个piso循环(每个timestep里有好几个piso循环,通过这样的循环得到正确的p,然后再进入下一个timestep)。这个变量是基于省略了压力项的动量方程算出来的速度预测值,用于更新变量。

可以看出,每个piso循环的最终结果,都应该满足这个2合1方程式,这样得到的预测结果物理上就说得通了,因为两个物理constraint都满足了嘛。

步骤:

1. 首先是获取第一个压力预测 p^{(1)}
这个值通过将第一步算出来的U^r带入压力泊松方程而求得:

泊松方程是一个数学上的常见形式,高中时候物理课上就学过这种公式,算电场用的。这种方程的特点就是方程的一边必须是已知数。这里因为U^r已知,所以 \nabla(HbyA^r) 也就已知了。这个方程就变成了一个经典的泊松方程。 p^{(1)} 就可以解出来了。

2. 然后是根据第一个压力预测修正速度

将动量方程改一下,带入第一个压力预测 p^{(1)}

可以看出,新的速度预测 U_p^{(1)} 是这个方程中的唯一未知量,那当然可以由这个式子求出来。OpenFoam里一般会化简一下,由HbyA表示:

3. 循环重复直到收敛

重复1和2的动作,直到新得到的速度严格满足动量方程。需要特别注意的是第二次带入时就要根据新预测的速度 U_p^{(1)} 来重新计算 \nabla(HbyA^r) 了。也就是用每次预测出来的速度替换U^r。

李东岳博士有另一种概括方法

应该是一个意思。


  • 【观察代码】

以上内容都是数学公式,下一步是观察代码,看看OpenFOAM是怎么弄的:

这张图里凝聚了我对icoFoam代码的理解,希望这张图片能够以清晰的方式上传。大家最好点大图观看。


到目前为止我大致了解了OpenFOAM里PISO算法的意思,但是还有几个小疑问:

  1. 看上去PISO loop里并没有改变UEqn这个变量,而HbyA这个变量又是由UEqn得到的,这么说每个PISO loop用的都是同样的HbyA变量咯?那为什么不把HbyA放到PISO循环外面减少计算量?HbyA应该是不变的吗?
  2. PISOloop不是只有一个循环吗?为什么代码里有两个循环?
  3. 看上去并没有图中验证“收敛”的步骤,而是通过固定的常数nCorr和nNonOrthCorr求出来的。这行吗?
  4. 似乎solve(UEqn == -fvc::grad(p));之后速度变量U就直接undate了?C++中可以又这样的操作吗?U又不是UEqn的子变量,这行代码也没有赋值给U,怎么就update了?还有其他类似情况吗?


  • 【数据流动】

我们现在跑一个1*2*3网格的case。就跑自带的cavity case。

cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity ~/Desktop/

把blockMeshDict改成1*2*3的:

blocks
(
    hex (0 1 2 3 4 5 6 7) (2 3 1) simpleGrading (1 1 1)
);

先跑一跑,看看网格形状:

blockMesh
icoFoam
paraFoam

可以看出,一共6个cell,id顺序从左往右,从下到上。

我只看两个timestep,从0到第一个timestep,和从第一个到第二个timestep。

1. 从0到第一个timestep

while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        #include "readPISOControls.H"
        #include "CourantNo.H"
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        solve(UEqn == -fvc::grad(p));

可以看出这段代码的目的就是把动量方程放进去,solve命令算出第一个预测速度U^r。

我的第一个疑问是就如刚才提出的问题4:似乎solve(UEqn == -fvc::grad(p));之后速度变量U就直接update了?C++中可以有这样的操作吗?U又不是UEqn的子变量,这行代码也没有赋值给U,怎么就update了?还有其他类似情况吗?

验证了一下,在solve语句的前后都输出了一次U,结果两次果然不同,这说明确实solve语句更新了速度变量U。

这是一个很大的问题,说明solve这个语句可以修改这条语句外的变量而不用特意声明,意味着我光看一个solve语句并不能知道他改变了哪些外面的变量,这太危险了,怎么会有这么危险的表达方式,万一solve里面出了个bug,把所有变量都改了,那不是乱了套了。

冷静了一下,我猜测这里可能用到了指针(pointer)或者引用(reference)中的某一个,从而导致了任何对U的修改都直接作用于其本身,而不用创建临时变量,从而提高代码效率。

问了些大佬,都说:“这种操作在C里很正常,C中一般都是在一个函数中传入一个量进行修改。反而直接返回,对于不了解内存管理的人比较危险”。唉,都怪我过去贪玩,不爱用C,就爱matlab和python这种傻瓜代码,现在遇到C就傻了。

这里复制一个reference在函数参数传递中的例子

例三: 
           class MyClass 
           { 
               public: 
                   int a; 
                   void method(); 
           }; 
           void fun1(MyClass mc) 
           { 
                 mc.a=40; 
                 mc.method(); 
           } 
           void fun2(MyClass* mc) 
           { 
                 mc->a=60; 
                 mc->method(); 
           } 
           void fun3(MyClass& mc) 
           { 
                 mc.a=80; 
                 mc.method(); 
           } 
           MyClass   myclass; 
           void fun1(MyClass); 
           void fun2(MyClass*); 
           void fun3(MyClass&); 
           fun1(myclass);      //执行完函数调用后,myclass.a=20不变。 
           fun2(&myclass);     //执行完函数调用后,myclass.a=60,改变了。
           fun3(myclass);      //执行完函数调用后,myclass.a=80,改变了。
           //注意fun1和fun3的实参,证明了:使用对象和使用对象的引用,在语法格式上是一样的 

icoFoam.C中没看到&和*的符号,所以我估计应该是用fun3那种方法传递参数的吧。

现在的问题是:是UEqn中引用的U,还是solve中引用的U?如何引用的?我需要看到代码,才能确定OpenFOAM确实是使用了引用/指针的方式update的U值。

参考介红恩博主和苏军伟博主的两篇博客:

OpenFOAM常用类的一些总结_Jay_见贤思齐_新浪博客blog.sina.com.cn
OpenFOAM中的神奇方程定义方式的背后_苏军伟_新浪博客blog.sina.com.cn

fvMatrix中包含两组数据,一组是A,一组是一部分b,比如明显fvm::ddt(U)会产生一部分需要放到右边的系数,那就是b的组成部分。所以我在图上用[1/2]和[2/2]表示fvMatrix由两部分组成,一部分构成了A,一部分构成了b。

看上去这个逻辑很简单,由UEqn得到A和部分b,再由fvc::grad(p)得到剩下的b,最后solve把两个b相加得到完整的b,再执行一下 x=a^{-1}b 求得x就行了。

感觉没有能引用U的地方啊?UEqn里只储存了A和b,尚未对U进行计算,就不可能对其赋值;solve函数会对U赋值,但solve函数的输入变量里又没有U。

想了一会儿,这个逻辑里肯定有问题,应该是“UEqn里只储存了A和b”这句话有错,UEqn里应该不止储存了A和b,还储存了U的地址。solve函数应该执行了这样一个操作(我猜的):

solve(fvVectorMatrix& fvM)
{
    A = fvM.A();
    b = fvM.b();
    fvM.refU = b/A;   //将算出来的x值赋给UEqn中的引用变量refU,也就等于赋给U了
}

而fvVectorMatrix的结构应该是这样的:

class fvVectorMatrix:
{
matrix A;
vector b;
volVectorField& refU;  //UEqun所在的fvVectorMatrix类中有这么一个引用变量refU
}

去代码里找找看是不是这样:

~/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

template<class Type>
class fvMatrix
:
    public refCount,
    public lduMatrix
{
    // Private data

        //- Const reference to GeometricField<Type, fvPatchField, volMesh>
        //  Converted into a non-const reference at the point of solution.
        const GeometricField<Type, fvPatchField, volMesh>& psi_;

在fvMatrix.H的定义中可以看到,他的子变量中确实有一个引用变量psi_,我认为这就是引用U的地方;

关于solve函数,我在这里写了怎么追踪solve函数的过程:OpenFOAM中的动态多态(Info,solve)

fvMatrixSolve.C

template<class Type>
  104 Foam::solverPerformance Foam::fvMatrix<Type>::solveSegregated
  105 (
  106     const dictionary& solverControls
  107 )
  108 {
...
  118     GeometricField<Type, fvPatchField, volMesh>& psi =
  119        const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
...
  211         psi.internalField().replace(cmpt, psiCmpt);
  212         diag() = saveDiag;

可以清楚得看出,solve文件里对场的赋值,就是直接赋值在原来的场上。

(solveSegregated是solve函数里最终干活的函数,他在118行建立了一个对psi_场的引用psi,这个被引用的psi_就是刚才我们在fvMatrix.H文件最后一行看到的;之后211行对psi的操作用了replace这个词,那意思很明显了,就是把原来的场替换掉。)


至此,我们基本证明了solve里面确实是使用了引用/指针的方式来update的U值。接下来看看矩阵是怎么变化的:

修改源文件,让其输出solve公式中的变量:

fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        Info<< "U = " << U << endl;
        Info<< "phi = " << phi << endl;
        Info<< "fvm::ddt(U) = " << fvm::ddt(U) << endl;
        Info<< "fvm::div(phi, U) = " << fvm::div(phi, U) << endl;
        Info<< "- fvm::laplacian(nu, U) = " << - fvm::laplacian(nu, U) << endl;
        Info<< "-fvc::grad(p) = " << -fvc::grad(p) << endl;
        Info<< "UEqn = " << UEqn << endl;
        Info<< "UEqn + fvc::grad(p) = " << UEqn + fvc::grad(p) << endl;

        solve(UEqn == -fvc::grad(p));

        Info<< "U = " << U << endl;

结果是:

U = dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }
    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    frontAndBack
    {
        type            empty;
    }
}

phi = dimensions      [0 3 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    movingWall
    {
        type            calculated;
        value           uniform 0;
    }
    fixedWalls
    {
        type            calculated;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
        value           nonuniform 0();
    }
}

fvm::ddt(U) = false true false 6(0.00333333 0.00333333 0.00333333 0.00333333 0.00333333 0.00333333)
fvm::div(phi, U) = true true true 7{-0}6{0}7{0}
- fvm::laplacian(nu, U) = false true true 6(0.000216667 0.000216667 0.000366667 0.000366667 0.000216667 0.000216667)7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
-fvc::grad(p) = dimensions      [0 1 -2 0 0 0 0];

internalField   uniform (-0 -0 -0);

boundaryField
{
    movingWall
    {
        type            zeroGradient;
    }
    fixedWalls
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

UEqn = true true true 7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)6(0.00355 0.00355 0.0037 0.0037 0.00355 0.00355)7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
[0 4 -2 0 0 0 0]
6{(0 0 0)}

3
(
2((0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
8((0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
0()
)


3
(
2((0.0003 0 0) (0.0003 0 0))
8{(0 0 0)}
0()
)


UEqn + fvc::grad(p) = true true true 7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)6(0.00355 0.00355 0.0037 0.0037 0.00355 0.00355)7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 4.42335e-06, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
U = dimensions      [0 1 -1 0 0 0 0];

internalField   nonuniform List<vector> 6((0.000117175 0 0) (0.000117161 0 0) (0.00305958 0 0) (0.0030594 0 0) (0.0767127 0 0) (0.0767127 0 0));

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }
    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    frontAndBack
    {
        type            empty;
    }
}

这里有3个疑问:

(1). fvm::ddt(U) = false true false这里我不太清楚前三个true和false是什么意思?

(2).为什么同样是fvMatrix类,fvm::ddt(U)不会输出fvm.dimensions_之类的信息,而UEqn会?

(3).ddt(U)应该要包涵A和部分b,为什么这里只输出6个数?这6个数是什么?

根据OpenFOAM中的动态多态(Info,solve)一文中的推理:

fvMatrix.C:

// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //

template<class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
{
    os  << static_cast<const lduMatrix&>(fvm) << nl
        << fvm.dimensions_ << nl
        << fvm.source_ << nl
        << fvm.internalCoeffs_ << nl
        << fvm.boundaryCoeffs_ << endl;

    os.check("Ostream& operator<<(Ostream&, fvMatrix<Type>&");

    return os;
}

IduMatrix.C:

Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum)
{
    Switch hasLow = ldum.hasLower();
    Switch hasDiag = ldum.hasDiag();
    Switch hasUp = ldum.hasUpper();

    os  << hasLow << token::SPACE << hasDiag << token::SPACE
        << hasUp << token::SPACE;

    if (hasLow)
    {
        os  << ldum.lower();
    }

    if (hasDiag)
    {
        os  << ldum.diag();
    }

    if (hasUp)
    {
        os  << ldum.upper();
    }

    os.check("Ostream& operator<<(Ostream&, const lduMatrix&");

    return os;
}

(1). false true false分别是确定这个矩阵有没有lower,diag,upper这三个矩阵的基本组成部分。然后如果有,就依次输出来。那么可以看出,fvm::ddt(U)形成的矩阵就只有对角线上有值,所以只有diag是true,只输出一列数据。

(2). 原因可能是Ueqn是被fvVectorMatrix正式定义的,定义的时候constructor会给其dimensions_之类的东西赋值,而直接用ddt返回的就不正式,不会给dimensions_赋值。其实这块我也搞不太清楚,我现在知道的是fvm::ddt(U)返回的并不是真正的fvVectorMatrix类,而是一个叫做:

class Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >

的类。在我的理解里Foam::fvMatrix<Foam::Vector<double> >就是fvVectorMatrix,应该是前面的Foam::tmp起到了作用。查了一下,doxygen说class Foam::tmp是:

A class for managing temporary objects.

以我的理解,可能因为是管理临时对象,所以功能不够全面,不包含dimensions_之类的东西,需要先用fvVectorMatrix把这个临时对象固定住,才能正确输出。

(3). 这里只显示A不显示b应该是因为输出方式不对:

fvVectorMatrix ddt_tmp(fvm::ddt(U));
Info<< "fvVectorMatrix ddt_tmp(fvm::ddt(U)); = " << ddt_tmp << endl;

这样先create一个fvVectorMatrix类的变量ddt_tmp,再输出,就全面了:

fvVectorMatrix ddt_tmp(fvm::ddt(U)); = false true false 6(0.00333333 0.00333333 0.00333333 0.00333333 0.00333333 0.00333333)
[0 4 -2 0 0 0 0]
6{(0 0 0)}

3
(
2{(0 0 0)}
8{(0 0 0)}
0()
)


3
(
2{(0 0 0)}
8{(0 0 0)}
0()
)

看来还是我的调用方式太粗暴,以后调用变量还是应该像UEqn那样,先包装一层外壳,再让他输出,遵守OpenFOAM的规矩。就像这样:

 fvVectorMatrix ddt_tmp(fvm::ddt(U));
        Info<< "fvVectorMatrix ddt_tmp(fvm::ddt(U)); = " << ddt_tmp << endl;

        fvVectorMatrix div_tmp(fvm::div(phi, U));
        Info<< "fvVectorMatrix div_tmp(fvm::div(phi, U)); = " << div_tmp << endl;

        fvVectorMatrix laplacian_tmp(- fvm::laplacian(nu, U));
        Info<< "fvVectorMatrix laplacian_tmp(- fvm::laplacian(nu, U)); = " << laplacian_tmp << endl;

        Info<< "-fvc::grad(p) = " << -fvc::grad(p) << endl;
        Info<< "UEqn = " << UEqn << endl;
        fvVectorMatrix UEqn_gradp
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + fvc::grad(p)
        );
        Info<< "UEqn_gradp = " << UEqn_gradp << endl;

        solve(UEqn == -fvc::grad(p));
        Info<< "U = " << U << endl;

解决完这些破事之后,我要分析一下fvVectorMatrix中的矩阵A和b分别是什么变量,怎么储存的。


首先看一下最后的UEqn_gradp变量,这个变量就是整个公式形成的矩阵,把这个公式里的A逆和b相除应该就能得到速度U的新值,正好跟新的U值对比一下。

以下是UEqn_gradp的输出结果:

UEqn_gradp = true true true 7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)6(0.00355 0.00355 0.0037 0.0037 0.00355 0.00355)7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
[0 4 -2 0 0 0 0]
6{(0 0 0)}

3
(
2((0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
8((0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
0()
)


3
(
2((0.0003 0 0) (0.0003 0 0))
8{(0 0 0)}
0()
)

以上信息肯定包涵矩阵A和b,只是不知道是哪一个。

翻了几篇文献:

OpenFOAM guide/Matrices in OpenFOAMopenfoamwiki.net

还是没看懂,图太少,很多东西讲得很含糊。而且我总觉得少了点什么,U的公式应该是3个公式,因为速度有三个方向,所以未知数应该有3*6=18个,但是UEqn_gradp输出的结构似乎是个6*6的矩阵而不是18*18的,这里不太理解。

https://www.slideshare.net/fumiyanozaki96/openfoam-32087641www.slideshare.net

这篇就清楚多了,日本人干事就是细心!幸亏我本科的时候修了二外日语,没想到在这种时候派上用场了!

矩阵所代表的空间关系
IduMatrix的结构-diag,upper,lower
Iduaddressing

上图翻译:upper,lower,diag三个变量只容纳了矩阵种的非零元素,因此只看这三个变量的话,是不知道这些数字是在矩阵种的哪行哪列的。

这个位置的信息,就储存在“IduAddressing”这个类里。

具体来说,内部的每个face都有其编号,这个face两侧的两个control volume(CV)也有编号。IduAddressing就是把这些对应关系保存下来的一个类。

比如说左图里,7号face连接了3号CV和5号CV,那么就定义:

lowerAddr[7]=3
upperAddr[7]=5
举个例子

看完这些我大概明白了:

UEqn_gradp = true true true 
7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
6(0.00355 0.00355 0.0037 0.0037 0.00355 0.00355)
7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)

一开始输出的这三个变量767是对应lower,diag,upper三部分,因为一共有6个cell所以diag有6个元素;因为6个cell在空间上有7个接触面,会形成7个cell之间的关联,所以lower和upper三角都是7个元素。

现在不理解的还是刚才的问题,U的公式应该是3个公式,因为速度有三个方向,所以未知数应该有3*6=18个,但是UEqn_gradp输出的结构似乎是个6*6的矩阵,这里不太理解。


////////////////////////////////////////////////////////////////////////////

之前的代码写法可能有问题,fvVectorMatrix UEqn_gradp这个函数不应该直接加上fvc的grad(p),因为fvc返回的结构都不是矩阵,而是geometricField,所以不能相加。


    fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

一个一个项看吧:

  1. 先看ddt项:

几何形状看blockMesh:

convertToMeters 0.1;

vertices
(
    (0 0 0)
    (1 0 0)
    (1 1 0)
    (0 1 0)
    (0 0 0.1)
    (1 0 0.1)
    (1 1 0.1)
    (0 1 0.1)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (2 3 1) simpleGrading (1 1 1)
);

按照这个公式,矩阵系数应该= \frac{V_P}{\Delta t}=\frac{0.1/2*0.1/3*0.01}{0.005}=0.00333333333

source系数应该=0,因为初始速度都为0.但是这样就不好验证了:

为了方便验证,我们令初始内部场为 U_x=0.1,U_y=0.2

internalField   uniform (0.1 0.2 0);

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }

这样source项就应该为: \frac{V_P}{\Delta t}*U_p^n=\frac{0.1/2*0.1/3*0.01}{0.005}*(0.1,0.2,0)=0.00333*(0.1,0.2,0)=*(0.000333,0.000666,0)

对比一下输出结果:

fvVectorMatrix ddt_tmp(fvm::ddt(U)); = false true false 
6(0.00333333 0.00333333 0.00333333 0.00333333 0.00333333 0.00333333)
[0 4 -2 0 0 0 0]
6((0.000333333 0.000666667 0) (0.000333333 0.000666667 0) (0.000333333 0.000666667 0) (0.000333333 0.000666667 0) (0.000333333 0.000666667 0) (0.000333333 0.000666667 0))

3
(
2{(0 0 0)}
8{(0 0 0)}
0()
)


3
(
2{(0 0 0)}
8{(0 0 0)}
0()
)

确实是一样的。

可以看出,这里矩阵是6*6的,但是source是6*3的,说明openfoam的处理方式是这样的:

(6*6)*(6*3)=(6*3) 一共18个方程。比较疑惑的是如何保证同一个矩阵A可以用三次,还是每次都要重新算A,需要接下来继续验证。

相关源代码如下:

2. 接着看div项:

+ fvm::div(phi, U)

其中F是这样推导出来的:

外积的计算方式是:

因为 U_f 不能直接带入方程计算,所以还要对 U_f 操作一下,因为这里用的是均匀mesh,所以我们可以认为用的是中心差分法:

带回去原式,原式变成:

按照这个公式,对流项的计算也用到了高斯定理,求和符号代表了相邻cell对本cell的贡献之和。那么对任何一个相邻的cell而言,其与本cell之间的矩阵系数应该是:

\frac{F_f^n}{2}=\frac{U^n_fS_f}{2}=\frac{\phi_f}{2}

公式中的phi代表面上的通量,也就是速度乘以面积向量:

phi = dimensions      [0 3 -1 0 0 0 0];

internalField   nonuniform List<scalar> 
7(3.33333e-05 0.0001 0.0001 3.33333e-05 0.0001 0.0001 3.33333e-05);

boundaryField
{
    movingWall
    {
        type            calculated;
        value           uniform 0;
    }
    fixedWalls
    {
        type            calculated;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
        value           nonuniform 0();
    }
}

手算了一下,没问题。

x方向的phi=0.033*0.01*0.1=0.000033

y方向的phi=0.05*0.01*0.2=0.0001

一致的。

对比一下输出结果:

fvVectorMatrix div_tmp(fvm::div(phi, U)); = true true true 
7(-1.66667e-05 -5e-05 -5e-05 -1.66667e-05 -5e-05 -5e-05 -1.66667e-05)
6(6.66667e-05 3.33333e-05 1.66667e-05 -1.66667e-05 -3.33333e-05 -6.66667e-05)
7(1.66667e-05 5e-05 5e-05 1.66667e-05 5e-05 5e-05 1.66667e-05)
[0 4 -2 0 0 0 0]
6{(0 0 0)}

3
(
2{(0 0 0)}
8{(0 0 0)}
0()
)


3
(
2{(-0 -0 -0)}
8{(-0 -0 -0)}
0()

矩阵元素是一致的。

但奇怪的是这里source term居然都为0?

想了想,source term不是 \frac{F_f^nU_p^n}{2} ,而是 \sum_{}^{}{\frac{F_f^nU_p^n}{2}} ,所以source=0的原因是两个cell互相的影响抵消了。因为每一对cell所产生的两个source大小一样,正负号相反,所以正好抵消。

internalCoeffs和boundaryCoeffs居然也都为0?

想了一下,两个系数为0代表边界条件没有贡献。这是因为边界上通量为0,所以为0是对的。

对比了一下离散公式:

dxdydz的部分应该是“隐藏”进了通量phi的计算中。所以这一项也是一个矩阵可以用3次的。

相关源代码如下:

3. 再看laplacian项,之前追踪过一维的:

陈与论:用GdbOF追踪laplacian函数在OpenFOAM中的定义过程zhuanlan.zhihu.com图标
- fvm::laplacian(nu, U)

按照这个公式,相邻cell之间的矩阵系数就是-\frac{\nu|S_f|}{|d|}

nu是0.01,d是0.05或0.033,Sf刚刚算过了,是0.00033或0.0005. 则系数是-0.00015或-0.000066。

source也是都抵消了,都是0.

对比一下输出结果:

fvVectorMatrix laplacian_tmp(- fvm::laplacian(nu, U)); = false true true 
6(0.000216667 0.000216667 0.000366667 0.000366667 0.000216667 0.000216667)
7(-6.66667e-05 -0.00015 -0.00015 -6.66667e-05 -0.00015 -0.00015 -6.66667e-05)
[0 4 -2 0 0 0 0]
6{(-0 -0 -0)}

3
(
2((0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
8((0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.000133333 0.000133333 0.000133333) (0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
0()
)


3
(
2((0.0003 0 0) (0.0003 0 0))
8{(0 0 0)}
0()
)

一致的。

边界条件系数凭经验看也是对的,跟矩阵系数应该是2倍的关系,原因是最靠近壁面的cell中心离壁面的具体是半个cell,所以乘以2。其他的我没有细算,应该是正确的,我懒得一个个对了。

两个固定的cell之间,|d|应该也是固定的。所以这一项也是一个矩阵可以用3次的。

相关源代码如下:

4. 合并

      fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        solve(UEqn == -fvc::grad(p));

刚才我们算的1~3项,就已经把所有需要算的矩阵算完了,现在跟右边的压力梯度项联合起来,就能联立方程组求预测速度 U_p^r 了。

压力梯度全是0,因为我设定边界条件都是zerogradient,所以甚至可以当这一项不存在。

-fvc::grad(p) = dimensions      [0 1 -2 0 0 0 0];

internalField   uniform (-0 -0 -0);

boundaryField
{
    movingWall
    {
        type            zeroGradient;
    }
    fixedWalls
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

我们把1-3算出来的矩阵相加,跟OpenFOAM输出的UEqn对比:

可以看出,相加得到的A和b,与最终的A和b是一样的。

但是到这里还不能松懈,因为还有internalCoeffs和boundaryCoeffs没有比较。OpenFOAM不是简单的拿A/b算,而是在算的时候再把boundary condition考虑进矩阵。具体到代码里就是internalCoeffs和boundaryCoeffs,区别是前者加到矩阵A里,后者加到source项b里。

其实也不用比了,三项中只有laplacian项有internalCoeffs和boundaryCoeffs,直接比较laplacian和UEqn就行了:

laplacian:

3
(
2((0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
8((0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.0003 0.0003 0.0003) 
(0.0003 0.0003 0.0003))
0()
)


3
(
2((0.0003 0 0) (0.0003 0 0))
8{(0 0 0)}
0()
)

UEqn:

3
(
2((0.0003 0.0003 0.0003) (0.0003 0.0003 0.0003))
8((0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.000133333 0.000133333 0.000133333) 
(0.0003 0.0003 0.0003) 
(0.0003 0.0003 0.0003))
0()
)


3
(
2((0.0003 0 0) (0.0003 0 0))
8{(0 0 0)}
0()
)

一模一样。

这样,基本就说明全过程都在我们掌握之中,没有地方动了其他手脚。

现在算A/b:

先把边界条件加进A和b里,得到最终的矩阵A_f和b_f:

b fvMatrixSolve.C:189
r
pfvmatrixfull this matrix.txt

可惜的是我只找到了OpenFOAM里算矩阵的位置,却没找到算source的位置。Ueqn_f的source是我手算的出来的。

现在我们用matlab算一算预测速度,看看跟OpenFOAM算的是否一样:

对比OpenFAOM结果:

U = dimensions      [0 1 -1 0 0 0 0];

internalField   nonuniform List<vector> 6(
(0.085768 0.171434 0) 
(0.0872187 0.174331 0) 
(0.0966623 0.189262 0) 
(0.0984789 0.192758 0) 
(0.167387 0.180687 0) 
(0.170292 0.1839 0)
);

一致的,说明就是这么算的。

///////////////////////////////////////////////////////////////////////////

到目前为止,我们只做了piso循环的第一步,就是预测速度,接下来我们需要走之后的步骤。在此之前,我们先回顾一下piso循环的基本流程:

说起来也很好理解,两个方程AB和两个未知数ab。先用上一个时间步长的未知数值ab,解方程A,得到未知数a的预测值,然后把a的预测值带入到方程B,得到未知数b的预测值,然后再把b带回A,得到a的新预测值。之后就是重复把新生成的预测值轮流带到方程里,直到预测值不再变化。这应该是一种很自然而然的想法,换一个时空,给足够的时间,也有人能想出来。

编辑于 2018-01-25

文章被以下专栏收录