diff --git a/bin/foamMonitor b/bin/foamMonitor
index 2d9fed9ff8ba0f0cf9dee4d824e231279a85c7a8..a15bd5b168ffc43c125cf7ca03c5082e77e47e0f 100755
--- a/bin/foamMonitor
+++ b/bin/foamMonitor
@@ -6,7 +6,8 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #-------------------------------------------------------------------------------
-#                           | Copyright (C) 2015 OpenFOAM Foundation
+#    Copyright (C) 2015 OpenFOAM Foundation
+#    Copyright (C) 2019 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM.
@@ -40,11 +41,12 @@ usage() {
 
 Usage: ${0##*/} [OPTION] <file>
 options:
-  -h | -help            print the usage
+  -h | -help            prints the usage
   -i | -idle <time>     stops if <file> unchanging for <time> sec (default = 60)
   -l | -logscale        plots data (y-axis) on log scale, e.g. for residuals
   -r | -refresh <time>  refreshes display every <time> sec (default = 10)
   -y | -yrange <range>  sets data (y-axis) <range>, format "[0:1]"
+  -g | -grid            draws grid lines on the plot
 
 Monitor data with Gnuplot from time-value(s) graphs written by OpenFOAM
 e.g. by functionObjects
@@ -62,6 +64,7 @@ plotFileHeader() {
 set term x11 1 font "helvetica,17" linewidth 1.5 persist noraise
 $LOGSCALE
 $YRANGE
+$GRID
 set title "Data Monitoring"
 set xlabel "$XLABEL"
 plot \\
@@ -82,6 +85,7 @@ IDLE=60
 REFRESH=10
 LOGSCALE=""
 YRANGE=""
+GRID=""
 GNUPLOT=$(which gnuplot)
 ! [ "x$GNUPLOT" = "x" ] || usage "Gnuplot not installed"
 
@@ -111,6 +115,10 @@ do
         YRANGE="set yrange $2"
         shift 2
         ;;
+    -g | -grid)
+        GRID="set grid"
+        shift 1
+        ;;
     -*)
         usage "unknown option: '$*'"
         ;;
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
index 1887e54180a7fd3204c1b1cf43740c28518b0f73..86a19cd1336dd7bdd01e7f53dc23342212fc2ad1 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
@@ -94,6 +94,9 @@ makeRASModel(LRR);
 #include "SSG.H"
 makeRASModel(SSG);
 
+#include "kEpsilonPhitF.H"
+makeRASModel(kEpsilonPhitF);
+
 
 // -------------------------------------------------------------------------- //
 // LES models
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
index 2d846cddc67a5e394b429bccf74fe49a14867b05..6004df24d32410bf4b3c93fd05e63609ab26b77a 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
@@ -90,6 +90,9 @@ makeRASModel(LRR);
 #include "SSG.H"
 makeRASModel(SSG);
 
+#include "kEpsilonPhitF.H"
+makeRASModel(kEpsilonPhitF);
+
 
 // -------------------------------------------------------------------------- //
 // LES models
diff --git a/src/TurbulenceModels/turbulenceModels/Make/files b/src/TurbulenceModels/turbulenceModels/Make/files
index 7a6ad66294741a48ff7c7138530747cf59b4348e..8a3c77f09e57ecf827cfbc3b4ba6210bf7803667 100644
--- a/src/TurbulenceModels/turbulenceModels/Make/files
+++ b/src/TurbulenceModels/turbulenceModels/Make/files
@@ -1,5 +1,6 @@
 turbulenceModel.C
 RAS/v2f/v2fBase.C
+RAS/kEpsilonPhitF/kEpsilonPhitFBase.C
 
 LESdelta = LES/LESdeltas
 
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C
new file mode 100644
index 0000000000000000000000000000000000000000..2da4eee3edfee698aa5029a104b7d9848ec1a773
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C
@@ -0,0 +1,513 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kEpsilonPhitF.H"
+#include "fvOptions.H"
+#include "bound.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void kEpsilonPhitF<BasicTurbulenceModel>::correctNut()
+{
+    // (LUU:p. 173)
+    this->nut_ = Cmu_*phit_*k_*T_;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+template<class BasicTurbulenceModel>
+volScalarField kEpsilonPhitF<BasicTurbulenceModel>::Ts() const
+{
+    // (LUU:Eq. 7)
+    return
+        max
+        (
+            k_/epsilon_,
+            CT_*sqrt
+            (
+                max
+                (
+                    this->nu(),
+                    dimensionedScalar(this->nu()().dimensions(), Zero)
+                )/epsilon_
+            )
+        );
+}
+
+
+template<class BasicTurbulenceModel>
+volScalarField kEpsilonPhitF<BasicTurbulenceModel>::Ls() const
+{
+    // (LUU:Eq. 7)
+    return
+        CL_*max
+        (
+            pow(k_, 1.5)/epsilon_,
+            Ceta_*pow025
+            (
+                pow3
+                (
+                    max
+                    (
+                        this->nu(),
+                        dimensionedScalar(this->nu()().dimensions(), Zero)
+                    )
+                )/epsilon_
+            )
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+kEpsilonPhitF<BasicTurbulenceModel>::kEpsilonPhitF
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+    kEpsilonPhitFBase(),
+
+    Cmu_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.22
+        )
+    ),
+    Ceps1a_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ceps1a",
+            this->coeffDict_,
+            1.4
+        )
+    ),
+    Ceps1b_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ceps1b",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    Ceps1c_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ceps1c",
+            this->coeffDict_,
+            0.05
+        )
+    ),
+    Ceps2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ceps2",
+            this->coeffDict_,
+            1.9
+        )
+    ),
+    Cf1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cf1",
+            this->coeffDict_,
+            1.4
+        )
+    ),
+    Cf2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cf2",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    CL_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "CL",
+            this->coeffDict_,
+            0.25
+        )
+    ),
+    Ceta_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ceta",
+            this->coeffDict_,
+            110.0
+        )
+    ),
+    CT_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "CT",
+            this->coeffDict_,
+            6.0
+        )
+    ),
+    sigmaK_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaK",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    sigmaEps_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaEps",
+            this->coeffDict_,
+            1.3
+        )
+    ),
+    sigmaPhit_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaPhit",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    phit_
+    (
+        IOobject
+        (
+            IOobject::groupName("phit", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    f_
+    (
+        IOobject
+        (
+            IOobject::groupName("f", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    T_
+    (
+        IOobject
+        (
+            "T",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh_,
+        dimensionedScalar(dimTime, Zero)
+    ),
+
+    phitMin_(dimensionedScalar("phitMin", phit_.dimensions(), SMALL)),
+    fMin_(dimensionedScalar("fMin", f_.dimensions(), SMALL)),
+    TMin_(dimensionedScalar("TMin", dimTime, SMALL)),
+    L2Min_(dimensionedScalar("L2Min", sqr(dimLength), SMALL))
+{
+    bound(k_, this->kMin_);
+    bound(epsilon_, this->epsilonMin_);
+    bound(phit_, phitMin_);
+    bound(f_, fMin_);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+
+    if
+    (
+        mag(sigmaK_.value()) < VSMALL
+     || mag(sigmaEps_.value()) < VSMALL
+     || mag(sigmaPhit_.value()) < VSMALL
+    )
+    {
+        FatalErrorInFunction
+            << "Non-zero values are required for the model constants:" << nl
+            << "sigmaK = " << sigmaK_ << nl
+            << "sigmaEps = " << sigmaEps_ << nl
+            << "sigmaPhit = " << sigmaPhit_ << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool kEpsilonPhitF<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        Ceps1a_.readIfPresent(this->coeffDict());
+        Ceps1b_.readIfPresent(this->coeffDict());
+        Ceps1c_.readIfPresent(this->coeffDict());
+        Ceps2_.readIfPresent(this->coeffDict());
+        Cf1_.readIfPresent(this->coeffDict());
+        Cf2_.readIfPresent(this->coeffDict());
+        CL_.readIfPresent(this->coeffDict());
+        Ceta_.readIfPresent(this->coeffDict());
+        CT_.readIfPresent(this->coeffDict());
+        sigmaK_.readIfPresent(this->coeffDict());
+        sigmaEps_.readIfPresent(this->coeffDict());
+        sigmaPhit_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicTurbulenceModel>
+void kEpsilonPhitF<BasicTurbulenceModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Construct local convenience references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+    const volVectorField& U = this->U_;
+    volScalarField& nut = this->nut_;
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
+
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))
+    );
+
+    tmp<volSymmTensorField> tS(symm(fvc::grad(U)));
+    volScalarField G(this->GName(), nut*(2.0*(dev(tS()) && tS())));
+    tS.clear();
+
+    T_ = Ts();
+    bound(T_, TMin_);
+
+    const volScalarField L2(type() + "L2", sqr(Ls()) + L2Min_);
+
+    const volScalarField::Internal Ceps1
+    (
+        "Ceps1",
+        Ceps1a_*(Ceps1b_ + Ceps1c_*sqrt(1.0/phit_()))
+    );
+
+
+    // Update epsilon and G at the wall
+    epsilon_.boundaryFieldRef().updateCoeffs();
+
+    // Turbulent kinetic energy dissipation rate equation (LUU:Eq. 4)
+    // k/T ~ epsilon
+    tmp<fvScalarMatrix> epsEqn
+    (
+        fvm::ddt(alpha, rho, epsilon_)
+      + fvm::div(alphaRhoPhi, epsilon_)
+      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
+      ==
+        alpha()*rho()*Ceps1*G()/T_()
+      - fvm::SuSp
+        (
+            (2.0/3.0*Ceps1)*(alpha()*rho()*divU),
+            epsilon_
+        )
+      - fvm::Sp(alpha()*rho()*Ceps2_/T_(), epsilon_)
+      + fvOptions(alpha, rho, epsilon_)
+    );
+
+    epsEqn.ref().relax();
+    fvOptions.constrain(epsEqn.ref());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
+    solve(epsEqn);
+    fvOptions.correct(epsilon_);
+    bound(epsilon_, this->epsilonMin_);
+
+
+    // Turbulent kinetic energy equation (LUU:Eq. 3)
+    // epsilon/k ~ 1/Ts
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaRhoPhi, k_)
+      - fvm::laplacian(alpha*rho*DkEff(), k_)
+      ==
+        alpha()*rho()*G()
+      - fvm::SuSp(2.0/3.0*alpha()*rho()*divU, k_)
+      - fvm::Sp(alpha()*rho()/T_(), k_)
+      + fvOptions(alpha, rho, k_)
+    );
+
+    kEqn.ref().relax();
+    fvOptions.constrain(kEqn.ref());
+    solve(kEqn);
+    fvOptions.correct(k_);
+    bound(k_, this->kMin_);
+
+
+    // Elliptic relaxation function equation (LUU:Eq. 18)
+    // All source terms are non-negative functions (LUU:p. 176)
+    tmp<fvScalarMatrix> fEqn
+    (
+      - fvm::laplacian(f_)
+      ==
+      - fvm::Sp(1.0/L2(), f_)
+      - (
+            (Cf1_ - 1.0)*(phit_() - 2.0/3.0)/T_()
+           -(Cf2_*G())/k_()
+           +(Cf2_*(2.0/3.0)*divU)
+           -(2.0*this->nu()*(fvc::grad(phit_) & fvc::grad(k_)))()/k_()
+           -(this->nu()*fvc::laplacian(phit_))()
+        )/L2()
+    );
+
+    fEqn.ref().relax();
+    solve(fEqn);
+    bound(f_, fMin_);
+
+
+    // Normalised wall-normal fluctuating velocity scale equation (LUU:Eq. 17)
+    // All source terms are non-negative functions (LUU:p. 176)
+    tmp<fvScalarMatrix> phitEqn
+    (
+        fvm::ddt(alpha, rho, phit_)
+      + fvm::div(alphaRhoPhi, phit_)
+      - fvm::laplacian(alpha*rho*DphitEff(), phit_)
+      ==
+        alpha()*rho()*f_()
+      - fvm::SuSp
+        (
+            alpha()*rho()*
+            (
+                G()/k_()
+              - (2.0/3.0)*divU
+              - (2.0*nut*(fvc::grad(phit_) & fvc::grad(k_)))()
+                /(k_()*sigmaPhit_*phit_())
+            )
+          , phit_
+        )
+      + fvOptions(alpha, rho, phit_)
+    );
+
+    phitEqn.ref().relax();
+    fvOptions.constrain(phitEqn.ref());
+    solve(phitEqn);
+    fvOptions.correct(phit_);
+    bound(phit_, phitMin_);
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.H b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.H
new file mode 100644
index 0000000000000000000000000000000000000000..18d07e829067478f248b62ff08cdb3259cefc6d2
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.H
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::RASModels::kEpsilonPhitF
+
+Group
+    grpRASTurbulence
+
+Description
+    The k-epsilon-phit-f turbulence closure model for incompressible and
+    compressible flows.
+
+    The model is a three-transport-equation linear-eddy-viscosity turbulence
+    closure model alongside an elliptic relaxation equation:
+      - Turbulent kinetic energy, \c k,
+      - Turbulent kinetic energy dissipation rate, \c epsilon,
+      - Normalised wall-normal fluctuating velocity scale, \c phit,
+      - Elliptic relaxation factor, \c f.
+
+    Reference:
+    \verbatim
+        Standard model (Tag:LUU):
+            Laurence, D. R., Uribe, J. C., & Utyuzhnikov, S. V. (2005).
+            A robust formulation of the v2−f model.
+            Flow, Turbulence and Combustion, 73(3-4), 169–185.
+            DOI:10.1007/s10494-005-1974-8
+    \endverbatim
+
+    The default model coefficients are (LUU:Eqs. 19-20):
+    \verbatim
+        kEpsilonPhitFCoeffs
+        {
+            Cmu         0.22,    // Turbulent viscosity constant
+            Ceps1a      1.4,     // Model constant for epsilon
+            Ceps1b      1.0,     // Model constant for epsilon
+            Ceps1c      0.05,    // Model constant for epsilon
+            Ceps2       1.9,     // Model constant for epsilon
+            Cf1         1.4,     // Model constant for f
+            Cf2         0.3,     // Model constant for f
+            CL          0.25,    // Model constant for L
+            Ceta        110.0,   // Model constant for L
+            CT          6.0,     // Model constant for T
+            sigmaK      1.0,     // Turbulent Prandtl number for k
+            sigmaEps    1.3,     // Turbulent Prandtl number for epsilon
+            sigmaPhit   1.0,     // Turbulent Prandtl number for phit = sigmaK
+        }
+    \endverbatim
+
+Note
+    The name of the original variable replacing 'v2' is 'phi' (LUU:Eq. 14).
+    However, the name 'phi' preexisted in OpenFOAM; therefore, this name was
+    replaced by 'phit' herein.
+
+SourceFiles
+    kEpsilonPhitF.C
+
+SeeAlso
+    v2f.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kEpsilonPhitF_H
+#define kEpsilonPhitF_H
+
+#include "kEpsilonPhitFBase.H"
+#include "RASModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class kEpsilonPhitF Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class kEpsilonPhitF
+:
+    public eddyViscosity<RASModel<BasicTurbulenceModel>>,
+    public kEpsilonPhitFBase
+{
+    // Private Member Functions
+
+        //- No copy construct
+        kEpsilonPhitF(const kEpsilonPhitF&) = delete;
+
+        //- No copy assignment
+        void operator=(const kEpsilonPhitF&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+            dimensionedScalar Ceps1a_;
+            dimensionedScalar Ceps1b_;
+            dimensionedScalar Ceps1c_;
+            dimensionedScalar Ceps2_;
+            dimensionedScalar Cf1_;
+            dimensionedScalar Cf2_;
+            dimensionedScalar CL_;
+            dimensionedScalar Ceta_;
+            dimensionedScalar CT_;
+            dimensionedScalar sigmaK_;
+            dimensionedScalar sigmaEps_;
+            dimensionedScalar sigmaPhit_;
+
+
+        // Fields
+
+            //- Turbulent kinetic energy [m2/s2]
+            volScalarField k_;
+
+            //- Turbulent kinetic energy dissipation rate [m2/s3]
+            volScalarField epsilon_;
+
+            //- Normalised wall-normal fluctuating velocity scale [-]
+            volScalarField phit_;
+
+            //- Elliptic relaxation factor [1/s]
+            volScalarField f_;
+
+            //- Turbulent time scale [s]
+            volScalarField T_;
+
+
+        // Bounding values
+
+            dimensionedScalar phitMin_;
+            dimensionedScalar fMin_;
+            dimensionedScalar TMin_;
+            dimensionedScalar L2Min_;
+
+
+    // Protected Member Functions
+
+        //- Update nut with the latest available k, phit, and T
+        virtual void correctNut();
+
+        //- Return the turbulent time scale, T
+        volScalarField Ts() const;
+
+        //- Return the turbulent length scale, L
+        volScalarField Ls() const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kEpsilonPhitF");
+
+
+    // Constructors
+
+        //- Construct from components
+        kEpsilonPhitF
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kEpsilonPhitF() = default;
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DkEff",
+                    this->nut_/sigmaK_ + this->nu()
+                )
+            );
+        }
+
+        //- Return the effective diffusivity for epsilon
+        tmp<volScalarField> DepsilonEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DepsilonEff",
+                    this->nut_/sigmaEps_ + this->nu()
+                )
+            );
+        }
+
+        //- Return the effective diffusivity for phit
+        tmp<volScalarField> DphitEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DphitEff",
+                    this->nut_/sigmaPhit_ + this->nu()
+                )
+            );
+        }
+
+        //- Return the turbulent kinetic energy field
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulent kinetic energy dissipation rate field
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
+        //- Return the normalised wall-normal fluctuating velocity scale field
+        virtual tmp<volScalarField> phit() const
+        {
+            return phit_;
+        }
+
+        //- Return the elliptic relaxation factor field
+        virtual tmp<volScalarField> f() const
+        {
+            return f_;
+        }
+
+        //- Solve the transport equations and correct the turbulent viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "kEpsilonPhitF.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
\ No newline at end of file
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..e4a56e54bbe7bad2f5838fdb21db63322134d334
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kEpsilonPhitFBase.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+    defineTypeNameAndDebug(kEpsilonPhitFBase, 0);
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.H b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..622bc3254cdd9d2f60dbcd85ce65a6bd80553535
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitFBase.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::RASModels::kEpsilonPhitFBase
+
+Group
+    grpRASTurbulence
+
+Description
+    Abstract base-class for the \c k-epsilon-phit-f model to provide boundary
+    condition access to the \c phit and \c f fields.
+
+See also
+    Foam::RASModels::kEpsilonPhitF
+
+SourceFiles
+    kEpsilonPhitFBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kEpsilonPhitFBase_H
+#define kEpsilonPhitFBase_H
+
+#include "RASModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class kEpsilonPhitFBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class kEpsilonPhitFBase
+{
+public:
+
+    //- Runtime type information
+    TypeName("kEpsilonPhitFBase");
+
+
+    // Constructors
+
+        kEpsilonPhitFBase()
+        {}
+
+
+    //- Destructor
+    virtual ~kEpsilonPhitFBase()
+    {}
+
+
+    // Member Functions
+
+        //- Return the normalised wall-normal fluctuating velocity scale field
+        virtual tmp<volScalarField> phit() const = 0;
+
+        //- Return the elliptic relaxation factor field
+        virtual tmp<volScalarField> f() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
index da3175bf5137fe0100ed928dd5a81918c379f901..3573bcd319026cf54a5961fea3b023b2489887dc 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
@@ -29,6 +29,7 @@ License
 #include "fWallFunctionFvPatchScalarField.H"
 #include "nutWallFunctionFvPatchScalarField.H"
 #include "v2f.H"
+#include "kEpsilonPhitF.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -111,7 +112,6 @@ void fWallFunctionFvPatchScalarField::updateCoeffs()
             internalField().group()
         )
     );
