Added stabilisation to the nut update.

This commit is contained in:
henry 2009-05-13 00:02:36 +01:00
parent 0ee817a060
commit 0d2fd6b252

View File

@ -49,21 +49,20 @@ scalar nutRoughWallFunctionFvPatchScalarField::fnRough
const scalar kappa
) const
{
// Set deltaB based on non-dimensional roughness height
scalar deltaB = 0.0;
// Return fn based on non-dimensional roughness height
if (KsPlus < 90.0)
{
deltaB =
1.0/kappa
*log((KsPlus - 2.25)/87.75 + Cs*KsPlus)
*sin(0.4258*(log(KsPlus) - 0.811));
return pow
(
(KsPlus - 2.25)/87.75 + Cs*KsPlus,
sin(0.4258*(log(KsPlus) - 0.811))
);
}
else
{
deltaB = 1.0/kappa*log(1.0 + Cs*KsPlus);
return (1.0 + Cs*KsPlus);
}
return exp(min(deltaB*kappa, 50.0));
}
@ -183,7 +182,7 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
const scalar Cmu25 = pow(Cmu, 0.25);
const scalar kappa = ras.kappa().value();
const scalar E = ras.E().value();
scalar yPlusLam = ras.yPlusLam();
const scalar yPlusLam = ras.yPlusLam();
const scalarField& y = ras.y()[patch().index()];
@ -199,17 +198,35 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
label faceCellI = patch().faceCells()[faceI];
scalar uStar = Cmu25*sqrt(k[faceCellI]);
scalar yPlus = uStar*y[faceI]/nuw[faceI];
scalar KsPlus = uStar*Ks_[faceI]/nuw[faceI];
scalar Edash = E;
scalar yPlusLamNew = yPlusLam;
if (KsPlus > 2.25)
{
Edash = E/fnRough(KsPlus, Cs_[faceI], kappa);
yPlusLam = ras.yPlusLam(kappa, Edash);
}
if (yPlus > yPlusLam)
{
scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
// To avoid oscillations limit the change in the wall viscosity
// which is particularly important if it temporarily becomes zero
nutw[faceI] =
max
(
min
(
nuw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1),
2*limitingNutw
), 0.5*limitingNutw
);
}
else
{
nutw[faceI] = 0.0;
}
if (debug)
@ -217,18 +234,9 @@ void nutRoughWallFunctionFvPatchScalarField::updateCoeffs()
Info<< "yPlus = " << yPlus
<< ", KsPlus = " << KsPlus
<< ", Edash = " << Edash
<< ", yPlusLam = " << yPlusLam
<< ", nutw = " << nutw[faceI]
<< endl;
}
if (yPlus > yPlusLamNew)
{
nutw[faceI] = nuw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1);
}
else
{
nutw[faceI] = 0.0;
}
}
}