Improved combustion solver robustness for coarse grids

Common/include/option_structure.hpp

@@ -1446,6 +1446,7 @@ struct FluidFlamelet_ParsedOptions {
   su2double* spark_reaction_rates; /*!< \brief Source terms for flamelet spark ignition option. */
   unsigned short nspark;           /*!< \brief Number of source terms for spark initialization. */
   bool preferential_diffusion = false;  /*!< \brief Preferential diffusion physics for flamelet solver.*/
+  bool thickenedflame_correction{true}; /*!< \brief Thickened flame correction. */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1429,6 +1429,9 @@ void CConfig::SetConfig_Options() {
   /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/
   addDoubleListOption("SPARK_REACTION_RATES", flamelet_ParsedOptions.nspark, flamelet_ParsedOptions.spark_reaction_rates);
 
+  /*!\brief THICKENED_FLAME_CORRECTION \n DESCRIPTION: Dampen source terms and enhance diffusion based on the flame length scale. \ingroup Config*/
+  addBoolOption("THICKENED_FLAME_CORRECTION", flamelet_ParsedOptions.thickenedflame_correction, true);
+
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
   /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/

SU2_CFD/include/solvers/CSolver.hpp

@@ -580,6 +580,12 @@ class CSolver {
    */
   inline virtual void SetPrimitive_Limiter(CGeometry *geometry, const CConfig *config) { }
 
+  /*!
+   * \brief A virtual member.
+   * \return flame thickness value.
+   */
+  virtual su2double GetFlameThickness() const {return 1.0;}
+
   /*!
    * \brief Compute the projection of a variable for MUSCL reconstruction.
    * \note The result should be halved when added to i (or subtracted from j).

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -37,6 +37,9 @@
  */
 class CSpeciesFlameletSolver final : public CSpeciesSolver {
  private:
+  const su2double default_flame_thickness{1.0};
+  su2double global_flame_thickness;
+  bool calc_flame_thickness{false};
   FluidFlamelet_ParsedOptions flamelet_config_options;
   bool include_mixture_fraction = false; /*!< \brief include mixture fraction as a controlling variable. */
   /*!
@@ -81,10 +84,11 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    * \param[in] iPoint - node ID.
    * \param[in] scalars - local scalar solution.
    * \param[in] table_source_names - variable names of scalar source terms.
+   * \param[in] F - flame thickness correction factor.
    * \return - within manifold bounds (0) or outside manifold bounds (1).
    */
   unsigned long SetScalarSources(const CConfig* config, CFluidModel* fluid_model_local, unsigned long iPoint,
-                                 const vector<su2double>& scalars);
+                                 const vector<su2double>& scalars, const su2double F=1.0);
 
   /*!
    * \brief Retrieve passive look-up data from manifold.
@@ -105,6 +109,22 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    */
   unsigned long SetPreferentialDiffusionScalars(CFluidModel* fluid_model_local,
                                                 unsigned long iPoint, const vector<su2double>& scalars);
+
+  /*!
+   * \brief Calculate correction factor for flame propagation on coarse grids.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] iPoint - node ID.
+   * \return - flame thickness correction factor.
+   */                                              
+  su2double ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const;   
+
+  /*!
+   * \brief Approximate the minimum flame thickness value used for the thickened flame model.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \return - approximate flame thickness value.
+   */
+  su2double GetOverallFlameThickness(CGeometry * geometry,  CSolver **solver_container) const;
 
  public:
   /*!
@@ -195,4 +215,10 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
    */
   void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics,
                         const CConfig* config) final;
+
+  /*!
+   * \brief Obtain the overall flame thickness value.
+   * \return flame thickness value.
+   */
+  virtual su2double GetFlameThickness() const {return global_flame_thickness;}
 };

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1069,6 +1069,8 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) {
         const auto& CV_name = flamelet_config_options.controlling_variable_names[iCV];
         AddHistoryOutput("RMS_"+CV_name, "rms["+CV_name+"]",ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean squared residual of " + CV_name + " controlling variable equation.", HistoryFieldType::RESIDUAL);
       }
+      if (flamelet_config_options.thickenedflame_correction)
+        AddHistoryOutput("THICKNESS","flamethickness",ScreenOutputFormat::FIXED, "RMS_RES", "Flame thickness used for thickened flame correction model.", HistoryFieldType::RESIDUAL);
 
       /*--- auxiliary species transport ---*/
       for (auto i_scalar = 0u; i_scalar < flamelet_config_options.n_user_scalars; i_scalar++){
@@ -1304,6 +1306,9 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co
         }
       }
 
