From 898924aa48a98eb785dd4aab6e83cf0fb0548f13 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 16 Feb 2022 11:12:26 +0000 Subject: [PATCH] multiphaseEulerFoam::fvModels::interfaceTurbulenceDamping: Free-surface phase turbulence damping function Implementation of the interFoam VoFTurbulenceDamping for multiphaseEulerFoam. In this implementation no distinction is made between a dispersed phase and the interface so it is formally only applicable when interface compression is used between the phase and the other phases. Special handling for dispersed phases may be added in the future. Description Free-surface phase turbulence damping function Adds an extra source term to the mixture or phase epsilon or omega equation to reduce turbulence generated near a free-surface. The implementation is based on Reference: \verbatim Frederix, E. M. A., Mathur, A., Dovizio, D., Geurts, B. J., & Komen, E. M. J. (2018). Reynolds-averaged modeling of turbulence damping near a large-scale interface in two-phase flow. Nuclear engineering and design, 333, 122-130. \endverbatim but with an improved formulation for the coefficient \c A appropriate for unstructured meshes including those with split-cell refinement patterns. However the dimensioned length-scale coefficient \c delta remains and must be set appropriatly for the case by performing test runs and comparing with known results. Clearly this model is far from general and more research is needed in order that \c delta can be obtained directly from the interface flow and turbulence conditions. Usage Example usage: \verbatim interfaceTurbulenceDamping { type interfaceTurbulenceDamping; libs ("libmultiphaseEulerFoamFvModels.so"); phase water; // Interface turbulence damping length scale // This is a required input as described in section 3.3 of the paper delta 1e-4; } \endverbatim --- .../interfaceTurbulenceDamping.C | 324 ++++++++++++++++++ .../interfaceTurbulenceDamping.H | 223 ++++++++++++ 2 files changed, 547 insertions(+) create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C b/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C new file mode 100644 index 0000000000..b18fd6659c --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C @@ -0,0 +1,324 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022 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 . + +\*---------------------------------------------------------------------------*/ + +#include "interfaceTurbulenceDamping.H" +#include "surfaceInterpolate.H" +#include "fvcGrad.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +namespace Foam +{ + namespace fv + { + defineTypeNameAndDebug(interfaceTurbulenceDamping, 0); + + addToRunTimeSelectionTable + ( + fvModel, + interfaceTurbulenceDamping, + dictionary + ); + } +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::tmp +Foam::fv::interfaceTurbulenceDamping::interfaceFraction +( + const volScalarField& alpha +) const +{ + const fvMesh& mesh = this->mesh(); + + tmp tA + ( + volScalarField::Internal::New + ( + "A", + mesh, + dimensionedScalar(dimless, 0) + ) + ); + volScalarField::Internal& A = tA.ref(); + + const surfaceVectorField& Sf = mesh.Sf(); + const labelUList& own = mesh.owner(); + const labelUList& nei = mesh.neighbour(); + + const surfaceScalarField alphaf(fvc::interpolate(alpha)); + + const volVectorField gradAlpha(fvc::grad(alpha)); + const volVectorField::Internal n + ( + gradAlpha()/(mag(gradAlpha()) + phase_.fluid().deltaN()) + ); + + const scalarField& ialpha = alpha; + const scalarField& ialphaf = alphaf; + scalarField sumnSf(mesh.nCells(), 0); + + forAll(own, facei) + { + { + const scalar nSf(mag(n[own[facei]] & Sf[facei])); + A[own[facei]] += nSf*(ialphaf[facei] - ialpha[own[facei]]); + sumnSf[own[facei]] += nSf; + } + { + const scalar nSf(mag(n[nei[facei]] & Sf[facei])); + A[nei[facei]] += nSf*(ialphaf[facei] - ialpha[nei[facei]]); + sumnSf[nei[facei]] += nSf; + } + } + + forAll(mesh.boundary(), patchi) + { + const labelUList& own = mesh.boundary()[patchi].faceCells(); + const fvsPatchScalarField& palphaf = alphaf.boundaryField()[patchi]; + + forAll(mesh.boundary()[patchi], facei) + { + const scalar nSf(mag(n[own[facei]] & Sf[facei])); + A[own[facei]] += nSf*(palphaf[facei] - ialpha[own[facei]]); + sumnSf[own[facei]] += nSf; + } + } + + scalarField& a = A.field(); + forAll(a, i) + { + if (sumnSf[i] > small) + { + a[i] = 2*mag(a[i])/sumnSf[i]; + } + else + { + a[i] = 0; + } + } + + return tA; +} + + +template +void Foam::fv::interfaceTurbulenceDamping::addRhoSup +( + const RhoType& rho, + fvMatrix& eqn, + const word& fieldName +) const +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + const phaseSystem::phaseModelPartialList& movingPhases = + phase_.fluid().movingPhases(); + + volScalarField::Internal aSqrnu + ( + movingPhases[0]*sqr(movingPhases[0].thermo().nu()()()) + ); + + for (label phasei=1; phasei(IOobject::groupName("alpha", phaseName_)) + ), + turbulence_ + ( + mesh.lookupObject + ( + IOobject::groupName + ( + momentumTransportModel::typeName, + phaseName_ + ) + ) + ), + C2_("C2", dimless, 0), + betaStar_("betaStar", dimless, 0), + beta_("beta", dimless, 0) +{ + const word epsilonName(IOobject::groupName("epsilon", phaseName_)); + const word omegaName(IOobject::groupName("omega", phaseName_)); + + if (mesh.foundObject(epsilonName)) + { + fieldName_ = epsilonName; + C2_.read(turbulence_.coeffDict()); + } + else if (mesh.foundObject(omegaName)) + { + fieldName_ = omegaName; + betaStar_.read(turbulence_.coeffDict()); + + // Read beta for k-omega models or beta1 for k-omega SST + if (turbulence_.coeffDict().found("beta")) + { + beta_.read(turbulence_.coeffDict()); + } + else + { + beta_ = + dimensionedScalar("beta1", dimless, turbulence_.coeffDict()); + } + } + else + { + FatalIOErrorInFunction(dict) + << "Cannot find either " << epsilonName << " or " << omegaName + << " field for fvModel " << typeName << exit(FatalIOError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::wordList Foam::fv::interfaceTurbulenceDamping::addSupFields() const +{ + return wordList(1, fieldName_); +} + + +void Foam::fv::interfaceTurbulenceDamping::addSup +( + fvMatrix& eqn, + const word& fieldName +) const +{ + addRhoSup(one(), eqn, fieldName); +} + + +void Foam::fv::interfaceTurbulenceDamping::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName +) const +{ + addRhoSup(rho(), eqn, fieldName); +} + + +void Foam::fv::interfaceTurbulenceDamping::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName +) const +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + volScalarField::Internal aSqrnu + ( + alpha*sqr(phase_.thermo().nu()()()) + ); + + if (fieldName == IOobject::groupName("epsilon", phaseName_)) + { + eqn += rho()*interfaceFraction(alpha) + *C2_*aSqrnu*turbulence_.k()()/pow4(delta_); + } + else if (fieldName == IOobject::groupName("omega", phaseName_)) + { + eqn += rho()*interfaceFraction(alpha) + *beta_*aSqrnu/(sqr(betaStar_)*pow4(delta_)); + } + else + { + FatalErrorInFunction + << "Support for field " << fieldName << " is not implemented" + << exit(FatalError); + } +} + + +void Foam::fv::interfaceTurbulenceDamping::updateMesh(const mapPolyMesh&) +{} + + +void Foam::fv::interfaceTurbulenceDamping::distribute +( + const mapDistributePolyMesh& +) +{} + + +bool Foam::fv::interfaceTurbulenceDamping::movePoints() +{ + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H b/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H new file mode 100644 index 0000000000..62479f1e4c --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H @@ -0,0 +1,223 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022 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 . + +Class + Foam::fv::interfaceTurbulenceDamping + +Description + Free-surface phase turbulence damping function + + Adds an extra source term to the mixture or phase epsilon or omega + equation to reduce turbulence generated near a free-surface. The + implementation is based on + + Reference: + \verbatim + Frederix, E. M. A., Mathur, A., Dovizio, D., Geurts, B. J., + & Komen, E. M. J. (2018). + Reynolds-averaged modeling of turbulence damping + near a large-scale interface in two-phase flow. + Nuclear engineering and design, 333, 122-130. + \endverbatim + + but with an improved formulation for the coefficient \c A appropriate for + unstructured meshes including those with split-cell refinement patterns. + However the dimensioned length-scale coefficient \c delta remains and must + be set appropriatly for the case by performing test runs and comparing with + known results. Clearly this model is far from general and more research is + needed in order that \c delta can be obtained directly from the interface + flow and turbulence conditions. + +Usage + Example usage: + \verbatim + interfaceTurbulenceDamping + { + type interfaceTurbulenceDamping; + + libs ("libmultiphaseEulerFoamFvModels.so"); + + phase water; + + // Interface turbulence damping length scale + // This is a required input as described in section 3.3 of the paper + delta 1e-4; + } + \endverbatim + +SourceFiles + interfaceTurbulenceDamping.C + +\*---------------------------------------------------------------------------*/ + +#ifndef interfaceTurbulenceDamping_H +#define interfaceTurbulenceDamping_H + +#include "fvModel.H" +#include "phaseSystem.H" +#include "phaseCompressibleMomentumTransportModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class interfaceTurbulenceDamping Declaration +\*---------------------------------------------------------------------------*/ + +class interfaceTurbulenceDamping +: + public fvModel +{ + // Private Data + + //- The name of the Lagrangian phase + word phaseName_; + + //- Field name + word fieldName_; + + //- Interface turbulence damping length scale + // This is a required input as described in section 3.3 of the paper + dimensionedScalar delta_; + + //- Reference to the phase + const phaseModel& phase_; + + //- Reference to the mixture turbulence model + const phaseCompressible::momentumTransportModel& turbulence_; + + // Turbulence model coefficients + + dimensionedScalar C2_; + dimensionedScalar betaStar_; + dimensionedScalar beta_; + + + // Private Member Functions + + //- Interface fraction in a cell + tmp interfaceFraction + ( + const volScalarField& alpha + ) const; + + //- Add explicit contribution to incompressible or compressible + // mixture epsilon or omega equation + template + void addRhoSup + ( + const RhoType& rho, + fvMatrix& eqn, + const word& fieldName + ) const; + + +public: + + //- Runtime type information + TypeName("interfaceTurbulenceDamping"); + + + // Constructors + + //- Construct from explicit source name and mesh + interfaceTurbulenceDamping + ( + const word& sourceName, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + //- Disallow default bitwise copy construction + interfaceTurbulenceDamping + ( + const interfaceTurbulenceDamping& + ) = delete; + + + // Member Functions + + //- Return the list of fields for which the option adds source term + // to the transport equation + virtual wordList addSupFields() const; + + //- Add explicit contribution to mixture epsilon or omega equation + virtual void addSup + ( + fvMatrix& eqn, + const word& fieldName + ) const; + + //- Add explicit contribution to compressible + // mixture epsilon or omega equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName + ) const; + + //- Add explicit contribution to phase epsilon or omega equation + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName + ) const; + + + // Mesh changes + + //- Update for mesh changes + virtual void updateMesh(const mapPolyMesh&); + + //- Update mesh corresponding to the given distribution map + virtual void distribute(const mapDistributePolyMesh&); + + //- Update for mesh motion + virtual bool movePoints(); + + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const interfaceTurbulenceDamping&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //