pEqn.H中phiHbyA计算修正项的问题

yhdthu

在pEqn.H中，计算phiHbyA时会有一个修正项，是利用上一时间步通量和速度向面上的插值之差来修正本时间步的值，这一项会增加额外的耗散，好处是计算鲁棒性好

surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);

以上是代码和修正通量的表达式，但是有个问题是，当计算压力时，phiHbyA作为泊松方程的源项，修正项phi_2是会对压力的计算结果产生影响，比如之前所发的一个帖子传送门

经过我的测试，如果去掉这项，会得出比较合理的压力分布，但是随着计算进行，会发生发散现象，所以不知道这项真正的意义是什么，为什么要用上一个时间步的物理量修正本时间，对压力求解的影响如何避免？

李东岳

虽然方程和代码是一致的，fvc::ddtCorr(U, phi)并不是求解NS方程必须的项，而是一个稳定项，类似/max(x, 1e-6)。我对fvc::ddtCorr(U, phi)的数值影响比较感兴趣，但是并没有时间去做研究。如果你也感兴趣，相对于研究多相fvc::ddtCorr(U, phi)数值的影响，我建议你从icoFoam中的fvc::ddtCorr(U, phi)入手，研究他的数值作用。

yhdthu

@李东岳好的，多谢前辈回复，我会研究一下关于计算稳定性的问题，如果有需要我会跟您交流～

程迪

@yhdthu
这玩意儿涉及到面流量场和体速度场的相容性的问题，即插值后的U和phi是否一致。由于在很多分离算法（SIMPLE, PISO等）中，再考虑在同位网格条件下，phi和U由于更新有先后顺序之分，可能会导致U和phi不一致，所以需要修正。

$$
\delta \phi=\phi^{o}-U^o\cdot S_f\
K_c = 1-\min{(1,\frac{|\delta\phi|}{|\phi^o|+\epsilon})}\
\phi_{c}=\delta \phi\cdot K_c\cdot \frac 1 A\cdot \frac 1 {\Delta t}
$$

参考资料：https://zh.scribd.com/doc/48195039/ddtPhiCorr

ddtCorr和ddtPhiCorr之类的都是相互调用的。而且of41和of2x相差也特别大。深坑绕行。

yhdthu

@程迪你好，谢谢你的回复，我有以下疑问：

1.为什么这个修正放在了phiHbyA里面呢？看起来，这个phiHbyA由两部分组成，预测速度矩阵rAU*U.H()计算的flux(HbyA)，和上一时间步通量/速度修正项目。难道必须要和上一时间步产生关联么？

2.理论上说，是不是在每一个时间步收敛后，时间滞后效应基本消除，所得速度基本可以同时满足动量方程和压力方程，这个修正项应该为0？

3.我做了算例测试，发现这个修正项对网格大小很敏感，即如果网格比较密，这项会比较小；也做了同样网格下不同时间步的测试，发现如果时间步越小flux(HbyA)越小，此时修正项对求解压力泊松方程的影响就很大了，这项一直都不是0，请问这正常么？

4.我计算多相流问题，压力场一直不稳定，详情在最上方的传送门中，现在怀疑是由这项引起的，这个修正项的影响相当于在泊松方程里加入了一个源项，导致压力计算爆掉了，但是计算单相流又没事，请问问题是不是出在这呢？

谢谢

程迪

@yhdthu 在 pEqn.H中phiHbyA计算修正项的问题中说：

3.我做了算例测试，发现这个修正项对网格大小很敏感，即如果网格比较密，这项会比较小；也做了同样网格下不同时间步的测试，发现如果时间步越小flux(HbyA)越小，此时修正项对求解压力泊松方程的影响就很大了，这项一直都不是0，请问这正常么？

我现在也不清楚了，你可以考虑按照Jasak的做法，直接去掉这个项。

#Jasak的搞法

Jasak似乎玩了另外的花样。在foam-extend 4.0中，Jasak的icoFoam和OpenFOAM.org以及OpenFOAMplus的icoFoam就不一样了。

