mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
settlingFoam: improved handling of the turbulence and wall-functions
This commit is contained in:
Henry
@ -41,12 +41,6 @@
|
||||
{
|
||||
label faceCelli = currPatch.faceCells()[facei];
|
||||
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]
|
||||
*::sqrt(k[faceCelli])
|
||||
/nub_;
|
||||
|
||||
|
||||
// For corner cells (with two boundary or more faces),
|
||||
// epsilon and G in the near-wall cell are calculated
|
||||
// as an average
|
||||
@ -57,8 +51,6 @@
|
||||
Cmu75*::pow(k[faceCelli], 1.5)
|
||||
/(kappa_*y[patchi][facei]);
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
G[faceCelli] +=
|
||||
(nutbw[facei] + nub_)*magFaceGradU[facei]
|
||||
*Cmu25*::sqrt(k[faceCelli])
|
||||
@ -66,7 +58,6 @@
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// perform the averaging
|
||||
@ -83,6 +74,7 @@
|
||||
|
||||
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
cellBoundaryFaceCount[faceCelli] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -20,16 +20,11 @@
|
||||
|
||||
// calculate yPlus
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]
|
||||
*::sqrt(k[faceCelli])
|
||||
/nub_;
|
||||
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nub_;
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
nutw[facei] =
|
||||
yPlus*nub_*kappa_
|
||||
/::log(E_*yPlus)
|
||||
- nub_;
|
||||
nutw[facei] = nub_*(yPlus*kappa_/::log(E_*yPlus) -1);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -41,12 +41,6 @@
|
||||
{
|
||||
label faceCelli = currPatch.faceCells()[facei];
|
||||
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]
|
||||
*::sqrt(k[faceCelli])
|
||||
/nu2_;
|
||||
|
||||
|
||||
// For corner cells (with two boundary or more faces),
|
||||
// epsilon and G in the near-wall cell are calculated
|
||||
// as an average
|
||||
@ -57,8 +51,6 @@
|
||||
Cmu75*::pow(k[faceCelli], 1.5)
|
||||
/(kappa_*y[patchi][facei]);
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
G[faceCelli] +=
|
||||
(nut2w[facei] + nu2_)*magFaceGradU[facei]
|
||||
*Cmu25*::sqrt(k[faceCelli])
|
||||
@ -66,7 +58,6 @@
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// perform the averaging
|
||||
@ -83,6 +74,7 @@
|
||||
|
||||
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
cellBoundaryFaceCount[faceCelli] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -20,16 +20,11 @@
|
||||
|
||||
// calculate yPlus
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]
|
||||
*::sqrt(k[faceCelli])
|
||||
/nu2_;
|
||||
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
nutw[facei] =
|
||||
yPlus*nu2_*kappa_
|
||||
/::log(E_*yPlus)
|
||||
- nu2_;
|
||||
nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -9,7 +9,9 @@
|
||||
(Alpha/(scalar(1.001) - Alpha))*(sqr(rhoc)/rho)*Vdj*Vdj,
|
||||
"div(phiVdj,Vdj)"
|
||||
)
|
||||
- fvm::laplacian(mu, U, "laplacian(muEff,U)")
|
||||
- fvm::laplacian(muEff, U, "laplacian(muEff,U)")
|
||||
- (fvc::grad(U) & fvc::grad(muEff))
|
||||
//- fvc::div(muEff*dev2(T(fvc::grad(U))))
|
||||
);
|
||||
|
||||
UEqn.relax();
|
||||
|
||||
@ -10,13 +10,16 @@
|
||||
phi + rhoc*(mesh.Sf() & fvc::interpolate(Vdj))
|
||||
);
|
||||
|
||||
solve
|
||||
fvScalarMatrix AlphaEqn
|
||||
(
|
||||
fvm::ddt(rho, Alpha)
|
||||
+ fvm::div(phiAlpha, Alpha)
|
||||
- fvm::laplacian(mut, Alpha)
|
||||
);
|
||||
|
||||
AlphaEqn.relax();
|
||||
AlphaEqn.solve();
|
||||
|
||||
Info<< "Solid phase fraction = "
|
||||
<< Alpha.weightedAverage(mesh.V()).value()
|
||||
<< " Min(Alpha) = " << min(Alpha).value()
|
||||
|
||||
@ -34,4 +34,6 @@
|
||||
)
|
||||
+ mul;
|
||||
}
|
||||
|
||||
mul = min(mul, muMax);
|
||||
}
|
||||
|
||||
@ -61,6 +61,7 @@
|
||||
dimensionedScalar rhod(transportProperties.lookup("rhod"));
|
||||
|
||||
dimensionedScalar muc(transportProperties.lookup("muc"));
|
||||
dimensionedScalar muMax(transportProperties.lookup("muMax"));
|
||||
|
||||
dimensionedScalar plasticViscosityCoeff
|
||||
(
|
||||
@ -328,12 +329,12 @@
|
||||
);
|
||||
|
||||
|
||||
Info<< "Calculating field mu\n" << endl;
|
||||
volScalarField mu
|
||||
Info<< "Calculating field muEff\n" << endl;
|
||||
volScalarField muEff
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"mu",
|
||||
"muEff",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
@ -345,7 +346,7 @@
|
||||
|
||||
Info<< "Calculating field (g.h)f\n" << endl;
|
||||
volScalarField gh("gh", g & mesh.C());
|
||||
surfaceScalarField ghf(surfaceScalarField("gh", g & mesh.Cf()));
|
||||
surfaceScalarField ghf("gh", g & mesh.Cf());
|
||||
|
||||
volScalarField p
|
||||
(
|
||||
|
||||
@ -5,7 +5,9 @@ if (turbulence)
|
||||
y.correct();
|
||||
}
|
||||
|
||||
dimensionedScalar k0("k0", k.dimensions(), 0);
|
||||
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
|
||||
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
|
||||
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
|
||||
|
||||
volScalarField divU(fvc::div(phi/fvc::interpolate(rho)));
|
||||
@ -28,13 +30,13 @@ if (turbulence)
|
||||
+ fvm::div(phi, epsilon)
|
||||
- fvm::laplacian
|
||||
(
|
||||
mut/sigmaEps + mul, epsilon,
|
||||
mut/sigmaEps + muc, epsilon,
|
||||
"laplacian(DepsilonEff,epsilon)"
|
||||
)
|
||||
==
|
||||
C1*G*epsilon/k
|
||||
C1*G*epsilon/(k + kMin)
|
||||
- fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
|
||||
- fvm::Sp(C2*rho*epsilon/k, epsilon)
|
||||
- fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
|
||||
);
|
||||
|
||||
#include "wallDissipation.H"
|
||||
@ -42,7 +44,7 @@ if (turbulence)
|
||||
epsEqn.relax();
|
||||
epsEqn.solve();
|
||||
|
||||
bound(epsilon, epsilonMin);
|
||||
bound(epsilon, epsilon0);
|
||||
|
||||
|
||||
// Turbulent kinetic energy equation
|
||||
@ -52,19 +54,19 @@ if (turbulence)
|
||||
+ fvm::div(phi, k)
|
||||
- fvm::laplacian
|
||||
(
|
||||
mut/sigmak + mul, k,
|
||||
mut/sigmak + muc, k,
|
||||
"laplacian(DkEff,k)"
|
||||
)
|
||||
==
|
||||
G
|
||||
- fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
|
||||
- fvm::Sp(rho*epsilon/k, k)
|
||||
- fvm::Sp(rho*epsilon/(k + kMin), k)
|
||||
);
|
||||
|
||||
kEqn.relax();
|
||||
kEqn.solve();
|
||||
|
||||
bound(k, kMin);
|
||||
bound(k, k0);
|
||||
|
||||
|
||||
//- Re-calculate viscosity
|
||||
@ -73,4 +75,4 @@ if (turbulence)
|
||||
#include "wallViscosity.H"
|
||||
}
|
||||
|
||||
mu = mut + mul;
|
||||
muEff = mut + mul;
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -57,11 +57,15 @@ int main(int argc, char *argv[])
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.loop())
|
||||
while (runTime.run())
|
||||
{
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
#include "readTimeControls.H"
|
||||
#include "compressibleCourantNo.H"
|
||||
#include "setDeltaT.H"
|
||||
|
||||
runTime++;
|
||||
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
#include "rhoEqn.H"
|
||||
|
||||
|
||||
@ -1,9 +1,10 @@
|
||||
{
|
||||
labelList cellBoundaryFaceCount(epsilon.size(), 0);
|
||||
|
||||
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
scalar Cmu75 = ::pow(Cmu.value(), 0.75);
|
||||
scalar kappa_ = kappa.value();
|
||||
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
|
||||
const scalar kappa_ = kappa.value();
|
||||
const scalar muc_ = muc.value();
|
||||
|
||||
const fvPatchList& patches = mesh.boundary();
|
||||
|
||||
@ -32,9 +33,6 @@
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
const scalarField& rhow = rho.boundaryField()[patchi];
|
||||
|
||||
const scalarField muw(mul.boundaryField()[patchi]);
|
||||
const scalarField& mutw = mut.boundaryField()[patchi];
|
||||
|
||||
scalarField magFaceGradU
|
||||
@ -46,10 +44,6 @@
|
||||
{
|
||||
label faceCelli = curPatch.faceCells()[facei];
|
||||
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
|
||||
/(muw[facei]/rhow[facei]);
|
||||
|
||||
// For corner cells (with two boundary or more faces),
|
||||
// epsilon and G in the near-wall cell are calculated
|
||||
// as an average
|
||||
@ -57,19 +51,17 @@
|
||||
cellBoundaryFaceCount[faceCelli]++;
|
||||
|
||||
epsilon[faceCelli] +=
|
||||
Cmu75*rho[faceCelli]*::pow(k[faceCelli], 1.5)
|
||||
Cmu75*::pow(k[faceCelli], 1.5)
|
||||
/(kappa_*y[patchi][facei]);
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
G[faceCelli] +=
|
||||
mutw[facei]*magFaceGradU[facei]
|
||||
(mutw[facei] + muc_)
|
||||
*magFaceGradU[facei]
|
||||
*Cmu25*::sqrt(k[faceCelli])
|
||||
/(kappa_*y[patchi][facei]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// perform the averaging
|
||||
@ -86,6 +78,7 @@
|
||||
|
||||
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
cellBoundaryFaceCount[faceCelli] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,7 +1,9 @@
|
||||
{
|
||||
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
scalar kappa_ = kappa.value();
|
||||
scalar E_ = E.value();
|
||||
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
const scalar kappa_ = kappa.value();
|
||||
const scalar E_ = E.value();
|
||||
const scalar muc_ = muc.value();
|
||||
const scalar nuc_ = muc_/rhoc.value();
|
||||
|
||||
const fvPatchList& patches = mesh.boundary();
|
||||
|
||||
@ -11,9 +13,6 @@
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
const scalarField& rhow = rho.boundaryField()[patchi];
|
||||
|
||||
const scalarField muw(mul.boundaryField()[patchi]);
|
||||
scalarField& mutw = mut.boundaryField()[patchi];
|
||||
|
||||
forAll(curPatch, facei)
|
||||
@ -22,13 +21,12 @@
|
||||
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
|
||||
/(muw[facei]/rhow[facei]);
|
||||
/nuc_;
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
mutw[facei] =
|
||||
muw[facei]
|
||||
*(yPlus*kappa_/::log(E_*yPlus) - 1);
|
||||
muc_*(yPlus*kappa_/::log(E_*yPlus) - 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -23,12 +23,13 @@ boundaryField
|
||||
inlet
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0.0191 0 0);
|
||||
value $internalField;
|
||||
}
|
||||
|
||||
outlet
|
||||
{
|
||||
type zeroGradient;
|
||||
type pressureInletOutletVelocity;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
bottomWall
|
||||
|
||||
@ -23,12 +23,13 @@ boundaryField
|
||||
inlet
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 0.001;
|
||||
value $internalField;
|
||||
}
|
||||
|
||||
outlet
|
||||
{
|
||||
type zeroGradient;
|
||||
type inletOutlet;
|
||||
inletValue $internalField;
|
||||
}
|
||||
|
||||
bottomWall
|
||||
|
||||
@ -23,12 +23,13 @@ boundaryField
|
||||
inlet
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 1.50919e-06;
|
||||
value $internalField;
|
||||
}
|
||||
|
||||
outlet
|
||||
{
|
||||
type zeroGradient;
|
||||
type inletOutlet;
|
||||
inletValue $internalField;
|
||||
}
|
||||
|
||||
bottomWall
|
||||
|
||||
@ -23,12 +23,13 @@ boundaryField
|
||||
inlet
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 0.00015;
|
||||
value $internalField;
|
||||
}
|
||||
|
||||
outlet
|
||||
{
|
||||
type zeroGradient;
|
||||
type inletOutlet;
|
||||
inletValue $internalField;
|
||||
}
|
||||
|
||||
bottomWall
|
||||
|
||||
@ -22,7 +22,8 @@ boundaryField
|
||||
{
|
||||
inlet
|
||||
{
|
||||
type zeroGradient;
|
||||
type buoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
outlet
|
||||
@ -39,7 +40,8 @@ boundaryField
|
||||
|
||||
endWall
|
||||
{
|
||||
type zeroGradient;
|
||||
type buoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
top
|
||||
|
||||
@ -63,7 +63,7 @@ boundary
|
||||
}
|
||||
bottomWall
|
||||
{
|
||||
type patch;
|
||||
type wall;
|
||||
faces
|
||||
(
|
||||
(0 1 7 6)
|
||||
|
||||
@ -31,7 +31,7 @@ FoamFile
|
||||
}
|
||||
bottomWall
|
||||
{
|
||||
type patch;
|
||||
type wall;
|
||||
nFaces 200;
|
||||
startFace 15804;
|
||||
}
|
||||
|
||||
@ -29,6 +29,8 @@ yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 1050.8;
|
||||
|
||||
yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
|
||||
|
||||
muMax muMax [ 1 -1 -1 0 0 0 0 ] 10.0;
|
||||
|
||||
rhoc rhoc [ 1 -3 0 0 0 0 0 ] 996;
|
||||
|
||||
rhod rhod [ 1 -3 0 0 0 0 0 ] 1996;
|
||||
|
||||
@ -27,7 +27,7 @@ endTime 6400;
|
||||
|
||||
deltaT 0.1;
|
||||
|
||||
writeControl runTime;
|
||||
writeControl adjustableRunTime;
|
||||
|
||||
writeInterval 20;
|
||||
|
||||
@ -45,4 +45,10 @@ timePrecision 6;
|
||||
|
||||
runTimeModifiable yes;
|
||||
|
||||
adjustTimeStep on;
|
||||
|
||||
maxCo 0.5;
|
||||
|
||||
maxDeltaT 1;
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -28,9 +28,9 @@ gradSchemes
|
||||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss limitedLinearV 1;
|
||||
div(phi,k) Gauss limitedLinear 1;
|
||||
div(phi,epsilon) Gauss limitedLinear 1;
|
||||
div(phi,U) Gauss linearUpwind grad(U);
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div(phiAlpha,Alpha) Gauss limitedLinear01 1;
|
||||
div(phiVdj,Vdj) Gauss linear;
|
||||
}
|
||||
|
||||
@ -36,14 +36,14 @@ solvers
|
||||
{
|
||||
solver PBiCG;
|
||||
preconditioner DILU;
|
||||
tolerance 1e-07;
|
||||
tolerance 1e-8;
|
||||
relTol 0.1;
|
||||
}
|
||||
|
||||
"(U|Alpha|k|epsilon)Final"
|
||||
{
|
||||
$k;
|
||||
tolerance 1e-07;
|
||||
tolerance 1e-8;
|
||||
relTol 0;
|
||||
}
|
||||
|
||||
@ -51,21 +51,21 @@ solvers
|
||||
{
|
||||
solver PCG;
|
||||
preconditioner DIC;
|
||||
tolerance 1e-07;
|
||||
tolerance 1e-8;
|
||||
relTol 0.1;
|
||||
}
|
||||
|
||||
rhoFinal
|
||||
{
|
||||
$rho;
|
||||
tolerance 1e-07;
|
||||
tolerance 1e-8;
|
||||
relTol 0;
|
||||
}
|
||||
}
|
||||
|
||||
PIMPLE
|
||||
{
|
||||
nCorrectors 2;
|
||||
nCorrectors 3;
|
||||
nNonOrthogonalCorrectors 0;
|
||||
}
|
||||
|
||||
@ -76,6 +76,7 @@ relaxationFactors
|
||||
}
|
||||
equations
|
||||
{
|
||||
"Alpha.*" 1;
|
||||
"U.*" 1;
|
||||
"k.*" 1;
|
||||
"epsilon.*" 1;
|
||||
|
||||
@ -29,6 +29,8 @@ yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 95.25;
|
||||
|
||||
yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
|
||||
|
||||
muMax muMax [ 1 -1 -1 0 0 0 0 ] 10.0;
|
||||
|
||||
rhoc rhoc [ 1 -3 0 0 0 0 0 ] 1000;
|
||||
|
||||
rhod rhod [ 1 -3 0 0 0 0 0 ] 1042;
|
||||
|
||||
@ -76,6 +76,7 @@ relaxationFactors
|
||||
}
|
||||
equations
|
||||
{
|
||||
"Alpha.*" 1;
|
||||
"U.*" 1;
|
||||
"k.*" 1;
|
||||
"epsilon.*" 1;
|
||||
|
||||
Reference in New Issue
Block a user