rhoCentralFoam求解器的问题

vivian

最近在看rhoCentralFoam求解器，有一些问题，请教各位：
在使用openfoam中rhoCentralFoam求解器时，我参照其中的例子结合自身算例，给定了边界条件及初始值，在计算4000步后，出现如下错误：

--> FOAM FATAL ERROR: 
Maximum number of iterations exceeded

    From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::perfectGas<Foam::specie> >; Type = Foam::sensibleInternalEnergy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>]
    in file /home/liyue/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-3.0.1/src/OSspecific/POSIX/printStack.C:218
#1  Foam::error::abort() at ~/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/error.C:249
#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/errorManip.H:85 (discriminator 4)
#3  Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::T(double, double, double, double (Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::*)(double, double) const, double (Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::*)(double, double) const, double (Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::*)(double) const) const at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/thermoI.H:66
#4  Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::TEs(double, double, double) const at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/thermoI.H:425
#5  Foam::sensibleInternalEnergy<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> >::THE(Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> const&, double, double, double) const at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/sensibleInternalEnergy.H:126
#6  Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::THE(double, double, double) const at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/thermoI.H:365
#7  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C:46 (discriminator 2)
#8  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ~/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C:143
#9  ? at ~/OpenFOAM/OpenFOAM-3.0.1/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C:232
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  ? at ??:?
Aborted (core dumped)

结合前面求解e方程残差无变化，我分析，这是由于e的求解修正出现了错误，为查找错误来源，我想在rhoCentralFoam的求解器中加断点，看求解过程中的数值传递与变化。

于是，我采用了debug模式在rhoCentralFoam.C的215行加入断点，想输出e,结果屏幕输出了大量信息，我应该如何输出才能得到想要的场信息呢？

如果我想梳理出程序运行的主线，搞清楚从场的离散点到通过离散方法，获得最终求解的代数方程矩阵，我认为最好的办法是从建立场开始，关注场中存储的离散点的信息变化。那么请问如何输出存储点压力，面通量信息的矩阵呢？

CFD中文网

@vivian
错误显示迭代步数超出了。

我采用了debug模式在rhoCentralFoam.C的215行加入断点，想输出e,结果屏幕输出了大量信息，我应该如何输出才能得到想要的场信息呢？

e是个场，你要输出场的话，就会输出所有的场的值，不清楚你要看什么，也可以单独输出场的最大值最小值看看越界否。如果一定要查看整个场的信息，输出到log再看，这样log文件会变得很大。

rhoCentralFoam &> log

如果我想梳理出程序运行的主线，搞清楚从场的离散点到通过离散方法，获得最终求解的代数方程矩阵，我认为最好的办法是从建立场开始，关注场中存储的离散点的信息变化。那么请问如何输出存储点压力，面通量信息的矩阵呢？

从离散的矩阵系数如fvScalarMatrix组建我们常规的稀疏矩阵不是一俩行代码可以搞定（估计不到100行？）因为openfoam里面存储的是稀疏矩阵，并用指针+指针。所以，可能得用点时间。

更简单的方法，是回到问题本身，很少是因为矩阵系统解不出来的原因，而是设置有问题导致发散。

Wayne

@cfd-china 厉害了我的哥

vivian

谢谢您的回复！:big_mouth:
结合tutorial中wedge15的算例，我的模拟模型为马赫数2的压缩拐角流动。入口设置了速度边界，出口为超声速出口。采用k-epsilon湍流模型，所需k,epsilon值根据经验公式给出。

压力场设定为：

internalField   uniform 100000;

boundaryField
{
    INLET
    
    {
        type            inletOutlet;
        inletValue      uniform 101325;
        value           uniform 101325;
    }
    

    OUTLET
    {
        type            waveTransmissive;
        field           p;
        phi             phi;
        rho             rho;
        psi             thermo:psi;
        gamma           1.3;
        fieldInf        100000;
        lInf            1;
        value           uniform 100000;
    }

    BOTTOM
    {
        type            zeroGradient;

    }
    TOP
    {
        type            zeroGradient;
    }
    frontAndBackPlanes
    {
        type            empty;
        
    }    
}

速度设定为：

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (600 0 0);

boundaryField
{
    INLET
    {
        type            supersonicFreestream;
        pInf            100000;
        TInf            300;
        UInf            (600 0 0);
        gamma           1.4;
        value           uniform (662 0 0);
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      uniform (662 0 0);
        value           uniform (662 0 0);
    }

    BOTTOM
    {
        type            zeroGradient;

    }
    TOP
    {
        type            symmetryPlane;
    }
    frontAndBackPlanes
    {
        type            empty;
        
    } 
    
}

温度场设置为：

internalField   uniform 300;

boundaryField
{
    INLET
    {
        type            inletOutlet;
        inletValue      uniform 300;
        value           uniform 300;
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      uniform 300;
        value           uniform 300;
    }

    BOTTOM
    {
        type            zeroGradient;

    }
    TOP
    {
        type            symmetryPlane;
    }
    frontAndBackPlanes
    {
        type            empty;
        
    }

求解器的设置是根据相似算例给出，这个的设置有什么原则吗？我设置的如下，还请多多指点：

solvers
{
    "rho.*"
    {
        solver          diagonal;
    }

    "p.*"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-12;
        relTol          0;
    }

    "(U|e).*"
    {
        $p;
        tolerance       1e-9;
    }

    "(nuTilda).*"
    {
        $p;
        tolerance       1e-10;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors      2;
    nNonOrthogonalCorrectors 0;
}

Q1：不知道我这样设置边界条件是否合理哪？
Q2：求解器的选择应该如何设置呢，有什么选择原则，请问有什么资料推荐码？
Q3：在外部参数设置找不出差错，又有发散产生时，请问如何debug呢？

liailei0627

你好，出现了跟你一样的问题，用rhoCentralFoam LES，请问问题解决了吗，可否指教下思路。

zhangdong

不知道大家有没有关注过thermophysicalProperties里面的设置。rhoCentralFoam 的wedge15的算例里面的mu设置为零。如果马赫数为2的时候，设置了湍流，那么thermophysicalProperties相关项也要发生变化。mu应该不是零，并且要随温度场发生变化。这种情况下rhoCentralFoam是否还合适不好说。我在sonicFoam里面看到有考虑湍流的，给出的入口速度是600。如果改为sonicFoam求解一下试试呢？

吴优

@zhangdong 如果mu随温度变化，transport就不选const而是sutherland，然后在transport中规定As，Ts取值

宝丁

@liailei0627 +1

zhangdong

@vivian 您能贴出thermoProperties的设置吗？我用OpenFOAM-6，rhoCentralFoam的tutorial里面wedge15Ma5的算例是马赫数5，该算例通过设置transport中Cp的选项，使得声速是1(https://openfoam.com/documentation/tutorial-guide/tutorialse6.php)，所以入口速度设置的是5，不知您是否把thermoProperties更改过来了？
另外，我用sonicFoam计算圆柱扰流，分别在Ma=1, Ma=3的情况下，进行计算，Re=1.0e+6，采用k-epsilon模型，可以计算。

Rachel0096

@vivian 我也遇到了一样的问题，请问您的问题解决了吗，可以指点一下吗？