+      if (flamelet_config_options.thickenedflame_correction)
+        SetHistoryOutputValue("THICKNESS", solver[SPECIES_SOL]->GetFlameThickness());
+
       SetHistoryOutputValue("LINSOL_ITER_FLAMELET", solver[SPECIES_SOL]->GetIterLinSolver());
       SetHistoryOutputValue("LINSOL_RESIDUAL_FLAMELET", log10(solver[SPECIES_SOL]->GetResLinSolver()));
     }

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -39,6 +39,8 @@
 
   /*--- Retrieve options from config. ---*/
   flamelet_config_options = config->GetFlameletParsedOptions();
+  global_flame_thickness = default_flame_thickness;
+  calc_flame_thickness = flamelet_config_options.thickenedflame_correction;
 
   /*--- Dimension of the problem. ---*/
   nVar = flamelet_config_options.n_scalars;
@@ -65,45 +67,79 @@
 
   /*--- Add the solver name. ---*/
   SolverName = "FLAMELET";
+
+  if (calc_flame_thickness && rank==MASTER_NODE) {
+    cout << "Applying thickened flame source and diffusion correction." << endl;
+  }
 }
 
 void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config,
                                            unsigned short iMesh, unsigned short iRKStep,
                                            unsigned short RunTime_EqSystem, bool Output) {
   unsigned long n_not_in_domain_local = 0, n_not_in_domain_global = 0;
   vector<su2double> scalars_vector(nVar);
+
   unsigned long spark_iter_start, spark_duration;
   bool ignition = false;
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
   /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/
+  unsigned long iter;
+  if (config->GetMultizone_Problem()) {
+    iter = config->GetOuterIter();
+  } else if (config->GetTime_Domain()) {
+    iter = config->GetTimeIter();
+  } else {
+    iter = config->GetInnerIter();
+  }
   if ((flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK)) {
     auto spark_init = flamelet_config_options.spark_init;
     spark_iter_start = ceil(spark_init[4]);
     spark_duration = ceil(spark_init[5]);
-    unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter();
+
     ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
-
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
+  /* Update global flame thickness value. */
+  if (calc_flame_thickness) {
+    su2double test_thickness = GetOverallFlameThickness(geometry, solver_container);
+    if (test_thickness < global_flame_thickness) {
+      global_flame_thickness = test_thickness;
+    } else {
+      global_flame_thickness = 0.95*global_flame_thickness + 0.05*test_thickness;
+    }
+  }
+
+  /* Flame thickness correction factors */
+  su2double F{1.0}, F_source{1.0};
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPoint; i_point++) {
     CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel();
     su2double* scalars = nodes->GetSolution(i_point);
+
+    /*--- Calculate correction factor for flame propagation on coarse grids. ---*/
+    if (calc_flame_thickness) {
+      F = ThickenedFlameCorrection(geometry, i_point);
+      F_source = 1.0 / F;
+    }
+
     for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar];
 
-    /*--- Compute total source terms from the production and consumption. ---*/
-    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector);
+    /*--- Only apply thickened flame correction factor to sources for steady problems. ---*/
+    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source);
 
     if (ignition) {
       /*--- Apply source terms within spark radius. ---*/
       su2double dist_from_center = 0,
                 spark_radius = flamelet_config_options.spark_init[3];
       dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data());
       if (dist_from_center < pow(spark_radius,2)) {
-        for (auto iVar = 0u; iVar < nVar; iVar++)
-          nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar]);
+        for (auto iVar = 0u; iVar < nVar; iVar++) {
+          su2double source_total = nodes->GetScalarSources(i_point)[iVar] + F_source * flamelet_config_options.spark_reaction_rates[iVar];
+          nodes->SetScalarSource(i_point, iVar,  source_total);
+        }
       }
     }
 
@@ -115,9 +151,9 @@
     /*--- Set mass diffusivity based on thermodynamic state. ---*/
     auto T = flowNodes->GetTemperature(i_point);
     fluid_model_local->SetTDState_T(T, scalars);
-    /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table ---*/
+    /*--- set the diffusivity in the fluid model to the diffusivity obtained from the lookup table, multiplied by flame thickness correction factor ---*/
     for (auto i_scalar = 0u; i_scalar < nVar; ++i_scalar) {
-      nodes->SetDiffusivity(i_point, fluid_model_local->GetMassDiffusivity(i_scalar), i_scalar);
+      nodes->SetDiffusivity(i_point, F * (fluid_model_local->GetMassDiffusivity(i_scalar)), i_scalar);
     }
 
     /*--- Obtain preferential diffusion scalar values. ---*/
