关于 phi 和 fvVectorMatrix 的两个问题

mengweilm425

@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，且需要迭代求解。

wwzhao

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

李东岳

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

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

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

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;
}

wwzhao

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

wwzhao

lduMatrix::diag()和fvMatrix<Type>::D()并不等价。相比lduMatrix::diag()，fvMatrix<Type>::D()考虑了边界对系数矩阵的影响，具体代码对应这句：addCmptAvBoundaryDiag(tdiag.ref());。

CFD中文网

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

wwzhao

@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()。

mengweilm425

@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]，还请各位大神指点。

再次感谢大家的帮助！

wwzhao

@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_。

mengweilm425

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

vivian

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

CFD中文网

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

vivian

@cfd-china 您好，谢谢您的回复。是替换了的，是我给他速度赋值赋错啦。。

mengweilm425

@vivian @cfd-china 看来我来晚了:happy:

Wayne

厉害了我的哥哥姐姐们

CFD中文网

关于 phi 和 fvVectorMatrix 的两个问题