Fixed values for turbulence quantities in upstream half-plane

Common/include/CConfig.hpp

@@ -117,6 +117,7 @@ class CConfig {
   InvDesign_Cp,             /*!< \brief Flag to know if the code is going to compute and plot the inverse design. */
   InvDesign_HeatFlux,       /*!< \brief Flag to know if the code is going to compute and plot the inverse design. */
   Wind_Gust,                /*!< \brief Flag to know if there is a wind gust. */
+  Turb_Fixed_Values,        /*!< \brief Flag to know if there are fixed values for turbulence quantities in one half-plane. */
   Aeroelastic_Simulation,   /*!< \brief Flag to know if there is an aeroelastic simulation. */
   Weakly_Coupled_Heat,      /*!< \brief Flag to know if a heat equation should be weakly coupled to the incompressible solver. */
   Rotating_Frame,           /*!< \brief Flag to know if there is a rotating frame. */
@@ -921,6 +922,8 @@ class CConfig {
   Gust_Ampl,                  /*!< \brief Gust amplitude. */
   Gust_Begin_Time,            /*!< \brief Time at which to begin the gust. */
   Gust_Begin_Loc;             /*!< \brief Location at which the gust begins. */
+  /*! \brief Maximal scalar product of the normed far-field velocity vector and a space coordinate where fixed turbulence quantities are set. */
+  su2double Turb_Fixed_Values_MaxScalarProd;
   long Visualize_CV;          /*!< \brief Node number for the CV to be visualized */
   bool ExtraOutput;           /*!< \brief Check if extra output need. */
   bool Wall_Functions;           /*!< \brief Use wall functions with the turbulence model */
@@ -8051,6 +8054,20 @@ class CConfig {
    */
   su2double GetGust_Begin_Loc(void) const { return Gust_Begin_Loc; }
 
+  /*!
+   * \brief Get whether fixed values for turbulence quantities are applied.
+   * \return <code>TRUE</code> if fixed values are applied; otherwise <code>FALSE</code>.
+   */
+  bool GetTurb_Fixed_Values(void) const { return Turb_Fixed_Values; }
+
+  /*!
+   * \brief Get shift of the upstream half-plane where fixed values for turbulence quantities are applied.
+   * \details This half-plane is given by the condition that the dot product between the
+   * coordinate vector and the normalized far-field velocity vector is less than what this
+   * function returns.
+   */
+  su2double GetTurb_Fixed_Values_MaxScalarProd(void) const { return Turb_Fixed_Values_MaxScalarProd; }
+
   /*!
    * \brief Get the number of iterations to evaluate the parametric coordinates.
    * \return Number of iterations to evaluate the parametric coordinates.

Common/src/CConfig.cpp

@@ -2105,6 +2105,12 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Direction of the gust X or Y dir */
   addEnumOption("GUST_DIR", Gust_Dir, Gust_Dir_Map, Y_DIR);
 
+  /* Fixed values for turbulence quantities to keep them at inflow conditions. */
+  /* DESCRIPTION: Fix turbulence quantities to far-field values inside an upstream half-space. */
+  addBoolOption("TURB_FIXED_VALUES", Turb_Fixed_Values, false);
+  /* DESCRIPTION: Shift of the fixed values half-space, in space units in the direction of far-field velocity. */
+  addDoubleOption("TURB_FIXED_VALUES_DOMAIN", Turb_Fixed_Values_MaxScalarProd, numeric_limits<su2double>::lowest());
+
   /* Harmonic Balance config */
   /* DESCRIPTION: Omega_HB = 2*PI*frequency - frequencies for Harmonic Balance method */
   addDoubleListOption("OMEGA_HB", nOmega_HB, Omega_HB);
@@ -4364,6 +4370,10 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
     SU2_MPI::Error("The LM transition model is under maintenance.", CURRENT_FUNCTION);
   }
 
+  if(Turb_Fixed_Values && Turb_Fixed_Values_MaxScalarProd==numeric_limits<su2double>::lowest()){
+    SU2_MPI::Error("TURB_FIXED_VALUES activated, but no domain set with TURB_FIXED_VALUES_DOMAIN.", CURRENT_FUNCTION);
+  }
+
   /*--- Check for constant lift mode. Initialize the update flag for
    the AoA with each iteration to false  ---*/
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -1363,6 +1363,14 @@ class CSolver {
                                               CNumerics *visc_numerics,
                                               CConfig *config,
                                               unsigned short val_marker) { }
+  /*!
+   * \brief Virtual function to apply something like a strong BC to the whole domain.
+   * \details Overridden in CTurbSolver to impose fixed values to turbulence quantities
+   * in a specified upstream half-plane.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void Impose_Fixed_Values(const CGeometry *geometry, const CConfig *config) { }
 
  /*!
    * \brief Get the outer state for fluid interface nodes.

SU2_CFD/include/solvers/CTurbSASolver.hpp

@@ -39,7 +39,7 @@
 
 class CTurbSASolver final : public CTurbSolver {
 private:
-  su2double nu_tilde_Inf, nu_tilde_Engine, nu_tilde_ActDisk;
+  su2double nu_tilde_Engine, nu_tilde_ActDisk;
 
   /*!
    * \brief A virtual member.
@@ -365,16 +365,6 @@ class CTurbSASolver final : public CTurbSolver {
                   unsigned short val_marker,
                   bool val_inlet_surface) override;
 
-  /*!
-   * \brief Set the solution using the Freestream values.
-   * \param[in] config - Definition of the particular problem.
-   */
-  inline void SetFreeStream_Solution(const CConfig *config) override {
-    SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++)
-      nodes->SetSolution(iPoint, 0, nu_tilde_Inf);
-  }
-
   /*!
    * \brief Store of a set of provided inlet profile values at a vertex.
    * \param[in] val_inlet - vector containing the inlet values for the current vertex.
@@ -416,6 +406,6 @@ class CTurbSASolver final : public CTurbSolver {
    * \brief Get the value of nu tilde at the far-field.
    * \return Value of nu tilde at the far-field.
    */
-  inline su2double GetNuTilde_Inf(void) const override { return nu_tilde_Inf; }
+  inline su2double GetNuTilde_Inf(void) const override { return Solution_Inf[0]; }
 
 };

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -38,9 +38,7 @@
 class CTurbSSTSolver final : public CTurbSolver {
 private:
   su2double
-  constants[10] = {0.0}, /*!< \brief Constants for the model. */
-  kine_Inf,              /*!< \brief Free-stream turbulent kinetic energy. */
-  omega_Inf;             /*!< \brief Free-stream specific dissipation. */
+  constants[10] = {0.0}; /*!< \brief Constants for the model. */
 
 public:
   /*!
@@ -235,18 +233,6 @@ class CTurbSSTSolver final : public CTurbSolver {
    */
   inline const su2double* GetConstants() const override { return constants; }
 
-  /*!
-   * \brief Set the solution using the Freestream values.
-   * \param[in] config - Definition of the particular problem.
-   */
-  inline void SetFreeStream_Solution(const CConfig *config) override {
-    SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
-      nodes->SetSolution(iPoint, 0, kine_Inf);
-      nodes->SetSolution(iPoint, 1, omega_Inf);
-    }
-  }
-
   /*!
    * \brief Store of a set of provided inlet profile values at a vertex.
    * \param[in] val_inlet - vector containing the inlet values for the current vertex.
@@ -287,12 +273,12 @@ class CTurbSSTSolver final : public CTurbSolver {
    * \brief Get the value of the turbulent kinetic energy.
    * \return Value of the turbulent kinetic energy.
    */
