OpenFOAM动网格的通量修正-卡核

OpenFOAM处理动网格的思路很简单，就是在网格变形（mesh.controledUpdate()）之后，对速度通量进行修正。其中包括：correctPhi.H，fvc::makeRelative(phi,U)。下面一个个进行学习并记录。

correctPhi.H

CorrectPhi
(
    U,
    phi,
    p,
    dimensionedScalar("rAUf", dimTime, 1),
    geometricZeroField(),
    pimple
);

#include "continuityErrs.H"

correctPhi在做动网格的时候都需要在fvSolution的PIMPLE下设置correctPhi true。静网格则不需要。它用来修正迭代求解之前因为网格变动而需要同时改变的面通量phi，以保持连续性。

if (correctPhi)
{
    // Calculate absolute flux
    // from the mapped surface velocity
    phi = mesh.Sf() & Uf();

    #include "correctPhi.H"

    // Make the flux relative to the mesh motion
    fvc::makeRelative(phi, U);
}

可以看到，在correctPhi之前，phi重新计算了一遍。这里Uf是上一步pEqn.H的最后进行修正的：

// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);

因为网格变动，所以面速度Uf需要根据网格运动来修正:

 void Foam::fvc::correctUf
 (
     autoPtr<surfaceVectorField>& Uf,
     const volVectorField& U,
     const surfaceScalarField& phi
 )
 {
     const fvMesh& mesh = U.mesh();
 
     if (mesh.dynamic())
     {
         Uf() = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
         Uf() += n*(phi/mesh.magSf() - (n & Uf()));
     }
 }

phi是网格变形之前的速度通量，而Uf是上一步求解之后的面上速度。其实有一点是不太确定的，那就是在这一步求解完之后，phi和Uf不应该是一致的吗？那么那一行Uf() += n*(phi/mesh.magSf() – (n & Uf()));又是做什么呢？我们知道动网格下网格变形后，网格面上的通量phi因为面位置变化而不可靠了，所以需要Uf来起到一个双保险的作用。

尔后在下一个时间步开始前，网格变形后，先用 phi = mesh.Sf() & Uf();来重新计算面上通量值，但因为这个时候Uf也是网格变形之前的值，所以，需要调用correctPhi进行修正。在这里，Uf是作为储存面上速度的作用，毕竟要变成phi还是得与面法向矢量进行内积。这个时候，mesh.Sf()就已经是网格变形之后的面法向矢量了，所以相当于修正了一半？另一半需要在correctPhi里对速度通量进行剩下的修正。correctPhi在src/finiteVolume/cfdtools/general的文件夹中，看到它调用的类构造函数：

template<class RAUfType, class DivUType>
void Foam::CorrectPhi
(
    volVectorField& U,
    surfaceScalarField& phi,
    const volScalarField& p,
    const RAUfType& rAUf,
    const DivUType& divU,
    pimpleControl& pimple
)
{
    const fvMesh& mesh = U.mesh();
    const Time& runTime = mesh.time();

    correctUphiBCs(U, phi);

    // Initialize BCs list for pcorr to zero-gradient
    wordList pcorrTypes
    (
        p.boundaryField().size(),
        zeroGradientFvPatchScalarField::typeName
    );

    // Set BCs of pcorr to fixed-value for patches at which p is fixed
    forAll(p.boundaryField(), patchi)
    {
        if (p.boundaryField()[patchi].fixesValue())
        {
            pcorrTypes[patchi] = fixedValueFvPatchScalarField::typeName;
        }
    }

    volScalarField pcorr
    (
        IOobject
        (
            "pcorr",
            runTime.timeName(),
            mesh
        ),
        mesh,
        dimensionedScalar(p.dimensions(), Zero),
        pcorrTypes
    );

    if (pcorr.needReference())
    {
        fvc::makeRelative(phi, U);
        adjustPhi(phi, U, pcorr);
        fvc::makeAbsolute(phi, U);
    }

    mesh.setFluxRequired(pcorr.name());

    while (pimple.correctNonOrthogonal())
    {
        // Solve for pcorr such that the divergence of the corrected flux
        // matches the divU provided (from previous iteration, time-step...)
        fvScalarMatrix pcorrEqn
        (
            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
        );

        pcorrEqn.setReference(0, 0);

        pcorrEqn.solve
        (
            mesh.solver(pcorr.select(pimple.finalNonOrthogonalIter()))
        );

        if (pimple.finalNonOrthogonalIter())
        {
            phi -= pcorrEqn.flux();
        }
    }
}

correctPhi理论上是为了满足压力方程的相容性条件，实际上调用了adjustPhi。在icoFoam等求解器adjustPhi是对p进行的，而在这里则是对pcorr进行的。通过求解pcorrEqn，强制满足压力方程的相容性条件，得到面通量修正值pcorrEqn.flux()，再在最后对phi进行修正。

最后是makeRelative函数：

 void Foam::fvc::makeRelative
 (
     surfaceScalarField& phi,
     const volVectorField& U
 )
 {
     if (phi.mesh().moving())
     {
         phi -= fvc::meshPhi(U);
     }
 }

其解释是“make the given flux relative to the mesh motion”，字面意思是使给定的通量与网格运动相关。查找meshPhi(U)，返回一个ddtScheme有关的通量。我想，这与时间格式相关。如二阶时间格式backward下的meshPhi(U)：

template<class Type>
 tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const scalar deltaT = deltaT_();
     const scalar deltaT0 = deltaT0_(vf);
 
     // Coefficient for t-3/2 (between times 0 and 00)
     const scalar coefft0_00 = deltaT/(deltaT + deltaT0);
 
     // Coefficient for t-1/2 (between times n and 0)
     const scalar coefftn_0 = 1 + coefft0_00;
 
     return surfaceScalarField::New
     (
         mesh().phi().name(),
         coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
     );
 }

可以看到，其返回的通量场，不仅包括当前时间步的通量，而且也包含上一时间步的通量，由特定的系数进行权重的匹配。