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
+
+// ************************************************************************* //