-  inline su2double GetTke_Inf(void) const override { return kine_Inf; }
+  inline su2double GetTke_Inf(void) const override { return Solution_Inf[0]; }
 
   /*!
    * \brief Get the value of the turbulent frequency.
    * \return Value of the turbulent frequency.
    */
-  inline su2double GetOmega_Inf(void) const override { return omega_Inf; }
+  inline su2double GetOmega_Inf(void) const override { return Solution_Inf[1]; }
 
 };

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -55,6 +55,8 @@ class CTurbSolver : public CSolver {
   Gamma_Minus_One;              /*!< \brief Fluids's Gamma - 1.0  . */
   vector<su2activematrix> Inlet_TurbVars;  /*!< \brief Turbulence variables at inlet profiles */
 
+  su2double Solution_Inf[MAXNVAR] = {0.0}; /*!< \brief Far-field solution. */
+
   /*--- Sliding meshes variables. ---*/
 
   vector<su2matrix<su2double*> > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone)
@@ -249,6 +251,24 @@ class CTurbSolver : public CSolver {
                           CNumerics *visc_numerics,
                           CConfig *config) final;
 
+  /*!
+   * \brief Set the solution using the Freestream values.
+   * \param[in] config - Definition of the particular problem.
+   */
+  inline void SetFreeStream_Solution(const CConfig *config) final {
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+      nodes->SetSolution(iPoint, Solution_Inf);
+    }
+  }
+
+  /*!
+   * \brief Impose fixed values to turbulence quantities.
+   * \details Turbulence quantities are set to far-field values in an upstream half-plane
+   * in order to keep them from decaying.
+   */
+  void Impose_Fixed_Values(const CGeometry *geometry, const CConfig *config) final;
+
   /*!
    * \brief Prepare an implicit iteration.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/src/integration/CIntegration.cpp

@@ -146,6 +146,9 @@ void CIntegration::Space_Integration(CGeometry *geometry,
     }
   }
 
+  /*--- Modification of the system on the whole domain, like for a strong BC. ---*/
+  solver_container[MainSolver]->Impose_Fixed_Values(geometry, config);
+
   /*--- Strong boundary conditions (Navier-Stokes and Dirichlet type BCs) ---*/
 
   for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++)

SU2_CFD/src/solvers/CTurbSASolver.cpp

@@ -113,11 +113,13 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
   /*--- Factor_nu_Inf in [3.0, 5.0] ---*/
 
   Factor_nu_Inf = config->GetNuFactor_FreeStream();
