# 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);
}``````

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

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

`````` 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来起到一个双保险的作用。

``````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(),
);

// 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);
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();
}
}
}``````

`````` void Foam::fvc::makeRelative
(
surfaceScalarField& phi,
const volVectorField& U
)
{
if (phi.mesh().moving())
{
phi -= fvc::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()
);
}``````

### 参考

OpenFOAM: src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C Source File