surfaceInterpolation::makeWeight: improve handling of very small faces

If division by the face area might cause an overflow use the cell-face distances
directly to calculate the weights.
This commit is contained in:
Henry Weller
2019-03-19 20:27:55 +00:00
parent aaf4a0bfcc
commit 5a978bbc65
3 changed files with 71 additions and 19 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -43,15 +43,30 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const
{
const cyclicFvPatch& nbrPatch = neighbFvPatch();
const scalarField deltas(nf()&coupledFvPatch::delta());
const scalarField nbrDeltas(nbrPatch.nf()&nbrPatch.coupledFvPatch::delta());
const vectorField delta(coupledFvPatch::delta());
const vectorField nbrDelta(nbrPatch.coupledFvPatch::delta());
forAll(deltas, facei)
const scalarField nfDelta(nf() & delta);
const scalarField nbrNfDelta(nbrPatch.nf() & nbrDelta);
forAll(delta, facei)
{
scalar di = deltas[facei];
scalar dni = nbrDeltas[facei];
const scalar ndoi = nfDelta[facei];
const scalar ndni = nbrNfDelta[facei];
const scalar ndi = ndoi + ndni;
w[facei] = dni/(di + dni);
if (ndni/vGreat < ndi)
{
w[facei] = ndni/ndi;
}
else
{
const scalar doi = mag(delta[facei]);
const scalar dni = mag(nbrDelta[facei]);
const scalar di = doi + dni;
w[facei] = dni/di;
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -42,20 +42,44 @@ void Foam::processorFvPatch::makeWeights(scalarField& w) const
{
if (Pstream::parRun())
{
const vectorField delta(coupledFvPatch::delta());
// The face normals point in the opposite direction on the other side
scalarField neighbFaceCentresCn
const vectorField nbrDelta
(
procPolyPatch_.neighbFaceCentres()
- procPolyPatch_.neighbFaceCellCentres()
);
const scalarField nfDelta(nf() & delta);
const scalarField nbrNfDelta
(
(
procPolyPatch_.neighbFaceAreas()
/(mag(procPolyPatch_.neighbFaceAreas()) + vSmall)
)
& (
procPolyPatch_.neighbFaceCentres()
- procPolyPatch_.neighbFaceCellCentres())
) & nbrDelta
);
w = neighbFaceCentresCn
/((nf()&coupledFvPatch::delta()) + neighbFaceCentresCn);
forAll(delta, facei)
{
const scalar ndoi = nfDelta[facei];
const scalar ndni = nbrNfDelta[facei];
const scalar ndi = ndoi + ndni;
if (ndni/vGreat < ndi)
{
w[facei] = ndni/ndi;
}
else
{
const scalar doi = mag(delta[facei]);
const scalar dni = mag(nbrDelta[facei]);
const scalar di = doi + dni;
w[facei] = dni/di;
}
}
}
else
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -178,9 +178,22 @@ void Foam::surfaceInterpolation::makeWeights() const
// 90 deg and the dot-product will be positive. For invalid
// meshes (d & s <= 0), this will stabilise the calculation
// but the result will be poor.
scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
w[facei] = SfdNei/(SfdOwn + SfdNei);
const scalar SfdOwn = mag(Sf[facei]&(Cf[facei] - C[owner[facei]]));
const scalar SfdNei = mag(Sf[facei]&(C[neighbour[facei]] - Cf[facei]));
const scalar SfdOwnNei = SfdOwn + SfdNei;
if (SfdNei/vGreat < SfdOwnNei)
{
w[facei] = SfdNei/SfdOwnNei;
}
else
{
const scalar dOwn = mag(Cf[facei] - C[owner[facei]]);
const scalar dNei = mag(C[neighbour[facei]] - Cf[facei]);
const scalar dOwnNei = dOwn + dNei;
w[facei] = dNei/dOwnNei;
}
}
surfaceScalarField::Boundary& wBf =