ENH-DOC: fix inconsistent Cmu lookup (fixes #647)

src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C

@@ -48,8 +48,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(p, iF),
+    kName_("k"),
     mixingLength_(0.0),
-    kName_("k")
+    Cmu_(0.0)
 {
     this->refValue() = 0.0;
     this->refGrad() = 0.0;
@@ -67,8 +68,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf, p, iF, mapper),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -81,8 +83,12 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(p, iF),
-    mixingLength_(dict.get<scalar>("mixingLength")),
-    kName_(dict.getOrDefault<word>("k", "k"))
+    kName_(dict.getOrDefault<word>("k", "k")),
+    mixingLength_
+    (
+        dict.getCheck<scalar>("mixingLength", scalarMinMax::ge(SMALL))
+    ),
+    Cmu_(dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL)))
 {
     this->phiName_ = dict.getOrDefault<word>("phi", "phi");
 
@@ -101,8 +107,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -114,8 +121,9 @@ turbulentMixingLengthDissipationRateInletFvPatchScalarField
 )
 :
     inletOutletFvPatchScalarField(ptf, iF),
+    kName_(ptf.kName_),
     mixingLength_(ptf.mixingLength_),
-    kName_(ptf.kName_)
+    Cmu_(ptf.Cmu_)
 {}
 
 
@@ -138,10 +146,8 @@ void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
         )
     );
 
-    const scalar Cmu =
-        turbModel.coeffDict().getOrDefault<scalar>("Cmu", 0.09);
-
-    const scalar Cmu75 = pow(Cmu, 0.75);
+    Cmu_ = (readIfPresent("Cmu", turbModel.coeffDict())).value();
+    const scalar Cmu75 = pow(Cmu_, 0.75);
 
     const fvPatchScalarField& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
@@ -149,7 +155,7 @@ void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
     const fvsPatchScalarField& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
 
-    this->refValue() = Cmu75*kp*sqrt(kp)/mixingLength_;
+    this->refValue() = (Cmu75/mixingLength_)*pow(kp, 1.5);
     this->valueFraction() = 1.0 - pos0(phip);
 
     inletOutletFvPatchScalarField::updateCoeffs();

src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,43 +31,63 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This boundary condition provides a turbulence dissipation, \f$\epsilon\f$
-    (epsilon) inlet condition based on a specified mixing length.  The patch
-    values are calculated using:
+    This boundary condition provides an inlet condition for turbulent kinetic
+    energy dissipation rate, i.e. \c epsilon, based on a specified mixing
+    length.  The patch values are calculated using:
 
         \f[
             \epsilon_p = \frac{C_{\mu}^{0.75} k^{1.5}}{L}
         \f]
 
     where
-
     \vartable
-        \epsilon_p | patch epsilon values
-        C_{\mu} | Model coefficient, set to 0.09
-        k       | turbulence kinetic energy
-        L       | length scale
+      \epsilon_p | Patch epsilon values     [m2/s3]
+      C_\mu      | Empirical model constant retrived from turbulence model
+      k          | Turbulent kinetic energy [m2/s2]
+      L          | Mixing length scale      [m]
     \endvartable
 
 Usage
-    \table
-        Property     | Description             | Required    | Default value
-        mixingLength | Length scale [m]        | yes         |
-        phi          | flux field name         | no          | phi
-        k            | turbulence kinetic energy field name | no | k
-    \endtable
-
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
+        // Mandatory entries (unmodifiable)
         type            turbulentMixingLengthDissipationRateInlet;
+
+        // Mandatory entries (runtime modifiable)
         mixingLength    0.005;
-        value           uniform 200;   // placeholder
+
+        // Optional entries (runtime modifiable)
+        Cmu             0.09;
+        k               k;
+        phi             phi;
+
+        // Placeholder
+        value           uniform 200;
     }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property     | Description                  | Type   | Req'd | Dflt
+      mixingLength | Mixing length scale [m]      | scalar |  yes  | -
+      Cmu          | Empirical model constant     | scalar |  no   | 0.09
+      phi          | Name of flux field           | word   |  no   | phi
+      k            | Name of turbulent kinetic energy field | word | no | k
+    \endtable
+
 Note
-    In the event of reverse flow, a zero-gradient condition is applied
+    - The boundary condition is derived from \c inletOutlet condition.
+      Therefore, in the event of reverse flow, a zero-gradient condition
+      is applied.
+    - The order of precedence to input the empirical model constant \c Cmu is:
+      turbulence model, boundary condition dictionary, and default value=0.09.
+    - The empirical model constant \c Cmu is not a spatiotemporal variant field.
+      Therefore, the use of the boundary condition may not be fully consistent
+      with the turbulence models where \c Cmu is a variant field, such as
+      \c realizableKE closure model in this respect. Nevertheless, workflow
+      observations suggest that the matter poses no importance.
 
 See also
     Foam::inletOutletFvPatchField
@@ -94,14 +115,17 @@ class turbulentMixingLengthDissipationRateInletFvPatchScalarField
 :
     public inletOutletFvPatchScalarField
 {
-    // Private data
-
-        //- Turbulent length scale
-        scalar mixingLength_;
+    // Private Data
 
         //- Name of the turbulent kinetic energy field
         word kName_;
 
+        //- Mixing length scale
+        scalar mixingLength_;
+
+        //- Empirical model constant
+        scalar Cmu_;
+
 
 public:
 
@@ -127,8 +151,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  turbulentMixingLengthDissipationRateInletFvPatchScalarField
-        //  onto a new patch
+        //- turbulentMixingLengthDissipationRateInletFvPatchScalarField
+        //- onto a new patch
         turbulentMixingLengthDissipationRateInletFvPatchScalarField
         (
             const turbulentMixingLengthDissipationRateInletFvPatchScalarField&,
@@ -179,7 +203,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();