关于icoFoam、pisoFoam等求解过程的连续性误差的疑惑

Shihang Chen

各位老师好，我关于piso算法进行速度修正有一点疑惑，想请教一下。

U = HbyA - rAU*fvc::grad(p);修正后的速度场准确满足连续性条件吗？
在速度压力耦合计算结束后，如何构建满足连续性条件的phi？

谢谢！

以下是具体问题：
这是icoFoam中最后一部分构建压力泊松方程并且迭代求解的代码：

            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector

                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();

下面是 #include "continuityErrs.H"


{
    volScalarField contErr(fvc::div(phi));

    scalar sumLocalContErr = runTime.deltaTValue()*
        mag(contErr)().weightedAverage(mesh.V()).value();

    scalar globalContErr = runTime.deltaTValue()*
        contErr.weightedAverage(mesh.V()).value();
    cumulativeContErr += globalContErr;

    Info<< "time step continuity errors : sum local = " << sumLocalContErr
        << ", global = " << globalContErr
        << ", cumulative = " << cumulativeContErr
        << endl;
}

使用cavity进行测试发现此时的
time step continuity errors : sum local = 9.66354e-09，在1e-8这个量级，满足连续性条件。

理论上最后根据（U = HbyA - rAUfvc::grad(p);）得到的流场应该满足连续性条件才对，那为什么是在U修正之前进行连续性的计算呢？如果在（U = HbyA - rAUfvc::grad(p);）之后再次构建phi，进行连续性判断会发现局部连续性误差显著增大。

...
U = HbyA - rAU*fvc::grad(p);

surfaceScalarField phi2 = fvc::flux(U);
            
volScalarField contErr(fvc::div(phi2));

scalar sumLocalContErr = runTime.deltaTValue()*
                mag(contErr)().weightedAverage(mesh.V()).value();

scalar globalContErr = runTime.deltaTValue()*
                contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;

Info<< "time step continuity errors : sum local = " << sumLocalContErr
                << ", global = " << globalContErr
                << ", cumulative = " << cumulativeContErr
                << endl;

使用cavity进行测试发现此时的
time step continuity errors : sum local = 00222664，在2e-3左右，误差较大。
因此，是否是这样的phi的构造方式有问题呢？如何得到速度修正之后，最终满足连续性条件的phi呢。

李东岳

U = HbyA - rAU*fvc::grad(p);修正后的速度场准确满足连续性条件吗？

不满足。满足连续性条件的是phi，而不是U。

在速度压力耦合计算结束后，如何构建满足连续性条件的phi？

phi = phiHbyA - pEqn.flux();这一行。这个讨论在论坛里面有很多讨论，你可以搜一下。我记得 @coolhhh 发了很多帖子讨论这个问题。

coolhhh

可以参考这个帖子 pisoFoam计算的U，求div(U)结果为何不是严格等于0

Shihang Chen

@李东岳
 @coolhhh
感谢两位老师，我提出这个问题的主要原因是在现有的CFD-DEM流固耦合计算中，一种基于FD/IBM的细网格方法是需要对速度场U进行改变，并针对连续性方程进行修正的。而这个连续性修正的过程就需要再次显式构建连续性条件：div(U)=0。当然这是理论上的条件。
$\nabla \cdot (\nabla \phi) = \nabla \cdot U$
而下面代码中的连续性修正因子的phiIB理论上只在颗粒存在的条件下才有意义，同样才有值。但是如果div(U)在颗粒不存在得情况就不为0了，那么这样的修正就是错误的。虽然当前这个方面的主要算法就是下面展示的这样。

目前看来在PISO外部，U和phi是分离的，但是CFD-DEM的修正是针对U的，这样操作应该是有问题的，误差大小有待考察。而我原本的目标是希望能够在PISO外，得到满足连续性的U或者phi，不过目前看来是很难实现的。

上述过程具体代码如下：

void Foam::cfdemCloud::calcVelocityCorrection
(
    volScalarField& p,
    volVectorField& U,
    volScalarField& phiIB,
    volScalarField& voidfraction//颗粒体积分数（颗粒内部为1
)
{

    
    void Foam::cfdemCloudIB::calcVelocityCorrection
    
    (
    
        volScalarField& p,
        volVectorField& U,
        volScalarField& phiIB,
        volScalarField& voidfraction
    
    )
    {
        setParticleVelocity(U);//改变颗粒所在区域流体的速度
        // 修改速度会使其不满足连续性，因此使用一个phiIB进行修正
    
        fvScalarMatrix phiIBEqn 
        (  
            fvm::laplacian(phiIB) == fvc::div(U) + fvc::ddt(voidfraction)
        );
    
        if(phiIB.needReference()) 
    
        {
    
             phiIBEqn.setReference(pRefCell_, pRefValue_);
    
        }
    
        phiIBEqn.solve();

        U=U-fvc::grad(phiIB);
    
        U.correctBoundaryConditions();

        // correct the pressure as well
    
        p=p+phiIB/U.mesh().time().deltaT(); 
    
        p.correctBoundaryConditions();  
    }

CFD中文网

关于icoFoam、pisoFoam等求解过程的连续性误差的疑惑