Implementing k-eps and k-omega RANS models

docs/references/references.bib

@@ -395,3 +395,21 @@ @misc{Adcock:2021
   howpublished = {APS DFD},
   year = {2021},
 }
+
+@article{chien1982predictions,
+  title={Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model},
+  author={Chien, Kuei-Yuan},
+  journal={AIAA journal},
+  volume={20},
+  number={1},
+  pages={33--38},
+  year={1982}
+}
+
+@book{wilcox1998turbulence,
+  title={Turbulence modeling for CFD},
+  author={Wilcox, David C and others},
+  volume={2},
+  year={1998},
+  publisher={DCW industries La Canada, CA}
+}

docs/source/theory/supportedEquationSet.rst

@@ -685,14 +685,22 @@ For simulations in which a buoyancy source term is desired, the code supports th
 
 .. _eqn_komega_sst:
 
-Shear Stress Transport (SST) RANS Model Suite
-+++++++++++++++++++++++++++++++++++++++++++++
-
-SST Formulation
-~~~~~~~~~~~~~~~
-
-Although Nalu-Wind is primarily expected to be a LES simulation tool, RANS
-modeling is supported through the activation of the SST equation set.
+RANS Model Suite
+++++++++++++++++
+
+Although Nalu-Wind is primarily expected to be a LES simulation tool,
+RANS modeling is supported through the activation of different
+two-equation RANS models: the Chien :math:`k-\epsilon` model
+:cite:`chien1982predictions`, the Wilcox 1998 :math:`k-\omega` model
+:cite:`wilcox1998turbulence`, and the SST model. For the first two
+models, the reader is referred to the reference papers and the NASA
+Turbulence Modeling Resource for the `Chien
+<https://turbmodels.larc.nasa.gov/ke-chien.html>`_ and `Wilcox
+<https://turbmodels.larc.nasa.gov/wilcox.html>`_ models. The SST model
+is explained in more details below.
+
+Shear Stress Transport (SST) Formulation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 It has been observed that standard 1998 :math:`k-\omega` models display
 a strong sensitivity to the free stream value of :math:`\omega` (see

docs/source/theory/turbulenceModeling.rst

@@ -166,14 +166,14 @@ One Equation :math:`k^{sgs}`
 
 See :math:`k^{sgs}` :ref:`PDE <theory_ksgs_les_model>` section.
 
-SST RANS Model
-++++++++++++++
+RANS Models
++++++++++++
 
-As noted, Nalu-Wind does support a SST RANS-based model (the reader is
-referred to the SST equation set description).
+As noted, Nalu-Wind supports several different RANS-based model (the
+reader is referred to the RANS model equation set description).
 
 Hybrid RANS/LES Models
-+++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++
 
 Nalu-Wind supports the Active Model Split (AMS) hybrid RANS/LES turbulence model :cite:`Haering-etal:2020`. 
 The reader is referred to the AMS equation set for more details.

include/ChienKEpsilonEquationSystem.h

@@ -0,0 +1,89 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+
+#ifndef ChienKEpsilonEquationSystem_h
+#define ChienKEpsilonEquationSystem_h
+
+#include <EquationSystem.h>
+#include <FieldTypeDef.h>
+#include <NaluParsedTypes.h>
+
+namespace stk {
+struct topology;
+namespace mesh {
+class Part;
+}
+} // namespace stk
+
+namespace sierra {
+namespace nalu {
+
+class EquationSystems;
+class AlgorithmDriver;
+class TurbKineticEnergyEquationSystem;
+class TotalDissipationRateEquationSystem;
+
+class ChienKEpsilonEquationSystem : public EquationSystem
+{
+
+public:
+  ChienKEpsilonEquationSystem(EquationSystems& equationSystems);
+  virtual ~ChienKEpsilonEquationSystem();
+
+  virtual void load(const YAML::Node&);
+
+  virtual void initialize();
+
+  virtual void register_nodal_fields(stk::mesh::Part* part);
+
+  virtual void register_wall_bc(
+    stk::mesh::Part* part,
+    const stk::topology& theTopo,
+    const WallBoundaryConditionData& wallBCData);
+
+  virtual void register_interior_algorithm(stk::mesh::Part* part);
+
+  virtual void solve_and_update();
+
+  void initial_work();
+  virtual void post_external_data_transfer_work();
+  virtual void post_iter_work();
+
+  void clip_min_distance_to_wall();
+  void compute_dplus_function();
+  void update_and_clip();
+  void clip_ke(
+    const stk::mesh::NgpMesh& ngpMesh,
+    const stk::mesh::Selector& sel,
+    stk::mesh::NgpField<double>& tke,
+    stk::mesh::NgpField<double>& tdr);
+
+  TurbKineticEnergyEquationSystem* tkeEqSys_;
+  TotalDissipationRateEquationSystem* tdrEqSys_;
+
+  ScalarFieldType* tke_;
+  ScalarFieldType* tdr_;
+  ScalarFieldType* minDistanceToWall_;
+  ScalarFieldType* dplus_;
+
+  bool isInit_;
+
+  // saved of mesh parts that are for wall bcs
+  std::vector<stk::mesh::Part*> wallBcPart_;
+
+  bool resetAMSAverages_;
+
+  const double tkeMinValue_{1.0e-8};
+  const double tdrMinValue_{1.0e-8};
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

include/Enums.h

@@ -74,15 +74,16 @@ enum EquationType {
   EQ_ENTHALPY = 6,
   EQ_MESH_DISPLACEMENT = 7,
   EQ_SPEC_DISS_RATE = 8,
-  EQ_MASS_FRACTION = 9,
-  EQ_PNG   = 10,
-  EQ_PNG_P = 11,
-  EQ_PNG_Z = 12,
-  EQ_PNG_H = 13,
-  EQ_PNG_U = 14,
-  EQ_PNG_TKE = 15, // FIXME... Last PNG managed like this..
-  EQ_WALL_DISTANCE = 16,
-  EQ_GAMMA_TRANS = 17,
+  EQ_TOT_DISS_RATE = 9,
+  EQ_MASS_FRACTION = 10,
+  EQ_PNG = 11,
+  EQ_PNG_P = 12,
+  EQ_PNG_Z = 13,
+  EQ_PNG_H = 14,
+  EQ_PNG_U = 15,
+  EQ_PNG_TKE = 16, // FIXME... Last PNG managed like this..
+  EQ_WALL_DISTANCE = 17,
+  EQ_GAMMA_TRANS = 18,
   EquationSystemType_END
 };
 
@@ -96,15 +97,15 @@ static const std::string EquationTypeMap[] = {
   "Enthalpy",
   "MeshVelocity",
   "Specific_Dissipation_Rate",
+  "Total_Dissipation_Rate",
   "Mass_Fraction",
   "PNG",
   "PNG_P",
   "PNG_Z",
   "PNG_H",
   "PNG_U",
   "PNG_TKE",
-  "Wall_Distance"
-};
+  "Wall_Distance"};
 
 enum UserDataType {
   CONSTANT_UD = 0,
@@ -166,19 +167,15 @@ enum class TurbulenceModel {
   SST_DES,
   SST_AMS,
   SST_IDDES,
+  KE,
+  KO,
   TurbulenceModel_END
 };
 
 // matching string name index into above enums (must match PERFECTLY)
 static const std::string TurbulenceModelNames[] = {
-  "laminar",
-  "ksgs",
-  "smagorinsky",
-  "wale",
-  "sst",
-  "sst_des",
-  "sst_ams",
-  "sst_iddes"};
+  "laminar", "ksgs",      "smagorinsky", "wale", "sst", "sst_des",
+  "sst_ams", "sst_iddes",   "ke",  "ko"};
 
 enum TurbulenceModelConstant {
   TM_cMu = 0,
@@ -234,15 +231,23 @@ enum TurbulenceModelConstant {
   TM_ams_peclet_offset = 50,
   TM_ams_peclet_slope = 51,
   TM_ams_peclet_scale = 52,
-  TM_tkeAmb = 53,
-  TM_sdrAmb = 54,
-  TM_caOne = 55,
-  TM_caTwo = 56,
-  TM_ceOne = 57,
-  TM_ceTwo = 58,
-  TM_c0t = 59,
-  TM_avgTimeCoeff = 60,
-  TM_END = 61
+  TM_fMuExp = 53,
+  TM_utau = 54,
+  TM_cEpsOne = 55,
+  TM_cEpsTwo = 56,
+  TM_fOne = 57,
+  TM_sigmaK = 58,
+  TM_sigmaEps = 59,
+  TM_tkeAmb = 60,
+  TM_sdrAmb = 61,
+  TM_avgTimeCoeff = 62,
+  TM_alphaInf = 63,
+  TM_caOne = 64,
+  TM_caTwo = 65,
+  TM_ceOne = 66,
+  TM_ceTwo = 67,
+  TM_c0t = 68,
+  TM_END = 69
 };
 
 static const std::string TurbulenceModelConstantNames[] = {
@@ -299,9 +304,22 @@ static const std::string TurbulenceModelConstantNames[] = {
   "ams_peclet_offset",
   "ams_peclet_slope",
   "ams_peclet_scale",
+  "fMuExp",
+  "utau",
+  "cEpsOne",
+  "cEpsTwo",
+  "fOne",
+  "sigmaK",
+  "sigmaEps",
   "tke_amb",
   "sdr_amb",
   "avgTimeCoeff",
+  "alphaInf",
+  "caOne",
+  "caTwo",
+  "ceOne",
+  "ceTwo",
+  "c0t",
   "END"};
 
 enum ActuatorType {

include/NaluParsedTypes.h

@@ -54,6 +54,12 @@ struct SpecDissRate {
   {}
 };
 
+struct TotDissRate
+{
+  double totDissRate_;
+  TotDissRate() : totDissRate_(0.0) {}
+};
+
 struct GammaInf {
   double gamma_;
   GammaInf()

include/NaluParsing.h

@@ -53,6 +53,7 @@ struct WallUserData : public UserData {
   Velocity u_;
   Velocity dx_;
   TurbKinEnergy tke_;
+  TotDissRate tdr_;
   MixtureFraction mixFrac_;
   MassFraction massFraction_;
   NormalHeatFlux q_;
@@ -96,19 +97,27 @@ struct InflowUserData : public UserData {
   Velocity u_;
   TurbKinEnergy tke_;
   SpecDissRate sdr_;
+  TotDissRate tdr_;
   MixtureFraction mixFrac_;
   MassFraction massFraction_;
   GammaInf gamma_;
 
   bool uSpec_;
   bool tkeSpec_;
   bool sdrSpec_;
+  bool tdrSpec_;
   bool mixFracSpec_;
   bool massFractionSpec_;
   bool gammaSpec_;
   InflowUserData()
     : UserData(),
-    uSpec_(false), tkeSpec_(false), sdrSpec_(false), mixFracSpec_(false), massFractionSpec_(false), gammaSpec_(false)
+      uSpec_(false),
+      tkeSpec_(false),
+      sdrSpec_(false),
+      tdrSpec_(false),
+      mixFracSpec_(false),
+      massFractionSpec_(false),
+      gammaSpec_(false)
   {}
 };
 
@@ -117,6 +126,7 @@ struct OpenUserData : public UserData {
   Pressure p_;
   TurbKinEnergy tke_;
   SpecDissRate sdr_;
+  TotDissRate tdr_;
   MixtureFraction mixFrac_;
   MassFraction massFraction_;
   GammaOpen gamma_;
@@ -125,6 +135,7 @@ struct OpenUserData : public UserData {
   bool pSpec_;
   bool tkeSpec_;
   bool sdrSpec_;
+  bool tdrSpec_;
   bool mixFracSpec_;
   bool massFractionSpec_;
   bool totalP_;
@@ -133,9 +144,16 @@ struct OpenUserData : public UserData {
 
   OpenUserData()
     : UserData(),
-      uSpec_(false), pSpec_(false), tkeSpec_(false), 
-      sdrSpec_(false), mixFracSpec_(false), massFractionSpec_(false), 
-      totalP_{false}, gammaSpec_(false), entrainMethod_{EntrainmentMethod::COMPUTED}
+      uSpec_(false),
+      pSpec_(false),
+      tkeSpec_(false),
+      sdrSpec_(false),
+      tdrSpec_(false),
+      mixFracSpec_(false),
+      massFractionSpec_(false),
+      totalP_{false},
+      gammaSpec_(false),
+      entrainMethod_{EntrainmentMethod::COMPUTED}
   {}
 };
 
@@ -450,6 +468,12 @@ template<> struct convert<sierra::nalu::SpecDissRate> {
   static bool decode(const Node& node, sierra::nalu::SpecDissRate& rhs) ;
 };
 
+template <>
+struct convert<sierra::nalu::TotDissRate>
+{
+  static bool decode(const Node& node, sierra::nalu::TotDissRate& rhs);
+};
+
 template<> struct convert<sierra::nalu::GammaInf> {
   static bool decode(const Node& node, sierra::nalu::GammaInf& rhs) ;
 };

include/ShearStressTransportEquationSystem.h

@@ -73,7 +73,7 @@ class ShearStressTransportEquationSystem : public EquationSystem
 
   TurbKineticEnergyEquationSystem* tkeEqSys_;
   SpecificDissipationRateEquationSystem* sdrEqSys_;
-  GammaEquationSystem *gammaEqSys_;
+  GammaEquationSystem* gammaEqSys_;
 
   ScalarFieldType* tke_;
   ScalarFieldType* sdr_;