Bugfix: MRF zone flux check

This commit is contained in:
Hrvoje Jasak 2017-01-24 16:44:22 +00:00
parent 58edb15fa9
commit 78a9142214
5 changed files with 86 additions and 60 deletions

View file

@ -139,13 +139,15 @@ void Foam::MRFZone::setMRFFaces()
labelList nIncludedFaces(patches.size(), 0);
labelList nExcludedFaces(patches.size(), 0);
register label faceI;
forAll (patches, patchi)
{
const polyPatch& pp = patches[patchi];
forAll (pp, patchFaceI)
{
label faceI = pp.start() + patchFaceI;
faceI = pp.start() + patchFaceI;
if (faceType[faceI] == 1)
{
@ -174,7 +176,7 @@ void Foam::MRFZone::setMRFFaces()
forAll (pp, patchFaceI)
{
label faceI = pp.start() + patchFaceI;
faceI = pp.start() + patchFaceI;
if (faceType[faceI] == 1)
{
@ -187,7 +189,6 @@ void Foam::MRFZone::setMRFFaces()
}
}
if (debug)
{
faceSet internalFaces
@ -489,42 +490,51 @@ void Foam::MRFZone::absoluteFlux
void Foam::MRFZone::meshPhi
(
surfaceVectorField& phi
surfaceScalarField& phi
) const
{
const surfaceVectorField& Cf = mesh_.Cf();
const surfaceVectorField& Sf = mesh_.Sf();
const vector& origin = origin_.value();
const vector rotVel = Omega();
// Internal faces
const vectorField& CfIn = Cf.internalField();
const vectorField& SfIn = Sf.internalField();
register label faceI, patchFaceI;
forAll (internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] = (rotVel ^ (Cf[facei] - origin));
faceI = internalFaces_[i];
phi[faceI] = SfIn[faceI] & (rotVel ^ (CfIn[faceI] - origin));
}
// Included patches
forAll (includedFaces_, patchi)
{
forAll (includedFaces_[patchi], i)
{
label patchFaceI = includedFaces_[patchi][i];
phi.boundaryField()[patchi][patchFaceI] =
(rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin));
forAll (includedFaces_, patchI)
{
forAll (includedFaces_[patchI], i)
{
patchFaceI = includedFaces_[patchI][i];
phi.boundaryField()[patchI][patchFaceI] =
Sf.boundaryField()[patchI][patchFaceI]
& (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin));
}
}
// Excluded patches
forAll (excludedFaces_, patchi)
forAll (excludedFaces_, patchI)
{
forAll (excludedFaces_[patchi], i)
forAll (excludedFaces_[patchI], i)
{
label patchFaceI = excludedFaces_[patchi][i];
patchFaceI = excludedFaces_[patchI][i];
phi.boundaryField()[patchi][patchFaceI] =
(rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin));
phi.boundaryField()[patchI][patchFaceI] =
Sf.boundaryField()[patchI][patchFaceI]
& (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin));
}
}
}
@ -534,21 +544,23 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
const vector& origin = origin_.value();
const vector rotVel = Omega();
register label patchFaceI;
// Included patches
forAll (includedFaces_, patchi)
forAll (includedFaces_, patchI)
{
const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
const vectorField& patchC = mesh_.Cf().boundaryField()[patchI];
vectorField pfld(U.boundaryField()[patchi]);
vectorField pfld(U.boundaryField()[patchI]);
forAll (includedFaces_[patchi], i)
forAll (includedFaces_[patchI], i)
{
label patchFaceI = includedFaces_[patchi][i];
patchFaceI = includedFaces_[patchI][i];
pfld[patchFaceI] = (rotVel ^ (patchC[patchFaceI] - origin));
}
U.boundaryField()[patchi] == pfld;
U.boundaryField()[patchI] == pfld;
}
}

View file

@ -221,7 +221,7 @@ public:
) const;
//- Compute the pseudo mesh phi of the MRF region
void meshPhi(surfaceVectorField& phi) const;
void meshPhi(surfaceScalarField& phi) const;
//- Correct the boundary velocity for the roation of the MRF region
void correctBoundaryVelocity(volVectorField& U) const;

View file

