diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createRASTurbulence.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createRASTurbulence.H new file mode 100644 index 0000000000..bacd835655 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createRASTurbulence.H @@ -0,0 +1,178 @@ + IOdictionary RASProperties + ( + IOobject + ( + "RASProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + + Switch turbulence + ( + RASProperties.lookup("turbulence") + ); + + dictionary kEpsilonDict + ( + RASProperties.subDictPtr("kEpsilonCoeffs") + ); + + dimensionedScalar Cmu + ( + dimensionedScalar::lookupOrAddToDict + ( + "Cmu", + kEpsilonDict, + 0.09 + ) + ); + + dimensionedScalar C1 + ( + dimensionedScalar::lookupOrAddToDict + ( + "C1", + kEpsilonDict, + 1.44 + ) + ); + + dimensionedScalar C2 + ( + dimensionedScalar::lookupOrAddToDict + ( + "C2", + kEpsilonDict, + 1.92 + ) + ); + + dimensionedScalar alpha1k + ( + dimensionedScalar::lookupOrAddToDict + ( + "alpha1k", + kEpsilonDict, + 1.0 + ) + ); + + dimensionedScalar alpha1Eps + ( + dimensionedScalar::lookupOrAddToDict + ( + "alpha1Eps", + kEpsilonDict, + 0.76923 + ) + ); + + dictionary wallFunctionDict + ( + RASProperties.subDictPtr("wallFunctionCoeffs") + ); + + dimensionedScalar kappa + ( + dimensionedScalar::lookupOrAddToDict + ( + "kappa", + wallFunctionDict, + 0.41 + ) + ); + + dimensionedScalar E + ( + dimensionedScalar::lookupOrAddToDict + ( + "E", + wallFunctionDict, + 9.8 + ) + ); + + if (RASProperties.lookupOrDefault("printCoeffs", false)) + { + Info<< "kEpsilonCoeffs" << kEpsilonDict << nl + << "wallFunctionCoeffs" << wallFunctionDict << endl; + } + + + nearWallDist y(mesh); + + + Info<< "Reading field k\n" << endl; + volScalarField k + ( + IOobject + ( + "k", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field epsilon\n" << endl; + volScalarField epsilon + ( + IOobject + ( + "epsilon", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Calculating field nut2\n" << endl; + volScalarField nut2 + ( + IOobject + ( + "nut2", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + Cmu*sqr(k)/epsilon + ); + + Info<< "Calculating field nuEff1\n" << endl; + volScalarField nuEff1 + ( + IOobject + ( + "nuEff1", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sqr(Ct)*nut2 + nu1 + ); + + Info<< "Calculating field nuEff2\n" << endl; + volScalarField nuEff2 + ( + IOobject + ( + "nuEff2", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + nut2 + nu2 + ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H new file mode 100644 index 0000000000..1340587ea3 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H @@ -0,0 +1,64 @@ +if (turbulence) +{ + if (mesh.changing()) + { + y.correct(); + } + + tmp tgradU2 = fvc::grad(U2); + volScalarField G(2*nut2*(tgradU2() && dev(symm(tgradU2())))); + tgradU2.clear(); + + #include "wallFunctions.H" + + // Dissipation equation + fvScalarMatrix epsEqn + ( + fvm::ddt(epsilon) + + fvm::div(phi2, epsilon) + - fvm::Sp(fvc::div(phi2), epsilon) + - fvm::laplacian + ( + alpha1Eps*nuEff2, epsilon, + "laplacian(DepsilonEff,epsilon)" + ) + == + C1*G*epsilon/k + - fvm::Sp(C2*epsilon/k, epsilon) + ); + + #include "wallDissipation.H" + + epsEqn.relax(); + epsEqn.solve(); + + epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15)); + + + // Turbulent kinetic energy equation + fvScalarMatrix kEqn + ( + fvm::ddt(k) + + fvm::div(phi2, k) + - fvm::Sp(fvc::div(phi2), k) + - fvm::laplacian + ( + alpha1k*nuEff2, k, + "laplacian(DkEff,k)" + ) + == + G + - fvm::Sp(epsilon/k, k) + ); + kEqn.relax(); + kEqn.solve(); + + k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8)); + + //- Re-calculate turbulence viscosity + nut2 = Cmu*sqr(k)/epsilon; + + #include "wallViscosity.H" +} + +nuEff2 = nut2 + nu2; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/wallDissipation.H b/applications/solvers/multiphase/twoPhaseEulerFoam/wallDissipation.H new file mode 100644 index 0000000000..14224fd646 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/wallDissipation.H @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + wallDissipation + +Description + Set wall dissipation in the epsilon matrix + +\*---------------------------------------------------------------------------*/ + +{ + const fvPatchList& patches = mesh.boundary(); + + forAll(patches, patchi) + { + const fvPatch& p = patches[patchi]; + + if (isA(p)) + { + epsEqn.setValues + ( + p.faceCells(), + epsilon.boundaryField()[patchi].patchInternalField() + ); + } + } +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/wallFunctions.H b/applications/solvers/multiphase/twoPhaseEulerFoam/wallFunctions.H new file mode 100644 index 0000000000..c91cce9a00 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/wallFunctions.H @@ -0,0 +1,81 @@ +{ + labelList cellBoundaryFaceCount(epsilon.size(), 0); + + scalar Cmu25 = ::pow(Cmu.value(), 0.25); + scalar Cmu75 = ::pow(Cmu.value(), 0.75); + scalar kappa_ = kappa.value(); + scalar nu2_ = nu2.value(); + + const fvPatchList& patches = mesh.boundary(); + + //- Initialise the near-wall P field to zero + forAll(patches, patchi) + { + const fvPatch& currPatch = patches[patchi]; + + if (isA(currPatch)) + { + forAll(currPatch, facei) + { + label faceCelli = currPatch.faceCells()[facei]; + + epsilon[faceCelli] = 0.0; + G[faceCelli] = 0.0; + } + } + } + + //- Accumulate the wall face contributions to epsilon and G + // Increment cellBoundaryFaceCount for each face for averaging + forAll(patches, patchi) + { + const fvPatch& currPatch = patches[patchi]; + + if (isA(currPatch)) + { + const scalarField& nut2w = nut2.boundaryField()[patchi]; + + scalarField magFaceGradU(mag(U2.boundaryField()[patchi].snGrad())); + + forAll(currPatch, facei) + { + label faceCelli = currPatch.faceCells()[facei]; + + // For corner cells (with two boundary or more faces), + // epsilon and G in the near-wall cell are calculated + // as an average + + cellBoundaryFaceCount[faceCelli]++; + + epsilon[faceCelli] += + Cmu75*::pow(k[faceCelli], 1.5) + /(kappa_*y[patchi][facei]); + + G[faceCelli] += + (nut2w[facei] + nu2_)*magFaceGradU[facei] + *Cmu25*::sqrt(k[faceCelli]) + /(kappa_*y[patchi][facei]); + } + } + } + + + // perform the averaging + + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isA(curPatch)) + { + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + G[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + cellBoundaryFaceCount[faceCelli] = 1; + } + } + } +} diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/wallViscosity.H b/applications/solvers/multiphase/twoPhaseEulerFoam/wallViscosity.H new file mode 100644 index 0000000000..bd51ed7dd1 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/wallViscosity.H @@ -0,0 +1,36 @@ +{ + scalar Cmu25 = ::pow(Cmu.value(), 0.25); + scalar kappa_ = kappa.value(); + scalar E_ = E.value(); + scalar nu2_ = nu2.value(); + + const fvPatchList& patches = mesh.boundary(); + + forAll(patches, patchi) + { + const fvPatch& currPatch = patches[patchi]; + + if (isA(currPatch)) + { + scalarField& nutw = nut2.boundaryField()[patchi]; + + forAll(currPatch, facei) + { + label faceCelli = currPatch.faceCells()[facei]; + + // calculate yPlus + scalar yPlus = + Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_; + + if (yPlus > 11.6) + { + nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) -1); + } + else + { + nutw[facei] = 0.0; + } + } + } + } +}