关于 phi 和 fvVectorMatrix 的两个问题



  • 刚开始学习openfoam,遇到两个问题, 前来请教。

    我在cavity的算例里面建立了4*4的网格,如下

    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) (4 4 1) simpleGrading (1 1 1)
    );
    

    然后输出了flux phi的值

    Info<< " phi " << phi << nl <<endl;
    
    Info<< " phi size " << phi.size() << nl <<endl;
    

    得到如下结果:

    phi dimensions [0 3 -1 0 0 0 0];
    
    internalField nonuniform List<scalar>
    24
    (
    -2.28233e-06
    2.28233e-06
    -3.2402e-06
    9.57866e-07
    -2.28434e-06
    -9.55866e-07
    -2.28433e-06
    -3.94416e-06
    6.22649e-06
    -5.34147e-06
    2.35517e-06
    -3.95016e-06
    -2.34716e-06
    -6.2345e-06
    -6.72257e-06
    1.29491e-05
    -7.99247e-06
    3.62507e-06
    -6.71657e-06
    -3.62306e-06
    -1.29511e-05
    1.29491e-05
    1.65741e-05
    1.29511e-05
    )
    ;
    
    boundaryField
    {
    movingWall
    {
    type calculated;
    value uniform 0;
    }
    fixedWalls
    {
    type calculated;
    value uniform 0;
    }
    frontAndBack
    {
    type empty;
    value nonuniform 0();
    }
    }
    
    
    phi size 24
    
    

    我的第一个问题是, 为什么phi只在internal face上面赋值了, 边界上的face不用计算吗?在 fvm::div(phi, U) 也用到了phi值, 边界上的flux是怎么考虑的?

    第二个问题是, 我声明了一个新的fvVectorMatrix

    fvVectorMatrix UEqn_time_term
    (
    fvm::ddt(U)
    );
    

    查看源项的值

    Info<< " UEqn_time_term.s " << UEqn_time_term.source() << nl << endl;
    

    输出结果:

    UEqn_time_term.s
    16
    (
    (-5.83396e-06 5.74984e-06 0)
    (-1.42259e-05 2.43986e-06 0)
    (-1.42171e-05 -2.34961e-06 0)
    (-5.92331e-06 -5.6699e-06 0)
    (-1.00527e-05 2.19123e-05 0)
    (-2.40783e-05 8.37173e-06 0)
    (-2.41399e-05 -8.28764e-06 0)
    (-1.01346e-05 -2.19111e-05 0)
    (-1.39378e-05 4.90818e-05 0)
    (-3.80295e-05 1.57254e-05 0)
    (-3.80131e-05 -1.56865e-05 0)
    (-1.39231e-05 -4.90753e-05 0)
    (9.96865e-05 3.28844e-05 0)
    (7.57897e-05 9.67425e-06 0)
    (7.57862e-05 -9.65121e-06 0)
    (9.97123e-05 -3.2868e-05 0)
    )
    

    这个fvmatrix的源项不是应该为零吗?为什么会有值呢?

    非常感谢!!



  • Hi,

    我的第一个问题是, 为什么phi只在internal face上面赋值了, 边界上的face不用计算吗?在 fvm::div(phi, U) 也用到了phi值, 边界上的flux是怎么考虑的?

    phi是通量值,只在面上有值,也就是说只有面有通量。看这个:http://www.dyfluid.com/flux.html 边界处的通量计算而来。要么给定,要么通过边界条件计算。

    这个fvmatrix的源项不是应该为零吗?为什么会有值呢?

    不能是0. 是零的话:AU=0。你的速度都是零。具体的在离散过程中:
    0_1474522774670_捕获.JPG
    对时间的离散调用了显性的值。



  • 东岳老师,非常感谢您的回复!
    对于第一问题,我又进行了尝试,下面是我的理解,如果有不对的地方,还请您指正。
    边界上的phi可以由边界条件得到,在cavity这个算例中,边界上没有法向的速度,因而phi全部为‘0’。而在internal face上的phi,可以由已知的速度场进行差值得到face上的速度,再乘以该face的面积和法向,就得到了phi。

    对于第二个问题,我又输出了U的值。

    Info<< "   U   " << U << nl << endl;
    

    得到

       U   dimensions      [0 1 -1 0 0 0 0];
    
    internalField   nonuniform List<vector>
    16
    (
    (-0.00109257 0.00108567 0)
    (-0.00260853 0.000420926 0)
    (-0.00260827 -0.000416558 0)
    (-0.0010966 -0.00108186 0)
    (-0.00181643 0.00400068 0)
    (-0.00421015 0.00141225 0)
    (-0.00421276 -0.00140849 0)
    (-0.00182014 -0.00400052 0)
    (-0.0033937 0.00937929 0)
    (-0.00774594 0.0025887 0)
    (-0.00774527 -0.00258687 0)
    (-0.00339316 -0.00937899 0)
    (0.0210076 0.00646146 0)
    (0.0145416 0.00159495 0)
    (0.0145415 -0.00159377 0)
    (0.0210087 -0.00646063 0)
    )
    ;
    
    boundaryField
    {
        movingWall
        {
            type            fixedValue;
            value           uniform (1 0 0);
        }
        fixedWalls
        {
            type            noSlip;
        }
        frontAndBack
        {
            type            empty;
        }
    }
    

    发现 UEqn_time_term.source() 里面的值对应为速度的值除以时间步长再乘以对应cell的体积。(时间步长为0.001,体积为0.01乘以0.025乘以0.025)



  • 0_1474590590150_捕获.JPG
    这一项吧?对应icoFoam解析公式6的第一项。
    这个研究思路很好。小伙子我看好你。:sunglasses: 哈哈

    边界上的phi可以由边界条件得到,在cavity这个算例中,边界上没有法向的速度,因而phi全部为‘0’。而在internal face上的phi,可以由已知的速度场进行差值得到face上的速度,再乘以该face的面积和法向,就得到了phi。

    是的



  • @李东岳

    东岳老师,我在研究完icoFoam中的 ‘fvm::ddt(U)’ 和 ‘fvm::div(phi, U)’ 之后又进而研究了‘fvm::laplacian(nu, U)’,有如下几个问题没有弄明白。

    为了看的清楚,附上网格编号的图。

    0_1474592981661_grid.JPG

    我分别定义了不同的 fvVectorMatrix

            fvVectorMatrix UEqn_time_term
            (
               fvm::ddt(U)
            );
    
            fvVectorMatrix UEqn_advection_term
            (
               fvm::div(phi, U)
            );
    
            fvVectorMatrix UEqn_diffusion
            (
               fvm::laplacian(nu, U)
            );
    

    对他们进行输出之后发现 UEqn_time_term 和 UEqn_advection_term 的diag() 和A() 每个元素都是相差一个cell 的体积

    UEqn.A() = UEqn.diag()/mesh.V()
    

    但是UEqn_diffusion的diag() 和A()只有在非边界的单元才满足上述关系。如下

      UEqn_diffusion.diag
    16
    (
    -0.0002
    -0.0003
    -0.0003
    -0.0002
    -0.0003
    -0.0004 \\这个满足上面的关系
    -0.0004 \\这个满足上面的关系
    -0.0003
    -0.0003
    -0.0004 \\这个满足上面的关系
    -0.0004 \\这个满足上面的关系
    -0.0003 
    -0.0002
    -0.0003
    -0.0003
    -0.0002
    )
    
    
      UEqn_diffusion.A  dimensions      [0 0 -1 0 0 0 0];
    
    internalField   nonuniform List<scalar>
    16
    (
    -96
    -80
    -80
    -96
    -80
    -64
    -64
    -80
    -80
    -64
    -64
    -80
    -96
    -80
    -80
    -96
    )
    ;
    
    boundaryField
    {
        movingWall
        {
            type            extrapolatedCalculated;
            value           nonuniform List<scalar> 4(-96 -80 -80 -96);
        }
        fixedWalls
        {
            type            extrapolatedCalculated;
            value           nonuniform List<scalar>
    12
    (
    -96
    -80
    -80
    -96
    -96
    -80
    -80
    -96
    -96
    -80
    -80
    -96
    )
    ;
        }
        frontAndBack
        {
            type            empty;
        }
    }
    

    我的第一个问题是 在UEqn_diffusion中.A()是怎么计算的?

    我在输出UEqn_diffusion的时候也发现它与其他的输出结果不同,在网上搜了半天也没有得到答案。

    Info<< "  UEqn_diffusion  " << UEqn_diffusion << nl << endl;
    

    输出结果如下:

      UEqn_diffusion  false true true
    
    16 \\这个应该是diag()
    (
    -0.0002
    -0.0003
    -0.0003
    -0.0002
    -0.0003
    -0.0004
    -0.0004
    -0.0003
    -0.0003
    -0.0004
    -0.0004
    -0.0003
    -0.0002
    -0.0003
    -0.0003
    -0.0002
    )
    
    24 \\这个应该是upper
    (
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    0.0001
    )
    
    [0 4 -2 0 0 0 0]
    16{(0 0 0)} \\我猜这个是source,这里面没有源项,4by4的网格,应该都是零。如果不对,请您指正。
    
    3 \\从这个开始往下我不知道是什么了,其他的matrix这些项都是零,请问这些值都代表了什么?
    (
    4((-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002))
    
    12
    (
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    (-0.0002 -0.0002 -0.0002)
    )
    
    0()
    )
    
    
    3
    (
    4((-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0))
    12{(-0 -0 -0)}
    0()
    )
    

    第二个问题 写在了代码的注释里,主要是不明白这个输出的结果里面除了diag(), upper() 和source()之外的那些项是什么?

    我在推倒UEqn_diffusion 的系数矩阵的时候发现,由边界条件得到边界上的速度值,那么最上面的滑动边界上的cell (编号为P)

    \begin{equation}
    S_{e}\nu_{e}\frac{u_{E}-u_{P}}{\Delta x} + S_{w}\nu_{w}\frac{u_{W}-u_{P}}{\Delta x} + S_{n}\nu_{n}\frac{u_{N}-u_{P}}{\Delta y} + S_{s}\nu_{s}\frac{u_{S}-u_{P}}{\Delta y}
    \end{equation}

    里面的$u_{N}$可以根据边界条件差值

    \begin{equation}
    u_{N} = 2 u_{n} - u_{P}
    \end{equation}

    其中$u_{n}$为上面滑动边界已知的速度。

    但是这个边界条件并没有在系数矩阵的diag()中体现出来,这也带来了我的第三个问题,这个边界条件在 UEqn_diffusion 是怎么体现出来的?是否与第二个问题有关?

    第四个问题就是也是又第三个问题引出来的,我试着求解了一下UEqn_diffusion,并输出了一下U的结果

    solve(UEqn_diffusion);
    
    Info<< "   U solution  " << U << nl << endl;
    

    但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同, 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。

    我的问题有点儿多,我也会继续查找相关的信息,期待您的解答!

    非常感谢!



  • @cfd-china 谢谢您!我是刚开始学Openfoam,希望能多多向前辈们请教!



  • @mengweilm425

    补充一下上面的帖子

    solve(UEqn_diffusion);
    
    Info<< "   U solution  " << U << nl << endl;
    

    结果如下:

    smoothSolver:  Solving for Ux, Initial residual = 0.921181, Final residual = 9.32208e-06, No Iterations 11
    smoothSolver:  Solving for Uy, Initial residual = 0.999778, Final residual = 8.66271e-06, No Iterations 8
       U solution  dimensions      [0 1 -1 0 0 0 0];
    
    internalField   nonuniform List<vector>
    16
    (
    (0.0189067 -1.2469e-08 0)
    (0.0430646 -3.74215e-08 0)
    (0.0430639 -4.90001e-08 0)
    (0.0189057 -2.83642e-08 0)
    (0.0703756 -3.73926e-08 0)
    (0.153354 -1.08185e-07 0)
    (0.153353 -1.33266e-07 0)
    (0.070374 -6.9139e-08 0)
    (0.179619 -4.88555e-08 0)
    (0.34663 -1.32978e-07 0)
    (0.34663 -1.48122e-07 0)
    (0.179618 -6.64227e-08 0)
    (0.481091 -2.81724e-08 0)
    (0.706929 -6.88449e-08 0)
    (0.706929 -6.63391e-08 0)
    (0.481091 -2.54916e-08 0)
    )
    ;
    
    boundaryField
    {
        movingWall
        {
            type            fixedValue;
            value           uniform (1 0 0);
        }
        fixedWalls
        {
            type            noSlip;
        }
        frontAndBack
        {
            type            empty;
        }
    }
    

    但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同, 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。感觉openfoam求解了这个方程,从UEqn_diffusion的系数和source()项为零来看,求解结果应该是零,但是显然是我想错了。如果是$u_{P}$为未知,$u_{W}$等相邻的cell的速度用的是上一个时间步的值,那么线性方程组只存在对角线的值,不需要迭代,直接就能出结果。在这个问题上没有想明白,还希望您能指点一下这个方程是怎么求解的。



  • 不好意思信息量太大。

    如果对时间离散fvm::ddt(T),其离散后的为对角阵。因此如果源项为0,AT=0则T也为零。并且不需要迭代求解。对于对流项fvm::div(phi,T),离散后非对角阵因此不一定T为0,且需要迭代求解。


  • 版主

    @李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b,而是指加在网格单元中心的力源项。



  • @wwzhao
    即使没有在方程等号右边添加源项的fvMatrix,比如:

        fvScalarMatrix TEqn
        (
            fvm::ddt(T) 
        );
    

    这个矩阵系统依然具有源项。显性离散的值全部进入了fvMatrixsource()

    P.S. 这个问题@mengweilm425 不是在一楼测试过了么?
    用一个网格单元测试:

    #include "fvCFD.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"
    
        #include "createTime.H"
        #include "createMesh.H"
    
        volScalarField T
        (
            IOobject
            (
                "T",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            1
        );
    
        fvScalarMatrix TEqn
        (
            fvm::ddt(T) 
        );
    
        Info<< "source" << TEqn.source() << endl;
    }
    


  • 针对@mengweilm425 提到的A()diag()。区别是A()=diag()/V;可测试:

    #include "fvCFD.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"
    
        #include "createTime.H"
        #include "createMesh.H"
    
        volScalarField T
        (
            IOobject
            (
                "T",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            1
        );
    
        fvScalarMatrix TEqn
        (
            fvm::ddt(T) 
        );
    
        Info<< "source" << TEqn.source() << endl;
        Info<< "source" << TEqn.A() << endl;
        Info<< "source" << TEqn.diag() << endl;//diag() equals TEqn.A() multiply mesh.V()
    }
    

    源代码:

    template<class Type>
    Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
    {
        tmp<volScalarField> tAphi
        (
            new volScalarField
            (
                IOobject
                (
                    "A("+psi_.name()+')',
                    psi_.instance(),
                    psi_.mesh(),
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                psi_.mesh(),
                dimensions_/psi_.dimensions()/dimVol,
                extrapolatedCalculatedFvPatchScalarField::typeName
            )
        );
    
        tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
        tAphi.ref().correctBoundaryConditions();
    
        return tAphi;
    }
    
    template<class Type>
    Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
    {
        tmp<scalarField> tdiag(new scalarField(diag()));
        addCmptAvBoundaryDiag(tdiag.ref());
        return tdiag;
    }
    

  • 版主

    @李东岳 我是这么理解的,fvm::ddt 反映的是上一时刻(或迭代步)对当前时刻(或迭代步)的物理量的影响,OpenFOAM把这种影响当作源项处理,因此将离散后的系数放在source()里面。


  • 版主

    @李东岳 纠正一下,lduMatrix::diag()fvMatrix<Type>::D()并不等价。相比lduMatrix::diag()fvMatrix<Type>::D()考虑了边界对系数矩阵的影响,具体代码对应这句:addCmptAvBoundaryDiag(tdiag.ref());



  • 有关diag()和D()严格来说是的。你可以测试一下然后给楼主提供更详细的信息。:sunglasses:


  • 版主

    @cfd-china 离散的具体实现可参阅相关源码。

    fvm::ddt(T)为例,调用的是

    template<class Type>
    tmp<fvMatrix<Type> >
    ddt
    (
        const GeometricField<Type, fvPatchField, volMesh>& vf
    )
    {
        return fv::ddtScheme<Type>::New
        (
            vf.mesh(),
            vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
        )().fvmDdt(vf);
    }
    

    这里调用了fv::ddtScheme<Type>::New,这个函数是一个runtime construction的实现,即在运行时读取字典文件中的离散格式并构造相应的对象:

    template<class Type>
    tmp<ddtScheme<Type> > ddtScheme<Type>::New
    (
        const fvMesh& mesh,
        Istream& schemeData
    )
    {
        if (fv::debug)
        {
            Info<< "ddtScheme<Type>::New(const fvMesh&, Istream&) : "
                   "constructing ddtScheme<Type>"
                << endl;
        }
    
        if (schemeData.eof())
        {
            FatalIOErrorIn
            (
                "ddtScheme<Type>::New(const fvMesh&, Istream&)",
                schemeData
            )   << "Ddt scheme not specified" << endl << endl
                << "Valid ddt schemes are :" << endl
                << IstreamConstructorTablePtr_->sortedToc()
                << exit(FatalIOError);
        }
    
        const word schemeName(schemeData);
    
        typename IstreamConstructorTable::iterator cstrIter =
            IstreamConstructorTablePtr_->find(schemeName);
    
        if (cstrIter == IstreamConstructorTablePtr_->end())
        {
            FatalIOErrorIn
            (
                "ddtScheme<Type>::New(const fvMesh&, Istream&)",
                schemeData
            )   << "Unknown ddt scheme " << schemeName << nl << nl
                << "Valid ddt schemes are :" << endl
                << IstreamConstructorTablePtr_->sortedToc()
                << exit(FatalIOError);
        }
    
        return cstrIter()(mesh, schemeData);
    }
    

    以常见的Euler格式为例,将会构造EulerDdtScheme类型的对象,并用其fvmDdt()方法返回的fvMatrix,看fvmDdt()的实现

    template<class Type>
    tmp<fvMatrix<Type> >
    EulerDdtScheme<Type>::fvmDdt
    (
        const GeometricField<Type, fvPatchField, volMesh>& vf
    )
    {
        tmp<fvMatrix<Type> > tfvm
        (
            new fvMatrix<Type>
            (
                vf,
                vf.dimensions()*dimVol/dimTime
            )
        );
    
        fvMatrix<Type>& fvm = tfvm();
    
        scalar rDeltaT = 1.0/mesh().time().deltaTValue();
    
        fvm.diag() = rDeltaT*mesh().Vsc();
    
        if (mesh().moving())
        {
            fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
        }
        else
        {
            fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
        }
    
        return tfvm;
    }
    

    可以看到fvm::ddt(T)产生的fvm.diag()为时间倒数乘以网格单元体积,其产生的源项fvm.source()为时间倒数乘以上一时刻的T乘以网格单元体积。

    对于fvm::div及fvm::laplacian可以做类似分析。

    插一句,fvm::div离散不会产生source(),fvm::laplacian会产生source()。



  • @wwzhao @李东岳 @cfd-china

    非常感谢大家的回复!我现在的阶段还看不好代码,但是我又对fvMatrix做了一些测试。

    还是如下的fvMatrix

        fvVectorMatrix UEqn_diffusion
        (
           fvm::laplacian(nu, U)
        );
    

    我发现在解这个方程的时候

    \begin{equation}
    [A][x]=[b]
    \end{equation}

    矩阵[A] 实际上是 UEqn.D() + UEqn.upper() + UEqn.lower(), 而没有使用UEqn.diag()。UEqn.D()在UEqn.diag()的基础上加上了边界条件对当前网格的影响。上面方程中的[b],也不是UEqn.source(),而是和边界条件中已知的速度相关的值。我直接求解得到的结果和openfoam算出来的结果是一样的。由于没有真正的去看代码,不知道Openfoam里面有没有相关的function可以直接输出上述$[A][x]=[b]$方程组中的[A]和[b],还请各位大神指点。

    再次感谢大家的帮助!


  • 版主

    @mengweilm425 $[A][x]=[b]$ 中的 [A] 可以看成是 internal field 离散后的系数矩阵和 boundary conditions 对系数矩阵的影响两部分之和,其中 internal field 离散之后的系数矩阵表示为 upper(),diag() 和 lower()。boundary condition 对系数矩阵的影响则放在 internalCoeffs_ 里面。[b] 也一样,可以分为 internal field 产生的源项 source(),及 boundary conditions 对源项的影响 boundaryCoeffs_。

    OpenFOAM中没有直接输出 [A] 和 [b] 的函数,关于如何将 boundary conditions 的影响加到 fvMatrix 里面,你可以参考 addBoundaryDiag,addCmptAvBoundaryDiag 及 addBoundarySource 这几个函数。

    @wwzhao关于 phi 和 fvVectorMatrix 的两个问题 中说:

    @李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b,而是指加在网格单元中心的力源项。

    上面的帖子描述不太准确,确切说 [b] 是 internal field 产生的源项 source() + boundary field 产生的源项 boundaryCoeffs_。



  • @wwzhao 谢谢您的回复,这就明朗多了。我再好好看看这部分代码。:happy:



  • @mengweilm425 您好 看了你的研究过程,觉得很好。我也想这样分步骤运行各个程序,以分清各个量含义。请问您是怎么做到的。我先编译了自己求解器,在文件中加入自己想输出的量,希望在输出中出现想看到的量。结果并没有看到想看的内容。请问您是如何。。



  • @vivian
    我猜测你的求解器名字没有替换?如果编译后的是icoFoam2,记得输入icoFoam2而不是原来的icoFoam…


登录后回复
 

与 CFD 中国 的连接断开,我们正在尝试重连,请耐心等待