interFoam-UEqn.H/pEqn.H

压力方程 “pEqn.H”

{

    volScalarField rAU("rAU", 1.0/UEqn.A());      // rAU:在速度方程的的最后一个解中，矩阵对角项系数的倒数
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));       //转换为表面标量场
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));      //定义HbyA
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );         //将上述体积矢量场转化为面心标量场（参考有限体积法）,并保证速度通量的全局守恒，以确保压力方程有解
    MRF.makeRelative(phiHbyA);    //将绝对通量转化为相对于面的通量(参考动网格技术)
    adjustPhi(phiHbyA, U, p_rgh);        //对方程进行修正，保证速度边界条件守恒
    surfaceScalarField phig     
    (
        (
            mixture.surfaceTensionForce()
          - ghf*fvc::snGrad(rho)
        )*rAUf*mesh.magSf()
    );
    phiHbyA += phig;              //更新面通量场
    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);       // Update the pressure BCs to ensure flux consistency
    while (pimple.correctNonOrthogonal())          // 非正交压力修正循环 ，根据fvSolution字典文件中设定值n,求解以下方程n-1次 
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );
        //只有求解区域所有的压力边界都为第二类边界条件时，上面的值才会用到。
        //如果有第一类边界条件，压力参考值为该点处边界值。
        //对于不可压缩流动压力值为相对值，上面的参考值的大小对结果无影响，除非存在边界压力。
        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));		//通过查询system/fvSolution下的PIMPLE更新压力参考网格并重新设定参考值
        if (pimple.finalNonOrthogonalIter())    //若外迭代次数为n,则压力场的松弛仅在n-1次外循环前进行；若外迭代次数为1，这里使用piso算法，压力不进行亚松驰
        {
            phi = phiHbyA - p_rghEqn.flux();           //在最后一次非正交内循环中，使用最新压力校正通量
            p_rgh.relax();
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);       //校正速度，满足边界条件(主要针对第二类边界条件）
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }
    #include "continuityErrs.H"         //计算连续性方程误差
    p == p_rgh + rho*gh;       // 计算总压=动压+静压
    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );   
        p_rgh = p - rho*gh;       //更新动压
    }
}

p.relax(); || UEqn().relax();
场松弛会改变压力的值，方程松弛会改变压力的变化量。

UEqn.H

MRF.correctBoundaryVelocity(U);   //MRF:基于旋转坐标系下的N-S方程修正边界速度,不需要可忽略
fvVectorMatrix UEqn
(
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)       //该项包含两部分：1.分子扩散项 2.雷诺应力的偏分量(dev)的散度。雷诺应力主分量（hyd）的散度归结到了压力项中，这是大多数雷诺时均和大涡模型实现的一贯做法。因此压力比层流模型中多了一项,那就是雷诺应力的主分量，通常被称为湍动压力（动压：p_rgh）                   
     ==
        fvOptions(rho, U)     //源项或约束，不需要可忽略
);
UEqn.relax();       //松弛迭代,加速收敛
fvOptions.constrain(UEqn);  //更正边界条件，维持通量守恒
//momentumPredictor：动量预测求解开关，对多相流以及低雷诺数一般设置为off；如果进行速度预测（on），则求解完整的动量方程得到预测速度；如果不进行速度预测(off），则预测速度直接取当前已知时间步的速度
if (pimple.momentumPredictor())
{
        solve
        (
            UEqn
         ==
            fvc::reconstruct   //显示重构网格中心面（光滑了整个分布情况），转换为体中心通量
            (
                (
                    mixture.surfaceTensionForce()     //表面张力项
                  - ghf*fvc::snGrad(rho)              //体积力项（垂直梯度）
                  - fvc::snGrad(p_rgh)                // `p_rgh:动压`;压力梯度项
                ) * mesh.magSf()                      //网格面矢量绝对值      
            )
        )；
        fvOptions.correct(U);   保证动量守恒
}