diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1361a12072d..dff478254dd 100644 --- a/Common/include/option_structure.hpp +++ b/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.*/ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18884a52a35..33eeaec75f6 100644 --- a/Common/src/CConfig.cpp +++ b/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()) { diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index df61813a3e2..26f44286691 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/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& scalars); + const vector& 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& 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: /*! diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 622a1a6b068..e2548c19158 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/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 scalars_vector(nVar); + const su2double F_default{1.0}; unsigned long spark_iter_start, spark_duration; bool ignition = false; auto* flowNodes = su2staticcast_p(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& scalars) { + unsigned long iPoint, const vector& scalars, const su2double F) { /*--- Compute total source terms from the production and consumption. ---*/ vector 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; +} \ No newline at end of file