@ -44,34 +44,42 @@ void Foam::MRFZone::relativeRhoFlux
const vector rotVel = Omega();
// Internal faces
const vectorField& CfIn = Cf.internalField();
const vectorField& SfIn = Sf.internalField();
register label faceI, patchFaceI;
forAll (internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] -= rho[facei]*(rotVel ^ (Cf[facei] - origin)) & Sf[facei];
faceI = internalFaces_[i];
phi[faceI] -=
rho[faceI]*(rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI];
}
// Included patches
forAll (includedFaces_, patchi)
// Included patches: reset the flux to exactly zero to avoid
// round-off issues
forAll (includedFaces_, patchI)
{
forAll (includedFaces_[patchi], i)
forAll (includedFaces_[patchI], i)
{
label patchFacei = includedFaces_[patchi][i];
patchFaceI = includedFaces_[patchI][i];
phi.boundaryField()[patchi][patchFacei] = 0.0;
phi.boundaryField()[patchI][patchFaceI] = 0.0;
}
}
// Excluded patches
forAll (excludedFaces_, patchi)
forAll (excludedFaces_, patchI)
{
forAll (excludedFaces_[patchi], i)
forAll (excludedFaces_[patchI], i)
{
label patchFacei = excludedFaces_[patchi][i];
patchFaceI = excludedFaces_[patchI][i];
phi.boundaryField()[patchi][patchFacei] -=
rho.boundaryField()[patchi][patchFacei]
*(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
& Sf.boundaryField()[patchi][patchFacei];
phi.boundaryField()[patchI][patchFaceI] -=
rho.boundaryField()[patchI][patchFaceI]
*(rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
& Sf.boundaryField()[patchI][patchFaceI];
}
}
}
@ -91,35 +99,41 @@ void Foam::MRFZone::absoluteRhoFlux
const vector rotVel = Omega();
// Internal faces
const vectorField& CfIn = Cf.internalField();
const vectorField& SfIn = Sf.internalField();
register label faceI, patchFaceI;
forAll (internalFaces_, i)
{
label facei = internalFaces_[i];
phi[facei] += (rotVel ^ (Cf[facei] - origin)) & Sf[facei];
faceI = internalFaces_[i];
phi[faceI] += (rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI];
}
// Included patches
forAll (includedFaces_, patchi)
forAll (includedFaces_, patchI)
{
forAll (includedFaces_[patchi], i)
forAll (includedFaces_[patchI], i)
{
label patchFacei = includedFaces_[patchi][i];
patchFaceI = includedFaces_[patchI][i];
phi.boundaryField()[patchi][patchFacei] +=
(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
& Sf.boundaryField()[patchi][patchFacei];
phi.boundaryField()[patchI][patchFaceI] +=
(rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
& Sf.boundaryField()[patchI][patchFaceI];
}
}
// Excluded patches
forAll (excludedFaces_, patchi)
forAll (excludedFaces_, patchI)
{
forAll (excludedFaces_[patchi], i)
forAll (excludedFaces_[patchI], i)
{
label patchFacei = excludedFaces_[patchi][i];
patchFaceI = excludedFaces_[patchI][i];
phi.boundaryField()[patchi][patchFacei] +=
(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
& Sf.boundaryField()[patchi][patchFacei];
phi.boundaryField()[patchI][patchFaceI] +=
(rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin))
& Sf.boundaryField()[patchI][patchFaceI];
}
}
}

View file

@ -86,11 +86,11 @@ Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::fluxCorrection() const
}
Foam::tmp<Foam::surfaceVectorField> Foam::MRFZones::meshPhi() const
Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::meshPhi() const
{
tmp<surfaceVectorField> tMRFZonesPhiCorr
tmp<surfaceScalarField> tMRFZonesPhiCorr
(
new surfaceVectorField
new surfaceScalarField
(
IOobject
(
@ -101,10 +101,10 @@ Foam::tmp<Foam::surfaceVectorField> Foam::MRFZones::meshPhi() const
IOobject::NO_WRITE
),
mesh_,
dimensionedVector("zero", dimVelocity, vector::zero)
dimensionedScalar("zero", dimVolume/dimTime, 0)
)
);
surfaceVectorField& MRFZonesPhiCorr = tMRFZonesPhiCorr();
surfaceScalarField& MRFZonesPhiCorr = tMRFZonesPhiCorr();
forAll (*this, i)
{

View file

@ -81,7 +81,7 @@ public:
tmp<surfaceScalarField> fluxCorrection() const;
//- Return pseudo mesh flux
tmp<surfaceVectorField> meshPhi() const;
tmp<surfaceScalarField> meshPhi() const;
// Incompressible MRF
@ -134,7 +134,7 @@ public:
void absoluteVelocity(volVectorField& Urel) const;
//- Correct the boundary velocity for the roation of the MRF region
//- Correct the boundary velocity for the rotation of the MRF region
void correctBoundaryVelocity(volVectorField& U) const;
//- Compute source term for volScalarFields