Improved combustion solver robustness for coarse grids

Common/include/option_structure.hpp

@@ -1446,6 +1446,8 @@ 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.*/
+
+  su2double flame_lengthscale{1e-3}; /*!< \brief Lengthscale above which chemical source terms are dampened.*/
 };
 
 /*!

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 FLAME_LENGTHSCALE \n DESCRIPTION: Lengthscale above which species source terms are dampened. \ingroup Config*/
+  addDoubleOption("FLAME_LENGTHSCALE", flamelet_ParsedOptions.flame_lengthscale, 1e-3);
+
   /*--- Options related to mass diffusivity and thereby the species solver. ---*/
 
   /*!\brief DIFFUSIVITY_MODEL\n DESCRIPTION: mass diffusivity model \n DEFAULT constant disffusivity \ingroup Config*/
@@ -5790,6 +5793,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION);
     /*--- We can have additional user defined transported scalars ---*/
     flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars;
+
+    if (flamelet_ParsedOptions.flame_lengthscale <= 0.0)
+      SU2_MPI::Error("Flame length scale value should be positive.", CURRENT_FUNCTION);
+
   }
 
   if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -81,10 +81,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 +106,14 @@ 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);   
 
  public:
   /*!

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -72,6 +72,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
                                            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);
+  const su2double F_default{1.0};
   unsigned long spark_iter_start, spark_duration;
   bool ignition = false;
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
@@ -81,7 +82,14 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     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();
+    unsigned long iter;
+    if (config->GetMultizone_Problem()) {
+      iter = config->GetOuterIter();
+    } else if (config->GetTime_Domain()) {
+      iter = config->GetTimeIter();
+    } else {
+      iter = config->GetInnerIter();
+    }
     ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
 
@@ -91,10 +99,17 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
   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. ---*/
+    su2double F = ThickenedFlameCorrection(geometry, i_point);
+
     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. ---*/
+    su2double F_source = config->GetTime_Domain() ? F_default : 1.0 / F;
+    unsigned long misses = SetScalarSources(config, fluid_model_local, i_point, scalars_vector, F_source);
 
     if (ignition) {
       /*--- Apply source terms within spark radius. ---*/
@@ -103,7 +118,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
       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]);
+          nodes->SetScalarSource(i_point, iVar,  nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar] / F);
       }
     }
 
@@ -115,9 +130,9 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
     /*--- 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,7 +228,6 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
 
     for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) {
       fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel();
-      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. ---*/
@@ -227,6 +241,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver**
 
         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;
@@ -382,7 +397,7 @@ void CSpeciesFlameletSolver::SetPreconditioner(CGeometry* geometry, CSolver** so
 
 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. ---*/
@@ -515,7 +530,7 @@ void CSpeciesFlameletSolver::BC_ConjugateHeat_Interface(CGeometry* geometry, CSo
 }
 
 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);
@@ -537,8 +552,9 @@ unsigned long CSpeciesFlameletSolver::SetScalarSources(const CConfig* config, CF
     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;
 }
 
@@ -803,3 +819,10 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo
   su2double pv_burnt = scalars[I_PROGVAR] - delta;
   return pv_burnt;
 }
+
+
+su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(CGeometry const * geometry, const unsigned long iPoint) {
+  su2double flame_vol = pow(flamelet_config_options.flame_lengthscale,nDim);
+  su2double F = max(1.0, geometry->nodes->GetVolume(iPoint) / flame_vol);
+  return F;
+}