diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 442d2124f36b..5ef3af41fc35 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -9,6 +9,7 @@ on: branches: - 'develop' - 'master' + - 'feature_flamelet' jobs: build: @@ -60,7 +61,7 @@ jobs: strategy: fail-fast: false matrix: - testscript: ['tutorials.py', 'parallel_regression.py', 'flamelet_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py'] + testscript: ['tutorials.py', 'parallel_regression.py', 'flamelet_regression.py', 'multicomp_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py'] include: - testscript: 'tutorials.py' tag: MPI @@ -68,6 +69,8 @@ jobs: tag: MPI - testscript: 'flamelet_regression.py' tag: MPI + - testscript: 'multicomp_regression.py' + tag: MPI - testscript: 'parallel_regression_AD.py' tag: MPI - testscript: 'serial_regression.py' @@ -92,7 +95,7 @@ jobs: uses: docker://su2code/test-su2:20200303 with: # -t -c - args: -b ${{github.ref}} -t develop -c feature_flamelet -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c feature_multicomp -s ${{matrix.testscript}} unit_tests: runs-on: ubuntu-latest name: Unit Tests diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index b13bb57aaa06..69305ede01b4 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -378,6 +378,7 @@ class CConfig { su2double *Surface_TotalPressure; /*!< \brief Total pressure at the boundaries. */ su2double *Surface_CO; /*!< \brief Mass fraction of CO at the boundaries. */ su2double *Surface_NOx; /*!< \brief Mass fraction of NO at the boundaries. */ + su2double *Surface_CH4; /*!< \brief Mass fraction of CH4 at the boundaries. */ //su2double **Surface_Scalar; su2double *Surface_PressureDrop; /*!< \brief Pressure drop between boundaries. */ su2double *Surface_DC60; /*!< \brief Specified surface DC60 for nacelle boundaries. */ @@ -792,8 +793,8 @@ class CConfig { Beta_Factor, /*!< \brief Value of the epsilon^2 multiplier for Beta for the incompressible preconditioner. */ Gas_Constant, /*!< \brief Specific gas constant. */ Gas_ConstantND, /*!< \brief Non-dimensional specific gas constant. */ - Molecular_Weight, /*!< \brief Molecular weight of an incompressible ideal gas (g/mol). */ - Specific_Heat_Cp, /*!< \brief Specific heat at constant pressure. */ + *Molecular_Weight, /*!< \brief Molecular weight of an incompressible ideal gas (g/mol). */ + *Specific_Heat_Cp, /*!< \brief Specific heat at constant pressure. */ Specific_Heat_CpND, /*!< \brief Non-dimensional specific heat at constant pressure. */ Specific_Heat_Cv, /*!< \brief Specific heat at constant volume. */ Specific_Heat_CvND, /*!< \brief Non-dimensional specific heat at constant volume. */ @@ -810,16 +811,16 @@ class CConfig { Pressure_Critical, /*!< \brief Critical Pressure for real fluid model. */ Density_Critical, /*!< \brief Critical Density for real fluid model. */ Acentric_Factor, /*!< \brief Acentric Factor for real fluid model. */ - Mu_Constant, /*!< \brief Constant viscosity for ConstantViscosity model. */ + *Mu_Constant, /*!< \brief Constant viscosity for ConstantViscosity model. */ Mu_ConstantND, /*!< \brief Non-dimensional constant viscosity for ConstantViscosity model. */ - Thermal_Conductivity_Constant, /*!< \brief Constant thermal conductivity for ConstantConductivity model. */ + *Thermal_Conductivity_Constant, /*!< \brief Constant thermal conductivity for ConstantConductivity model. */ Thermal_Conductivity_ConstantND, /*!< \brief Non-dimensional constant thermal conductivity for ConstantConductivity model. */ *Scalar_Init, /*!< \brief Initial uniform value for scalar transport. */ - Mu_Ref, /*!< \brief Reference viscosity for Sutherland model. */ + *Mu_Ref, /*!< \brief Reference viscosity for Sutherland model. */ Mu_RefND, /*!< \brief Non-dimensional reference viscosity for Sutherland model. */ - Mu_Temperature_Ref, /*!< \brief Reference temperature for Sutherland model. */ + *Mu_Temperature_Ref, /*!< \brief Reference temperature for Sutherland model. */ Mu_Temperature_RefND, /*!< \brief Non-dimensional reference temperature for Sutherland model. */ - Mu_S, /*!< \brief Reference S for Sutherland model. */ + *Mu_S, /*!< \brief Reference S for Sutherland model. */ Mu_SND; /*!< \brief Non-dimensional reference S for Sutherland model. */ array CpPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for specific heat Cp. */ array MuPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */ @@ -855,8 +856,8 @@ class CConfig { wallModel_B, /*!< \brief constant B for turbulence wall modeling */ wallModel_RelFac, /*!< \brief relaxation factor for the Newton method used in the wall model */ wallModel_MinYplus; /*!< \brief minimum Y+ value, below which the wall model is not used anymore */ - su2double Prandtl_Lam, /*!< \brief Laminar Prandtl number for the gas. */ - Prandtl_Turb, /*!< \brief Turbulent Prandtl number for the gas. */ + su2double *Prandtl_Lam, /*!< \brief Laminar Prandtl number for the gas. */ + *Prandtl_Turb, /*!< \brief Turbulent Prandtl number for the gas. */ Length_Ref, /*!< \brief Reference length for non-dimensionalization. */ Pressure_Ref, /*!< \brief Reference pressure for non-dimensionalization. */ Temperature_Ref, /*!< \brief Reference temperature for non-dimensionalization.*/ @@ -1147,11 +1148,21 @@ class CConfig { unsigned short nScreenOutput, /*!< \brief Number of screen output variables (max: 6). */ nHistoryOutput, nVolumeOutput; /*!< \brief Number of variables printed to the history file. */ bool Multizone_Residual; /*!< \brief Determines if memory should be allocated for the multizone residual. */ - - unsigned short n_scalars; - unsigned short n_lookups; - unsigned short n_table_sources; /* the number of transported scalars for combustion */ - + + unsigned short n_scalars, /*!< \brief Number of transported scalars. */ + n_species, /*!< \brief Number of species in multicomponent flow. */ + nMolecular_Weight, /*!< \brief Number of species molecular weights. */ + nSpecific_Heat_Cp, /*!< \brief Number of species specific heat constants at constant pressure. */ + nMu_Constant, /*!< \brief Number of species constant viscosities. */ + nMu_Ref, /*!< \brief Number of species reference constants for Sutherland model. */ + nMu_Temperature_Ref, /*!< \brief Number of species reference temperature for Sutherland model. */ + nMu_S, /*!< \brief Number of species reference S for Sutherland model. */ + nThermal_Conductiviy_Constant, /*!< \brief Number of species constant thermal conductivity. */ + nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */ + nPrandtl_Turb, /*!< \brief Number of species turbulent Prandtl number. */ + n_lookups, /*!< \brief Number of transported scalars for combustion. */ + n_table_sources; /*!< \brief Number of transported scalars for combustion. */ + vector table_scalar_names; /*!< \brief vector to store names of scalar variables. */ vector table_source_names; /*!< \brief vector to store names of scalar source variables. */ string* table_lookup_names; /*!< \brief vector to store names of look up variables. */ @@ -1180,7 +1191,7 @@ class CConfig { POD_KIND POD_Basis_Gen; /*!< \brief Type of POD basis generation (static or incremental). */ unsigned short maxBasisDim, /*!< \brief Maximum number of POD basis dimensions. */ rom_save_freq; /*!< \brief Frequency of unsteady time steps to save. */ - + /* other NEMO configure options*/ unsigned short nSpecies, /*!< \brief No of species present in flow */ iWall_Catalytic, @@ -1619,13 +1630,13 @@ class CConfig { * \brief Get the value of the molecular weight for an incompressible ideal gas (g/mol). * \return Value of the molecular weight for an incompressible ideal gas (g/mol). */ - su2double GetMolecular_Weight(void) const { return Molecular_Weight; } + su2double GetMolecular_Weight(unsigned short val_index = 0) const { return Molecular_Weight [val_index]; } /*! * \brief Get the value of specific heat at constant pressure. * \return Value of the constant: Cp */ - su2double GetSpecific_Heat_Cp(void) const { return Specific_Heat_Cp; } + su2double GetSpecific_Heat_Cp(unsigned short val_index = 0) const { return Specific_Heat_Cp [val_index]; } /*! * \brief Get the non-dimensional value of specific heat at constant pressure. @@ -1714,13 +1725,13 @@ class CConfig { * \brief Get the value of the laminar Prandtl number. * \return Laminar Prandtl number. */ - su2double GetPrandtl_Lam(void) const { return Prandtl_Lam; } + su2double GetPrandtl_Lam(unsigned short val_index = 0) const { return Prandtl_Lam [val_index]; } /*! * \brief Get the value of the turbulent Prandtl number. * \return Turbulent Prandtl number. */ - su2double GetPrandtl_Turb(void) const { return Prandtl_Turb; } + su2double GetPrandtl_Turb(unsigned short val_index = 0) const { return Prandtl_Turb [val_index]; } /*! * \brief Get the value of the von Karman constant kappa for turbulence wall modeling. @@ -2062,21 +2073,21 @@ class CConfig { * \return Initial temperature for incompressible flows. */ su2double GetInc_Temperature_Init(void) const { return Inc_Temperature_Init; } - + /*! * \brief Get the flag for activating scalar transport clipping * \return Flag for scalar clipping */ bool GetScalar_Clipping(void) { return Scalar_Clipping; } - + bool GetEnableRemeshing(void) { return enable_remeshing; } bool GetUseWeakScalarBC(void) { return use_weak_scalar_bc; } - + su2double *GetFFDBounds(void) { return ffd_bounds; } - + /*! * \brief Get the flame offset for flamelet model initialization * \return flame offset for flamelet model initialization @@ -2203,7 +2214,7 @@ class CConfig { /*! * \brief Get the file name of the look up table - * \return File name of the look up table + * \return File name of the look up table */ string GetFileNameLUT(void){ return file_name_lut; }; @@ -2457,7 +2468,8 @@ class CConfig { * \brief Set the value of the specific heat at constant pressure (incompressible fluids with energy equation). * \param[in] val_specific_heat_cp - specific heat at constant pressure. */ - void SetSpecific_Heat_Cp(su2double val_specific_heat_cp) { Specific_Heat_Cp = val_specific_heat_cp; } + void SetSpecific_Heat_Cp(su2double val_specific_heat_cp) { Specific_Heat_Cp[0] = val_specific_heat_cp; } + void SetSpecific_Heat_Cp(su2double val_specific_heat_cp, unsigned short val_index) { Specific_Heat_Cp[val_index] = val_specific_heat_cp; } /*! * \brief Set the non-dimensional value of the specific heat at constant pressure (incompressible fluids with energy equation). @@ -3826,12 +3838,12 @@ class CConfig { * \return Mass diffusivity model. */ DIFFUSIVITYMODEL GetKind_DiffusivityModel(void) const { return Kind_DiffusivityModel; } - + /*! * \brief Get the value of the constant viscosity. * \return Constant viscosity. */ - su2double GetMu_Constant(void) const { return Mu_Constant; } + su2double GetMu_Constant(unsigned short val_index = 0) const { return Mu_Constant [val_index]; } /*! * \brief Get the value of the non-dimensional constant viscosity. @@ -3843,44 +3855,43 @@ class CConfig { * \brief Get the value of the thermal conductivity. * \return Thermal conductivity. */ - su2double GetThermal_Conductivity_Constant(void) const { return Thermal_Conductivity_Constant; } + su2double GetThermal_Conductivity_Constant(unsigned short val_index = 0) const { return Thermal_Conductivity_Constant [val_index]; } /*! * \brief Get the value of the non-dimensional thermal conductivity. * \return Non-dimensional thermal conductivity. */ - su2double GetThermal_Conductivity_ConstantND(void) const { return Thermal_Conductivity_ConstantND; } - + su2double GetThermal_Conductivity_ConstantND(void) const { return Thermal_Conductivity_ConstantND; } /*! * \brief Get the value of the constant mass diffusivity for scalar transport. * \return Constant mass diffusivity. */ su2double GetDiffusivity_Constant(void) const { return Diffusivity_Constant; } - + /*! * \brief Get the value of the non-dimensional constant mass diffusivity. * \return Non-dimensional constant mass diffusivity. */ su2double GetDiffusivity_ConstantND(void) const { return Diffusivity_ConstantND; } - + /*! * \brief Get the value of the laminar Schmidt number for scalar transport. * \return Laminar Schmidt number for scalar transport. */ su2double GetSchmidt_Lam(void) const { return Schmidt_Lam; } - + /*! * \brief Get the value of the turbulent Schmidt number for scalar transport. * \return Turbulent Schmidt number for scalar transport. */ su2double GetSchmidt_Turb(void) const { return Schmidt_Turb; } - + /*! * \brief Get the value of the reference viscosity for Sutherland model. * \return The reference viscosity. */ - su2double GetMu_Ref(void) const { return Mu_Ref; } + su2double GetMu_Ref(unsigned short val_index = 0) const { return Mu_Ref [val_index]; } /*! * \brief Get the value of the non-dimensional reference viscosity for Sutherland model. @@ -3892,7 +3903,7 @@ class CConfig { * \brief Get the value of the reference temperature for Sutherland model. * \return The reference temperature. */ - su2double GetMu_Temperature_Ref(void) const { return Mu_Temperature_Ref; } + su2double GetMu_Temperature_Ref(unsigned short val_index = 0) const { return Mu_Temperature_Ref [val_index]; } /*! * \brief Get the value of the non-dimensional reference temperature for Sutherland model. @@ -3904,7 +3915,7 @@ class CConfig { * \brief Get the value of the reference S for Sutherland model. * \return The reference S. */ - su2double GetMu_S(void) const { return Mu_S; } + su2double GetMu_S(unsigned short val_index = 0) const { return Mu_S [val_index]; } /*! * \brief Get the value of the non-dimensional reference S for Sutherland model. @@ -3986,12 +3997,12 @@ class CConfig { * \brief Set the value of the non-dimensional constant mass diffusivity. */ void SetDiffusivity_ConstantND(su2double diffusivity_const) { Diffusivity_ConstantND = diffusivity_const; } - + /*! * \brief Set the value of the reference mass diffusivity. */ void SetDiffusivity_Ref(su2double diffusivity_ref); - + /*! * \brief Set the value of the non-dimensional reference viscosity for Sutherland model. */ @@ -4135,7 +4146,7 @@ class CConfig { * \return relaxation coefficient of the linear solver for the implicit formulation. */ su2double GetRelaxation_Factor_Scalar(void) { return Relaxation_Factor_Scalar; } - + /*! * \brief Get the relaxation coefficient of the CHT coupling. * \return relaxation coefficient of the CHT coupling. @@ -4344,8 +4355,8 @@ class CConfig { * \return Kind of the scalar transport model. */ unsigned short GetKind_Scalar_Model(void) const { return Kind_Scalar_Model; }; - - + + /*! * \brief Get the kind of time integration method. * \note This is the information that the code will use, the method will @@ -4435,7 +4446,7 @@ class CConfig { * \return MUSCL scheme. */ bool GetMUSCL_Scalar(void) { return MUSCL_Scalar; } - + /*! * \brief Get if the upwind scheme used MUSCL or not. * \note This is the information that the code will use, the method will @@ -4619,7 +4630,7 @@ class CConfig { * \return Method for limiting the spatial gradients solving the scalar transport equations. */ unsigned short GetKind_SlopeLimit_Scalar(void) { return Kind_SlopeLimit_Scalar; } - + /*! * \brief Get the method for limiting the spatial gradients. * \return Method for limiting the spatial gradients solving the adjoint turbulent equation. @@ -4735,7 +4746,7 @@ class CConfig { * \return Kind of integration scheme for the scalar transport equations. */ unsigned short GetKind_TimeIntScheme_Scalar(void) const { return Kind_TimeIntScheme_Scalar; } - + /*! * \brief Get the kind of convective numerical scheme for the scalar transport equations (upwind). * \note This value is obtained from the config file, and it is constant @@ -4743,7 +4754,7 @@ class CConfig { * \return Kind of convective numerical scheme for the scalar transport equations. */ unsigned short GetKind_ConvNumScheme_Scalar(void) const { return Kind_ConvNumScheme_Scalar; } - + /*! * \brief Get the kind of center convective numerical scheme for the scalar transport equations. * \note This value is obtained from the config file, and it is constant @@ -4751,7 +4762,7 @@ class CConfig { * \return Kind of center convective numerical scheme for the scalar transport equations. */ unsigned short GetKind_Centered_Scalar(void) const { return Kind_Centered_Scalar; } - + /*! * \brief Get the kind of upwind convective numerical scheme for the scalar transport equations. * \note This value is obtained from the config file, and it is constant @@ -4759,7 +4770,7 @@ class CConfig { * \return Kind of upwind convective numerical scheme for the scalar transport equations. */ unsigned short GetKind_Upwind_Scalar(void) const { return Kind_Upwind_Scalar; } - + /*! * \brief Get the kind of integration scheme (implicit) * for the turbulence equations. @@ -6601,11 +6612,11 @@ class CConfig { /*! * \brief Get the scalar values at an inlet boundary * \param[in] val_index - Index corresponding to the inlet boundary. - * \return The inlet scalar values. + * \return The inlet scalar values. */ // nijso: TODO we do not need inlet enthalpy, it is computed from temperature! su2double* GetInlet_ScalarVal(string val_index) const; - + /*! * \brief Get the temperature at a supersonic inlet boundary. * \param[in] val_index - Index corresponding to the inlet boundary. @@ -6666,7 +6677,7 @@ class CConfig { * \return Value of the CFL reduction for scalar transport equations. */ su2double GetCFLRedCoeff_Scalar(void) { return CFLRedCoeff_Scalar; } - + /*! * \brief Get the flow direction unit vector at an inlet boundary. * \param[in] val_index - Index corresponding to the inlet boundary. @@ -7658,6 +7669,13 @@ class CConfig { */ void SetSurface_NOx(unsigned short val_imarker, su2double val_surface_nox){ Surface_NOx[val_imarker] = val_surface_nox; }; + /*! + * \brief Set the CH4 mass fraction at the surface. + * \param[in] val_imarker - Index corresponding to the outlet boundary. + * \param[in] val_surface_no - Value of the CH4 mass fraction. + */ + void SetSurface_CH4(unsigned short val_imarker, su2double val_surface_ch4){ Surface_CH4[val_imarker] = val_surface_ch4; }; + /*! * \brief Set the pressure drop between two surfaces. * \param[in] val_marker - Index corresponding to the outlet boundary. @@ -7930,14 +7948,14 @@ class CConfig { * \return The CO mass fraction. */ su2double GetSurface_CO(unsigned short val_imarker) const { return Surface_CO[val_imarker]; } - + /*! * \brief Get the scalar mass fraction at an outlet boundary. * \param[in] val_index - Index corresponding to the outlet boundary. * \return The scalar mass fraction. */ //su2double GetSurface_Scalar(unsigned short val_imarker, unsigned short val_i_scalar) const { return Surface_Scalar[val_imarker][val_i_scalar]; } - + /*! * \brief Get the NOx mass fraction at an outlet boundary. * \param[in] val_index - Index corresponding to the outlet boundary. @@ -7945,6 +7963,13 @@ class CConfig { */ su2double GetSurface_NOx(unsigned short val_imarker) const { return Surface_NOx[val_imarker]; }; + /*! + * \brief Get the CH4 mass fraction at an outlet boundary. + * \param[in] val_index - Index corresponding to the outlet boundary. + * \return The CH4 mass fraction. + */ + su2double GetSurface_CH4(unsigned short val_imarker) const { return Surface_CH4[val_imarker]; }; + /*! * \brief Get the pressure drop between two surfaces. * \param[in] val_index - Index corresponding to the outlet boundary. @@ -9399,19 +9424,26 @@ class CConfig { void SetNScalars(unsigned short n_scalars) { this->n_scalars = n_scalars; } /*! - * \brief Get the number of transported scalars for combustion + * \brief Get the number of transported scalars for combustion */ unsigned short GetNScalars(void) const { return n_scalars; } /*! - * \brief Get the number of transported scalars for combustion + * \brief Get the number of transported scalars required for initialisation + */ + unsigned short GetNScalarsInit(void) const { return nScalar_Init; } + + void SetNScalarsInit(unsigned short nScalar_Init) { this->nScalar_Init = nScalar_Init; } + + /*! + * \brief Get the number of transported scalars for combustion */ unsigned short GetNLookups(void) const { return n_lookups; } void SetNTableSources(unsigned short n_table_sources) { this->n_table_sources = n_table_sources; } /*! - * \brief Get the number of transported scalars source terms for combustion + * \brief Get the number of transported scalars source terms for combustion */ unsigned short GetNTableSources(void) const { return n_table_sources; } @@ -9446,13 +9478,13 @@ class CConfig { * \brief Get the scalar source term name i_source */ string GetTableSourceName(unsigned short i_source) const { return table_source_names.at(i_source); } - + /*! * \brief Get the maximum bound for scalar transport clipping * \return Maximum value for scalar clipping */ su2double *GetScalar_Clipping_Max(void) { return Scalar_Clipping_Max; } - + /*! * \brief Get the minimum bound for scalar transport clipping * \return Minimum value for scalar clipping @@ -9463,7 +9495,7 @@ class CConfig { * \return Maximum value for scalar clipping */ su2double GetScalar_Clipping_Max(unsigned short iVal) { return Scalar_Clipping_Max[iVal]; } - + /*! * \brief Get the minimum bound for scalar transport clipping * \return Minimum value for scalar clipping @@ -9475,7 +9507,7 @@ class CConfig { * \return Minimum value for scalar clipping */ su2double GetScalar_Init(unsigned short ival) { return Scalar_Init[ival]; } - + /*! * \brief Get the convergence fields for monitoring * \param[in] iField - Index of the field @@ -9651,25 +9683,25 @@ class CConfig { * \return True if specified in config file. */ bool GetSave_libROM(void) const {return libROM; } - + /*! * \brief Get the name of the file for libROM to save. * \return Filename prefix for libROM to save to (default: "su2"). */ string GetlibROMbase_FileName(void) const { return libROMbase_FileName; } - + /*! * \brief Static or incremental toggle for POD basis generation type. * \return Type of POD generation type */ POD_KIND GetKind_PODBasis(void) const { return POD_Basis_Gen; } - + /*! * \brief Get maximum number of POD basis dimensions (default: 100). * \return Maximum number of POD basis vectors. */ unsigned short GetMax_BasisDim(void) const { return maxBasisDim; } - + /*! * \brief Get frequency of unsteady time steps to save (default: 1). * \return Save frequency for unsteady time steps. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 71bbc9131ec1..befb5f7d521f 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -564,6 +564,7 @@ enum ENUM_FLUIDMODEL { MUTATIONPP = 7, /*!< \brief Mutation++ gas model for nonequilibrium flow. */ SU2_NONEQ = 8, /*!< \brief User defined gas model for nonequilibrium flow. */ FLAMELET_FLUID_MODEL = 9, /*!, \brief Flamelet model */ + MIXTURE_FLUID_MODEL = 10,/*!, \brief Mixture model */ }; static const MapType FluidModel_Map = { MakePair("STANDARD_AIR", STANDARD_AIR) @@ -575,6 +576,7 @@ static const MapType FluidModel_Map = { MakePair("INC_IDEAL_GAS_POLY", INC_IDEAL_GAS_POLY) MakePair("MUTATIONPP", MUTATIONPP) MakePair("FLAMELET_FLUID_MODEL", FLAMELET_FLUID_MODEL) + MakePair("MIXTURE_FLUID_MODEL", MIXTURE_FLUID_MODEL) MakePair("SU2_NONEQ", SU2_NONEQ) }; @@ -712,12 +714,14 @@ enum class DIFFUSIVITYMODEL { CONSTANT_DIFFUSIVITY, /*!< \brief Constant mass diffusivity for scalar transport. */ CONSTANT_SCHMIDT, /*!< \brief Constant Schmidt number for mass diffusion in scalar transport. */ FLAMELET, /*!< \brief flamelet model */ + UNITY_LEWIS, /*!< \brief Unity Lewis model */ }; static const MapType DiffusivityModel_Map = { MakePair("CONSTANT_DIFFUSIVITY", DIFFUSIVITYMODEL::CONSTANT_DIFFUSIVITY) MakePair("CONSTANT_SCHMIDT", DIFFUSIVITYMODEL::CONSTANT_SCHMIDT) MakePair("FLAMELET", DIFFUSIVITYMODEL::FLAMELET) + MakePair("UNITY_LEWIS", DIFFUSIVITYMODEL::UNITY_LEWIS) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 6acdc54fd4f1..ae3138319b25 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -909,7 +909,7 @@ void CConfig::SetPointersNull(void) { Surface_Temperature = nullptr; Surface_Pressure = nullptr; Surface_Density = nullptr; Surface_Enthalpy = nullptr; Surface_NormalVelocity = nullptr; Surface_TotalTemperature = nullptr; Surface_TotalPressure = nullptr; Surface_PressureDrop = nullptr; Surface_DC60 = nullptr; Surface_IDC = nullptr; - Surface_CO = nullptr; Surface_NOx = nullptr; + Surface_CO = nullptr; Surface_NOx = nullptr; Surface_CH4 = nullptr; //Surface_Scalar = nullptr; Outlet_MassFlow = nullptr; Outlet_Density = nullptr; Outlet_Area = nullptr; @@ -1016,7 +1016,7 @@ void CConfig::SetPointersNull(void) { Scalar_Init = nullptr; Scalar_Clipping_Min = nullptr; Scalar_Clipping_Max = nullptr; - + /*--- Variable initialization ---*/ TimeIter = 0; @@ -1048,6 +1048,15 @@ void CConfig::SetPointersNull(void) { Kind_TimeNumScheme = EULER_IMPLICIT; Gas_Composition = nullptr; + Molecular_Weight = nullptr; + Mu_Constant = nullptr; + Mu_Ref = nullptr; + Mu_Temperature_Ref = nullptr; + Mu_S = nullptr; + Specific_Heat_Cp = nullptr; + Thermal_Conductivity_Constant = nullptr; + Prandtl_Lam = nullptr; + Prandtl_Turb = nullptr; } @@ -1156,13 +1165,13 @@ void CConfig::SetConfig_Options() { /*!\brief GAMMA_VALUE \n DESCRIPTION: Ratio of specific heats (1.4 (air), only for compressible flows) \ingroup Config*/ addDoubleOption("GAMMA_VALUE", Gamma, 1.4); /*!\brief CP_VALUE \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/ - addDoubleOption("SPECIFIC_HEAT_CP", Specific_Heat_Cp, 1004.703); + addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp); /*!\brief CP_VALUE \n DESCRIPTION: Specific heat at constant volume, Cp (717.645 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/ addDoubleOption("SPECIFIC_HEAT_CV", Specific_Heat_Cv, 717.645); /*!\brief THERMAL_EXPANSION_COEFF \n DESCRIPTION: Thermal expansion coefficient (0.00347 K^-1 (air), used for Boussinesq approximation for liquids/non-ideal gases) \ingroup Config*/ addDoubleOption("THERMAL_EXPANSION_COEFF", Thermal_Expansion_Coeff, 0.00347); /*!\brief MOLECULAR_WEIGHT \n DESCRIPTION: Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) \ingroup Config*/ - addDoubleOption("MOLECULAR_WEIGHT", Molecular_Weight, 28.96); + addDoubleListOption("MOLECULAR_WEIGHT", nMolecular_Weight, Molecular_Weight); ///* DESCRIPTION: Specify if Mutation++ library is used */ /*--- Reading gas model as string or integer depending on TC library used. ---*/ @@ -1203,16 +1212,16 @@ void CConfig::SetConfig_Options() { /*--- Options related to Constant Viscosity Model ---*/ /* DESCRIPTION: default value for AIR */ - addDoubleOption("MU_CONSTANT", Mu_Constant , 1.716E-5); + addDoubleListOption("MU_CONSTANT", nMu_Constant, Mu_Constant); /*--- Options related to Sutherland Viscosity Model ---*/ /* DESCRIPTION: Sutherland Viscosity Ref default value for AIR SI */ - addDoubleOption("MU_REF", Mu_Ref, 1.716E-5); + addDoubleListOption("MU_REF", nMu_Ref, Mu_Ref); /* DESCRIPTION: Sutherland Temperature Ref, default value for AIR SI */ - addDoubleOption("MU_T_REF", Mu_Temperature_Ref, 273.15); + addDoubleListOption("MU_T_REF", nMu_Temperature_Ref, Mu_Temperature_Ref); /* DESCRIPTION: Sutherland constant, default value for AIR SI */ - addDoubleOption("SUTHERLAND_CONSTANT", Mu_S, 110.4); + addDoubleListOption("SUTHERLAND_CONSTANT", nMu_S, Mu_S); /*--- Options related to Thermal Conductivity Model ---*/ @@ -1224,7 +1233,7 @@ void CConfig::SetConfig_Options() { /*--- Options related to Constant Thermal Conductivity Model ---*/ /* DESCRIPTION: default value for AIR */ - addDoubleOption("THERMAL_CONDUCTIVITY_CONSTANT", Thermal_Conductivity_Constant , 0.0257); + addDoubleListOption("THERMAL_CONDUCTIVITY_CONSTANT", nThermal_Conductiviy_Constant, Thermal_Conductivity_Constant); /*--- Options related to temperature polynomial coefficients for fluid models. ---*/ @@ -1236,7 +1245,7 @@ void CConfig::SetConfig_Options() { addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, kt_polycoeffs.data()); /*--- Options related to mass diffusivity ---*/ - + addEnumOption("DIFFUSIVITY_MODEL", Kind_DiffusivityModel, DiffusivityModel_Map, DIFFUSIVITYMODEL::CONSTANT_DIFFUSIVITY); /* DESCRIPTION: default value for AIR */ addDoubleOption("DIFFUSIVITY_CONSTANT", Diffusivity_Constant , 0.001); @@ -1244,15 +1253,15 @@ void CConfig::SetConfig_Options() { addDoubleOption("SCHMIDT_LAM", Schmidt_Lam, 1.0); /*!\brief SCHMIDT_TURB \n DESCRIPTION: Turbulent Schmidt number of mass diffusion \n DEFAULT 0.90 \ingroup Config*/ addDoubleOption("SCHMIDT_TURB", Schmidt_Turb, 0.7); - + /*!\brief REYNOLDS_NUMBER \n DESCRIPTION: Reynolds number (non-dimensional, based on the free-stream values). Needed for viscous solvers. For incompressible solvers the Reynolds length will always be 1.0 \n DEFAULT: 0.0 \ingroup Config */ addDoubleOption("REYNOLDS_NUMBER", Reynolds, 0.0); /*!\brief REYNOLDS_LENGTH \n DESCRIPTION: Reynolds length (1 m by default). Used for compressible solver: incompressible solver will use 1.0. \ingroup Config */ addDoubleOption("REYNOLDS_LENGTH", Length_Reynolds, 1.0); /*!\brief PRANDTL_LAM \n DESCRIPTION: Laminar Prandtl number (0.72 (air), only for compressible flows) \n DEFAULT: 0.72 \ingroup Config*/ - addDoubleOption("PRANDTL_LAM", Prandtl_Lam, 0.72); + addDoubleListOption("PRANDTL_LAM", nPrandtl_Lam, Prandtl_Lam); /*!\brief PRANDTL_TURB \n DESCRIPTION: Turbulent Prandtl number (0.9 (air), only for compressible flows) \n DEFAULT 0.90 \ingroup Config*/ - addDoubleOption("PRANDTL_TURB", Prandtl_Turb, 0.90); + addDoubleListOption("PRANDTL_TURB", nPrandtl_Turb, Prandtl_Turb); /*--- Options related to wall models. ---*/ @@ -1315,7 +1324,7 @@ void CConfig::SetConfig_Options() { //addDoubleOption("SCALAR_INIT", Scalar_Init, 0.0); addDoubleListOption("SCALAR_INIT", nScalar_Init, Scalar_Init); - + /*!\brief SCALAR_CLIPPING \n DESCRIPTION: Activate clipping for scalar transport equations \ingroup Config*/ addBoolOption("SCALAR_CLIPPING", Scalar_Clipping, false); @@ -1325,7 +1334,7 @@ void CConfig::SetConfig_Options() { /*!\brief SCALAR_CLIPPING_MAX \n DESCRIPTION: Maximum value for scalar clipping \ingroup Config*/ addDoubleListOption("SCALAR_CLIPPING_MAX", nScalar_Clipping_Max, Scalar_Clipping_Max); - + /*!\brief SCALAR_CLIPPING_MIN \n DESCRIPTION: Minimum value for scalar clipping \ingroup Config*/ addDoubleListOption("SCALAR_CLIPPING_MIN", nScalar_Clipping_Min, Scalar_Clipping_Min); @@ -1345,13 +1354,13 @@ void CConfig::SetConfig_Options() { /*!\brief FLAME_THICKNESS \n DESCRIPTION: Thickness for flame initialization using the flamelet model \ingroup Config*/ addDoubleOption("FLAME_THICKNESS", flame_thickness, 0.5e-3); - + /*!\brief FLAME_NORMAL \n DESCRIPTION: Normal for flame initialization using the flamelet model \ingroup Config*/ flame_normal[0] = 1.0; flame_normal[1] = 0.0; flame_normal[2] = 0.0; addDoubleArrayOption("FLAME_NORMAL", 3, flame_normal); - + addDoubleOption("BURNT_THICKNESS", burnt_thickness, 1); /*!\brief INC_INLET_DAMPING \n DESCRIPTION: Damping factor applied to the iterative updates to the velocity at a pressure inlet in incompressible flow (0.1 by default). \ingroup Config*/ @@ -1941,7 +1950,7 @@ void CConfig::SetConfig_Options() { /*!\brief CONV_NUM_METHOD_SCALAR * \n DESCRIPTION: Convective numerical method \ingroup Config*/ addConvectOption("CONV_NUM_METHOD_SCALAR", Kind_ConvNumScheme_Scalar, Kind_Centered_Scalar, Kind_Upwind_Scalar); - + /*!\brief MUSCL_FLOW \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/ addBoolOption("MUSCL_ADJTURB", MUSCL_AdjTurb, false); /*!\brief SLOPE_LIMITER_ADJTURB @@ -2878,25 +2887,25 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Size of the edge groups colored for thread parallel edge loops (0 forces the reducer strategy). */ addUnsignedLongOption("EDGE_COLORING_GROUP_SIZE", edgeColorGroupSize, 512); - + /*--- options that are used for libROM ---*/ /*!\par CONFIG_CATEGORY:libROM options \ingroup Config*/ - + /*!\brief SAVE_LIBROM \n DESCRIPTION: Flag for saving data with libROM. */ addBoolOption("SAVE_LIBROM", libROM, false); - + /*!\brief LIBROM_BASE_FILENAME \n DESCRIPTION: Output base file name for libROM \ingroup Config*/ addStringOption("LIBROM_BASE_FILENAME", libROMbase_FileName, string("su2")); - + /*!\brief BASIS_GENERATION \n DESCRIPTION: Flag for saving data with libROM. */ addEnumOption("BASIS_GENERATION", POD_Basis_Gen, POD_Map, POD_KIND::STATIC); - + /*!\brief MAX_BASIS_DIM \n DESCRIPTION: Maximum number of basis vectors.*/ addUnsignedShortOption("MAX_BASIS_DIM", maxBasisDim, 100); - + /*!\brief MAX_BASIS_DIM \n DESCRIPTION: Maximum number of basis vectors.*/ addUnsignedShortOption("ROM_SAVE_FREQ", rom_save_freq, 1); - + /* END_CONFIG_OPTIONS */ } @@ -3319,7 +3328,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i (Kind_FluidModel == INC_IDEAL_GAS) || (Kind_FluidModel == INC_IDEAL_GAS_POLY) || (Kind_FluidModel == CONSTANT_DENSITY) || - (Kind_FluidModel == FLAMELET_FLUID_MODEL)); + (Kind_FluidModel == FLAMELET_FLUID_MODEL) || + (Kind_FluidModel == MIXTURE_FLUID_MODEL)); bool noneq_gas = ((Kind_FluidModel == MUTATIONPP) || (Kind_FluidModel == SU2_NONEQ)); bool standard_air = ((Kind_FluidModel == STANDARD_AIR)); @@ -3702,20 +3712,145 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + /*--- Set default values for various fluid properties. ---*/ + + static const su2double Molecular_Weight_Default = 28.96; + static const su2double Mu_Constant_Default = 1.716E-5; + static const su2double Mu_Ref_Default = Mu_Constant_Default; + static const su2double Mu_Temperature_Ref_Default = 273.15; + static const su2double Mu_S_Default = 110.4; + static const su2double Specific_Heat_Cp_Default = 1004.703; + static const su2double Thermal_Conductivity_Constant_Default = 2.57E-2; + static const su2double Prandtl_Lam_Default = 0.72; + static const su2double Prandtl_Turb_Default = 0.9; + + if (Molecular_Weight == nullptr){ + Molecular_Weight = new su2double[1]; + Molecular_Weight[0] = Molecular_Weight_Default; + nMolecular_Weight = 1; + } + + if (Mu_Constant == nullptr){ + Mu_Constant = new su2double[1]; + Mu_Constant[0] = Mu_Constant_Default; + nMu_Constant = 1; + } + + if (Mu_Ref == nullptr && Mu_Temperature_Ref == nullptr && Mu_S == nullptr){ + Mu_Ref = new su2double[1]; + Mu_Temperature_Ref = new su2double[1]; + Mu_S = new su2double[1]; + Mu_Ref[0] = Mu_Ref_Default; + Mu_Temperature_Ref[0] = Mu_Temperature_Ref_Default; + Mu_S[0] = Mu_S_Default; + nMu_Ref = 1; + nMu_Temperature_Ref = 1; + nMu_S = 1; + } + + if (Specific_Heat_Cp == nullptr){ + Specific_Heat_Cp = new su2double[1]; + Specific_Heat_Cp[0] = Specific_Heat_Cp_Default; + nSpecific_Heat_Cp = 1; + } + + if (Thermal_Conductivity_Constant == nullptr){ + Thermal_Conductivity_Constant = new su2double[1]; + Thermal_Conductivity_Constant[0] = Thermal_Conductivity_Constant_Default; + nThermal_Conductiviy_Constant = 1; + } + + if (Prandtl_Lam == nullptr){ + Prandtl_Lam = new su2double[1]; + Prandtl_Lam[0] = Prandtl_Lam_Default; + nPrandtl_Lam = 1; + } + + if (Prandtl_Turb == nullptr){ + Prandtl_Turb = new su2double[1]; + Prandtl_Turb[0] = Prandtl_Turb_Default; + nPrandtl_Turb = 1; + } + + /*--- Check whether inputs for MIXTURE_FLUID_MODEL are correctly specified. ---*/ + unsigned short n_species = 1; //TODO TK:: make it static? + + if (Kind_FluidModel == MIXTURE_FLUID_MODEL) + n_species = nScalar_Init + 1; + + /*--- Check whether the number of entries of each specified fluid property equals the number of transported scalar equations solved + 1. + * nMolecular_Weight and nSpecific_Heat_Cp are used because they are required for the fluid mixing models. + * Cp is required in case of MIXTURE_FLUID_MODEL because the energy equation needs to be active. --- */ + if ((nMolecular_Weight != n_species) || (nSpecific_Heat_Cp != n_species)) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for MOLECULAR_WEIGHT and SPECIFIC_HEAT_CP,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + + switch (Kind_ViscosityModel) { + case VISCOSITYMODEL::CONSTANT: + if (nMu_Constant != n_species) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for MU_CONSTANT,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + break; + case VISCOSITYMODEL::SUTHERLAND: + if ((nMu_Ref != n_species) || (nMu_Temperature_Ref != n_species) || (nMu_S != n_species)) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for MU_REF, MU_T_REF and SUTHERLAND_CONSTANT,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + break; + default: + if (n_species != 1) SU2_MPI::Error("Viscosity model not available.", CURRENT_FUNCTION); + break; + } + + switch (Kind_ConductivityModel) { + case CONDUCTIVITYMODEL::CONSTANT: + if (Kind_ConductivityModel_Turb == CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL) { + if ((nThermal_Conductiviy_Constant != n_species) || (nPrandtl_Turb != n_species)) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for THERMAL_CONDUCTIVITY_CONSTANT and PRANDTL_TURB,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + } else { + if (nThermal_Conductiviy_Constant != n_species) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for THERMAL_CONDUCTIVITY_CONSTANT,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + } + break; + case CONDUCTIVITYMODEL::CONSTANT_PRANDTL: + if (Kind_ConductivityModel_Turb == CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL) { + if ((nPrandtl_Lam != n_species) || (nPrandtl_Turb != n_species)) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for PRANDTL_LAM and PRANDTL_TURB,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + } else { + if (nPrandtl_Lam != n_species) { + SU2_MPI::Error("The use of MIXTURE_FLUID_MODEL requires the number of entries for PRANDTL_LAM,\n" + "to be equal to the number of entries of SCALAR_INIT + 1", CURRENT_FUNCTION); + } + } + break; + default: + if (n_species != 1) SU2_MPI::Error("Conductivity model not available.", CURRENT_FUNCTION); + break; + } + /*--- Overrule the default values for viscosity if the US measurement system is used. ---*/ if (SystemMeasurements == US) { - /* Correct the viscosities, if they contain the default SI values. */ - if(fabs(Mu_Constant-1.716E-5) < 1.0E-15) Mu_Constant /= 47.88025898; - if(fabs(Mu_Ref-1.716E-5) < 1.0E-15) Mu_Ref /= 47.88025898; + for(unsigned short iVar = 0; iVar < n_species; iVar++){ + if(fabs(Mu_Constant[iVar]-Mu_Constant_Default) < 1.0E-15) Mu_Constant[iVar] /= 47.88025898; + if(fabs(Mu_Ref[iVar]-Mu_Constant_Default) < 1.0E-15) Mu_Ref[iVar] /= 47.88025898; - /* Correct the values with temperature dimension, if they contain the default SI values. */ - if(fabs(Mu_Temperature_Ref-273.15) < 1.0E-8) Mu_Temperature_Ref *= 1.8; - if(fabs(Mu_S-110.4) < 1.0E-8) Mu_S *= 1.8; + /* Correct the values with temperature dimension, if they contain the default SI values. */ + if(fabs(Mu_Temperature_Ref[iVar]-Mu_Temperature_Ref_Default) < 1.0E-8) Mu_Temperature_Ref[iVar] *= 1.8; + if(fabs(Mu_S[iVar]-Mu_S_Default) < 1.0E-8) Mu_S[iVar] *= 1.8; - /* Correct the thermal conductivity, if it contains the default SI value. */ - if(fabs(Thermal_Conductivity_Constant-0.0257) < 1.0E-10) Thermal_Conductivity_Constant *= 0.577789317; + /* Correct the thermal conductivity, if it contains the default SI value. */ + if(fabs(Thermal_Conductivity_Constant[iVar]-Thermal_Conductivity_Constant_Default) < 1.0E-10) Thermal_Conductivity_Constant[iVar] *= 0.577789317; + } } /*--- Check for Measurement System ---*/ @@ -4706,7 +4841,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } if (Kind_DensityModel == INC_DENSITYMODEL::VARIABLE) { - if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLAMELET_FLUID_MODEL) { + if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLAMELET_FLUID_MODEL && Kind_FluidModel != MIXTURE_FLUID_MODEL) { SU2_MPI::Error("Variable density incompressible solver limited to ideal gases.\n Check the fluid model options (use INC_IDEAL_GAS, INC_IDEAL_GAS_POLY).", CURRENT_FUNCTION); } } @@ -4719,8 +4854,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) { if (Kind_ViscosityModel == VISCOSITYMODEL::SUTHERLAND) { - if ((Kind_FluidModel != INC_IDEAL_GAS) && (Kind_FluidModel != INC_IDEAL_GAS_POLY) && (Kind_FluidModel != FLAMELET_FLUID_MODEL)) { - SU2_MPI::Error("Sutherland's law only valid for ideal gases in incompressible flows.\n Must use VISCOSITY_MODEL=CONSTANT_VISCOSITY and set viscosity with\n MU_CONSTANT, or use DENSITY_MODEL= VARIABLE with FLUID_MODEL= INC_IDEAL_GAS or INC_IDEAL_GAS_POLY, or FLAMELET_FLUID_MODEL for VISCOSITY_MODEL=SUTHERLAND.\n NOTE: FREESTREAM_VISCOSITY is no longer used for incompressible flows!", CURRENT_FUNCTION); + if ((Kind_FluidModel != INC_IDEAL_GAS) && (Kind_FluidModel != INC_IDEAL_GAS_POLY) && (Kind_FluidModel != FLAMELET_FLUID_MODEL) && (Kind_FluidModel != MIXTURE_FLUID_MODEL)) { + SU2_MPI::Error("Sutherland's law only valid for ideal gases in incompressible flows.\n Must use VISCOSITY_MODEL=CONSTANT_VISCOSITY and set viscosity with\n MU_CONSTANT, or use DENSITY_MODEL= VARIABLE with FLUID_MODEL= INC_IDEAL_GAS or INC_IDEAL_GAS_POLY, or FLAMELET_FLUID_MODEL or MIXTURE_FLUID_MODEL for VISCOSITY_MODEL=SUTHERLAND.\n NOTE: FREESTREAM_VISCOSITY is no longer used for incompressible flows!", CURRENT_FUNCTION); } } } @@ -4808,15 +4943,15 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("Must list two markers for the pressure drop objective function.\n Expected format: MARKER_ANALYZE= (outlet_name, inlet_name).", CURRENT_FUNCTION); } } - + /*--- Disable any scalar model until they are implemented. ---*/ - + if ((Kind_Scalar_Model != NO_SCALAR_MODEL) && (Kind_Scalar_Model != PASSIVE_SCALAR) && (Kind_Scalar_Model != PROGRESS_VARIABLE)) { SU2_MPI::Error(string("Selected scalar model not yet implemented.") , CURRENT_FUNCTION); } - + /*--- Check feasibility for Streamwise Periodic flow ---*/ if (Kind_Streamwise_Periodic != ENUM_STREAMWISE_PERIODIC::NONE) { if (Kind_Regime != ENUM_REGIME::INCOMPRESSIBLE) @@ -5314,6 +5449,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Surface_IDR = new su2double [nMarker_Analyze] (); Surface_CO = new su2double [nMarker_Analyze] (); Surface_NOx = new su2double [nMarker_Analyze] (); + Surface_CH4 = new su2double [nMarker_Analyze] (); //Surface_Scalar = new su2double*[nMarker_Analyze] (); //for (int i_scalar=0; i_scalar < nMarker_Analyze; ++i_scalar) @@ -7851,6 +7987,7 @@ CConfig::~CConfig(void) { delete[] Surface_IDR; delete[] Surface_CO; delete[] Surface_NOx; + delete[] Surface_CH4; //delete[] Surface_Scalar; delete[] Inlet_Ttotal; @@ -8026,6 +8163,16 @@ CConfig::~CConfig(void) { delete [] Scalar_Clipping_Max; delete [] Scalar_Init; + delete [] Molecular_Weight; + delete [] Mu_Constant; + delete [] Mu_Ref; + delete [] Mu_Temperature_Ref; + delete [] Mu_S; + delete [] Specific_Heat_Cp; + delete [] Thermal_Conductivity_Constant; + delete [] Prandtl_Lam; + delete [] Prandtl_Turb; + } string CConfig::GetFilename(string filename, string ext, int Iter) const { diff --git a/SU2_CFD/include/fluid/CConstantDiffusivity.hpp b/SU2_CFD/include/fluid/CConstantDiffusivity.hpp index 221b284771eb..e0e8be937e74 100644 --- a/SU2_CFD/include/fluid/CConstantDiffusivity.hpp +++ b/SU2_CFD/include/fluid/CConstantDiffusivity.hpp @@ -50,7 +50,8 @@ class CConstantDiffusivity final : public CDiffusivityModel { su2double rho, su2double mu_lam, su2double mu_turb, - su2double cp) override { } + su2double cp, + su2double kt) override { } private: su2double diff_{0.0}; diff --git a/SU2_CFD/include/fluid/CConstantSchmidt.hpp b/SU2_CFD/include/fluid/CConstantSchmidt.hpp index fc5b07d7ad2d..a18c1d1562f3 100644 --- a/SU2_CFD/include/fluid/CConstantSchmidt.hpp +++ b/SU2_CFD/include/fluid/CConstantSchmidt.hpp @@ -62,7 +62,7 @@ class CConstantSchmidt final: public CDiffusivityModel { /*! * \brief Set diffusivity. */ - void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override { + void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override { diff_ = mu_lam / (rho*sc_lam_); } diff --git a/SU2_CFD/include/fluid/CConstantSchmidtRANS.hpp b/SU2_CFD/include/fluid/CConstantSchmidtRANS.hpp index 1383fd61ec8e..67cd39fa3d95 100644 --- a/SU2_CFD/include/fluid/CConstantSchmidtRANS.hpp +++ b/SU2_CFD/include/fluid/CConstantSchmidtRANS.hpp @@ -57,7 +57,7 @@ class CConstantSchmidtRANS : public CDiffusivityModel { /*! * \brief Set mass diffusivity. */ - void SetDiffusivity(su2double t, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override { + void SetDiffusivity(su2double t, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override { diff_ = diff_lam_const_ + mu_turb/(rho*sc_turb_); } diff --git a/SU2_CFD/include/fluid/CDiffusivityModel.hpp b/SU2_CFD/include/fluid/CDiffusivityModel.hpp index d1e1523874a3..03de839b40c8 100644 --- a/SU2_CFD/include/fluid/CDiffusivityModel.hpp +++ b/SU2_CFD/include/fluid/CDiffusivityModel.hpp @@ -47,5 +47,5 @@ class CDiffusivityModel { /*! * \brief Set mass diffusivity */ - virtual void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) = 0; + virtual void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) = 0; }; diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index f440eeec28ae..66cd0886283d 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -125,6 +125,36 @@ class CFluidModel { */ su2double GetCv() const { return Cv; } + /*! + * \brief Get fluid mean specific heat capacity at constant pressure. + */ + su2double ComputeMeanSpecificHeatCp(const std::vector specific_heat_cp, const su2double * const val_scalars) { + su2double val_scalars_sum = 0.0; + unsigned short n_scalars = specific_heat_cp.size() - 1; + + for (int i_scalar = 0; i_scalar < n_scalars; i_scalar++){ + Cp += specific_heat_cp[i_scalar] * val_scalars[i_scalar]; + val_scalars_sum += val_scalars[i_scalar]; + } + return Cp += specific_heat_cp[n_scalars]*(1 - val_scalars_sum); + } + + /*! + * \brief Get fluid mean molecular weight. + */ + static su2double ComputeMeanMolecularWeight(const std::vector molar_masses, const su2double * const val_scalars) { + su2double OneOverMeanMolecularWeight = 0.0; + su2double val_scalars_sum = 0.0; + unsigned short n_scalars = molar_masses.size() - 1; + + for (int i_scalar = 0; i_scalar < n_scalars; i_scalar++){ + OneOverMeanMolecularWeight += val_scalars[i_scalar]/(molar_masses[i_scalar]/1000); + val_scalars_sum += val_scalars[i_scalar]; + } + OneOverMeanMolecularWeight += (1 - val_scalars_sum)/(molar_masses[n_scalars]/1000); + return 1/OneOverMeanMolecularWeight; + } + /*! * \brief Get the source term of the transported scalar */ @@ -164,11 +194,11 @@ class CFluidModel { * \brief Get fluid dynamic viscosity. */ virtual inline su2double GetLaminarViscosity() { - LaminarViscosity->SetViscosity(Temperature, Density); - Mu = LaminarViscosity->GetViscosity(); - LaminarViscosity->SetDerViscosity(Temperature, Density); - dmudrho_T = LaminarViscosity->Getdmudrho_T(); - dmudT_rho = LaminarViscosity->GetdmudT_rho(); + LaminarViscosity->SetViscosity(Temperature, Density); + Mu = LaminarViscosity->GetViscosity(); + LaminarViscosity->SetDerViscosity(Temperature, Density); + dmudrho_T = LaminarViscosity->Getdmudrho_T(); + dmudT_rho = LaminarViscosity->GetdmudT_rho(); return Mu; } @@ -188,7 +218,7 @@ class CFluidModel { * \brief Get fluid mass diffusivity. */ virtual inline su2double GetMassDiffusivity() { - MassDiffusivity->SetDiffusivity(Temperature, Density, Mu, Mu_Turb, Cp); + MassDiffusivity->SetDiffusivity(Temperature, Density, Mu, Mu_Turb, Cp, Kt); mass_diffusivity = MassDiffusivity->GetDiffusivity(); //ThermalConductivity->SetDerConductivity(Temperature, Density, dmudrho_T, dmudT_rho, Cp); //dktdrho_T = ThermalConductivity->Getdktdrho_T(); @@ -264,12 +294,12 @@ class CFluidModel { /*! * \brief Set viscosity model. */ - void SetLaminarViscosityModel(const CConfig* config); + virtual void SetLaminarViscosityModel(const CConfig* config); /*! * \brief Set thermal conductivity model. */ - void SetThermalConductivityModel(const CConfig* config); + virtual void SetThermalConductivityModel(const CConfig* config); /*! * \brief Set mass diffusivity model. diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp new file mode 100644 index 000000000000..096e2d6df9ed --- /dev/null +++ b/SU2_CFD/include/fluid/CFluidScalar.hpp @@ -0,0 +1,86 @@ +#pragma once + +#include +#include + +#include "CFluidModel.hpp" + +class CFluidScalar final : public CFluidModel { +private: + unsigned short n_scalars = 0; /*!< \brief Number of transported scalars. */ + unsigned short n_species_mixture = 0; /*!< \brief Number of species in mixture. */ + su2double Gas_Constant = 0.0; /*!< \brief Specific gas constant. */ + su2double Gamma = 0.0; /*!< \brief Ratio of specific heats of the gas. */ + su2double Pressure_Thermodynamic = 0.0; /*!< \brief Constant pressure thermodynamic. */ + + bool wilke; + bool davidson; + + std::vector massFractions; /*!< \brief Mass fractions of all species. */ + std::vector moleFractions; /*!< \brief Mole fractions of all species. */ + std::vector molarMasses; /*!< \brief Molar masses of all species. */ + std::vector laminarViscosity; /*!< \brief Laminar viscosity of all species. */ + std::vector specificHeat; /*!< \brief Specific heat of all species. */ + std::vector laminarthermalConductivity; /*!< \brief Laminar thermal conductivity of all species. */ + + static const int ARRAYSIZE = 100; + std::unique_ptr LaminarViscosityPointers[ARRAYSIZE]; + std::unique_ptr ThermalConductivityPointers[ARRAYSIZE]; + + /*! + * \brief Convert mass fractions to mole fractions. + * \param[in] val_scalars - Scalar mass fraction. + */ + std::vector& massToMoleFractions(const su2double * const val_scalars); + + /*! + * \brief Wilke mixing law for mixture viscosity. + * \param[in] val_scalars - Scalar mass fraction. + */ + su2double wilkeViscosity(const su2double * const val_scalars); + + /*! + * \brief Davidson mixing law for mixture viscosity. + * \param[in] val_scalars - Scalar mass fraction. + */ + su2double davidsonViscosity(const su2double * const val_scalars); + + /*! + * \brief Wilke mixing law for mixture thermal conductivity. + * \param[in] val_scalars - Scalar mass fraction. + */ + su2double wilkeConductivity(const su2double * const val_scalars); + +public: + /*! + * \brief Constructor of the class. + */ + CFluidScalar(CConfig *config, const su2double value_pressure_operating); + + /*! + * \brief Set viscosity model. + */ + void SetLaminarViscosityModel(const CConfig* config) override; + + /*! + * \brief Set thermal conductivity model. + */ + void SetThermalConductivityModel(const CConfig* config) override; + + /*! + * \brief Set the Dimensionless State using Temperature. + * \param[in] val_temperature - Temperature value at the point. + * \param[in] val_scalars - Scalar mass fraction. + */ + unsigned long SetTDState_T(const su2double val_temperature, su2double * const val_scalars) override; + + /*! + * \brief Get fluid dynamic viscosity. + */ + inline su2double GetLaminarViscosity() override {return Mu; } + + /*! + * \brief Get fluid thermal conductivity. + */ + inline su2double GetThermalConductivity() override { return Kt; } +}; diff --git a/SU2_CFD/include/fluid/CUnityLewisDiffusivity.hpp b/SU2_CFD/include/fluid/CUnityLewisDiffusivity.hpp new file mode 100644 index 000000000000..0d4e06307aa4 --- /dev/null +++ b/SU2_CFD/include/fluid/CUnityLewisDiffusivity.hpp @@ -0,0 +1,57 @@ +/*! + * \file CUnityLewisDiffusivity.hpp + * \brief Defines unity Lewis mass diffusivity. + * \author M.Heimgartner + * \version 7.0.6 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CDiffusivityModel.hpp" + +/*! + * \class CUnityLewisDiffusivity + * \brief Defines a unity Lewis mass diffusivity model for species equations. + * \author M.Heimgartner + */ +class CUnityLewisDiffusivity final : public CDiffusivityModel { +public: + /*! + * \brief Constructor of the class. + */ + CUnityLewisDiffusivity() {} + + su2double GetDiffusivity() const override {return diff_;} + + /*! + * \brief Set diffusivity. + */ + void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override { + diff_ = kt / (Lewis * rho * cp); + } + + private: + su2double diff_{0.0}; + su2double kt_{0.0}; + su2double Lewis{1}; +}; diff --git a/SU2_CFD/include/numerics/flamelet/scalarLegacy_sources.hpp b/SU2_CFD/include/numerics/flamelet/scalarLegacy_sources.hpp index 8fe85fe19985..8235a16aa1e0 100644 --- a/SU2_CFD/include/numerics/flamelet/scalarLegacy_sources.hpp +++ b/SU2_CFD/include/numerics/flamelet/scalarLegacy_sources.hpp @@ -104,7 +104,7 @@ class CSourcePieceWise_transportedScalar_general final : public CNumerics { if (flame) Residual[iVar] += yinv*Volume*(Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1]; else - Residual[iVar] += yinv*Volume*Density_i*(Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1]; + Residual[iVar] += yinv*Volume*(Density_i*Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1]; } } diff --git a/SU2_CFD/include/solvers/CPassiveScalarSolver.hpp b/SU2_CFD/include/solvers/CPassiveScalarSolver.hpp index 7e959e59d95d..da1eb52820ea 100644 --- a/SU2_CFD/include/solvers/CPassiveScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CPassiveScalarSolver.hpp @@ -88,6 +88,15 @@ class CPassiveScalarSolver: public CScalarLegacySolver { CConfig *config, unsigned short iMesh); + /*! + * \brief Compute the primitive variables (diffusivities) + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] Output - boolean to determine whether to print output. + * \return - The number of non-physical points. + */ + unsigned long SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output); + /*! * \brief Set the initial condition for the scalar transport problem. diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 66d533ae24a0..c02f349c0774 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -98,7 +98,7 @@ unsigned long CFluidFlamelet::SetScalarSources(su2double *val_scalars){ delete [] table_sources; - /*--- not used at the moment ---*/ + /*--- not used at the moment ---*/ return exit_code; } diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp index 7b4f64dcf84d..2fb9de5c5273 100644 --- a/SU2_CFD/src/fluid/CFluidModel.cpp +++ b/SU2_CFD/src/fluid/CFluidModel.cpp @@ -37,9 +37,11 @@ #include "../../include/fluid/CPolynomialConductivityRANS.hpp" #include "../../include/fluid/CPolynomialViscosity.hpp" #include "../../include/fluid/CSutherland.hpp" +#include "../../include/fluid/CFluidScalar.hpp" #include "../../include/fluid/CConstantDiffusivity.hpp" #include "../../include/fluid/CConstantSchmidtRANS.hpp" #include "../../include/fluid/CConstantSchmidt.hpp" +#include "../../include/fluid/CUnityLewisDiffusivity.hpp" void CFluidModel::SetLaminarViscosityModel(const CConfig* config) { switch (config->GetKind_ViscosityModel()) { @@ -114,6 +116,9 @@ void CFluidModel::SetMassDiffusivityModel (const CConfig* config) { case DIFFUSIVITYMODEL::FLAMELET: /* do nothing. Diffusivity is obtained from the table and set in setTDState_T */ break; + case DIFFUSIVITYMODEL::UNITY_LEWIS: + MassDiffusivity = unique_ptr(new CUnityLewisDiffusivity()); + break; default: SU2_MPI::Error("Diffusivity model not available.", CURRENT_FUNCTION); break; diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp new file mode 100644 index 000000000000..761c4d162c27 --- /dev/null +++ b/SU2_CFD/src/fluid/CFluidScalar.cpp @@ -0,0 +1,249 @@ +#include +#include +#include +#include +#include + +#include "../../include/fluid/CFluidScalar.hpp" +#include "../../include/fluid/CConstantViscosity.hpp" +#include "../../include/fluid/CSutherland.hpp" +#include "../../include/fluid/CPolynomialViscosity.hpp" +#include "../../include/fluid/CConstantConductivity.hpp" +#include "../../include/fluid/CConstantConductivityRANS.hpp" +#include "../../include/fluid/CConstantPrandtl.hpp" +#include "../../include/fluid/CConstantPrandtlRANS.hpp" +#include "../../include/fluid/CPolynomialConductivity.hpp" +#include "../../include/fluid/CPolynomialConductivityRANS.hpp" +#include "../../include/fluid/CIncIdealGas.hpp" + +CFluidScalar::CFluidScalar(CConfig *config, const su2double value_pressure_operating) : CFluidModel() { + n_scalars = config->GetNScalarsInit(); + config->SetNScalarsInit(n_scalars); + n_species_mixture = n_scalars + 1; + + specificHeat.resize(n_species_mixture); + molarMasses.resize(n_species_mixture); + massFractions.resize(n_species_mixture); + moleFractions.resize(n_species_mixture); + laminarViscosity.resize(n_species_mixture); + laminarthermalConductivity.resize(n_species_mixture); + + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + molarMasses.at(iVar) = config->GetMolecular_Weight(iVar); + specificHeat.at(iVar) = config->GetSpecific_Heat_Cp(iVar); + } + + wilke = false; + davidson = true; + + Pressure_Thermodynamic = value_pressure_operating; + Gas_Constant = config->GetGas_Constant(); + Gamma = 1.0; + + SetLaminarViscosityModel(config); + SetThermalConductivityModel(config); +} + +void CFluidScalar::SetLaminarViscosityModel(const CConfig* config) { + switch (config->GetKind_ViscosityModel()) { + case VISCOSITYMODEL::CONSTANT: + /* Build a list of LaminarViscosity pointers to be used in e.g. wilkeViscosity to get the species viscosities. */ + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + LaminarViscosityPointers[iVar] = std::unique_ptr(new CConstantViscosity(config->GetMu_Constant(iVar))); + } + break; + case VISCOSITYMODEL::SUTHERLAND: + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + LaminarViscosityPointers[iVar] = std::unique_ptr(new CSutherland(config->GetMu_Ref(iVar), config->GetMu_Temperature_Ref(iVar), config->GetMu_S(iVar))); + } + break; + case VISCOSITYMODEL::POLYNOMIAL: + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + LaminarViscosityPointers[iVar] = std::unique_ptr>( + new CPolynomialViscosity(config->GetMu_PolyCoeffND())); + } + break; + case VISCOSITYMODEL::FLAMELET: + /* Do nothing. Viscosity is obtained from the table and set in SetTDState_T */ + break; + default: + SU2_MPI::Error("Viscosity model not available.", CURRENT_FUNCTION); + break; + } +} + +void CFluidScalar::SetThermalConductivityModel(const CConfig* config) { + switch (config->GetKind_ConductivityModel()) { + case CONDUCTIVITYMODEL::CONSTANT: + for(int iVar = 0; iVar < n_species_mixture; iVar++){ + if (config->GetKind_ConductivityModel_Turb() == CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL) { + ThermalConductivityPointers[iVar] = std::unique_ptr( + new CConstantConductivityRANS(config->GetThermal_Conductivity_Constant(iVar), config->GetPrandtl_Turb(iVar))); + } else { + ThermalConductivityPointers[iVar] = std::unique_ptr(new CConstantConductivity(config->GetThermal_Conductivity_Constant(iVar))); + } + } + break; + case CONDUCTIVITYMODEL::CONSTANT_PRANDTL: + for(int iVar = 0; iVar < n_species_mixture; iVar++){ + if (config->GetKind_ConductivityModel_Turb() == CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL) { + ThermalConductivityPointers[iVar] = std::unique_ptr( + new CConstantPrandtlRANS(config->GetPrandtl_Lam(iVar), config->GetPrandtl_Turb(iVar))); + } else { + ThermalConductivityPointers[iVar] = std::unique_ptr(new CConstantPrandtl(config->GetPrandtl_Lam(iVar))); + } + } + break; + case CONDUCTIVITYMODEL::POLYNOMIAL: + for(int iVar = 0; iVar < n_species_mixture; iVar++){ + if (config->GetKind_ConductivityModel_Turb() == CONDUCTIVITYMODEL_TURB::CONSTANT_PRANDTL) { + ThermalConductivity = std::unique_ptr>( + new CPolynomialConductivityRANS(config->GetKt_PolyCoeffND(), config->GetPrandtl_Turb())); + } else { + ThermalConductivity = std::unique_ptr>( + new CPolynomialConductivity(config->GetKt_PolyCoeffND())); + } + } + break; + case CONDUCTIVITYMODEL::FLAMELET: + /* Do nothing. Conductivity is obtained from the table and set in SetTDState_T */ + break; + default: + SU2_MPI::Error("Conductivity model not available.", CURRENT_FUNCTION); + break; + } +} + +std::vector& CFluidScalar::massToMoleFractions(const su2double * const val_scalars){ + su2double mixtureMolarMass {0.0}; + su2double val_scalars_sum {0.0}; + + for(int i_scalar = 0; i_scalar < n_scalars; i_scalar++){ + massFractions[i_scalar] = val_scalars[i_scalar]; + val_scalars_sum += val_scalars[i_scalar]; + } + massFractions[n_scalars] = 1 - val_scalars_sum; + + for(int iVar = 0; iVar < n_species_mixture; iVar++){ + mixtureMolarMass += massFractions[iVar] / molarMasses[iVar]; + } + + for(int iVar = 0; iVar < n_species_mixture; iVar++){ + moleFractions.at(iVar) = (massFractions[iVar] / molarMasses[iVar]) / mixtureMolarMass; + } + + return moleFractions; +} + +su2double CFluidScalar::wilkeViscosity(const su2double * const val_scalars){ + std::vector phi; + std::vector wilkeNumerator; + std::vector wilkeDenumeratorSum; + su2double wilkeDenumerator = 0.0; + wilkeDenumeratorSum.clear(); + wilkeNumerator.clear(); + su2double viscosityMixture = 0.0; + + /* Fill laminarViscosity with n_species_mixture viscosity values. */ + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + LaminarViscosityPointers[iVar]->SetViscosity(Temperature, Density); + laminarViscosity.at(iVar) = LaminarViscosityPointers[iVar]->GetViscosity(); + } + + for(int i = 0; i < n_species_mixture; i++){ + for(int j = 0; j < n_species_mixture; j++){ + phi.push_back(pow(1 + sqrt(laminarViscosity[i] / laminarViscosity[j]) * pow(molarMasses[j] / molarMasses[i], 0.25), 2) / sqrt(8 * (1 + molarMasses[i] / molarMasses[j]))); + wilkeDenumerator += moleFractions[j] * phi[j]; + } + wilkeDenumeratorSum.push_back(wilkeDenumerator); + wilkeDenumerator = 0.0; + phi.clear(); + wilkeNumerator.push_back(moleFractions[i] * laminarViscosity[i]); + viscosityMixture += wilkeNumerator[i] / wilkeDenumeratorSum[i]; + } + return viscosityMixture; +} + +su2double CFluidScalar::davidsonViscosity(const su2double * const val_scalars){ + su2double viscosityMixture = 0.0; + su2double fluidity = 0.0; + su2double E = 0.0; + su2double mixtureFractionDenumerator = 0.0; + const su2double A = 0.375; + std::vector mixtureFractions; + mixtureFractions.clear(); + + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + LaminarViscosityPointers[iVar]->SetViscosity(Temperature, Density); + laminarViscosity.at(iVar) = LaminarViscosityPointers[iVar]->GetViscosity(); + } + + for(int i = 0; i < n_species_mixture; i++){ + mixtureFractionDenumerator += moleFractions[i] * sqrt(molarMasses[i]); + } + + for(int j = 0; j < n_species_mixture; j++){ + mixtureFractions.push_back((moleFractions[j] * sqrt(molarMasses[j])) / mixtureFractionDenumerator); + } + + for(int i = 0; i < n_species_mixture; i++){ + for(int j = 0; j < n_species_mixture; j++){ + E = (2*sqrt(molarMasses[i]) * sqrt(molarMasses[j])) / (molarMasses[i] + molarMasses[j]); + fluidity += ((mixtureFractions[i] * mixtureFractions[j]) / (sqrt(laminarViscosity[i]) * sqrt(laminarViscosity[j]))) * pow(E, A); + } + } + return viscosityMixture = 1 / fluidity; +} + +su2double CFluidScalar::wilkeConductivity(const su2double * const val_scalars){ + std::vector phi; + std::vector wilkeNumerator; + std::vector wilkeDenumeratorSum; + su2double wilkeDenumerator = 0.0; + su2double conductivityMixture = 0.0; + wilkeDenumeratorSum.clear(); + wilkeNumerator.clear(); + + for (int iVar = 0; iVar < n_species_mixture; iVar++){ + ThermalConductivityPointers[iVar]->SetConductivity(Temperature, Density, Mu, Mu_Turb, Cp); + laminarthermalConductivity.at(iVar) = ThermalConductivityPointers[iVar]->GetConductivity(); + } + + for(int i = 0; i < n_species_mixture; i++){ + for(int j = 0; j < n_species_mixture; j++){ + phi.push_back(pow(1 + sqrt((laminarViscosity[i]) / (laminarViscosity[j])) * pow(molarMasses[j] / molarMasses[i], 0.25), 2) / sqrt(8 * (1 + molarMasses[i] / molarMasses[j]))); + wilkeDenumerator += moleFractions[j] * phi[j]; + } + wilkeDenumeratorSum.push_back(wilkeDenumerator); + wilkeDenumerator = 0.0; + phi.clear(); + wilkeNumerator.push_back(moleFractions[i] * laminarthermalConductivity[i]); + conductivityMixture += wilkeNumerator[i] / wilkeDenumeratorSum[i]; + } + return conductivityMixture; +} + +unsigned long CFluidScalar::SetTDState_T(const su2double val_temperature, su2double * const val_scalars){ + const su2double MeanMolecularWeight = ComputeMeanMolecularWeight(molarMasses, val_scalars); + + // CFluidModel::ComputeMeanSpecificHeatCp(specificHeat, val_scalars); + const su2double CpAir300Kelvin = 1009.39; + const su2double RatioSpecificHeatsAir = 1.4; + Cp = CpAir300Kelvin; + Cv = Cp/RatioSpecificHeatsAir; + Temperature = val_temperature; + Density = Pressure_Thermodynamic / ((Temperature * UNIVERSAL_GAS_CONSTANT) / MeanMolecularWeight); + + massToMoleFractions(val_scalars); + + if(wilke){ + Mu = wilkeViscosity(val_scalars); + } + else if(davidson){ + Mu = davidsonViscosity(val_scalars); + } + + Kt = wilkeConductivity(val_scalars); + + return 0; +} diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index ce9883e63515..7b6f4266e286 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -10,7 +10,8 @@ su2_cfd_src += files(['fluid/CFluidModel.cpp', 'fluid/CNEMOGas.cpp', 'fluid/CMutationTCLib.cpp', 'fluid/CSU2TCLib.cpp', - 'fluid/CFluidFlamelet.cpp']) + 'fluid/CFluidFlamelet.cpp', + 'fluid/CFluidScalar.cpp']) su2_cfd_src += files(['output/COutputFactory.cpp', 'output/CAdjElasticityOutput.cpp', diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index a287d9713db8..2a964f3c5097 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -251,7 +251,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) { - su2double yinv, Velocity_i[3]; + su2double yinv, Velocity_i[3], total_viscosity; unsigned short iDim, iVar, jVar; if (Coord_i[1] > EPS) { @@ -313,9 +313,11 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo Eddy_Viscosity_i = V_i[nDim+5]; Thermal_Conductivity_i = V_i[nDim+6]; - su2double total_viscosity; + // nijso: FIXME: the viscous terms have serious convergence issues... + total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i); + // su2double total_conductivity = Thermal_Conductivity_i + Cp_i*Eddy_Viscosity_i/Prandtl_Turb; /*--- The full stress tensor is needed for variable density ---*/ ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity); @@ -323,12 +325,20 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo /*--- Viscous terms. ---*/ residual[0] -= 0.0; - residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0][0]); - residual[2] -= Volume*(yinv*2.0*total_viscosity*PrimVar_Grad_i[2][1] - - yinv* yinv*2.0*total_viscosity*Velocity_i[1] - - TWO3*AuxVar_Grad_i[0][1]); + //residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0][0]); + //residual[2] -= Volume*(yinv*2.0*total_viscosity*PrimVar_Grad_i[2][1] - + // yinv* yinv*2.0*total_viscosity*Velocity_i[1] - + // TWO3*AuxVar_Grad_i[0][1]); + //residual[3] -= Volume*yinv*Thermal_Conductivity_i*PrimVar_Grad_i[nDim+1][1]; + su2double tau_tt = 2.0*total_viscosity * yinv*Velocity_i[1]; + su2double tau_rr = 2.0*total_viscosity * PrimVar_Grad_i[2][1]; + su2double tau_rz = tau[0][1]; + + residual[1] -= Volume*(yinv*tau_rz); + residual[2] -= Volume*(yinv*tau_rr - yinv*tau_tt); residual[3] -= Volume*yinv*Thermal_Conductivity_i*PrimVar_Grad_i[nDim+1][1]; - + // residual[3] -= 0.0; + if (implicit) jacobian[2][2] = Volume*yinv*2.0*total_viscosity*yinv; } } else { diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 185617e12a7b..790f8de840a0 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -128,7 +128,7 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ break; case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: AddHistoryOutput("RMS_PASSIVE_SCALAR", "rms[c]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean squared residual of the passive scalar equation.", HistoryFieldType::RESIDUAL); @@ -171,7 +171,7 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ break; case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: AddHistoryOutput("MAX_PASSIVE_SCALAR", "max[c]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the passive scalar equation.", HistoryFieldType::RESIDUAL); @@ -215,7 +215,7 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ break; case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: AddHistoryOutput("BGS_PASSIVE_SCALAR", "bgs[c]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the passive scalar equation.", HistoryFieldType::RESIDUAL); @@ -321,7 +321,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("RMS_RAD_ENERGY", log10(rad_solver->GetRes_RMS(0))); - + switch(scalar_model){ case PASSIVE_SCALAR: SetHistoryOutputValue("RMS_PASSIVE_SCALAR", log10(scalar_solver->GetRes_RMS(0))); @@ -337,7 +337,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("LINSOL_RESIDUAL_SCALAR", log10(scalar_solver->GetResLinSolver())); break; break; } - + SetHistoryOutputValue("MAX_PRESSURE", log10(flow_solver->GetRes_Max(0))); SetHistoryOutputValue("MAX_VELOCITY-X", log10(flow_solver->GetRes_Max(1))); SetHistoryOutputValue("MAX_VELOCITY-Y", log10(flow_solver->GetRes_Max(2))); @@ -353,7 +353,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv break; case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: SetHistoryOutputValue("MAX_PASSIVE_SCALAR", log10(scalar_solver->GetRes_Max(0))); @@ -365,7 +365,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("MAX_Y_NOX" , log10(scalar_solver->GetRes_Max(I_NOX) )); break; } - + if (multiZone){ SetHistoryOutputValue("BGS_PRESSURE", log10(flow_solver->GetRes_BGS(0))); SetHistoryOutputValue("BGS_VELOCITY-X", log10(flow_solver->GetRes_BGS(1))); @@ -386,7 +386,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv if (config->AddRadiation()) SetHistoryOutputValue("BGS_RAD_ENERGY", log10(rad_solver->GetRes_BGS(0))); - + switch(scalar_model){ case PASSIVE_SCALAR: SetHistoryOutputValue("BGS_PASSIVE_SCALAR", log10(scalar_solver->GetRes_BGS(0))); @@ -395,7 +395,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("BGS_PROGRESS_VARIABLE", log10(scalar_solver->GetRes_BGS(I_PROG_VAR))); SetHistoryOutputValue("BGS_ENTHALPY" , log10(scalar_solver->GetRes_BGS(I_ENTHALPY))); SetHistoryOutputValue("BGS_Y_CO" , log10(scalar_solver->GetRes_BGS(I_CO) )); - SetHistoryOutputValue("BGS_Y_NOX" , log10(scalar_solver->GetRes_BGS(I_NOX) )); + SetHistoryOutputValue("BGS_Y_NOX" , log10(scalar_solver->GetRes_BGS(I_NOX) )); break; } } @@ -492,10 +492,14 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: AddVolumeOutput("PASSIVE_SCALAR", "Passive_Scalar", "SOLUTION", "Passive scalar solution"); + AddVolumeOutput("DIFFUSIVITY" , "Diffusivity" , "SOLUTION", "Passive scalar diffusivity"); + AddVolumeOutput("SPECIFIC_HEAT_CP" , "Specific_Heat_Cp" , "SOLUTION", "Mixture specific heat cp"); + AddVolumeOutput("CONDUCTIVITY" , "Conductivity" , "SOLUTION", "Mixture conductivity"); + AddVolumeOutput("MEAN_MOLECULAR_WEIGHT" , "Mean_Molecular_Weight" , "SOLUTION", "Mixture molecular weight"); break; case PROGRESS_VARIABLE: AddVolumeOutput("PROGRESS_VARIABLE", "Progress_Variable", "SOLUTION", "Progress variable solution"); @@ -504,13 +508,15 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("Y_NOX" , "Y_NOx" , "SOLUTION", "NOx Mass fraction solution"); AddVolumeOutput("DIFFUSIVITY_PV" , "Diffusivity_PV" , "SOLUTION", "Diffusivity of the progress variable"); AddVolumeOutput("DIFFUSIVITY_ENTH" , "Diffusivity_Enth" , "SOLUTION", "Diffusivity of the enthalpy"); + AddVolumeOutput("SPECIFIC_HEAT_CP" , "Specific_Heat_Cp" , "SOLUTION", "Mixture specific heat cp"); + for (int i_lookup = 0; i_lookup < config->GetNLookups(); ++i_lookup) - if (config->GetLookupName(i_lookup)!="NULL"){ + if (config->GetLookupName(i_lookup)!="NULL"){ string strname1="lookup_"+config->GetLookupName(i_lookup); AddVolumeOutput(config->GetLookupName(i_lookup),strname1,"LOOKUP",config->GetLookupName(i_lookup)); } AddVolumeOutput("TABLE_MISSES" , "Table_misses" , "SOLUTION", "Lookup table misses"); - + break; case NO_SCALAR_MODEL: @@ -588,7 +594,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ case TURB_MODEL::NONE: break; } - + switch(scalar_model){ case PASSIVE_SCALAR: AddVolumeOutput("RES_PASSIVE_SCALAR", "Residual_Passive_Scalar", "RESIDUAL", "Residual of passive scalar equation"); @@ -603,7 +609,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ case NO_SCALAR_MODEL: break; } - + if (config->GetKind_SlopeLimit_Flow() != NO_LIMITER && config->GetKind_SlopeLimit_Flow() != VAN_ALBADA_EDGE) { AddVolumeOutput("LIMITER_PRESSURE", "Limiter_Pressure", "LIMITER", "Limiter value of the pressure"); AddVolumeOutput("LIMITER_VELOCITY-X", "Limiter_Velocity_x", "LIMITER", "Limiter value of the x-velocity"); @@ -628,7 +634,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ break; } } - + switch(scalar_model){ case PASSIVE_SCALAR: AddVolumeOutput("LIMITER_PASSIVE_SCALAR", "Limiter_Passive_Scalar", "LIMITER", "Limiter value for the passive scalar"); @@ -642,7 +648,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ case NO_SCALAR_MODEL: break; } - + // Hybrid RANS-LES if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES){ AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value"); @@ -726,13 +732,23 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve case TURB_MODEL::NONE: break; } - + // Solution data su2double *scalars; unsigned long table_misses; + std::vector molar_masses(config->GetNScalarsInit() + 1); + for (unsigned int iVar = 0; iVar < molar_masses.size(); iVar++) { + molar_masses.at(iVar) = config->GetMolecular_Weight(iVar); + } switch(scalar_model){ case PASSIVE_SCALAR: SetVolumeOutputValue("PASSIVE_SCALAR", iPoint, Node_Scalar->GetSolution(iPoint, 0)); + SetVolumeOutputValue("DIFFUSIVITY" , iPoint, Node_Scalar->GetDiffusivity(iPoint, 0)); + // Commented because currently Cp from FluidModel cannot show large variations (residuals blow up). Probably a bug somewhere. + // SetVolumeOutputValue("SPECIFIC_HEAT_CP" , iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); + SetVolumeOutputValue("SPECIFIC_HEAT_CP" , iPoint, 2224.43 * Node_Scalar->GetSolution(iPoint)[0] + 1009.39 * (1- Node_Scalar->GetSolution(iPoint)[0])); + SetVolumeOutputValue("CONDUCTIVITY" , iPoint, Node_Flow->GetThermalConductivity(iPoint)); + SetVolumeOutputValue("MEAN_MOLECULAR_WEIGHT" , iPoint, solver[FLOW_SOL]->GetFluidModel()->ComputeMeanMolecularWeight(molar_masses, Node_Scalar->GetSolution(iPoint))); break; case PROGRESS_VARIABLE: SetVolumeOutputValue("PROGRESS_VARIABLE", iPoint, Node_Scalar->GetSolution(iPoint, I_PROG_VAR)); @@ -741,6 +757,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("Y_NOX" , iPoint, Node_Scalar->GetSolution(iPoint, I_NOX )); SetVolumeOutputValue("DIFFUSIVITY_PV" , iPoint, Node_Scalar->GetDiffusivity(iPoint, I_PROG_VAR)); SetVolumeOutputValue("DIFFUSIVITY_ENTH" , iPoint, Node_Scalar->GetDiffusivity(iPoint, I_ENTHALPY)); + SetVolumeOutputValue("SPECIFIC_HEAT_CP" , iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); // update lookup scalars = Node_Scalar->GetSolution(iPoint); @@ -834,11 +851,11 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("RES_ENTHALPY" , iPoint, solver[SCALAR_SOL]->LinSysRes(iPoint, I_ENTHALPY)); SetVolumeOutputValue("RES_Y_CO" , iPoint, solver[SCALAR_SOL]->LinSysRes(iPoint, I_CO )); SetVolumeOutputValue("RES_Y_NOX" , iPoint, solver[SCALAR_SOL]->LinSysRes(iPoint, I_NOX )); - break; + break; case NO_SCALAR_MODEL: break; } - + if (config->GetKind_SlopeLimit_Flow() != NO_LIMITER && config->GetKind_SlopeLimit_Flow() != VAN_ALBADA_EDGE) { SetVolumeOutputValue("LIMITER_PRESSURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 0)); @@ -865,7 +882,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve break; } } - + switch(scalar_model){ case PASSIVE_SCALAR: SetVolumeOutputValue("LIMITER_PASSIVE_SCALAR", iPoint, Node_Scalar->GetLimiter(iPoint, 0)); @@ -875,11 +892,11 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("LIMITER_ENTHALPY" , iPoint, Node_Scalar->GetLimiter(iPoint, I_ENTHALPY)); SetVolumeOutputValue("LIMITER_Y_CO" , iPoint, Node_Scalar->GetLimiter(iPoint, I_CO )); SetVolumeOutputValue("LIMITER_Y_NOX" , iPoint, Node_Scalar->GetLimiter(iPoint, I_NOX )); - break; + break; case NO_SCALAR_MODEL: break; } - + if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES){ SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint)); SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); @@ -903,7 +920,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve if(config->GetKind_TimeIntScheme_Flow()==EULER_IMPLICIT){ SetVolumeOutputValue("TIMESTEP", iPoint, Node_Flow->GetDelta_Time(iPoint)); } - + // Streamwise Periodicity if(streamwisePeriodic) { SetVolumeOutputValue("RECOVERED_PRESSURE", iPoint, Node_Flow->GetStreamwise_Periodic_RecoveredPressure(iPoint)); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index b9b653e7c983..7218d32cac43 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -73,6 +73,8 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ AddHistoryOutput("AVG_CO", "Avg_CO", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total average mass fraction of CO on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); /// DESCRIPTION: Average mass fraction of NO AddHistoryOutput("AVG_NOX", "Avg_NOx", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total average mass fraction of NO on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); + /// DESCRIPTION: Average mass fraction of CH4 + AddHistoryOutput("AVG_CH4", "Avg_CH4", ScreenOutputFormat::SCIENTIFIC, "SPECIES_COEFF", "Total average mass fraction of CH4 on all markers set in MARKER_ANALYZE", HistoryFieldType::COEFFICIENT); /// END_GROUP /// BEGIN_GROUP: AERO_COEFF_SURF, DESCRIPTION: Surface values on non-solid markers. @@ -107,10 +109,12 @@ void CFlowOutput::AddAnalyzeSurfaceOutput(const CConfig *config){ AddHistoryOutputPerSurface("SURFACE_TOTAL_TEMPERATURE","Avg_TotalTemp", ScreenOutputFormat::SCIENTIFIC, "FLOW_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); /// DESCRIPTION: Average total pressure AddHistoryOutputPerSurface("SURFACE_TOTAL_PRESSURE", "Avg_TotalPress", ScreenOutputFormat::SCIENTIFIC, "FLOW_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); - /// DESCRIPTION: Average mass fraction of CO + /// DESCRIPTION: Average mass fraction of CO AddHistoryOutputPerSurface("AVG_CO", "Avg_CO", ScreenOutputFormat::SCIENTIFIC, "FLOW_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); - /// DESCRIPTION: Average mass fraction of NO + /// DESCRIPTION: Average mass fraction of NO AddHistoryOutputPerSurface("AVG_NOX", "Avg_NOx", ScreenOutputFormat::SCIENTIFIC, "FLOW_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); + /// DESCRIPTION: Average mass fraction of CH4 + AddHistoryOutputPerSurface("AVG_CH4", "Avg_CH4", ScreenOutputFormat::SCIENTIFIC, "FLOW_COEFF_SURF", Marker_Analyze, HistoryFieldType::COEFFICIENT); /// END_GROUP } @@ -157,6 +161,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, vector Surface_MassFlow_Abs (nMarker,0.0); vector Surface_CO (nMarker,0.0); vector Surface_NOx (nMarker,0.0); + vector Surface_CH4 (nMarker,0.0); su2double Tot_Surface_MassFlow = 0.0; su2double Tot_Surface_Mach = 0.0; @@ -173,6 +178,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, su2double Tot_SecondOverUniformity = 0.0; su2double Tot_Surface_CO = 0.0; su2double Tot_Surface_NOx = 0.0; + su2double Tot_Surface_CH4 = 0.0; //su2double Tot_Surface_Scalar[n_scalars]; //for (int i_scalar = 0; i_scalar < n_scalars; ++i_scalar) // Tot_Surface_Scalar[i_scalar] = 0.0; @@ -279,6 +285,10 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, Surface_NOx[iMarker] += scalar_solver->GetNodes()->GetSolution(iPoint, I_NOX) * Weight; } + if (config->GetKind_FluidModel() == MIXTURE_FLUID_MODEL) { + Surface_CH4[iMarker] += scalar_solver->GetNodes()->GetSolution(iPoint, 0) * Weight; + } + /*--- For now, always used the area to weight the uniformities. ---*/ Weight = abs(Area); @@ -306,9 +316,10 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, vector Surface_TotalPressure_Local (nMarker_Analyze,0.0); vector Surface_Area_Local (nMarker_Analyze,0.0); vector Surface_MassFlow_Abs_Local (nMarker_Analyze,0.0); - vector Surface_CO_Local (nMarker_Analyze,0.0); + vector Surface_CO_Local (nMarker_Analyze,0.0); vector Surface_NOx_Local (nMarker_Analyze,0.0); - + vector Surface_CH4_Local (nMarker_Analyze,0.0); + vector Surface_MassFlow_Total (nMarker_Analyze,0.0); vector Surface_Mach_Total (nMarker_Analyze,0.0); vector Surface_Temperature_Total (nMarker_Analyze,0.0); @@ -324,6 +335,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, vector Surface_MassFlow_Abs_Total (nMarker_Analyze,0.0); vector Surface_CO_Total (nMarker_Analyze,0.0); vector Surface_NOx_Total (nMarker_Analyze,0.0); + vector Surface_CH4_Total (nMarker_Analyze,0.0); vector Surface_MomentumDistortion_Total (nMarker_Analyze,0.0); @@ -353,6 +365,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, Surface_MassFlow_Abs_Local[iMarker_Analyze] += Surface_MassFlow_Abs[iMarker]; Surface_CO_Local[iMarker_Analyze] += Surface_CO[iMarker]; Surface_NOx_Local[iMarker_Analyze] += Surface_NOx[iMarker]; + Surface_CH4_Local[iMarker_Analyze] += Surface_CH4[iMarker]; } } @@ -380,6 +393,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, Allreduce(Surface_MassFlow_Abs_Local, Surface_MassFlow_Abs_Total); Allreduce(Surface_CO_Local, Surface_CO_Total); Allreduce(Surface_NOx_Local,Surface_NOx_Total); + Allreduce(Surface_CH4_Local, Surface_CH4_Total); /*--- Compute the value of Surface_Area_Total, and Surface_Pressure_Total, and set the value in the config structure for future use ---*/ @@ -400,6 +414,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, Surface_TotalPressure_Total[iMarker_Analyze] /= Weight; Surface_CO_Total[iMarker_Analyze] /= Weight; Surface_NOx_Total[iMarker_Analyze] /= Weight; + Surface_CH4_Total[iMarker_Analyze] /= Weight; } else { Surface_Mach_Total[iMarker_Analyze] = 0.0; @@ -412,6 +427,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, Surface_TotalPressure_Total[iMarker_Analyze] = 0.0; Surface_CO_Total[iMarker_Analyze] = 0.0; Surface_NOx_Total[iMarker_Analyze] = 0.0; + Surface_CH4_Total[iMarker_Analyze] = 0.0; } /*--- Compute flow uniformity parameters separately (always area for now). ---*/ @@ -506,11 +522,16 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, SetHistoryOutputPerSurfaceValue("AVG_CO", y_CO, iMarker_Analyze); Tot_Surface_CO += y_CO; config->SetSurface_CO(iMarker_Analyze, y_CO); - + su2double y_NOx = Surface_NOx_Total[iMarker_Analyze]; SetHistoryOutputPerSurfaceValue("AVG_NOX", y_NOx, iMarker_Analyze); Tot_Surface_NOx += y_NOx; config->SetSurface_NOx(iMarker_Analyze, y_NOx); + + su2double y_CH4 = Surface_CH4_Total[iMarker_Analyze]; + SetHistoryOutputPerSurfaceValue("AVG_CH4", y_CH4, iMarker_Analyze); + Tot_Surface_CH4 += y_CH4; + config->SetSurface_CH4(iMarker_Analyze, y_CH4); } /*--- Compute the average static pressure drop between two surfaces. Note @@ -542,6 +563,7 @@ void CFlowOutput::SetAnalyzeSurface(CSolver **solver, const CGeometry *geometry, SetHistoryOutputValue("SURFACE_TOTAL_PRESSURE", Tot_Surface_TotalPressure); SetHistoryOutputValue("AVG_CO", Tot_Surface_CO); SetHistoryOutputValue("AVG_NOX", Tot_Surface_NOx); + SetHistoryOutputValue("AVG_CH4", Tot_Surface_CH4); if ((rank == MASTER_NODE) && !config->GetDiscrete_Adjoint() && output) { diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 014ca789b8a3..8fc7d74071e7 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -32,6 +32,7 @@ #include "../../include/fluid/CIncIdealGasPolynomial.hpp" #include "../../include/variables/CIncNSVariable.hpp" #include "../include/fluid/CFluidFlamelet.hpp" +#include "../include/fluid/CFluidScalar.hpp" #include "../../include/limiters/CLimiterDetails.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" @@ -314,6 +315,19 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i config->SetPressure_Thermodynamic(Pressure_Thermodynamic); break; + case MIXTURE_FLUID_MODEL: + n_scalars = config->GetNScalarsInit(); + config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT/(config->GetMolecular_Weight(n_scalars)/1000.0)); + // Note: Density_FreeStream = INC_DENSITY_INT and Temperature_FreeStream = INC_TEMPERATURE_INIT. + // Mixture density is computed based on the MeanMolecularWeight in CFluidScalar.cpp + Pressure_Thermodynamic = Density_FreeStream*Temperature_FreeStream*config->GetGas_Constant(); + auxFluidModel = new CFluidScalar(config, Pressure_Thermodynamic); + dummy_scalar = new su2double[n_scalars](); + dummy_scalar[n_scalars-1] = 1; + auxFluidModel->SetTDState_T(Temperature_FreeStream, dummy_scalar); + config->SetPressure_Thermodynamic(Pressure_Thermodynamic); + break; + default: SU2_MPI::Error("Fluid model not implemented for incompressible solver.", CURRENT_FUNCTION); @@ -487,6 +501,14 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i fluidModel->SetTDState_T(Temperature_FreeStream,dummy_scalar); delete[] dummy_scalar; break; + + case MIXTURE_FLUID_MODEL: + fluidModel = new CFluidScalar(config, Pressure_Thermodynamic); + n_scalars = fluidModel->GetNScalars(); + dummy_scalar = new su2double[n_scalars](); + fluidModel->SetTDState_T(Temperature_FreeStream, dummy_scalar); + delete[] dummy_scalar; + break; } if (viscous) { diff --git a/SU2_CFD/src/solvers/CPassiveScalarSolver.cpp b/SU2_CFD/src/solvers/CPassiveScalarSolver.cpp index 390e922ad905..b43a99b731ad 100644 --- a/SU2_CFD/src/solvers/CPassiveScalarSolver.cpp +++ b/SU2_CFD/src/solvers/CPassiveScalarSolver.cpp @@ -31,6 +31,8 @@ #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include + CPassiveScalarSolver::CPassiveScalarSolver(void) : CScalarLegacySolver() {} CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, @@ -48,9 +50,9 @@ CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, /*--- Dimension of the problem --> passive scalar will only ever have a single equation. Other child classes of CScalarLegacySolver can have variable numbers of equations. ---*/ - - nVar = 1; - nPrimVar = 1; + + nVar = config->GetNScalarsInit(); + nPrimVar = nVar; nPoint = geometry->GetnPoint(); nPointDomain = geometry->GetnPointDomain(); @@ -62,11 +64,11 @@ CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, /*--- Define geometry constants in the solver structure ---*/ nDim = geometry->GetnDim(); - + /*--- Fluid model pointer initialization ---*/ - + FluidModel = nullptr; - + /*--- Single grid simulation ---*/ @@ -153,9 +155,11 @@ CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, nodes = new CPassiveScalarVariable(Scalar_Inf, nPoint, nDim, nVar, config); + //TODO MH: for MIXTURE_FLUID_MODEL a nan is created due to all inputs of SetDiffusivity being zero /*--- initialize the mass diffusivity ---*/ for (auto iVar = 0u; iVar < nVar; iVar++){ - auto M = FluidModel->GetMassDiffusivity(); // returns a su2double, note that for more species this should be a vector + // auto M = FluidModel->GetMassDiffusivity(); // returns a su2double, note that for more species this should be a vector + auto M = config->GetDiffusivity_Constant(); // Right now DIFFUSIVITY_CONSTANT needs to be specified next to UNITY_LEWIS in order to initialise su2double M // loop over all points and set diffusivity // why construct the entire diffusivity matrix? for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) @@ -185,7 +189,7 @@ CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, /*-- Allocation of inlets has to happen in derived classes (not CScalarLegacySolver), due to arbitrary number of scalar variables. First, we also set the column index for any inlet profiles. ---*/ - + Inlet_Position = nDim*2+2; if (turbulent) { if (turb_SA) Inlet_Position += 1; @@ -214,7 +218,7 @@ CPassiveScalarSolver::CPassiveScalarSolver(CGeometry *geometry, SetImplicitPeriodic(true); /*--- Store the initial CFL number for all grid points. ---*/ - + const su2double CFL = config->GetCFL(MGLevel)*config->GetCFLRedCoeff_Scalar(); for (auto iPoint = 0u; iPoint < nPoint; iPoint++) { nodes->SetLocalCFL(iPoint, CFL); @@ -237,9 +241,49 @@ void CPassiveScalarSolver::Preprocessing(CGeometry *geometry, CSolver **solver_c unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool muscl = config->GetMUSCL_Scalar(); - const bool limiter = (config->GetKind_SlopeLimit_Scalar() != NO_LIMITER) && - (config->GetInnerIter() <= config->GetLimiterIter()); + const bool muscl = config->GetMUSCL_Scalar(); + const bool limiter = (config->GetKind_SlopeLimit_Scalar() != NO_LIMITER) && + (config->GetInnerIter() <= config->GetLimiterIter()); + + su2double * scalars = new su2double[nVar]; + + for (auto i_point = 0u; i_point < nPoint; i_point++) { + + // Still needed later on to access GetMassDiffusivity() + // CFluidModel * fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + + for(int iVar = 0; iVar < nVar; iVar++){ + scalars[iVar] = nodes->GetSolution(i_point, iVar); + } + + // Temporary solution: compute mass diffusivity here with variable Cp. + std::vector SpecificHeatCp; + unsigned short n_species_mixture = nVar + 1; + SpecificHeatCp.resize(n_species_mixture); + + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + SpecificHeatCp.at(iVar) = config->GetSpecific_Heat_Cp(iVar); + } + + su2double cp = SpecificHeatCp[0] * scalars[0] + SpecificHeatCp[1] * (1- scalars[0]); + su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(i_point); + su2double thermal_conductivity = solver_container[FLOW_SOL]->GetNodes()->GetThermalConductivity(i_point); + su2double Lewis = 1; + su2double mass_diffusivity = thermal_conductivity / (Lewis * density * cp); + + for(auto iVar = 0u; iVar < nVar; iVar++){ + nodes->SetDiffusivity(i_point, mass_diffusivity, iVar); + } + + // for(auto iVar = 0u; iVar < nVar; iVar++){ + // nodes->SetDiffusivity(i_point, fluid_model_local->GetMassDiffusivity(), iVar); + // } + + if (!Output) LinSysRes.SetBlock_Zero(i_point); + + } + + delete[] scalars; /*--- Clear residual and system matrix, not needed for * reducer strategy as we write over the entire matrix. ---*/ @@ -282,27 +326,27 @@ void CPassiveScalarSolver::SetInitialCondition(CGeometry **geometry, CConfig *config, unsigned long ExtIter) { bool Restart = (config->GetRestart() || config->GetRestart_Flow()); - - + + if ((!Restart) && ExtIter == 0) { if (rank == MASTER_NODE){ cout << "Initializing passive scalar (initial condition)." << endl; cout << "initialization = " << nVar << " " << config->GetScalar_Init(0)<GetnMGLevels(); i_mesh++) { for (unsigned long i_point = 0; i_point < geometry[i_mesh]->GetnPoint(); i_point++) { - + for (unsigned long i_var = 0; i_var < nVar; i_var++) Solution[i_var] = 0.0; - - for(int i_scalar = 0;i_scalar < nVar;i_scalar++){ - scalar_init[i_scalar] = config->GetScalar_Init(i_scalar); + + for(int iVar = 0; iVar < nVar; iVar++){ + scalar_init[iVar] = config->GetScalar_Init(iVar); } - + solver_container[i_mesh][SCALAR_SOL]->GetNodes()->SetSolution(i_point, scalar_init); } @@ -314,7 +358,7 @@ void CPassiveScalarSolver::SetInitialCondition(CGeometry **geometry, solver_container[i_mesh][FLOW_SOL]->CompleteComms(geometry[i_mesh], config, SOLUTION); solver_container[i_mesh][FLOW_SOL]->Preprocessing( geometry[i_mesh], solver_container[i_mesh], config, i_mesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); - + } @@ -326,51 +370,55 @@ void CPassiveScalarSolver::SetInitialCondition(CGeometry **geometry, void CPassiveScalarSolver::SetPreconditioner(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - + unsigned short iVar; unsigned long iPoint, total_index; - + su2double BetaInc2, Density, dRhodT, dRhodC, Temperature, Delta; - + bool variable_density = (config->GetKind_DensityModel() == INC_DENSITYMODEL::VARIABLE); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - + /*--- Access the primitive variables at this node. ---*/ - + Density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); BetaInc2 = solver_container[FLOW_SOL]->GetNodes()->GetBetaInc2(iPoint); Temperature = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint); - + unsigned short nVar_Flow = solver_container[FLOW_SOL]->GetnVar(); - + su2double SolP = solver_container[FLOW_SOL]->LinSysSol[iPoint*nVar_Flow+0]; su2double SolT = solver_container[FLOW_SOL]->LinSysSol[iPoint*nVar_Flow+nDim+1]; - + /*--- We need the derivative of the equation of state to build the preconditioning matrix. For now, the only option is the ideal gas law, but in the future, dRhodT should be in the fluid model. ---*/ - + if (variable_density) { dRhodT = -Density/Temperature; } else { dRhodT = 0.0; } - - /*--- Passive scalars have no impact on the density. ---*/ - - dRhodC = 0.0; - + + std::vector MolecularWeight; + unsigned short n_species_mixture = nVar + 1; + MolecularWeight.resize(n_species_mixture); + + for (int iVar = 0; iVar < n_species_mixture; iVar++) { + MolecularWeight.at(iVar) = config->GetMolecular_Weight(iVar); + } + /*--- Modify matrix diagonal with term including volume and time step. ---*/ - + su2double Vol = geometry->nodes->GetVolume(iPoint); Delta = Vol / (config->GetCFLRedCoeff_Scalar()* solver_container[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint)); - + /*--- Calculating the inverse of the preconditioning matrix that multiplies the time derivative during time integration. ---*/ - + if (implicit) { // nijso: do we need to wipe the entire jacobian for preconditioning? //for (int i_var = 0; i_var < nVar; i_var++) { @@ -378,33 +426,36 @@ void CPassiveScalarSolver::SetPreconditioner(CGeometry *geometry, CSolver **solv // Jacobian_i[i_var][j_var] = 0.0; // } //} - + for (iVar = 0; iVar < nVar; iVar++) { - + total_index = iPoint*nVar+iVar; - + su2double c = nodes->GetSolution(iPoint,iVar); - + dRhodC = -Density*((MolecularWeight[1]-MolecularWeight[0])/(MolecularWeight[0]+(MolecularWeight[1]-MolecularWeight[0])*c)); + /*--- Compute the lag terms for the decoupled linear system from the mean flow equations and add to the residual for the scalar. In short, we are effectively making these terms explicit. ---*/ - - su2double artcompc1 = SolP * c/(Density*BetaInc2); - su2double artcompc2 = SolT * dRhodT * c/(Density); - + + su2double artcompc1 = (SolP * c)/(BetaInc2*(dRhodC*c + Density)); + su2double artcompc2 = (SolT * dRhodT * c)/(dRhodC*c + Density); + // su2double artcompc1 = SolP * c/(Density*BetaInc2); + // su2double artcompc2 = SolT * dRhodT * c/(Density); + LinSysRes[total_index] += artcompc1 + artcompc2; - + /*--- Add the extra Jacobian term to the scalar system. ---*/ - + su2double Jaccomp = c * dRhodC + Density; //This is Gamma su2double JacTerm = Jaccomp*Delta; - + Jacobian.AddVal2Diag(iPoint, JacTerm); - + } - + } - + } } @@ -438,9 +489,9 @@ void CPassiveScalarSolver::Source_Residual(CGeometry *geometry, CSolver **solver /*--- Gradient of the primitive and conservative variables ---*/ numerics->SetPrimVarGradient(flowNodes->GetGradient_Primitive(iPoint), nullptr); - + /*--- Scalar variables w/o reconstruction, and its gradient ---*/ - + numerics->SetScalarVar(nodes->GetSolution(iPoint), nullptr); numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nullptr); @@ -532,14 +583,14 @@ void CPassiveScalarSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_conta /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ if (geometry->nodes->GetDomain(iPoint)) { - + /*--- Allocate the value at the outlet ---*/ - auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - + auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] = nodes->GetSolution(Point_Normal, iVar); nodes->SetSolution_Old(iPoint, Solution); - + LinSysRes.SetBlock_Zero(iPoint); for (auto iVar = 0u; iVar < nVar; iVar++){ @@ -560,7 +611,7 @@ void CPassiveScalarSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_conta void CPassiveScalarSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - + } diff --git a/SU2_CFD/src/solvers/CScalarLegacySolver.cpp b/SU2_CFD/src/solvers/CScalarLegacySolver.cpp index 0df7333086ea..62846352d5fc 100644 --- a/SU2_CFD/src/solvers/CScalarLegacySolver.cpp +++ b/SU2_CFD/src/solvers/CScalarLegacySolver.cpp @@ -361,6 +361,11 @@ void CScalarLegacySolver::PrepareImplicitIteration(CGeometry *geometry, CSolver* const auto flowNodes = solver_container[FLOW_SOL]->GetNodes(); +/*--- use preconditioner for scalar transport equations ---*/ +// const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); +// if (incompressible) SetPreconditioner(geometry, solver_container, config); + + /*--- Set shared residual variables to 0 and declare * local ones for current thread to work on. ---*/ diff --git a/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_axisymmetric_velocityinlets.cfg b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_axisymmetric_velocityinlets.cfg new file mode 100755 index 000000000000..7ab5b9e05b61 --- /dev/null +++ b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_axisymmetric_velocityinlets.cfg @@ -0,0 +1,313 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Composition dependent incompressible flow with variable fluid properties through axisymmetric venturi. % +% Mixing of methane and air. % +% Author: Mark Heimgartner % +% Institution: Bosch Thermotechniek Deventer % +% Date: 08/28/2021 % +% File Version 7.1.1 "Blackbird", branch: feature_multicomp % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER = INC_RANS + +% Specify turbulent model (NONE, SA, SA_NEG, SST) +KIND_TURB_MODEL = SST + +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM = DIRECT + +% Restart solution (NO, YES) +RESTART_SOL = NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL = VARIABLE + +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES + +% Initial density for incompressible flows +% (1.2886 kg/m^3 by default (air), 998.2 Kg/m^3 (water)) +INC_DENSITY_INIT = 1.1766 + +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT = (1.00, 0.0, 0.0 ) + +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT = 300.0 + +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM = DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, MIXTURE_FLUID_MODEL) +FLUID_MODEL = MIXTURE_FLUID_MODEL + +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +MOLECULAR_WEIGHT= 16.043, 28.965 + +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. (CONSTANT_DENSITY, INC_IDEAL_GAS) and the heat equation. +SPECIFIC_HEAT_CP = 2224.43, 1009.39 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Laminar Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL, +% POLYNOMIAL_CONDUCTIVITY). +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY + +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0258 + +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72, 0.72 + +% Definition of the turbulent thermal conductivity model for RANS +% (CONSTANT_PRANDTL_TURB by default, NONE). +TURBULENT_CONDUCTIVITY_MODEL= NONE + +% Turbulent Prandtl number (0.9 (air) by default) +PRANDTL_TURB= 0.90, 0.90 + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_SCALAR= SCALAR_UPWIND + +% Specify scalar transport model (NONE, PASSIVE_SCALAR, PROGRESS_VARIABLE) +KIND_SCALAR_MODEL= PASSIVE_SCALAR + +% Specify mass diffusivity model (CONSTANT_DIFFUSIVITY, CONSTANT_SCHMIDT, FLAMELET, UNITY_LEWIS) +%DIFFUSIVITY_MODEL = UNITY_LEWIS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the scalar transport equations. +MUSCL_SCALAR= NO + +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_SCALAR = NONE + +% Time discretization for scalar transport problems (EULER_IMPLICIT) +TIME_DISCRE_SCALAR= EULER_IMPLICIT + +% Scalar clipping +SCALAR_CLIPPING= NO +SCALAR_INIT = 0 +SCALAR_CLIPPING_MIN= 0.0 +SCALAR_CLIPPING_MAX= 1.0 + +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_TURB= NO + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY, POLYNOMIAL_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY + +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.1102e-05, 1.8551e-05 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Navier-Stokes (no-slip), constant heat flux wall marker(s) (NONE = no marker) +% Format: ( marker name, constant heat flux (J/m^2), ... ) +MARKER_HEATFLUX= ( wall, 0.0 ) + +% Symmetry boundary marker(s) (NONE = no marker) +% Implementation identical to MARKER_EULER. +MARKER_SYM= ( axis ) + +% Axisymmetric simulation, only compressible flows (NO, YES) +AXISYMMETRIC= YES + +% List of inlet types for incompressible flows. List length must +% match number of inlet markers. Options: VELOCITY_INLET, PRESSURE_INLET. +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET VELOCITY_INLET + +% Inlet inlet scalar: (inlet marker, scalar initialisation value) +MARKER_INLET_SCALAR = (gas_inlet, 1.0, air_axial_inlet, 0.0, air_radial_inlet, 0.0) + +% Inc. Velocity: (inlet marker, temperature, velocity magnitude, flow_direction_x, +% flow_direction_y, flow_direction_z, ... ) where flow_direction is +% a unit vector. +MARKER_INLET= ( gas_inlet, 300, 12.3, 1.0, 0.0, 0.0, air_axial_inlet, 300, 2.2, 1.0, 0.0, 0.0, air_radial_inlet, 300, 4.9, 0.0, -1.0, 0.0 ) + +% List of outlet types for incompressible flows. List length must +% match number of outlet markers. Options: PRESSURE_OUTLET, MASS_FLOW_OUTLET +INC_OUTLET_TYPE= PRESSURE_OUTLET + +% Outlet boundary marker(s) (NONE = no marker) +% Compressible: ( outlet marker, back pressure (static thermodynamic), ... ) +% Inc. Pressure: ( outlet marker, back pressure (static gauge in Pa), ... ) +% Inc. Mass Flow: ( outlet marker, mass flow target (kg/s), ... ) +MARKER_OUTLET= ( outlet, 0.0) + +% Marker(s) of the surface that is going to be analyzed in detail (massflow, average pressure, distortion, etc) +MARKER_ANALYZE = ( outlet ) + +% Method to compute the average value in MARKER_ANALYZE (AREA, MASSFLUX). +MARKER_ANALYZE_AVERAGE = MASSFLUX + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES + +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER = 100 + +% Reduction factor of the CFL coefficient in the scalar transport problem +CFL_REDUCTION_SCALAR = 1.0 + +% Reduction factor of the CFL coefficient in the turbulence problem +CFL_REDUCTION_TURB = 1.0 + +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO + +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) + +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% Number of total iterations +ITER= 50000 + +% Writing solution file frequency +OUTPUT_WRT_FREQ=1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver for implicit formulations (BCGSTAB, FGMRES) +LINEAR_SOLVER= FGMRES + +% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS) +LINEAR_SOLVER_PREC= ILU + +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-3 + +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES + +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_FLOW = NONE + +% Coefficient for the limiter (smooth regions) +VENKAT_LIMITER_COEFF= 10.0 + +% 2nd and 4th order artificial dissipation coefficients +JST_SENSOR_COEFF= ( 0.5, 0.02 ) + +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence criteria (CAUCHY, RESIDUAL) +CONV_CRITERIA= RESIDUAL + +CONV_FIELD = (RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE) + +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL = -10 + +% Start convergence criteria at iteration number +CONV_STARTITER= 10 + +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 + +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +SCREEN_OUTPUT = INNER_ITER WALL_TIME RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_PASSIVE_SCALAR LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_SCALAR LINSOL_RESIDUAL_SCALAR + +HISTORY_OUTPUT = RMS_RES MAX_RES FLOW_COEFF FLOW_COEFF_SURF + +VOLUME_OUTPUT = DENSITY SOLUTION PRIMITIVE SOURCE RESIDUAL + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Mesh input file +MESH_FILENAME = C6_pneumatic_venturi_planar_mesh.cgns + +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= CGNS + +% Mesh output file +MESH_OUT_FILENAME= mesh_out.su2 + +% Restart flow input file +SOLUTION_FILENAME= initial_solution.dat + +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat + +% Output file format (PARAVIEW, PARAVIEW_BINARY, TECPLOT, STL) +OUTPUT_FILES = (RESTART, PARAVIEW, SURFACE_PARAVIEW) +TABULAR_FORMAT = CSV + +% Output file convergence history (w/o extension) +CONV_FILENAME= history + +% Output file restart flow +RESTART_FILENAME= restart.dat + +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat + +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow + +% Overwrite solution and visualization files or not +WRT_SOL_OVERWRITE = NO + +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint + +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat + +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow + +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint + +% Write the performance summary at the end of a calculation. +WRT_PERFORMANCE = YES \ No newline at end of file diff --git a/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_outlettargetmassflowrate.cfg b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_outlettargetmassflowrate.cfg new file mode 100755 index 000000000000..1c77d986b7ef --- /dev/null +++ b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_outlettargetmassflowrate.cfg @@ -0,0 +1,315 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Composition dependent incompressible flow with variable fluid properties through planar venturi % +% Author: Mark Heimgartner % +% Institution: Bosch Thermotechniek Deventer % +% Date: 07/07/2021 % +% File Version 7.1.1 "Blackbird", branch: feature_multicomp % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER = INC_RANS + +% Specify turbulent model (NONE, SA, SA_NEG, SST) +KIND_TURB_MODEL = SST + +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM = DIRECT + +% Restart solution (NO, YES) +RESTART_SOL = NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL = VARIABLE + +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES + +% Initial density for incompressible flows +% (1.2886 kg/m^3 by default (air), 998.2 Kg/m^3 (water)) +INC_DENSITY_INIT = 1.1766 + +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT = (1.00, 0.0, 0.0 ) + +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT = 300.0 + +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM = DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, MIXTURE_FLUID_MODEL) +FLUID_MODEL = MIXTURE_FLUID_MODEL + +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +MOLECULAR_WEIGHT= 16.043, 28.965 + +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. (CONSTANT_DENSITY, INC_IDEAL_GAS) and the heat equation. +SPECIFIC_HEAT_CP = 2224.43, 1009.39 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Laminar Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL, +% POLYNOMIAL_CONDUCTIVITY). +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY + +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0258 + +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72, 0.72 + +% Definition of the turbulent thermal conductivity model for RANS +% (CONSTANT_PRANDTL_TURB by default, NONE). +TURBULENT_CONDUCTIVITY_MODEL= NONE + +% Turbulent Prandtl number (0.9 (air) by default) +PRANDTL_TURB= 0.90, 0.90 + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_SCALAR= SCALAR_UPWIND + +% Specify scalar transport model (NONE, PASSIVE_SCALAR, PROGRESS_VARIABLE) +KIND_SCALAR_MODEL= PASSIVE_SCALAR + +% Specify mass diffusivity model (CONSTANT_DIFFUSIVITY, CONSTANT_SCHMIDT, FLAMELET, UNITY_LEWIS) +DIFFUSIVITY_MODEL = UNITY_LEWIS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the scalar transport equations. +MUSCL_SCALAR= YES + +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_SCALAR = NONE + +% Time discretization for scalar transport problems (EULER_IMPLICIT) +TIME_DISCRE_SCALAR= EULER_IMPLICIT + +% Scalar clipping +SCALAR_CLIPPING= YES +SCALAR_INIT = 1 +SCALAR_CLIPPING_MIN= 0.0 +SCALAR_CLIPPING_MAX= 1.0 + +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_TURB= NO + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY, POLYNOMIAL_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY + +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.1102e-05, 1.8551e-05 + +% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) +MU_REF= 1.118e-05, 1.795e-05 + +% Sutherland Temperature Ref (273.15 K default value for AIR SI) +MU_T_REF= 300, 300 + +% Sutherland constant (110.4 default value for AIR SI) +SUTHERLAND_CONSTANT= 187.78, 129.4 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Navier-Stokes (no-slip), constant heat flux wall marker(s) (NONE = no marker) +% Format: ( marker name, constant heat flux (J/m^2), ... ) +MARKER_HEATFLUX= ( wall, 0.0 ) + +% Symmetry boundary marker(s) (NONE = no marker) +% Implementation identical to MARKER_EULER. +MARKER_SYM= ( axis ) + +% Axisymmetric simulation, only compressible flows (NO, YES) +AXISYMMETRIC= NO + +% List of inlet types for incompressible flows. List length must +% match number of inlet markers. Options: VELOCITY_INLET, PRESSURE_INLET. +INC_INLET_TYPE= PRESSURE_INLET PRESSURE_INLET PRESSURE_INLET + +% Inlet inlet scalar: (inlet marker, scalar initialisation value) +MARKER_INLET_SCALAR = (gas_inlet, 1.0, air_axial_inlet, 0.0, air_radial_inlet, 0.0) + +% Inc. Pressure: (inlet marker, temperature, total pressure, flow_direction_x, +% flow_direction_y, flow_direction_z, ... ) where flow_direction is +% a unit vector. +MARKER_INLET= ( gas_inlet, 300, -11, 1.0, 0.0, 0.0, air_axial_inlet, 300, 0, 1.0, 0.0, 0.0, air_radial_inlet, 300, 0, 0.0, -1.0, 0.0 ) + +% List of outlet types for incompressible flows. List length must +% match number of outlet markers. Options: PRESSURE_OUTLET, MASS_FLOW_OUTLET +INC_OUTLET_TYPE= MASS_FLOW_OUTLET + +% Outlet boundary marker(s) (NONE = no marker) +% Compressible: ( outlet marker, back pressure (static thermodynamic), ... ) +% Inc. Pressure: ( outlet marker, back pressure (static gauge in Pa), ... ) +% Inc. Mass Flow: ( outlet marker, mass flow target (kg/s), ... ) +MARKER_OUTLET= ( outlet, 0.184348) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES + +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER = 60 + +% Reduction factor of the CFL coefficient in the scalar transport problem +CFL_REDUCTION_SCALAR = 1.0 + +% Reduction factor of the CFL coefficient in the turbulence problem +CFL_REDUCTION_TURB = 1.0 + +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO + +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) + +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% Number of total iterations +ITER= 50000 + +% Writing solution file frequency +OUTPUT_WRT_FREQ=1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver for implicit formulations (BCGSTAB, FGMRES) +LINEAR_SOLVER= FGMRES + +% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS) +LINEAR_SOLVER_PREC= ILU + +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-8 + +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 50 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES + +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_FLOW = NONE + +% Coefficient for the limiter (smooth regions) +VENKAT_LIMITER_COEFF= 10.0 + +% 2nd and 4th order artificial dissipation coefficients +JST_SENSOR_COEFF= ( 0.5, 0.02 ) + +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence criteria (CAUCHY, RESIDUAL) +CONV_CRITERIA= RESIDUAL + +CONV_FIELD = (RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE) + +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL = -10 + +% Start convergence criteria at iteration number +CONV_STARTITER= 10 + +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 + +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +SCREEN_OUTPUT = INNER_ITER WALL_TIME RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_PASSIVE_SCALAR LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_SCALAR LINSOL_RESIDUAL_SCALAR + +HISTORY_OUTPUT = RMS_RES MAX_RES FLOW_COEFF + +VOLUME_OUTPUT = DENSITY SOLUTION PRIMITIVE SOURCE RESIDUAL LOOKUP + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Mesh input file +MESH_FILENAME = C6_pneumatic_venturi_planar_mesh.cgns + +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= CGNS + +% Mesh output file +MESH_OUT_FILENAME= mesh_out.su2 + +% Restart flow input file +SOLUTION_FILENAME= initial_solution.dat + +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat + +% Output file format (PARAVIEW, PARAVIEW_BINARY, TECPLOT, STL) +OUTPUT_FILES = (RESTART, PARAVIEW, SURFACE_PARAVIEW) +TABULAR_FORMAT = CSV + +% Output file convergence history (w/o extension) +CONV_FILENAME= history + +% Output file restart flow +RESTART_FILENAME= restart.dat + +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat + +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow + +% Overwrite solution and visualization files or not +WRT_SOL_OVERWRITE = NO + +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint + +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat + +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow + +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint + +% Write the performance summary at the end of a calculation. +WRT_PERFORMANCE = YES \ No newline at end of file diff --git a/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets.cfg b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets.cfg new file mode 100755 index 000000000000..f375cc894d84 --- /dev/null +++ b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets.cfg @@ -0,0 +1,315 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Composition dependent incompressible flow with variable fluid properties through planar venturi % +% Author: Mark Heimgartner % +% Institution: Bosch Thermotechniek Deventer % +% Date: 07/07/2021 % +% File Version 7.1.1 "Blackbird", branch: feature_multicomp % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER = INC_RANS + +% Specify turbulent model (NONE, SA, SA_NEG, SST) +KIND_TURB_MODEL = SST + +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM = DIRECT + +% Restart solution (NO, YES) +RESTART_SOL = NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL = VARIABLE + +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES + +% Initial density for incompressible flows +% (1.2886 kg/m^3 by default (air), 998.2 Kg/m^3 (water)) +INC_DENSITY_INIT = 1.1766 + +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT = (1.00, 0.0, 0.0 ) + +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT = 300.0 + +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM = DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, MIXTURE_FLUID_MODEL) +FLUID_MODEL = MIXTURE_FLUID_MODEL + +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +MOLECULAR_WEIGHT= 16.043, 28.965 + +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. (CONSTANT_DENSITY, INC_IDEAL_GAS) and the heat equation. +SPECIFIC_HEAT_CP = 2224.43, 1009.39 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Laminar Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL, +% POLYNOMIAL_CONDUCTIVITY). +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY + +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0258 + +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72, 0.72 + +% Definition of the turbulent thermal conductivity model for RANS +% (CONSTANT_PRANDTL_TURB by default, NONE). +TURBULENT_CONDUCTIVITY_MODEL= NONE + +% Turbulent Prandtl number (0.9 (air) by default) +PRANDTL_TURB= 0.90, 0.90 + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_SCALAR= SCALAR_UPWIND + +% Specify scalar transport model (NONE, PASSIVE_SCALAR, PROGRESS_VARIABLE) +KIND_SCALAR_MODEL= PASSIVE_SCALAR + +% Specify mass diffusivity model (CONSTANT_DIFFUSIVITY, CONSTANT_SCHMIDT, FLAMELET, UNITY_LEWIS) +DIFFUSIVITY_MODEL = UNITY_LEWIS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the scalar transport equations. +MUSCL_SCALAR= YES + +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_SCALAR = NONE + +% Time discretization for scalar transport problems (EULER_IMPLICIT) +TIME_DISCRE_SCALAR= EULER_IMPLICIT + +% Scalar clipping +SCALAR_CLIPPING= YES +SCALAR_INIT = 1 +SCALAR_CLIPPING_MIN= 0.0 +SCALAR_CLIPPING_MAX= 1.0 + +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_TURB= NO + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY, POLYNOMIAL_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY + +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.1102e-05, 1.8551e-05 + +% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) +MU_REF= 1.118e-05, 1.795e-05 + +% Sutherland Temperature Ref (273.15 K default value for AIR SI) +MU_T_REF= 300, 300 + +% Sutherland constant (110.4 default value for AIR SI) +SUTHERLAND_CONSTANT= 187.78, 129.4 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Navier-Stokes (no-slip), constant heat flux wall marker(s) (NONE = no marker) +% Format: ( marker name, constant heat flux (J/m^2), ... ) +MARKER_HEATFLUX= ( wall, 0.0 ) + +% Symmetry boundary marker(s) (NONE = no marker) +% Implementation identical to MARKER_EULER. +MARKER_SYM= ( axis ) + +% Axisymmetric simulation, only compressible flows (NO, YES) +AXISYMMETRIC= NO + +% List of inlet types for incompressible flows. List length must +% match number of inlet markers. Options: VELOCITY_INLET, PRESSURE_INLET. +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET VELOCITY_INLET + +% Inlet inlet scalar: (inlet marker, scalar initialisation value) +MARKER_INLET_SCALAR = (gas_inlet, 1.0, air_axial_inlet, 0.0, air_radial_inlet, 0.0) + +% Inc. Velocity: (inlet marker, temperature, velocity magnitude, flow_direction_x, +% flow_direction_y, flow_direction_z, ... ) where flow_direction is +% a unit vector. +MARKER_INLET= ( gas_inlet, 300, 12.353484, 1.0, 0.0, 0.0, air_axial_inlet, 300, 1.7857917, 1.0, 0.0, 0.0, air_radial_inlet, 300, 6.9766378, 0.0, -1.0, 0.0 ) + +% List of outlet types for incompressible flows. List length must +% match number of outlet markers. Options: PRESSURE_OUTLET, MASS_FLOW_OUTLET +INC_OUTLET_TYPE= PRESSURE_OUTLET + +% Outlet boundary marker(s) (NONE = no marker) +% Compressible: ( outlet marker, back pressure (static thermodynamic), ... ) +% Inc. Pressure: ( outlet marker, back pressure (static gauge in Pa), ... ) +% Inc. Mass Flow: ( outlet marker, mass flow target (kg/s), ... ) +MARKER_OUTLET= ( outlet, 0.0) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES + +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER = 60 + +% Reduction factor of the CFL coefficient in the scalar transport problem +CFL_REDUCTION_SCALAR = 1.0 + +% Reduction factor of the CFL coefficient in the turbulence problem +CFL_REDUCTION_TURB = 1.0 + +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO + +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) + +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% Number of total iterations +ITER= 50000 + +% Writing solution file frequency +OUTPUT_WRT_FREQ=1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver for implicit formulations (BCGSTAB, FGMRES) +LINEAR_SOLVER= FGMRES + +% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS) +LINEAR_SOLVER_PREC= ILU + +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-8 + +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 50 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES + +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_FLOW = NONE + +% Coefficient for the limiter (smooth regions) +VENKAT_LIMITER_COEFF= 10.0 + +% 2nd and 4th order artificial dissipation coefficients +JST_SENSOR_COEFF= ( 0.5, 0.02 ) + +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence criteria (CAUCHY, RESIDUAL) +CONV_CRITERIA= RESIDUAL + +CONV_FIELD = (RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE) + +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL = -10 + +% Start convergence criteria at iteration number +CONV_STARTITER= 10 + +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 + +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +SCREEN_OUTPUT = INNER_ITER WALL_TIME RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_PASSIVE_SCALAR LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_SCALAR LINSOL_RESIDUAL_SCALAR + +HISTORY_OUTPUT = RMS_RES MAX_RES FLOW_COEFF + +VOLUME_OUTPUT = DENSITY SOLUTION PRIMITIVE SOURCE RESIDUAL LOOKUP + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Mesh input file +MESH_FILENAME = C6_pneumatic_venturi_planar_mesh.cgns + +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= CGNS + +% Mesh output file +MESH_OUT_FILENAME= mesh_out.su2 + +% Restart flow input file +SOLUTION_FILENAME= initial_solution.dat + +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat + +% Output file format (PARAVIEW, PARAVIEW_BINARY, TECPLOT, STL) +OUTPUT_FILES = (RESTART, PARAVIEW, SURFACE_PARAVIEW) +TABULAR_FORMAT = CSV + +% Output file convergence history (w/o extension) +CONV_FILENAME= history + +% Output file restart flow +RESTART_FILENAME= restart.dat + +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat + +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow + +% Overwrite solution and visualization files or not +WRT_SOL_OVERWRITE = NO + +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint + +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat + +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow + +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint + +% Write the performance summary at the end of a calculation. +WRT_PERFORMANCE = YES \ No newline at end of file diff --git a/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets_restarted.cfg b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets_restarted.cfg new file mode 100755 index 000000000000..01b53d8851fc --- /dev/null +++ b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/C6_pneumatic_venturi_planar_velocityinlets_restarted.cfg @@ -0,0 +1,315 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Composition dependent incompressible flow with variable fluid properties through planar venturi % +% Author: Mark Heimgartner % +% Institution: Bosch Thermotechniek Deventer % +% Date: 07/07/2021 % +% File Version 7.1.1 "Blackbird", branch: feature_multicomp % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) +SOLVER = INC_RANS + +% Specify turbulent model (NONE, SA, SA_NEG, SST) +KIND_TURB_MODEL = SST + +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) +MATH_PROBLEM = DIRECT + +% Restart solution (NO, YES) +RESTART_SOL = YES + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL = VARIABLE + +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = YES + +% Initial density for incompressible flows +% (1.2886 kg/m^3 by default (air), 998.2 Kg/m^3 (water)) +INC_DENSITY_INIT = 1.1766 + +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT = (1.00, 0.0, 0.0 ) + +% Initial temperature for incompressible flows that include the +% energy equation (288.15 K by default). Value is ignored if +% INC_ENERGY_EQUATION is false. +INC_TEMPERATURE_INIT = 300.0 + +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +INC_NONDIM = DIMENSIONAL + +% -------------------- FLUID MODEL --------------------------------------- % +% +% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS, +% CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY, MUTATIONPP, SU2_NONEQ, MIXTURE_FLUID_MODEL) +FLUID_MODEL = MIXTURE_FLUID_MODEL + +% Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) +MOLECULAR_WEIGHT= 16.043, 28.965 + +% Specific heat at constant pressure, Cp (1004.703 J/kg*K (air)). +% Incompressible fluids with energy eqn. (CONSTANT_DENSITY, INC_IDEAL_GAS) and the heat equation. +SPECIFIC_HEAT_CP = 2224.43, 1009.39 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +% Laminar Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL, +% POLYNOMIAL_CONDUCTIVITY). +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY + +% Molecular Thermal Conductivity that would be constant (0.0257 by default) +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0258 + +% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL) +PRANDTL_LAM= 0.72, 0.72 + +% Definition of the turbulent thermal conductivity model for RANS +% (CONSTANT_PRANDTL_TURB by default, NONE). +TURBULENT_CONDUCTIVITY_MODEL= NONE + +% Turbulent Prandtl number (0.9 (air) by default) +PRANDTL_TURB= 0.90, 0.90 + +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_SCALAR= SCALAR_UPWIND + +% Specify scalar transport model (NONE, PASSIVE_SCALAR, PROGRESS_VARIABLE) +KIND_SCALAR_MODEL= PASSIVE_SCALAR + +% Specify mass diffusivity model (CONSTANT_DIFFUSIVITY, CONSTANT_SCHMIDT, FLAMELET, UNITY_LEWIS) +DIFFUSIVITY_MODEL = UNITY_LEWIS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the scalar transport equations. +MUSCL_SCALAR= YES + +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, +% BARTH_JESPERSEN, VAN_ALBADA_EDGE) +SLOPE_LIMITER_SCALAR = NONE + +% Time discretization for scalar transport problems (EULER_IMPLICIT) +TIME_DISCRE_SCALAR= EULER_IMPLICIT + +% Scalar clipping +SCALAR_CLIPPING= YES +SCALAR_INIT = 1 +SCALAR_CLIPPING_MIN= 0.0 +SCALAR_CLIPPING_MAX= 1.0 + +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_TURB= NO + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY, POLYNOMIAL_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY + +% Molecular Viscosity that would be constant (1.716E-5 by default) +MU_CONSTANT= 1.1102e-05, 1.8551e-05 + +% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) +MU_REF= 1.118e-05, 1.795e-05 + +% Sutherland Temperature Ref (273.15 K default value for AIR SI) +MU_T_REF= 300, 300 + +% Sutherland constant (110.4 default value for AIR SI) +SUTHERLAND_CONSTANT= 187.78, 129.4 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +% Navier-Stokes (no-slip), constant heat flux wall marker(s) (NONE = no marker) +% Format: ( marker name, constant heat flux (J/m^2), ... ) +MARKER_HEATFLUX= ( wall, 0.0 ) + +% Symmetry boundary marker(s) (NONE = no marker) +% Implementation identical to MARKER_EULER. +MARKER_SYM= ( axis ) + +% Axisymmetric simulation, only compressible flows (NO, YES) +AXISYMMETRIC= NO + +% List of inlet types for incompressible flows. List length must +% match number of inlet markers. Options: VELOCITY_INLET, PRESSURE_INLET. +INC_INLET_TYPE= VELOCITY_INLET VELOCITY_INLET VELOCITY_INLET + +% Inlet inlet scalar: (inlet marker, scalar initialisation value) +MARKER_INLET_SCALAR = (gas_inlet, 1.0, air_axial_inlet, 0.0, air_radial_inlet, 0.0) + +% Inc. Velocity: (inlet marker, temperature, velocity magnitude, flow_direction_x, +% flow_direction_y, flow_direction_z, ... ) where flow_direction is +% a unit vector. +MARKER_INLET= ( gas_inlet, 300, 12.353484, 1.0, 0.0, 0.0, air_axial_inlet, 300, 1.7857917, 1.0, 0.0, 0.0, air_radial_inlet, 300, 6.9766378, 0.0, -1.0, 0.0 ) + +% List of outlet types for incompressible flows. List length must +% match number of outlet markers. Options: PRESSURE_OUTLET, MASS_FLOW_OUTLET +INC_OUTLET_TYPE= PRESSURE_OUTLET + +% Outlet boundary marker(s) (NONE = no marker) +% Compressible: ( outlet marker, back pressure (static thermodynamic), ... ) +% Inc. Pressure: ( outlet marker, back pressure (static gauge in Pa), ... ) +% Inc. Mass Flow: ( outlet marker, mass flow target (kg/s), ... ) +MARKER_OUTLET= ( outlet, 0.0) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES + +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER = 60 + +% Reduction factor of the CFL coefficient in the scalar transport problem +CFL_REDUCTION_SCALAR = 1.0 + +% Reduction factor of the CFL coefficient in the turbulence problem +CFL_REDUCTION_TURB = 1.0 + +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO + +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) + +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% Number of total iterations +ITER= 50000 + +% Writing solution file frequency +OUTPUT_WRT_FREQ=1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver for implicit formulations (BCGSTAB, FGMRES) +LINEAR_SOLVER= FGMRES + +% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS) +LINEAR_SOLVER_PREC= ILU + +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-8 + +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 50 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, +% TURKEL_PREC, MSW) +CONV_NUM_METHOD_FLOW= FDS + +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES + +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_FLOW = NONE + +% Coefficient for the limiter (smooth regions) +VENKAT_LIMITER_COEFF= 10.0 + +% 2nd and 4th order artificial dissipation coefficients +JST_SENSOR_COEFF= ( 0.5, 0.02 ) + +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +% Convergence criteria (CAUCHY, RESIDUAL) +CONV_CRITERIA= RESIDUAL + +CONV_FIELD = (RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE) + +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL = -10 + +% Start convergence criteria at iteration number +CONV_STARTITER= 10 + +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 100 + +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 + +SCREEN_OUTPUT = INNER_ITER WALL_TIME RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_PASSIVE_SCALAR LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_SCALAR LINSOL_RESIDUAL_SCALAR + +HISTORY_OUTPUT = RMS_RES MAX_RES FLOW_COEFF + +VOLUME_OUTPUT = DENSITY SOLUTION PRIMITIVE SOURCE RESIDUAL LOOKUP + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +% Mesh input file +MESH_FILENAME = C6_pneumatic_venturi_planar_mesh.cgns + +% Mesh input file format (SU2, CGNS, NETCDF_ASCII) +MESH_FORMAT= CGNS + +% Mesh output file +MESH_OUT_FILENAME= mesh_out.su2 + +% Restart flow input file +SOLUTION_FILENAME= restart_2iter.dat + +% Restart adjoint input file +SOLUTION_ADJ_FILENAME= solution_adj.dat + +% Output file format (PARAVIEW, PARAVIEW_BINARY, TECPLOT, STL) +OUTPUT_FILES = (RESTART, PARAVIEW, SURFACE_PARAVIEW) +TABULAR_FORMAT = CSV + +% Output file convergence history (w/o extension) +CONV_FILENAME= history + +% Output file restart flow +%RESTART_FILENAME= restart.dat + +% Output file restart adjoint +RESTART_ADJ_FILENAME= restart_adj.dat + +% Output file flow (w/o extension) variables +VOLUME_FILENAME= flow + +% Overwrite solution and visualization files or not +WRT_SOL_OVERWRITE = NO + +% Output file adjoint (w/o extension) variables +VOLUME_ADJ_FILENAME= adjoint + +% Output objective function gradient (using continuous adjoint) +GRAD_OBJFUNC_FILENAME= of_grad.dat + +% Output file surface flow coefficient (w/o extension) +SURFACE_FILENAME= surface_flow + +% Output file surface adjoint coefficient (w/o extension) +SURFACE_ADJ_FILENAME= surface_adjoint + +% Write the performance summary at the end of a calculation. +WRT_PERFORMANCE = YES \ No newline at end of file diff --git a/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/README.md b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/README.md new file mode 100755 index 000000000000..3d7c4787c6b4 --- /dev/null +++ b/TestCases/incomp_rans/multicomponentflow_variablefluidproperties/README.md @@ -0,0 +1,38 @@ +# Test cases for multicomponent flow with composition dependent fluid properties + +## Introduction +Mixing laws for various fluid properties (viscosity, thermal conductivity, etc.) have been implemented in the SU2 software suite. This work is an extension of earlier work where equations for passive scalars and tabulated chemistry were implemented in the feature_flamelet branch. +The work in feature_multicomp deals with a transported scalar, i.e. there is a coupling between the scalar and flow solution. +Both test cases concern the non-reacting mixing flow of methane and air through a planar (2D) venturi. The mass fraction of methane, which is the input to the mixing laws, is obtained by solving an additional transported scalar equation. + +The test cases concern: +- velocity inlets + pressure outlet +- pressure inlets + target mass flow rate + +Except for the boundary markers all other configuration options are the same. +The mesh file is also the same for both test cases. + +## Files for this test case +Below is a list of files for this test case and explanation for each file. +- SU2 repository + - C6_pneumatic_venturi_planar_velocityinlets.cfg + - Configuration file. + - C6_pneumatic_venturi_planar_outlettargetmassflowrate.cfg + - Configuration file. +- TestCases repository + - C6_pneumatic_venturi_planar_mesh.cgns + - Mesh file. The mesh was made with Fluent meshing. + +## How to use composition dependent fluid properties framework + +### Config file +Enable composition dependent fluid properties using: +- FLUID_MODEL = MIXTURE_FLUID_MODEL + +For each fluid property specify an array of species’ fluid property e.g.: +- MOLECULAR_WEIGHT = 16.043, 28.965 + - Similar for specific heat, viscosity and thermal conductivity. + +Define mixture diffusivity: +- DIFFUSIVITY_MODEL = UNITY_LEWIS + - Mass diffusivity is computed based on unity Lewis assumption using mixture density, specific heat capacity (Cp) and thermal conductivity. \ No newline at end of file diff --git a/TestCases/multicomp_regression.py b/TestCases/multicomp_regression.py new file mode 100644 index 000000000000..49ed4aac9ecf --- /dev/null +++ b/TestCases/multicomp_regression.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python + +## \file parallel_regression.py +# \brief Python script for automated regression testing of SU2 examples +# \author A. Aranake, A. Campos, T. Economon, T. Lukaczyk, S. Padron +# \version 7.1.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2021, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +# make print(*args) function available in PY2.6+, does'nt work on PY < 2.6 +from __future__ import print_function + +import sys +from TestCase import TestCase + +def main(): + '''This program runs SU2 and ensures that the output matches specified values. + This will be used to do checks when code is pushed to github + to make sure nothing is broken. ''' + + test_list = [] + + ######################### + ## NEMO solver ### + ######################### + + # Multicomponent flow with variable fluid properties + multicompflow_variableprop_venturi_planar_velocityinlets_restarted = TestCase('multicompflow_variableprop_venturi_planar_velocityinlets_restarted') + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.cfg_dir = "incomp_rans/multicomponentflow_variablefluidproperties" + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.cfg_file = "C6_pneumatic_venturi_planar_velocityinlets_restarted.cfg" + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.test_iter = 0 + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.test_vals = [-5.202775, -4.104739, -3.882902, -5.143972, 50, -6.815594, 16, -8.370789] + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.su2_exec = "mpirun -n 2 SU2_CFD" + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.timeout = 1600 + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.new_output = True + multicompflow_variableprop_venturi_planar_velocityinlets_restarted.tol = 0.00001 + test_list.append(multicompflow_variableprop_venturi_planar_velocityinlets_restarted) + + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate = TestCase('multicompflow_variableprop_venturi_planar_outlettargetmassflowrate') + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.cfg_dir = "incomp_rans/multicomponentflow_variablefluidproperties" + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.cfg_file = "C6_pneumatic_venturi_planar_outlettargetmassflowrate.cfg" + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.test_iter = 5 + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.test_vals = [-5.880398, -5.129997, -5.188506, -6.137455, 50.000000, -7.907198, 18.000000, -8.406553] + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.su2_exec = "mpirun -n 2 SU2_CFD" + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.timeout = 1600 + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.new_output = True + multicompflow_variableprop_venturi_planar_outlettargetmassflowrate.tol = 0.00001 + test_list.append(multicompflow_variableprop_venturi_planar_outlettargetmassflowrate) + + multicompflow_variableprop_venturi_axisymmetric_velocityinlets = TestCase('multicompflow_variableprop_venturi_axisymmetric_velocityinlets') + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.cfg_dir = "incomp_rans/multicomponentflow_variablefluidproperties" + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.cfg_file = "C6_pneumatic_venturi_axisymmetric_velocityinlets.cfg" + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.test_iter = 5 + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.test_vals = [-5.658917, -4.475284, -4.442651, -5.811413, 20, -1.853076, 4, -3.547922] + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.su2_exec = "mpirun -n 2 SU2_CFD" + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.timeout = 1600 + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.new_output = True + multicompflow_variableprop_venturi_axisymmetric_velocityinlets.tol = 0.00001 + test_list.append(multicompflow_variableprop_venturi_axisymmetric_velocityinlets) + + ###################################### + ### RUN TESTS ### + ###################################### + + pass_list = [ test.run_test() for test in test_list ] + + # Tests summary + print('==================================================================') + print('Summary of the parallel tests') + print('python version:', sys.version) + for i, test in enumerate(test_list): + if (pass_list[i]): + print(' passed - %s'%test.tag) + else: + print('* FAILED - %s'%test.tag) + + if all(pass_list): + sys.exit(0) + else: + sys.exit(1) + # done + +if __name__ == '__main__': + main()