Jasak的icoFoam把UEqn分成了HUEqn和ddtUEqn两部分：
- HUEqn：对流和扩散项
- ddtUEqn：非定常项；
动量预测步：HUEqn + ddtUEqn == -grad(p)
PISO循环中：H从HUEqn.H()中产生！与时间步完全无关。
速度更新：

//https://github.com/Unofficial-Extend-Project-Mirror/foam-extend-foam-extend-4.0/blob/master/applications/solvers/incompressible/icoFoam/icoFoam.C +113
// Note: cannot call H(U) here because the velocity is not complete
// HJ, 22/Jan/2016
U = 1.0/(aU + ddtUEqn.A())*
(
    U*aU - fvc::grad(p) + ddtUEqn.H()
);

Jasak 认为这使得代码更加清晰了（时间导数算子中完全去掉了ddtPhiCorr和ddtCorr）。而且完全消除了时间步和松弛因子对结果的影响。结果对比参考文献[^1]

[^1] Numerics Improvements and Validation Results: FOAM-Extend Update on Work-in-Progress

还有一种是类似 F. Moukalled. L. Mangani. M. Darwish. The Finite. Volume Method in Computational. Fluid Dynamics. An Advanced Introduction with. OpenFOAM® and Matlab® 书中代码的搞法，分别加各种修正。

程迪

关于fvc::ddtCorr()和fvc::ddtPhiCorr()

对于非稳态的求解器，如pimpleFoam, pisoFoam和icoFoam中，phiHbyA在构造时除了用到HbyA 向的插值，还会用到fvc::ddtCorr()，这一项主要是为了修正非定常和松弛因子对解的影响。而且这一项还是和OpenFOAM版本历史相关的！

//OpenFOAM 5.x, 4.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +77
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA) //这里用了fvc::flux函数，而不是原来的interpolate 和 & 算符
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //用的是ddtCorr!
            );

//OpenFOAM 2.3.x, 2.4.x
//applications/solvers/incompressible/icoFoam/icoFoam.C  +73
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //用的是ddtCorr
            );
//OpenFOAM 2.2.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +73
            surfaceScalarField phiHbyA//用了新的变量名!
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf()) 
              + fvc::ddtPhiCorr(rAU, U, phi) //用的是ddtPhiCorr!
            );


//OpenFOAM 2.1.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +72
            phi = (fvc::interpolate(U) & mesh.Sf())
                + fvc::ddtPhiCorr(rAU, U, phi);

//foam-extend 4.0
//applications/solvers/incompressible/icoFoam/icoFoam.C +86
            U = HUEqn.H()/aU; //注意这里的HUEqn不包含非定常项
            phi = (fvc::interpolate(U) & mesh.Sf()); //所以这里也不包含。

太老的就不管了，因为根据OpenFOAM.org官网的更新说明，2.3.0版本时已经把ddtPhiCorr换成了ddtCorr。

首先找OpenFOAM 5.x的代码，看看ddtCorr()到底是干啥的：

//src/finiteVolume/finiteVolume/fvc/fvcDdt.C +183
template<class Type>
tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh>>
ddtCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField
    <
        typename flux<Type>::type,
        fvsPatchField,
        surfaceMesh
    >& phi
)
{
    return fv::ddtScheme<Type>::New //根据变量名返回tmp<ddtScheme<Type>>的临时对象
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    ).ref().fvcDdtPhiCorr(U, phi);
}

//这里涉及到RTS机制，
//为了简化，参考steadyState和经典的Euler格式是怎么做的

//src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C +305
//steadyState返回的是0，无修正，这很正常。
//返回值的量纲是：phi.dimensions()/dimTime
//对于不可压缩流应该是L^3/T^2
template<class Type>
tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
steadyStateDdtScheme<Type>::fvcDdtPhiCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const fluxFieldType& phi
)
{
    return tmp<fluxFieldType>
    (
        new fluxFieldType
        (
            IOobject
            (
                "ddtCorr(" + U.name() + ',' + phi.name() + ')',
                mesh().time().timeName(),
                mesh()
            ),
            mesh(),
            dimensioned<typename flux<Type>::type>
            (
                "0",
                phi.dimensions()/dimTime,
                Zero
            )
        )
    );
}

