diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1361a12072d..b524a8b583b 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.*/ + bool thickenedflame_correction{true}; /*!< \brief Thickened flame correction. */ + su2double ignition_temperature = 1100.0; /*!< \brief Temperature above which the solution is considered ignited. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 79fe6d6aab1..911bd6e3b20 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1429,6 +1429,12 @@ 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: Coarse grid correction for source terms and diffusive fluxes in reacting flows. \ingroup Config*/ + addBoolOption("THICKENED_FLAME_CORRECTION", flamelet_ParsedOptions.thickenedflame_correction, true); + + /*!\brief IGNITION_TEMPERATURE \n DESCRIPTION: Temperature above which the solution is considered ignited. */ + addDoubleOption("IGNITION_TEMPERATURE", flamelet_ParsedOptions.ignition_temperature, 1100.0); + /*--- 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..a792bba6e80 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(const CGeometry* geometry, 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. + */ + su2double GetFlameThickness() const override {return global_flame_thickness;} }; diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 446db73dc13..5c0d44ef5f6 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, "FLAME_THICKNESS", "Flame thickness used for thickened flame correction model.", HistoryFieldType::COEFFICIENT); /*--- 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..6e121e6710c 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,58 @@ 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 calc_thickness = GetOverallFlameThickness(geometry, solver_container); + su2double test_thickness = min(default_flame_thickness, calc_thickness); + 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 +137,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 +152,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 +250,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 +263,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 +419,6 @@ 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; } @@ -798,9 +834,76 @@ su2double CSpeciesFlameletSolver::GetBurntProgressVariable(CFluidModel* fluid_mo bool outside = false; while (!outside) { fluid_model->SetTDState_T(300, scalars); - if (fluid_model->GetExtrapolation() == 1 || (fluid_model->GetTemperature()>1000.)) outside = true; + if (fluid_model->GetExtrapolation() == 1 || fluid_model->GetTemperature()>flamelet_config_options.ignition_temperature) outside = true; scalars[I_PROGVAR] += delta; } su2double pv_burnt = scalars[I_PROGVAR] - delta; return pv_burnt; } + + +su2double CSpeciesFlameletSolver::ThickenedFlameCorrection(const CGeometry* geometry, unsigned long iPoint) const { + su2double F{1.0}; + if (fabs(global_flame_thickness - default_flame_thickness) > EPS * max(1.0, fabs(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 { + const CFlowVariable* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + su2double pvmax_local{-1e3}, pvmin_local{1e3}, pvmax_global{0.0},pvmin_global{1e3}, 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 < nPointDomain; 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]={0}; + 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) + EPS); + } + + /* Parallel projection of temperature gradient against projected progress variable gradient */ + su2double gradT[MAXNDIM]={0}; + for (auto iDim=0u; iDim < nDim; iDim++) + gradT[iDim] = flowNodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim); + + su2double proj_grad_T_u = GeometryToolbox::DotProduct(nDim, gradT, proj_grad_pv_u); + su2double mag_gradT = GeometryToolbox::Norm(nDim, gradT); + + proj_grad_T_u /= (pow(mag_gradT,2) + EPS); + 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 + su2double* MyFlameThickness = new su2double[3]; + su2double* TotalFlameThickness = new su2double[3]; + MyFlameThickness[0] = gradpv_local; + MyFlameThickness[1] = pvmax_local; + MyFlameThickness[2] = Tmax_local; + + SU2_MPI::Allreduce(MyFlameThickness, TotalFlameThickness, 3, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&pvmin_local, &pvmin_global, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); + gradpv_global = TotalFlameThickness[0]; + pvmax_global = TotalFlameThickness[1]; + Tmax_global = TotalFlameThickness[2]; + + delete [] MyFlameThickness; + delete [] TotalFlameThickness; + + /* Update flame thickness value. */ + su2double flame_thickness{default_flame_thickness}; + if (Tmax_global > flamelet_config_options.ignition_temperature) flame_thickness = (pvmax_global - pvmin_global) / (gradpv_global+EPS); + + return flame_thickness; +} \ No newline at end of file diff --git a/config_template.cfg b/config_template.cfg index 0dd67764bdf..0ecf0d5e82d 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -934,7 +934,7 @@ CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL) % Method used to ignite the solution: % FLAME_FRONT : the flame is initialized using a plane, defined by a point and a normal. On one side, the solution is initialized % using 'burnt' conditions and on the other side 'unburnt' conditions. The normal points in the direction of the 'burnt' -% condition. +% condition. The 'burnt' condition is defined as the value of the progress variable for which the temperature exceeds the IGNITION_TEMPERATURE. % SPARK : the solution is ignited through application of a set of source terms within a specified region for a set number % of solver iterations. This artificial spark can be applied after a certain number of iterations has passed, allowing for % the flow to evolve before igniting the mixture. @@ -974,6 +974,15 @@ 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 + +% Flame ignition temperature +% Flame thickness calculations are performed when the maximum temperature in the flow solution exceeds this value. +IGNITION_TEMPERATURE= 1100.0 + % % --------------------- INVERSE DESIGN SIMULATION -----------------------------% %