-  nu_tilde_Inf  = Factor_nu_Inf*Viscosity_Inf/Density_Inf;
+  su2double nu_tilde_Inf  = Factor_nu_Inf*Viscosity_Inf/Density_Inf;
   if (config->GetKind_Trans_Model() == BC) {
     nu_tilde_Inf  = 0.005*Factor_nu_Inf*Viscosity_Inf/Density_Inf;
   }
 
+  Solution_Inf[0] = nu_tilde_Inf;
+
   /*--- Factor_nu_Engine ---*/
   Factor_nu_Engine = config->GetNuFactor_Engine();
   nu_tilde_Engine  = Factor_nu_Engine*Viscosity_Inf/Density_Inf;
@@ -513,7 +515,7 @@ void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container
 
       /*--- Set turbulent variable at the wall, and at infinity ---*/
 
-      conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde_Inf);
+      conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), Solution_Inf);
 
       /*--- Set Normal (it is necessary to change the sign) ---*/
 
@@ -1957,7 +1959,7 @@ su2double CTurbSASolver::GetInletAtVertex(su2double *val_inlet,
 void CTurbSASolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) {
 
   for(unsigned long iVertex=0; iVertex < nVertex[iMarker]; iVertex++){
-    Inlet_TurbVars[iMarker][iVertex][0] = nu_tilde_Inf;
+    Inlet_TurbVars[iMarker][iVertex][0] = GetNuTilde_Inf();
   }
 
 }

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -130,8 +130,11 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh
 
   su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, VelInf);
 
-  kine_Inf  = 3.0/2.0*(VelMag2*Intensity*Intensity);
-  omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio);
+  su2double kine_Inf  = 3.0/2.0*(VelMag2*Intensity*Intensity);
+  su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio);
+
+  Solution_Inf[0] = kine_Inf;
+  Solution_Inf[1] = omega_Inf;
 
   /*--- Eddy viscosity, initialized without stress limiter at the infinity ---*/
   muT_Inf = rhoInf*kine_Inf/omega_Inf;
@@ -473,9 +476,7 @@ void CTurbSSTSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_containe
 
       /*--- Set turbulent variable at the wall, and at infinity ---*/
 
-      su2double solution_j[] = {kine_Inf, omega_Inf};
-
-      conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j);
+      conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), Solution_Inf);
 
       /*--- Set Normal (it is necessary to change the sign) ---*/
 
@@ -959,8 +960,8 @@ su2double CTurbSSTSolver::GetInletAtVertex(su2double *val_inlet,
 void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) {
 
   for(unsigned long iVertex=0; iVertex < nVertex[iMarker]; iVertex++){
-    Inlet_TurbVars[iMarker][iVertex][0] = kine_Inf;
-    Inlet_TurbVars[iMarker][iVertex][1] = omega_Inf;
+    Inlet_TurbVars[iMarker][iVertex][0] = GetTke_Inf();
+    Inlet_TurbVars[iMarker][iVertex][1] = GetOmega_Inf();
   }
 
 }

SU2_CFD/src/solvers/CTurbSolver.cpp

@@ -80,6 +80,7 @@ CTurbSolver::~CTurbSolver(void) {
   }
 
   delete nodes;
+
 }
 
 void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container,
@@ -508,6 +509,41 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta
 
 }
 
+void CTurbSolver::Impose_Fixed_Values(const CGeometry *geometry, const CConfig *config){
+
+  /*--- Check whether turbulence quantities are fixed to far-field values on a half-plane. ---*/
+  if(config->GetTurb_Fixed_Values()){
+
+    const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+
+    /*--- Form normalized far-field velocity ---*/
+    const su2double* velocity_inf = config->GetVelocity_FreeStreamND();
+    su2double velmag_inf = GeometryToolbox::Norm(nDim, velocity_inf);
+    if(velmag_inf==0)
+      SU2_MPI::Error("Far-field velocity is zero, cannot fix turbulence quantities to inflow values.", CURRENT_FUNCTION);
+    su2double unit_velocity_inf[MAXNDIM];
+    for(unsigned short iDim=0; iDim<nDim; iDim++)
+      unit_velocity_inf[iDim] = velocity_inf[iDim] / velmag_inf;
+
+    SU2_OMP_FOR_DYN(omp_chunk_size)
+    for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+      if( GeometryToolbox::DotProduct(nDim, geometry->nodes->GetCoord(iPoint), unit_velocity_inf)
+        < config->GetTurb_Fixed_Values_MaxScalarProd() ) {
+        /*--- Set the solution values and zero the residual ---*/
+        nodes->SetSolution_Old(iPoint, Solution_Inf);
+        nodes->SetSolution(iPoint, Solution_Inf);
+        LinSysRes.SetBlock_Zero(iPoint);
+        if (implicit) {
+          /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/
+          for(unsigned long iVar=0; iVar<nVar; iVar++)
+            Jacobian.DeleteValsRowi(iPoint*nVar+iVar);
+        }
+      }
+    }
+  }
+
+}
+
 void CTurbSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver** solver_container, CConfig *config) {
 
   const auto flowNodes = solver_container[FLOW_SOL]->GetNodes();