@@ -213,8 +249,6 @@
 
     for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) {
       fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel();
-      if (flame_front_ignition) prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
-
       for (auto iVar = 0u; iVar < nVar; iVar++) scalar_init[iVar] = config->GetSpecies_Init()[iVar];
 
       /*--- Set enthalpy based on initial temperature and scalars. ---*/
@@ -228,6 +262,7 @@
 
         if (flame_front_ignition) {
 
+          prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init);
           /*--- Determine if point is above or below the plane, assuming the normal
             is pointing towards the burned region. ---*/
           point_loc = 0.0;
@@ -383,7 +418,7 @@
 
 void CSpeciesFlameletSolver::Source_Residual(CGeometry* geometry, CSolver** solver_container,
                                              CNumerics** numerics_container, CConfig* config, unsigned short iMesh) {
-
+   
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto i_point = 0u; i_point < nPointDomain; i_point++) {
     /*--- Add source terms from the lookup table directly to the residual. ---*/
@@ -516,7 +551,7 @@
 }
 
 unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CFluidModel* fluid_model_local,
-                                                       unsigned long iPoint, const vector<su2double>& scalars) {
+                                                       unsigned long iPoint, const vector<su2double>& scalars, const su2double F) {
   /*--- Compute total source terms from the production and consumption. ---*/
 
   vector<su2double> table_sources(flamelet_config_options.n_control_vars + 2 * flamelet_config_options.n_user_scalars);
@@ -538,8 +573,9 @@
     su2double source_cons = table_sources[flamelet_config_options.n_control_vars + 2 * i_aux + 1];
     source_scalar[flamelet_config_options.n_control_vars + i_aux] = source_prod + source_cons * y_aux;
   }
+  /*--- Source term is divided by flame thickness correction factor to improve stability on coarse grids. ---*/
   for (auto i_scalar = 0u; i_scalar < nVar; i_scalar++)
-    nodes->SetScalarSource(iPoint, i_scalar, source_scalar[i_scalar]);
+    nodes->SetScalarSource(iPoint, i_scalar,  F*source_scalar[i_scalar]);
   return misses;
 }
 
@@ -804,3 +840,61 @@
   su2double pv_burnt = scalars[I_PROGVAR] - delta;
   return pv_burnt;
 }
+
+
+su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) const {
+  su2double F{1.0};
+  if (global_flame_thickness != default_flame_thickness) {
+    su2double max_flame_vol = pow(global_flame_thickness, nDim);
+    F = max(1.0, geometry->nodes->GetVolume(iPoint) / max_flame_vol);
+  }
+  return F;
+}
+
+su2double CSpeciesFlameletSolver::GetOverallFlameThickness(CGeometry * geometry, CSolver **solver_container) const {
+  CFlowVariable const * flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+  su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{0.0}, gradpv_local{0.0}, gradpv_global{0.0}, Tmax_local{-1e6},Tmax_global{0.0};
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
+      su2double pv_local = nodes->GetSolution(iPoint, I_PROGVAR);
+      su2double T_local = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
+
+      /* Parallel projection of progress variable gradient against velocity */
+      su2double proj_grad_pv_u[MAXNDIM];
+      for (auto iDim=0u; iDim < nDim; iDim++) {
+        su2double gradpv = nodes->GetGradient(iPoint, I_PROGVAR, iDim);
+        su2double val_u = flowNodes->GetVelocity(iPoint, iDim);
+        proj_grad_pv_u[iDim] = gradpv * val_u * val_u / (flowNodes->GetVelocity2(iPoint) + 1e-5);
+      }
+
+      /* Parallel projection of temperature gradient against projected progress variable gradient */
+      su2double proj_grad_T_u = 0, gradT2{0}, mag_gradT{0};
+      for (auto iDim=0u; iDim < nDim; iDim++) {
+        su2double gradT = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
+        proj_grad_T_u += gradT * proj_grad_pv_u[iDim];
+        gradT2 += gradT*gradT;
+      }
+      mag_gradT = sqrt(gradT2);
+      proj_grad_T_u /= (gradT2 +  1e-5);
+      proj_grad_T_u *= mag_gradT;
+
+      /* Update minimum and maximum values. */
+      gradpv_local = max(gradpv_local, proj_grad_T_u);
+      pvmax_local = max(pvmax_local, pv_local);
+      pvmin_local = min(pvmin_local, pv_local);
+      Tmax_local = max(Tmax_local, T_local);
+  }
+  END_SU2_OMP_FOR
+  SU2_MPI::Barrier(SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&gradpv_local, &gradpv_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&pvmax_local, &pvmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&Tmax_local, &Tmax_global, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
+
+  /* Update flame thickness value. */
+  su2double flame_thickness{default_flame_thickness};
+  if (Tmax_global > 800.0) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+1e-5);
+
+  return flame_thickness;
+}

config_template.cfg

@@ -974,6 +974,11 @@ SPARK_INIT= (0.004, 0.0, 0.0, 1e-4, 100, 10)
 % The number of terms should equate the total number of species in the flamelet problem.
 SPARK_REACTION_RATES= (1000, 0, 0)
 
+% Thickened flame correction
+% Approximate the thickness of the flame and suppress source terms and enhance diffusivity in cells with a lengthscale higher
+% than the flame thickness.
+THICKENED_FLAME_CORRECTION= YES
+
 %
 % --------------------- INVERSE DESIGN SIMULATION -----------------------------%
 %