diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1361a12072d..7263b98acc3 100644 --- a/Common/include/option_structure.hpp +++ b/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. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 79fe6d6aab1..de803a1640b 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 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*/ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 67a8bb8112c..da5244c9298 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/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). diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index df61813a3e2..a882068cd59 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/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& scalars); + const vector& 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& 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;} }; diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 446db73dc13..2ca35e74c7c 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/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())); } diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index afc8ed5e2e5..015eac3daf6 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -39,6 +39,8 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con /*--- 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,6 +67,10 @@ CSpeciesFlameletSolver::CSpeciesFlameletSolver(CGeometry* geometry, CConfig* con /*--- 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, @@ -72,29 +78,57 @@ 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); + unsigned long spark_iter_start, spark_duration; bool ignition = false; auto* flowNodes = su2staticcast_p(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. ---*/ @@ -102,8 +136,10 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver 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 @@ 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,8 +249,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(); - 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 @@ 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; @@ -383,7 +418,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. ---*/ @@ -516,7 +551,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); @@ -538,8 +573,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; } @@ -804,3 +840,61 @@ 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) 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(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; +} \ No newline at end of file diff --git a/config_template.cfg b/config_template.cfg index 0dd67764bdf..62ca6b4e464 100644 --- a/config_template.cfg +++ b/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 -----------------------------% %