//Euler格式的修正
//src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +534
template<class Type>
tmp<typename EulerDdtScheme<Type>::fluxFieldType>
EulerDdtScheme<Type>::fvcDdtPhiCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const fluxFieldType& phi
)
{
    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); //时间增量的倒数

    fluxFieldType phiCorr
    (
        phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
    );

    return tmp<fluxFieldType>
    (
        new fluxFieldType
        (
            IOobject
            (
                "ddtCorr(" + U.name() + ',' + phi.name() + ')',
                mesh().time().timeName(),
                mesh()
            ),
            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
           *rDeltaT*phiCorr
        )
    );
}

//涉及到fvcDdtPhiCoeff系数，继续查找，发现这个函数位于ddtScheme中。
//src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C +151
template<class Type>
tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const fluxFieldType& phi,
    const fluxFieldType& phiCorr
)
{
    tmp<surfaceScalarField> tddtCouplingCoeff = scalar(1)
      - min
        (
            mag(phiCorr)
           /(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
            scalar(1)
        );

    surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();

    surfaceScalarField::Boundary& ccbf =
        ddtCouplingCoeff.boundaryFieldRef();

    forAll(U.boundaryField(), patchi)
    {
        if
        (
            U.boundaryField()[patchi].fixesValue()
         || isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
        )
        {
            ccbf[patchi] = 0.0;
        }
    }

    if (debug > 1)
    {
        InfoInFunction
            << "ddtCouplingCoeff mean max min = "
            << gAverage(ddtCouplingCoeff.primitiveField())
            << " " << gMax(ddtCouplingCoeff.primitiveField())
            << " " << gMin(ddtCouplingCoeff.primitiveField())
            << endl;
    }

    return tddtCouplingCoeff;
}

所以OpenFOAM的ddtCorr的算法如下：

首先用过去时间的通量$\phi^{n-1}$和速度$U^{n-1}$估算通量的修正$\phi^{c}$。
$$
\phi^{c} = \phi^{n-1} - U_f(U^{n-1})\cdot \mathbf S_f
$$
计算ddt耦合系数$K_c$：
$$
K_c = 1 - \min\left( \frac {|\phi^c|}{|\phi^{n-1}|+\epsilon}, 1\right)
$$
同时，应该注意到，ddt耦合系数$K_c$是一个大于0的无量纲界面张量场！它的边界条件也被特殊地处理了，对于固定速度边界或cyclicAMI边界，它的值是0.
输出修正量$\phi^k = K_c\phi^c/\Delta t$

但是为什么这么计算耦合系数$K_c$和通量$\phi^c$，我也不知道！

cfd-online上一个叫eugene的id提到ddtCorr实质上是文献[^19]中的所谓Choi Correction，对应的原始文献是[^20] ，但发帖者认为fvcDdtPhiCoeff = 1时是完整的Choi Correction，但此时OpenFOAM会不稳定，所以用了一个经验性的fvcDdtPhiCoeff()函数来限制这个Choi Correction。

Jasak的搞法

Jasak似乎玩了另外的花样。在foam-extend 4.0中，Jasak的icoFoam和OpenFOAM.org以及OpenFOAMplus的icoFoam就不一样了。

Jasak的icoFoam把UEqn分成了HUEqn和ddtUEqn两部分：
- HUEqn：对流和扩散项
- ddtUEqn：非定常项；
动量预测步：HUEqn + ddtUEqn == -grad(p)
PISO循环中：H从HUEqn.H()中产生！与时间步完全无关。
速度更新：

//https://github.com/Unofficial-Extend-Project-Mirror/foam-extend-foam-extend-4.0/blob/master/applications/solvers/incompressible/icoFoam/icoFoam.C +113
// Note: cannot call H(U) here because the velocity is not complete
// HJ, 22/Jan/2016
U = 1.0/(aU + ddtUEqn.A())*
(
    U*aU - fvc::grad(p) + ddtUEqn.H()
);

Jasak 认为这使得代码更加清晰了（时间导数算子中完全去掉了ddtPhiCorr和ddtCorr）。而且完全消除了时间步和松弛因子对结果的影响。结果对比参考文献[^17]

投影法

另外一种选择是用投影法，参考文献[^21]。投影法在同位网格上可以避免Rhie-Chow插值，从而避免这个ddtCorr()的问题。他们也做了数值试验来观察这个质量通量修正项（mass flux correction term）对PISO算法的效果。发现扔掉这一项，PISO算法的耗散大大降低。

OpenFOAMplus的搞法

OpenFOAMplus在v1706中的ddtScheme加入了一个scalar ddtPhiCoeff_的定义。从而使ddtCorr()的实现更加复杂化。

代码如下：

//https://develop.openfoam.com/Development/OpenFOAM-plus/blob/master/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
template<class Type>
class ddtScheme
:
    public tmp<ddtScheme<Type>>::refCount
{

protected:

    // Protected data

        const fvMesh& mesh_;

        //- Input for fvcDdtPhiCoeff (-1 default)
        scalar ddtPhiCoeff_; //ebf654f2这个集成rhoPimpleAdiabaticFoam的提交中才出现的。

//...

而此时更改为了：

// https://develop.openfoam.com/Development/OpenFOAM-plus/blob/master/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C +151
template<class Type>
tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const fluxFieldType& phi,
    const fluxFieldType& phiCorr
)
{
    tmp<surfaceScalarField> tddtCouplingCoeff
    (
        new surfaceScalarField
        (
            IOobject
            (
                "ddtCouplingCoeff",
                U.mesh().time().timeName(),
                U.mesh()
            ),
            U.mesh(),
            dimensionedScalar("one", dimless, 1.0)
        )
    );

    surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();

    if (ddtPhiCoeff_ < 0)
    {
        ddtCouplingCoeff -= min
        (
            mag(phiCorr)
           /(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
            scalar(1)
        );
    }
    else
    {
        ddtCouplingCoeff =
            dimensionedScalar("ddtPhiCoeff", dimless, ddtPhiCoeff_);
    }

    surfaceScalarField::Boundary& ccbf = ddtCouplingCoeff.boundaryFieldRef();

    forAll(U.boundaryField(), patchi)
    {
        if
        (
            U.boundaryField()[patchi].fixesValue()
         || isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
        )
        {
            ccbf[patchi] = 0.0;
        }
    }

    if (debug > 1)
    {
        InfoInFunction
            << "ddtCouplingCoeff mean max min = "
            << gAverage(ddtCouplingCoeff.primitiveField())
            << " " << gMax(ddtCouplingCoeff.primitiveField())
            << " " << gMin(ddtCouplingCoeff.primitiveField())
            << endl;
    }

    return tddtCouplingCoeff;
}

这样使得可以这个系数可以手工调整。

他们所引用的参考文献[^22] 暂时没找到全文

传送门：OpenFOAM的不可压缩流算法

李东岳

非常感谢 @程迪的讨论，我重新回顾一下看看能不能分享一些看法。

有关ddtPhi()这个函数的形式大家已经明白了，更重要的是明白为什么要这么做。MULES限制器和ddtPhi都是Henry Weller发明的，MULES限制器的植入以及原因目前都清楚了。但是ddtPhi的原因尚不明了。

@dyj19901127
宇老师大作被翻出来了，哈哈

CFD中文网

pEqn.H中phiHbyA计算修正项的问题

关于fvc::ddtCorr()和fvc::ddtPhiCorr()

Jasak的搞法

投影法

OpenFOAMplus的搞法