-    const v2fBase& v2fModel = refCast<const v2fBase>(turbModel);
 
     const nutWallFunctionFvPatchScalarField& nutw =
         nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
@@ -121,43 +121,60 @@ void fWallFunctionFvPatchScalarField::updateCoeffs()
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    const tmp<volScalarField> tk = turbModel.k();
-    const volScalarField& k = tk();
+    scalarField& f = *this;
 
-    const tmp<volScalarField> tepsilon = turbModel.epsilon();
-    const volScalarField& epsilon = tepsilon();
+    if (isA<v2fBase>(turbModel))
+    {
+        const v2fBase& v2fModel = refCast<const v2fBase>(turbModel);
 
-    const tmp<volScalarField> tv2 = v2fModel.v2();
-    const volScalarField& v2 = tv2();
+        const tmp<volScalarField> tk = turbModel.k();
+        const volScalarField& k = tk();
 
-    const scalar Cmu25 = pow025(nutw.Cmu());
-    const scalar N = 6.0;
+        const tmp<volScalarField> tepsilon = turbModel.epsilon();
+        const volScalarField& epsilon = tepsilon();
 
-    scalarField& f = *this;
+        const tmp<volScalarField> tv2 = v2fModel.v2();
+        const volScalarField& v2 = tv2();
 
-    // Set f wall values
-    forAll(f, facei)
-    {
-        const label celli = patch().faceCells()[facei];
+        const scalar Cmu25 = pow025(nutw.Cmu());
+        const scalar N = 6.0;
 
-        const scalar uTau = Cmu25*sqrt(k[celli]);
+        // Set f wall values
+        forAll(f, facei)
+        {
+            const label celli = patch().faceCells()[facei];
 
-        const scalar yPlus = uTau*y[facei]/nuw[facei];
+            const scalar uTau = Cmu25*sqrt(k[celli]);
 
-        if (nutw.yPlusLam() < yPlus)
-        {
-            const scalar v2c = v2[celli];
-            const scalar epsc = epsilon[celli];
-            const scalar kc = k[celli];
+            const scalar yPlus = uTau*y[facei]/nuw[facei];
 
-            f[facei] = N*v2c*epsc/(sqr(kc) + ROOTVSMALL);
-            f[facei] /= sqr(uTau) + ROOTVSMALL;
-        }
-        else
-        {
-            f[facei] = 0.0;
+            if (nutw.yPlusLam() < yPlus)
+            {
+                const scalar v2c = v2[celli];
+                const scalar epsc = epsilon[celli];
+                const scalar kc = k[celli];
+
+                f[facei] = N*v2c*epsc/(sqr(kc) + ROOTVSMALL);
+                f[facei] /= sqr(uTau) + ROOTVSMALL;
+            }
+            else
+            {
+                f[facei] = 0.0;
+            }
         }
     }
+    else if (isA<kEpsilonPhitFBase>(turbModel))
+    {
+        // (LUU:p. 176)
+        f = 0.0;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "The RAS model is neither the v2f nor kEpsilonPhitF model. "
+            << "Therefore, fWallFunction is not usable." << nl
+            << exit(FatalError);
+    }
 
     fixedValueFvPatchField<scalar>::updateCoeffs();
 
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
index 2462d8226fdd69c8587ebb45e0048f9f017c40c1..88fca25ff1ec16ed5c8fe0218d7a31f2b9f15908 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
@@ -32,24 +32,24 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the elliptic
-    relaxation factor, \c f, which is executed in the \c v2-f eddy viscosity
-    turbulence model. The condition is applicable for low- and high-Reynolds
-    number turbulent flow cases.
-
-    For \c f, the viscous sublayer and log-law region blending approaches are
-    claimed to be inviable (Popovac and Hanjalić (2007), p. 194). Therefore,
-    the only boundary condition blending mode is the stepwise mode
-    where the viscous sublayer and log-law region contributions switch over
-    \c y+ value derived from the \c kappa and \c E.
+    relaxation factor, \c f, which is executed in the \c v2-f based
+    turbulence closure models for low- and high-Reynolds number
+    turbulent flow cases.
 
     Reference:
     \verbatim
-        Remark on the blending approach:
-            Popovac, M., & Hanjalić, K. (2007).
-            Compound wall treatment for RANS computation of complex
-            turbulent flows and heat transfer.
-            Flow, turbulence and combustion, 78(2), 177-202.
-            doi:10.1007/s10494-006-9067-x
+    Remark on the blending approach for f (tag:PH):
+        Popovac, M., & Hanjalić, K. (2007).
+        Compound wall treatment for RANS computation of complex
+        turbulent flows and heat transfer.
+        Flow, turbulence and combustion, 78(2), 177-202.
+        DOI:10.1007/s10494-006-9067-x
+
+    Wall-boundary expression for f in kEpsilonPhitF model (tag:LUU):
+        Laurence, D. R., Uribe, J. C., & Utyuzhnikov, S. V. (2005).
+        A robust formulation of the v2−f model.
+        Flow, Turbulence and Combustion, 73(3-4), 169–185.
+        DOI:10.1007/s10494-005-1974-8
     \endverbatim
 
 Usage
@@ -69,6 +69,12 @@ Note
     the specified \c nutWallFunction in order to ensure that each patch
     possesses the same set of values for these coefficients.
 
+    For \c f, the viscous and inertial sublayer blending approaches were
+    claimed to be inviable (PH:p. 194). Therefore, the only blending mode
+    for the v2-f model is the stepwise mode where the viscous and inertial
+    sublayer contributions switch over a sublayer-intersection value of
+    \c y+ estimated from the \c kappa and \c E.
+
 See also
     Foam::fixedValueFvPatchField
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
index 698a3277d2239af4bc8d76345480074262a4c375..996575e5e2e56abf7426957a05db3f7f3e515278 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.C
@@ -57,7 +57,7 @@ Foam::turbulentInletFvPatchField<Type>::turbulentInletFvPatchField
     ranGen_(label(0)),
     fluctuationScale_(dict.get<Type>("fluctuationScale")),
     referenceField_("referenceField", dict, p.size()),
-    alpha_(dict.lookupOrDefault<scalar>("alpha", 0.1)),
+    alpha_(dict.getOrDefault<scalar>("alpha", 0.1)),
     curTimeIndex_(-1)
 {
     if (dict.found("value"))
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H
index 018354b17ffec89045eedef3acb8c42c7f68baf4..9119b7e2db06525421f37efde99d52b9b98c4a1e 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H
@@ -30,44 +30,68 @@ Group
     grpInletBoundaryConditions
 
 Description
-    This boundary condition generates a fluctuating inlet condition by adding
-    a random component to a reference (mean) field.
+    This boundary condition produces spatiotemporal-variant field by summing
+    a set of pseudo-random numbers and a given spatiotemporal-invariant mean
+    field. The field can be any type, e.g. scalarField. At a single point and
+    time, all components are summed by the same random number, e.g. velocity
+    components (u, v, w) are summed by the same random number, p; thus, output
+    is (u+p, v+p, w+p).
+
+    The pseudo-random number generator obeys the probability density function
+    of the uniform distribution constrained by the range [0:1]. The seed for
+    the random number generator is hard-coded; therefore, it will produce the
+    same sequence of random numbers at every execution.
 
     \f[
-        x_p = (1 - \alpha) x_p^{n-1} + \alpha (x_{ref} + s C_{RMS} x_{ref})
+        x_p = (1 - \alpha) x_p^{n - 1} + \alpha (x_{ref} + c s R |x_{ref}|)
     \f]
-
     where
 
     \vartable
-        x_p     | patch values
-        x_{ref} | reference patch values
+        x_p     | patch field
+        x_{ref} | spatiotemporal-invariant patch scalar
         n       | time level
-        \alpha  | fraction of new random component added to previous time value
-        C_{RMS} | RMS coefficient
-        s       | fluctuation scale
+        \alpha  | a scalar attempting to build two-temporal-point correlations
+                  by heuristically adding a fraction of the new random component
+                  to the previous time patch field
+        c       | a heuristic automatically calculated correction term
+                  to compensate energy level losses due to the alpha scalar
+        R       | pseudo-random number [HARD-CODED seed]
+        s       | fluctuation scale (proportional to the xRef)
     \endvartable
 
 Usage
     \table
-        Property     | Description             | Required    | Default value
+        Property         | Description             | Required    | Default value
         fluctuationScale | RMS fluctuation scale (fraction of mean) | yes |
-        referenceField | reference (mean) field | yes        |
-        alpha | fraction of new random component added to previous| no| 0.1
+        referenceField   | reference (mean) field  | yes         |
+        alpha | fraction of new random component added to previous  | no | 0.1
     \endtable
 
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
-        type            turbulentInlet;
-        fluctuationScale 0.1;
-        referenceField  uniform 10;
-        alpha           0.1;
+        // Mandatory entries
+        type             turbulentInlet;
+        fluctuationScale 0.1;               // the term `s` above
+        referenceField   uniform 10;        // the term `xRef` above
+
+        // Optional entries
+        alpha            0.1;               // the term `alpha` above
     }
     \endverbatim
 
-See also
+Note
+    This boundary condition should not be used for DES or LES computations as a
+    turbulent velocity inflow condition, because the BC will not produce
+    turbulence-alike time-series, and will decay almost immediately downstream
+    of the inlet boundary although its historical name suggests the opposite.
+
+    Nevertheless, the BC may be still used for other applications, e.g. as a
+    uniform-random noise source in aeroacoustics.
+
+SeeAlso
     Foam::fixedValueFvPatchField
 
 SourceFiles
@@ -95,7 +119,7 @@ class turbulentInletFvPatchField
 :
     public fixedValueFvPatchField<Type>
 {
-    // Private data
+    // Private Data
 
         //- Random number generator
         Random ranGen_;
@@ -137,7 +161,7 @@ public:
         );
 
         //- Construct by mapping given turbulentInletFvPatchField
-        //  onto a new patch
+        //- onto a new patch
         turbulentInletFvPatchField
         (
             const turbulentInletFvPatchField<Type>&,
@@ -181,7 +205,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         // Access
 
diff --git a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
index 9d63bfa70026a46a67f0ecda4fdbd6471c6bf904..a9eeea31b3592f78a7dbd9508958534d7ce56ed7 100644
--- a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
+++ b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
@@ -55,7 +55,7 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
     // Note: we want to add
     //      deltaU/deltaT
     // where deltaT is a local time scale:
-    //  U/(cbrt of volume)
+    //      U/(cbrt of volume)
     // Since directly manipulating the diagonal we multiply by volume.
 
     const scalarField& vol = mesh_.V();
@@ -66,10 +66,10 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
 
     forAll(U, cellI)
     {
-        scalar magU = mag(U[cellI]);
+        const scalar magU = mag(U[cellI]);
         if (magU > UMax_)
         {
-            scalar scale = sqr(Foam::cbrt(vol[cellI]));
+            const scalar scale = sqr(Foam::cbrt(vol[cellI]));
 
             diag[cellI] += scale*(magU-UMax_);
 
@@ -129,7 +129,7 @@ bool Foam::fv::velocityDampingConstraint::read(const dictionary& dict)
         if (!coeffs_.readIfPresent("UNames", fieldNames_))
         {
             fieldNames_.resize(1);
-            fieldNames_.first() = coeffs_.lookupOrDefault<word>("U", "U");
+            fieldNames_.first() = coeffs_.getOrDefault<word>("U", "U");
         }
 
         applied_.setSize(fieldNames_.size(), false);
diff --git a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H
index 83c65f29d822d59c7b467a02ed1db55701a8b2fe..cbbe294b861535f884246b545a63f06a86c6b786 100644
--- a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H
+++ b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H
@@ -31,28 +31,38 @@ Group
     grpFvOptionsConstraints
 
 Description
-    Constraint for velocity to dampen velocity fluctuations in
-    steady simulations
-
-    In case of velocity exceeding limit applies
-    a source sink to remove the excess velocity by adding a term
-    that drives the velocity to 0 (note that we cannot use the old
-    term velocity since it most likely is already excessive). This
-    additional constraint is scaled with (U-UMax)/dt
-    where dt is a local time scale (= velocity/cellsize) so it switches off
-    if the velocity is below UMax.
-
-    Constraint described by:
-
-        velocityDampingConstraintCoeffs
-        {
-            UMax        100;
-
-            // Optional: name of velocity field (default: U)
-            //U         U;
-            //UNames    (U);
-        }
-
+    Constraint on velocity field to dampen velocity fluctuations.
+
+    This constraint is primarily used to dampen velocity fluctuations in
+    the start-up phase of simulations. When the local velocity magnitude
+    exceeds the user-supplied maximum value a sink term is activated in
+    the affected region to lower the velocity to the limiting value.
+
+Usage
+    Example of functionality specification in \c fvOptions:
+    \verbatim
+    velocityDamper
+    {
+        // Mandatory entries
+        type            velocityDampingConstraint;
+        UMax            200;
+
+        // Optional entries
+        selectionMode   all;
+        UNames          (U);    // names of given velocity fields
+        U               U;      // name of given velocity field if
+                                //`UNames` is not present (default: U)
+
+        // Optional fvOptions entries
+        active          true;
+    }
+    \endverbatim
+
+Note
+    When active, this constraint manipulates the system of equations.
+    Users should ensure that it is not active when the case is converged
+    (steady-state) or during the period of interest (transient) to ensure
+    that its presence does not pollute the results.
 
 SourceFiles
     velocityDampingConstraint.C
@@ -82,9 +92,9 @@ class velocityDampingConstraint
 
 protected:
 
-    // Protected data
+    // Protected Data
 
-        //- Maximum velocity
+        //- Maximum velocity magnitude
         scalar UMax_;
 
 
diff --git a/tutorials/incompressible/simpleFoam/bump2D/0/epsilon b/tutorials/incompressible/simpleFoam/bump2D/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..8726020aa2dbc582736db94ebe3e0a1232fba889
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/bump2D/0/epsilon
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 469.872; // computed from omega and k
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    "(symUp|symDown)"
+    {
+        type            symmetryPlane;
+    }
+
+    bump
+    {
+        type            epsilonWallFunction;
+        value           $internalField;
+        lowReCorrection on;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/bump2D/0/f b/tutorials/incompressible/simpleFoam/bump2D/0/f
new file mode 100644
index 0000000000000000000000000000000000000000..60bc4e10d1a092025b13af83a47e57345994fb52
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/bump2D/0/f
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      f;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 -1 0 0 0 0 ];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    "(symUp|symDown)"
+    {
+        type            symmetryPlane;
+    }
+
+    bump
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/bump2D/0/nuTilda b/tutorials/incompressible/simpleFoam/bump2D/0/nuTilda
index 0e454300301ff96520845e7941bc75e926be6122..35a8d6c028436c1313fca3fa7ab5812315764041 100644
--- a/tutorials/incompressible/simpleFoam/bump2D/0/nuTilda
+++ b/tutorials/incompressible/simpleFoam/bump2D/0/nuTilda
@@ -45,7 +45,7 @@ boundaryField
 
     bump
     {
-        type            fixedValue; 
+        type            fixedValue;
         value           uniform 0;
     }
 
diff --git a/tutorials/incompressible/simpleFoam/bump2D/0/p b/tutorials/incompressible/simpleFoam/bump2D/0/p
index 14d7280e5b473ad6c9e220e656500ed1475a2a54..e50954ce2ed595f39dc9efda6cbe96d67ad03c09 100644
--- a/tutorials/incompressible/simpleFoam/bump2D/0/p
+++ b/tutorials/incompressible/simpleFoam/bump2D/0/p
@@ -43,7 +43,7 @@ boundaryField
 
     bump
     {
-        type            zeroGradient; 
+        type            zeroGradient;
     }
 
     frontAndBack
diff --git a/tutorials/incompressible/simpleFoam/bump2D/0/phit b/tutorials/incompressible/simpleFoam/bump2D/0/phit
new file mode 100644
index 0000000000000000000000000000000000000000..b15266813730155bcb9d685ec0dfbd7536f2c0e8
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/bump2D/0/phit
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      phit;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 0 0 0 0 ];
+
+internalField   uniform 0.66;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    "(symUp|symDown)"
+    {
+        type            symmetryPlane;
+    }
+
+    bump
+    {
+        type            fixedValue;
+        value           uniform 1e-10;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/bump2D/Allrun b/tutorials/incompressible/simpleFoam/bump2D/Allrun
index 32ad6c294786a5233a03f291a11928c27247afdf..962a15de4043ca8794c7dfe9f37f4e4069bd1a60 100755
--- a/tutorials/incompressible/simpleFoam/bump2D/Allrun
+++ b/tutorials/incompressible/simpleFoam/bump2D/Allrun
@@ -10,6 +10,7 @@ runApplication blockMesh
 models="
 kOmegaSST
 SpalartAllmaras
+kEpsilonPhitF
 "
 
 # Extract a value (Eg, from boundaryField/bump/value)
diff --git a/tutorials/incompressible/simpleFoam/bump2D/plot b/tutorials/incompressible/simpleFoam/bump2D/plot
index a400d328e33c391aae942e7d10087144fbc43bdd..52f8966b2755c78912bf97221468b397ca0a6554 100755
--- a/tutorials/incompressible/simpleFoam/bump2D/plot
+++ b/tutorials/incompressible/simpleFoam/bump2D/plot
@@ -4,7 +4,7 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 # Turbulence closure models
-models=('kOmegaSST' 'SpalartAllmaras')
+models=('kOmegaSST' 'SpalartAllmaras' 'kEpsilonPhitF')
 
 # Note: CFL3D data available from:
 # https://turbmodels.larc.nasa.gov/bump_sa.html
diff --git a/tutorials/incompressible/simpleFoam/bump2D/system/fvSchemes b/tutorials/incompressible/simpleFoam/bump2D/system/fvSchemes
index dfd865261c08e26445f290c270584a454404326f..a9cd9a7c3b3e005b7c5104bb2fb1f5638d71b7dc 100644
--- a/tutorials/incompressible/simpleFoam/bump2D/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/bump2D/system/fvSchemes
@@ -27,12 +27,14 @@ gradSchemes
 
 divSchemes
 {
-    default         none;
-    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    default          none;
+    div(phi,U)       bounded Gauss linearUpwind grad(U);
     div(phi,epsilon) bounded Gauss upwind;
-    div(phi,omega) bounded Gauss upwind;
-    div(phi,k)     bounded Gauss upwind;
+    div(phi,omega)   bounded Gauss upwind;
+    div(phi,k)       bounded Gauss upwind;
     div(phi,nuTilda) bounded Gauss upwind;
+    div(phi,phit)    bounded Gauss upwind;
+    div(phi,f)    bounded Gauss upwind;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/incompressible/simpleFoam/bump2D/system/fvSolution b/tutorials/incompressible/simpleFoam/bump2D/system/fvSolution
index 643f69a99099932aae3307bab41ed88e6f3505f6..eecd3ada6ab30e2eada26b13fbe21521e0282739 100644
--- a/tutorials/incompressible/simpleFoam/bump2D/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/bump2D/system/fvSolution
@@ -25,13 +25,21 @@ solvers
         relTol          0.1;
     }
 
-    "(U|k|epsilon|omega|nuTilda)"
+    "(U|k|epsilon|omega|nuTilda|phit)"
     {
         solver          PBiCGStab;
         preconditioner  DILU;
         tolerance       1e-8;
         relTol          0;
     }
+
+    f
+    {
+        solver          PBiCGStab;
+        preconditioner  DIC;
+        tolerance       1e-8;
+        relTol          0;
+    }
 }
 
 SIMPLE
@@ -42,7 +50,7 @@ SIMPLE
     {
         p               1e-5;
         U               1e-5;
-        "(k|epsilon|omega|nuTilda)" 1e-5;
+        "(k|epsilon|omega|nuTilda|phit|f)" 1e-5;
     }
 }
 
@@ -51,7 +59,7 @@ relaxationFactors
     equations
     {
         U               0.9;
-        "(k|epsilon|omega|nuTilda)" 0.7;
+        "(k|epsilon|omega|nuTilda|phit|f)" 0.7;
     }
 }
 
@@ -59,6 +67,8 @@ relaxationFactors
 cache
 {
     grad(U);
+    grad(k);
+    grad(phit);
 }
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
index fa0f0c2fcd8b5e5a1df62434f931b95fc48d9a79..e62a35eff37b1a34c001d7a0c0d6741600d989a2 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
@@ -12,6 +12,5 @@ zDir                 (0 0 1);
 flowDir              (1 0 0);
 z0                   uniform 0.1;
 zGround              uniform 935.0;
-value                $internalField;
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/Allclean b/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
index a1e5d09eb5072f67a786046673bc2a227ec847a0..bd93db45a94bda7c8f4636f50efa29350662900d 100755
--- a/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
@@ -5,7 +5,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 
 cleanCase0
 
-# Remove decomposeParDict
 rm -f system/decomposeParDict
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
index 0319d4b62b1be7f5523fcd3ac4394e348b776f7d..a4566f2fbc132f34ba02e5e57405708e7f136f7c 100755
--- a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
@@ -4,18 +4,19 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 # Make dummy 0 directory
+rm -rf 0
 mkdir 0
 
 runApplication blockMesh
-\cp system/decomposeParDict.hierarchical system/decomposeParDict
+cp system/decomposeParDict.hierarchical system/decomposeParDict
 runApplication decomposePar
 
-\cp system/decomposeParDict.ptscotch system/decomposeParDict
+cp system/decomposeParDict.ptscotch system/decomposeParDict
 runParallel snappyHexMesh -overwrite
 
 find . -iname '*level*' -type f -delete
 
-# - Set the initial fields
+# Set the initial fields
 restore0Dir -processor
 
 runParallel topoSet
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/constant/fvOptions b/tutorials/incompressible/simpleFoam/turbineSiting/constant/fvOptions
index 4048cd73a273f2d7d8ba13445cbd1086b113fdd5..9167d9ca5094c89b45aafc48c2c513929258aa8c 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/constant/fvOptions
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/constant/fvOptions
@@ -19,7 +19,7 @@ disk1
 {
     type            actuationDiskSource;
 
-    fields      (U);
+    fields          (U);
 
     selectionMode   cellSet;
     cellSet         actuationDisk1;
@@ -34,7 +34,7 @@ disk2
 {
     type            actuationDiskSource;
 
-    fields      (U);
+    fields          (U);
 
     selectionMode   cellSet;
     cellSet         actuationDisk2;