diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp index 30766e7b543..15dc1a83cc7 100644 --- a/SU2_CFD/include/fluid/CFluidModel.hpp +++ b/SU2_CFD/include/fluid/CFluidModel.hpp @@ -47,6 +47,7 @@ class CLookUpTable; class CFluidModel { protected: su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */ + su2double Enthalpy{0.0}; /*!<*\brief Enthalpy. */ su2double Entropy{0.0}; /*!< \brief Entropy. */ su2double Density{0.0}; /*!< \brief Density. */ su2double Pressure{0.0}; /*!< \brief Pressure. */ @@ -113,6 +114,11 @@ class CFluidModel { */ su2double GetStaticEnergy() const { return StaticEnergy; } + /*! + * \brief Get fluid enthalpy. + */ + su2double GetEnthalpy() const { return Enthalpy; } + /*! * \brief Get fluid density. */ @@ -186,6 +192,16 @@ class CFluidModel { return mass_diffusivity; } + /*! + * \brief Get heat diffusivity terms. + */ + virtual void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions = nullptr) {} + + /*! + * \brief Get gradient heat diffusivity terms. + */ + virtual void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions = nullptr) {} + /*! * \brief Get fluid pressure partial derivative. */ @@ -339,6 +355,13 @@ class CFluidModel { */ virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {} + /*! + * \brief Virtual member. + * \param[in] val_enthalpy - Enthalpy value at the point. + */ + virtual void ComputeTempFromEnthalpy(su2double val_enthalpy, su2double* val_temperature, + const su2double* val_scalars = nullptr) {} + /*! * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity. */ diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp index 206de692ae0..1db6307cd46 100644 --- a/SU2_CFD/include/fluid/CFluidScalar.hpp +++ b/SU2_CFD/include/fluid/CFluidScalar.hpp @@ -91,6 +91,11 @@ class CFluidScalar final : public CFluidModel { */ su2double ComputeMeanSpecificHeatCp(const su2double* val_scalars); + /*! + * \brief Compute Enthalpy given the temperature and scalars. + */ + su2double ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars); + /*! * \brief Compute gas constant for mixture. */ @@ -137,6 +142,22 @@ class CFluidScalar final : public CFluidModel { */ inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[ivar]; } + /*! + * \brief Get enthalpy diffusivity terms. + */ + void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) override; + + /*! + * \brief Get gradient enthalpy diffusivity terms. + */ + void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) override; + + /*! + * \brief Compute Temperature from Enthalpy and scalars. + */ + void ComputeTempFromEnthalpy(const su2double val_enthalpy, su2double* val_temperature, + const su2double* val_scalars) override; + /*! * \brief Set the Dimensionless State using Temperature. * \param[in] t - Temperature value at the point. diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 2e45916edb1..cec666bc72d 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -130,6 +130,10 @@ class CNumerics { const su2double *ScalarVar_i, /*!< \brief Vector of scalar variables at point i. */ *ScalarVar_j; /*!< \brief Vector of scalar variables at point j. */ + su2double + HeatFluxDiffusion; /*!< \brief Heat flux due to enthalpy diffusion for multicomponent. */ + su2double + Jac_HeatFluxDiffusion; /*!< \brief Heat flux jacobian due to enthalpy diffusion for multicomponent. */ const su2double *TransVar_i, /*!< \brief Vector of turbulent variables at point i. */ *TransVar_j; /*!< \brief Vector of turbulent variables at point j. */ @@ -188,6 +192,8 @@ class CNumerics { bool nemo; /*!< \brief Flag for NEMO problems */ + bool energy_multicomponent = false; /*!< \brief Flag for multicomponent and reacting flow */ + bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ public: @@ -751,6 +757,20 @@ class CNumerics { Diffusion_Coeff_j = val_diffusioncoeff_j; } + /*! + * \brief Set the heat flux due to enthalpy diffusion + * \param[in] val_heatfluxdiffusion - Value of the heat flux due to enthalpy diffusion. + */ + inline void SetHeatFluxDiffusion(su2double val_heatfluxdiffusion) { HeatFluxDiffusion = val_heatfluxdiffusion; } + + /*! + * \brief Set Jacobian of the heat flux due to enthalpy diffusion + * \param[in] val_jacheatfluxdiffusion - Value of the heat flux jacobian due to enthalpy diffusion. + */ + inline void SetJacHeatFluxDiffusion(su2double val_jac_heatfluxdiffusion) { + Jac_HeatFluxDiffusion = val_jac_heatfluxdiffusion; + } + /*! * \brief Set the laminar viscosity. * \param[in] val_laminar_viscosity_i - Value of the laminar viscosity at point i. diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index 654bc91e6ac..f3b997a9701 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -51,6 +51,8 @@ class CAvgGrad_Base : public CNumerics { *heat_flux_jac_i = nullptr, /*!< \brief Jacobian of the molecular + turbulent heat flux vector, projected onto the normal vector. */ **tau_jacobian_i = nullptr; /*!< \brief Jacobian of the viscous + turbulent stress tensor, projected onto the normal vector. */ su2double *Mean_PrimVar = nullptr; /*!< \brief Mean primitive variables. */ + su2double Mean_HeatFluxDiffusion; /*!< \brief Mean heat flux due to enthalpy diffusion for multicomponent flows. */ + su2double Mean_JacHeatFluxDiffusion; /*!< \brief Mean Jacobian heat flux due to enthalpy diffusion for multicomponent flows. */ const su2double *PrimVar_i = nullptr, *PrimVar_j = nullptr; /*!< \brief Primitives variables at point i and j. */ @@ -287,6 +289,7 @@ class CAvgGrad_Flow final : public CAvgGrad_Base { class CAvgGradInc_Flow final : public CAvgGrad_Base { private: su2double Mean_Thermal_Conductivity; /*!< \brief Mean value of the effective thermal conductivity. */ + su2double Mean_Heat_Capacity; /*!< \brief Mean value of the heat capacity. */ bool energy; /*!< \brief computation with the energy equation. */ /*! @@ -298,10 +301,10 @@ class CAvgGradInc_Flow final : public CAvgGrad_Base { * \param[in] val_gradprimvar - Gradient of the primitive variables. * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. * \param[in] val_thermal_conductivity - Thermal conductivity. + * \param[in] val_heatDiffusion - Heat diffusion */ - void GetViscousIncProjFlux(const su2double* const *val_gradprimvar, - const su2double *val_normal, - su2double val_thermal_conductivity); + void GetViscousIncProjFlux(const su2double* const* val_gradprimvar, const su2double* val_normal, + su2double val_thermal_conductivity, su2double val_heatFluxDiffusion); /*! * \brief Compute the projection of the viscous Jacobian matrices. diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index e4cb64e4b54..5fb9ef05be9 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -126,7 +126,8 @@ class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow { class CSourceIncAxisymmetric_Flow final : public CSourceBase_Flow { bool implicit, /*!< \brief Implicit calculation. */ viscous, /*!< \brief Viscous incompressible flows. */ - energy; /*!< \brief computation with the energy equation. */ + energy, /*!< \brief computation with the energy equation. */ + multicomponent; /*!< \brief multicomponent incompressible flows. */ public: /*! diff --git a/SU2_CFD/include/output/CFlowIncOutput.hpp b/SU2_CFD/include/output/CFlowIncOutput.hpp index fe7d7874427..f9c18d3d776 100644 --- a/SU2_CFD/include/output/CFlowIncOutput.hpp +++ b/SU2_CFD/include/output/CFlowIncOutput.hpp @@ -40,6 +40,7 @@ class CFlowIncOutput final: public CFlowOutput { private: TURB_MODEL turb_model; /*!< \brief The kind of turbulence model*/ bool heat; /*!< \brief Boolean indicating whether have a heat problem*/ + bool multicomponent; /*!< \brief Boolean indicating whether have a multicomponent problem*/ bool weakly_coupled_heat; /*!< \brief Boolean indicating whether have a weakly coupled heat equation*/ bool flamelet; /*!< \brief Boolean indicating whether we solve the flamelet equations */ unsigned short streamwisePeriodic; /*!< \brief Boolean indicating whether it is a streamwise periodic simulation. */ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 2aa31880bae..6f548e3ad31 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -69,7 +69,8 @@ class CFVMFlowSolverBase : public CSolver { su2double Mach_Inf = 0.0; /*!< \brief Mach number at the infinity. */ su2double Density_Inf = 0.0; /*!< \brief Density at the infinity. */ su2double Energy_Inf = 0.0; /*!< \brief Energy at the infinity. */ - su2double Temperature_Inf = 0.0; /*!< \brief Energy at the infinity. */ + su2double Temperature_Inf = 0.0; /*!< \brief Temperature at the infinity. */ + su2double Enthalpy_Inf = 0.0; /*!< \brief Enthalpy at the infinity. */ su2double Pressure_Inf = 0.0; /*!< \brief Pressure at the infinity. */ su2double* Velocity_Inf = nullptr; /*!< \brief Flow Velocity vector at the infinity. */ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index d7a634e8227..a4e2f884931 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -434,7 +434,7 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); - CVariable* turbNodes = nullptr; + CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); /*--- Points, coordinates and normal vector in edge ---*/ diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 990c9b7095d..50c52d4654a 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -191,6 +191,19 @@ class CIncEulerSolver : public CFVMFlowSolverBase struct CIndices { @@ -59,11 +59,11 @@ class CIncEulerVariable : public CFlowVariable { inline IndexType ThermalConductivity() const { return nDim+6; } inline IndexType CpTotal() const { return nDim+7; } inline IndexType CvTotal() const { return nDim+8; } + inline IndexType Enthalpy() const { return nDim + 9; } /*--- For compatible interface with NEMO. ---*/ inline IndexType SpeciesDensities() const { return std::numeric_limits::max(); } inline IndexType Temperature_ve() const { return std::numeric_limits::max(); } - inline IndexType Enthalpy() const { return std::numeric_limits::max(); } }; protected: @@ -121,6 +121,14 @@ class CIncEulerVariable : public CFlowVariable { return val_temperature <= 0.0; } + /*! + * \brief Set the value of the enthalpy for multicomponent incompressible flows with energy equation. + * \param[in] iPoint - Point index. + */ + inline void SetEnthalpy(unsigned long iPoint, su2double val_enthalpy) { + Primitive(iPoint, indices.Enthalpy()) = val_enthalpy; + } + /*! * \brief Set the value of the beta coeffient for incompressible flows. * \param[in] iPoint - Point index. @@ -153,6 +161,12 @@ class CIncEulerVariable : public CFlowVariable { */ inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint, indices.Temperature()); } + /*! + * \brief Get the enthalpy of the flow. + * \return Value of the enthalpy of the flow. + */ + inline su2double GetEnthalpy(unsigned long iPoint) const final { return Primitive(iPoint, indices.Enthalpy()); } + /*! * \brief Get the velocity of the flow. * \param[in] iDim - Index of the dimension. diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index 138026140b5..14d5e6bd7ac 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/SU2_CFD/include/variables/CIncNSVariable.hpp @@ -40,6 +40,7 @@ class CIncNSVariable final : public CIncEulerVariable { private: VectorType Tau_Wall; /*!< \brief Magnitude of the wall shear stress from a wall function. */ VectorType DES_LengthScale; + bool Energy_Multicomponent = false; public: /*! diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp index a57ddc7797e..b52c3b0b0b5 100644 --- a/SU2_CFD/src/fluid/CFluidScalar.cpp +++ b/SU2_CFD/src/fluid/CFluidScalar.cpp @@ -212,6 +212,35 @@ su2double CFluidScalar::ComputeMeanSpecificHeatCp(const su2double* val_scalars) return mean_cp; } +su2double CFluidScalar::ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars){ + su2double val_Enthalpy = Cp * (val_temperature - 298.15); + return val_Enthalpy; +} + +void CFluidScalar::ComputeTempFromEnthalpy(const su2double val_enthalpy, su2double* val_temperature, + const su2double* val_scalars) { + MassToMoleFractions(val_scalars); + su2double val_cp = ComputeMeanSpecificHeatCp(val_scalars); + *val_temperature = val_enthalpy / val_cp + 298.15; +} + +void CFluidScalar::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) { + for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) { + enthalpy_diffusions[iVar] = Density * + (specificHeat[iVar] * massDiffusivity[iVar] - + specificHeat[n_species_mixture - 1] * massDiffusivity[n_species_mixture - 1]) * + (Temperature - 298.15); + } +} + +void CFluidScalar::GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions){ + for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) { + grad_enthalpy_diffusions[iVar] = Density * + (specificHeat[iVar] * massDiffusivity[iVar] - + specificHeat[n_species_mixture - 1] * massDiffusivity[n_species_mixture - 1]); + } +} + void CFluidScalar::SetTDState_T(const su2double val_temperature, const su2double* val_scalars) { MassToMoleFractions(val_scalars); ComputeGasConstant(); @@ -219,6 +248,7 @@ void CFluidScalar::SetTDState_T(const su2double val_temperature, const su2double Density = Pressure_Thermodynamic / (Temperature * Gas_Constant); Cp = ComputeMeanSpecificHeatCp(val_scalars); Cv = Cp - Gas_Constant; + Enthalpy = ComputeEnthalpyFromT(Temperature, val_scalars); if (wilke) { Mu = WilkeViscosity(val_scalars); diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 2c36556f808..f78e4b60288 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -54,6 +54,7 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Prandtl_Lam = config->GetPrandtl_Lam(); Prandtl_Turb = config->GetPrandtl_Turb(); Gas_Constant = config->GetGas_ConstantND(); + energy_multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) && (config->GetEnergy_Equation())); tau = new su2double* [nDim]; for (iDim = 0; iDim < nDim; iDim++) @@ -278,10 +279,20 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou val_Proj_Jac_Tensor[2][2] = val_scale*((*val_density)*(proj_vel + val_normal[1]*val_velocity[1])); val_Proj_Jac_Tensor[2][3] = val_scale*((*val_dRhodT)*val_velocity[1]*proj_vel); - val_Proj_Jac_Tensor[3][0] = val_scale*((*val_cp)*(*val_temperature)*proj_vel/(*val_betainc2)); - val_Proj_Jac_Tensor[3][1] = val_scale*((*val_cp)*(*val_temperature)*val_normal[0]*(*val_density)); - val_Proj_Jac_Tensor[3][2] = val_scale*((*val_cp)*(*val_temperature)*val_normal[1]*(*val_density)); - val_Proj_Jac_Tensor[3][3] = val_scale*((*val_cp)*((*val_temperature)*(*val_dRhodT) + (*val_density))*proj_vel); + if (energy_multicomponent) { + val_Proj_Jac_Tensor[3][0] = val_scale * ((*val_temperature) * proj_vel / (*val_betainc2)); + val_Proj_Jac_Tensor[3][1] = val_scale * ((*val_temperature) * val_normal[0] * (*val_density)); + val_Proj_Jac_Tensor[3][2] = val_scale * ((*val_temperature) * val_normal[1] * (*val_density)); + val_Proj_Jac_Tensor[3][3] = + val_scale * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel; + + } else { + val_Proj_Jac_Tensor[3][0] = val_scale * ((*val_cp) * (*val_temperature) * proj_vel / (*val_betainc2)); + val_Proj_Jac_Tensor[3][1] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[0] * (*val_density)); + val_Proj_Jac_Tensor[3][2] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[1] * (*val_density)); + val_Proj_Jac_Tensor[3][3] = + val_scale * ((*val_cp) * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel); + } } else { @@ -309,12 +320,21 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou val_Proj_Jac_Tensor[3][3] = val_scale*((*val_density)*(proj_vel + val_normal[2]*val_velocity[2])); val_Proj_Jac_Tensor[3][4] = val_scale*((*val_dRhodT)*val_velocity[2]*proj_vel); - val_Proj_Jac_Tensor[4][0] = val_scale*((*val_cp)*(*val_temperature)*proj_vel/(*val_betainc2)); - val_Proj_Jac_Tensor[4][1] = val_scale*((*val_cp)*(*val_temperature)*val_normal[0]*(*val_density)); - val_Proj_Jac_Tensor[4][2] = val_scale*((*val_cp)*(*val_temperature)*val_normal[1]*(*val_density)); - val_Proj_Jac_Tensor[4][3] = val_scale*((*val_cp)*(*val_temperature)*val_normal[2]*(*val_density)); - val_Proj_Jac_Tensor[4][4] = val_scale*((*val_cp)*((*val_temperature)*(*val_dRhodT) + (*val_density))*proj_vel); - + if (energy_multicomponent) { + val_Proj_Jac_Tensor[4][0] = val_scale * ((*val_temperature) * proj_vel / (*val_betainc2)); + val_Proj_Jac_Tensor[4][1] = val_scale * ((*val_temperature) * val_normal[0] * (*val_density)); + val_Proj_Jac_Tensor[4][2] = val_scale * ((*val_temperature) * val_normal[1] * (*val_density)); + val_Proj_Jac_Tensor[4][3] = val_scale * ((*val_temperature) * val_normal[2] * (*val_density)); + val_Proj_Jac_Tensor[4][4] = + val_scale * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel; + } else { + val_Proj_Jac_Tensor[4][0] = val_scale * ((*val_cp) * (*val_temperature) * proj_vel / (*val_betainc2)); + val_Proj_Jac_Tensor[4][1] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[0] * (*val_density)); + val_Proj_Jac_Tensor[4][2] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[1] * (*val_density)); + val_Proj_Jac_Tensor[4][3] = val_scale * ((*val_cp) * (*val_temperature) * val_normal[2] * (*val_density)); + val_Proj_Jac_Tensor[4][4] = + val_scale * ((*val_cp) * ((*val_temperature) * (*val_dRhodT) + (*val_density)) * proj_vel); + } } AD::EndPassive(wasActive); } @@ -328,7 +348,11 @@ void CNumerics::GetPreconditioner(const su2double *val_density, const su2double val_Precon[0][0] = 1.0/(*val_betainc2); for (iDim = 0; iDim < nDim; iDim++) val_Precon[iDim+1][0] = val_velocity[iDim]/(*val_betainc2); - val_Precon[nDim+1][0] = (*val_cp)*(*val_temperature)/(*val_betainc2); + if (energy_multicomponent){ + val_Precon[nDim+1][0] = (*val_temperature)/(*val_betainc2); + }else{ + val_Precon[nDim+1][0] = (*val_cp)*(*val_temperature)/(*val_betainc2); + } for (jDim = 0; jDim < nDim; jDim++) { val_Precon[0][jDim+1] = 0.0; @@ -342,7 +366,11 @@ void CNumerics::GetPreconditioner(const su2double *val_density, const su2double val_Precon[0][nDim+1] = (*val_drhodt); for (iDim = 0; iDim < nDim; iDim++) val_Precon[iDim+1][nDim+1] = val_velocity[iDim]*(*val_drhodt); - val_Precon[nDim+1][nDim+1] = (*val_cp)*((*val_drhodt)*(*val_temperature) + (*val_density)); + if (energy_multicomponent){ + val_Precon[nDim+1][nDim+1] = (*val_drhodt)*(*val_temperature) + (*val_density); + }else{ + val_Precon[nDim+1][nDim+1] = (*val_cp)*((*val_drhodt)*(*val_temperature) + (*val_density)); + } } diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp index ae48e59f336..115a966896b 100644 --- a/SU2_CFD/src/numerics/flow/convection/fds.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp @@ -94,7 +94,7 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config su2double ProjGridVel = 0.0; AD::StartPreacc(); - AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(V_i, nDim+10); AD::SetPreaccIn(V_j, nDim+10); AD::SetPreaccIn(Normal, nDim); if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); @@ -118,7 +118,13 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config DensityInc_i = V_i[nDim+2]; DensityInc_j = V_j[nDim+2]; BetaInc2_i = V_i[nDim+3]; BetaInc2_j = V_j[nDim+3]; Cp_i = V_i[nDim+7]; Cp_j = V_j[nDim+7]; - Enthalpy_i = Cp_i*Temperature_i; Enthalpy_j = Cp_j*Temperature_j; + if (energy_multicomponent) { + Enthalpy_i = V_i[nDim + 9]; + Enthalpy_j = V_j[nDim + 9]; + } else { + Enthalpy_i = Cp_i * Temperature_i; + Enthalpy_j = Cp_j * Temperature_j; + } ProjVelocity = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -157,9 +163,15 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config MeandRhodT = 0.0; dRhodT_i = 0.0; dRhodT_j = 0.0; if (variable_density) { - MeandRhodT = -MeanDensity/MeanTemperature; - dRhodT_i = -DensityInc_i/Temperature_i; - dRhodT_j = -DensityInc_j/Temperature_j; + if (energy_multicomponent) { + MeandRhodT = -MeanDensity / (MeanCp * MeanTemperature); + dRhodT_i = -DensityInc_i / (Cp_i * Temperature_i); + dRhodT_j = -DensityInc_j / (Cp_j * Temperature_j); + } else { + MeandRhodT = -MeanDensity / MeanTemperature; + dRhodT_i = -DensityInc_i / Temperature_i; + dRhodT_j = -DensityInc_j / Temperature_j; + } } /*--- Compute ProjFlux_i ---*/ @@ -192,8 +204,11 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Lambda[iVar] = fabs(Lambda[iVar]); /*--- Build the preconditioning matrix using mean values ---*/ - - GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanTemperature, &MeandRhodT, Precon); + if (energy_multicomponent) { + GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanEnthalpy, &MeandRhodT, Precon); + } else { + GetPreconditioner(&MeanDensity, MeanVelocity, &MeanBetaInc2, &MeanCp, &MeanTemperature, &MeandRhodT, Precon); + } /*--- Build the absolute value of the preconditioned Jacobian, i.e., |A_precon| = P x |Lambda| x inv(P), where P diagonalizes the matrix @@ -206,13 +221,26 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Diff_V[0] = Pressure_j - Pressure_i; for (iDim = 0; iDim < nDim; iDim++) Diff_V[iDim+1] = Velocity_j[iDim] - Velocity_i[iDim]; - Diff_V[nDim+1] = Temperature_j - Temperature_i; + if (energy_multicomponent){ + Diff_V[nDim+1] = Enthalpy_j - Enthalpy_i; + }else{ + Diff_V[nDim+1] = Temperature_j - Temperature_i; + } /*--- Build the inviscid Jacobian w.r.t. the primitive variables ---*/ if (implicit) { - GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &Temperature_i, &dRhodT_i, Normal, 0.5, Jacobian_i); - GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &Temperature_j, &dRhodT_j, Normal, 0.5, Jacobian_j); + if (energy_multicomponent) { + GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &Enthalpy_i, &dRhodT_i, Normal, 0.5, + Jacobian_i); + GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &Enthalpy_j, &dRhodT_j, Normal, 0.5, + Jacobian_j); + } else { + GetInviscidIncProjJac(&DensityInc_i, Velocity_i, &BetaInc2_i, &Cp_i, &Temperature_i, &dRhodT_i, Normal, 0.5, + Jacobian_i); + GetInviscidIncProjJac(&DensityInc_j, Velocity_j, &BetaInc2_j, &Cp_j, &Temperature_j, &dRhodT_j, Normal, 0.5, + Jacobian_j); + } } /*--- Compute dissipation as Precon x |A_precon| x dV. If implicit, @@ -258,8 +286,13 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; } - Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; - Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; + if (energy_multicomponent){ + Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i; + Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j; + } else { + Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; + Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; + } } } } diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index a2cdb893533..ddef652b8db 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -583,6 +583,10 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Laminar_Viscosity_i = V_i[nDim+4]; Laminar_Viscosity_j = V_j[nDim+4]; Eddy_Viscosity_i = V_i[nDim+5]; Eddy_Viscosity_j = V_j[nDim+5]; Thermal_Conductivity_i = V_i[nDim+6]; Thermal_Conductivity_j = V_j[nDim+6]; + if (energy_multicomponent) { + Cp_i = V_i[nDim + 7]; + Cp_j = V_j[nDim + 7]; + } /*--- Mean transport properties ---*/ @@ -590,6 +594,9 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); + if (energy_multicomponent) { + Mean_Heat_Capacity = 0.5 * (Cp_i + Cp_j); + } /*--- Mean gradient approximation ---*/ @@ -623,7 +630,12 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); - GetViscousIncProjFlux(Mean_GradPrimVar, Normal, Mean_Thermal_Conductivity); + if (energy_multicomponent) { + Mean_HeatFluxDiffusion = HeatFluxDiffusion; + Mean_JacHeatFluxDiffusion = Jac_HeatFluxDiffusion; + } + + GetViscousIncProjFlux(Mean_GradPrimVar, Normal, Mean_Thermal_Conductivity, Mean_HeatFluxDiffusion); /*--- Implicit part ---*/ @@ -649,8 +661,15 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi proj_vector_ij += (Coord_j[iDim]-Coord_i[iDim])*Normal[iDim]; } proj_vector_ij = proj_vector_ij/dist_ij_2; - Jacobian_i[nDim+1][nDim+1] = -Mean_Thermal_Conductivity*proj_vector_ij; - Jacobian_j[nDim+1][nDim+1] = Mean_Thermal_Conductivity*proj_vector_ij; + if (energy_multicomponent){ + Jacobian_i[nDim + 1][nDim + 1] = + -(Mean_Thermal_Conductivity * proj_vector_ij + Mean_JacHeatFluxDiffusion) / Mean_Heat_Capacity; + Jacobian_j[nDim + 1][nDim + 1] = + (Mean_Thermal_Conductivity * proj_vector_ij + Mean_JacHeatFluxDiffusion) / Mean_Heat_Capacity; + } else { + Jacobian_i[nDim + 1][nDim + 1] = -Mean_Thermal_Conductivity * proj_vector_ij; + Jacobian_j[nDim + 1][nDim + 1] = Mean_Thermal_Conductivity * proj_vector_ij; + } } } @@ -677,7 +696,8 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi void CAvgGradInc_Flow::GetViscousIncProjFlux(const su2double* const *val_gradprimvar, const su2double *val_normal, - su2double val_thermal_conductivity) { + su2double val_thermal_conductivity, + su2double val_heatDiffusion) { /*--- Gradient of primitive variables -> [Pressure vel_x vel_y vel_z Temperature] ---*/ @@ -721,7 +741,9 @@ void CAvgGradInc_Flow::GetViscousIncProjFlux(const su2double* const *val_gradpri for (unsigned short iDim = 0; iDim < nDim; iDim++) Proj_Flux_Tensor[iVar] += Flux_Tensor[iVar][iDim] * val_normal[iDim]; } - + if (energy_multicomponent) { + Proj_Flux_Tensor[nVar - 1] += val_heatDiffusion; + } } void CAvgGradInc_Flow::GetViscousIncProjJacs(su2double val_dS, diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 6f8359a5f58..4aa072be466 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -245,7 +245,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); energy = config->GetEnergy_Equation(); viscous = config->GetViscous(); - + multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); } CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) { @@ -264,7 +264,11 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo DensityInc_i = V_i[nDim+2]; BetaInc2_i = V_i[nDim+3]; Cp_i = V_i[nDim+7]; - Enthalpy_i = Cp_i*Temp_i; + if (multicomponent && energy) { + Enthalpy_i = V_i[nDim + 9]; + } else { + Enthalpy_i = Cp_i * Temp_i; + } for (iDim = 0; iDim < nDim; iDim++) Velocity_i[iDim] = V_i[iDim+1]; @@ -296,7 +300,11 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo jacobian[3][0] = 0.0; jacobian[3][1] = 0.0; jacobian[3][2] = Enthalpy_i; - jacobian[3][3] = Cp_i*Velocity_i[1]; + if (multicomponent && energy) { + jacobian[3][3] = (1.0 - Enthalpy_i / (Cp_i * Temp_i)) * Velocity_i[1]; + } else { + jacobian[3][3] = Cp_i * Velocity_i[1]; + } for (iVar=0; iVar < nVar; iVar++) for (jVar=0; jVar < nVar; jVar++) @@ -326,8 +334,11 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo 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]; - + residual[3] -= Volume * yinv * Thermal_Conductivity_i * PrimVar_Grad_i[nDim + 1][1]; + if (multicomponent && energy) { + residual[3] -= Volume * yinv * HeatFluxDiffusion; + if (implicit) jacobian[3][3] -= Volume * yinv * Jac_HeatFluxDiffusion / Cp_i; + } } } else { diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 1f4c0809d1b..b7e05371bd0 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -37,6 +37,8 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp heat = config->GetEnergy_Equation(); + multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); + weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); flamelet = (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET); @@ -108,7 +110,13 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Root-mean square residual of the velocity z-component. if (nDim == 3) AddHistoryOutput("RMS_VELOCITY-Z", "rms[W]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the velocity z-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - if (heat || weakly_coupled_heat) AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat){ + if (multicomponent){ + AddHistoryOutput("RMS_ENTHALPY", "rms[h]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the enthalpy.", HistoryFieldType::RESIDUAL); + } else { + AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the temperature.", HistoryFieldType::RESIDUAL); + } + } AddHistoryOutputFields_ScalarRMS_RES(config); @@ -127,8 +135,13 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ if (nDim == 3) AddHistoryOutput("MAX_VELOCITY-Z", "max[W]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the velocity z-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - if (heat || weakly_coupled_heat) - AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Root-mean square residual of the temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat) { + if (multicomponent){ + AddHistoryOutput("MAX_ENTHALPY", "max[h]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the enthalpy.", HistoryFieldType::RESIDUAL); + } else { + AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature.", HistoryFieldType::RESIDUAL); + } + } AddHistoryOutputFields_ScalarMAX_RES(config); /// END_GROUP @@ -144,8 +157,13 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ if (nDim == 3) AddHistoryOutput("BGS_VELOCITY-Z", "bgs[W]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the velocity z-component.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the temperature. - if (heat || weakly_coupled_heat) - AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the temperature.", HistoryFieldType::RESIDUAL); + if (heat || weakly_coupled_heat) { + if (multicomponent){ + AddHistoryOutput("BGS_ENTHALPY", "bgs[h]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the enthalpy.", HistoryFieldType::RESIDUAL); + } else { + AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the temperature.", HistoryFieldType::RESIDUAL); + } + } AddHistoryOutputFields_ScalarBGS_RES(config); @@ -232,10 +250,18 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv if (heat) { SetHeatCoefficients(config, flow_solver); SetHistoryOutputValue("AVG_TEMPERATURE", flow_solver->GetTotal_AvgTemperature()); - SetHistoryOutputValue("RMS_TEMPERATURE", log10(flow_solver->GetRes_RMS(nDim + 1))); - SetHistoryOutputValue("MAX_TEMPERATURE", log10(flow_solver->GetRes_Max(nDim + 1))); - if (multiZone) { - SetHistoryOutputValue("BGS_TEMPERATURE", log10(flow_solver->GetRes_BGS(nDim + 1))); + if (multicomponent){ + SetHistoryOutputValue("RMS_ENTHALPY", log10(flow_solver->GetRes_RMS(nDim + 1))); + SetHistoryOutputValue("MAX_ENTHALPY", log10(flow_solver->GetRes_Max(nDim + 1))); + if (multiZone) { + SetHistoryOutputValue("BGS_ENTHALPY", log10(flow_solver->GetRes_BGS(nDim + 1))); + } + } else { + SetHistoryOutputValue("RMS_TEMPERATURE", log10(flow_solver->GetRes_RMS(nDim + 1))); + SetHistoryOutputValue("MAX_TEMPERATURE", log10(flow_solver->GetRes_Max(nDim + 1))); + if (multiZone) { + SetHistoryOutputValue("BGS_TEMPERATURE", log10(flow_solver->GetRes_BGS(nDim + 1))); + } } } @@ -295,8 +321,10 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector"); - if (heat || weakly_coupled_heat || flamelet) - AddVolumeOutput("TEMPERATURE", "Temperature","SOLUTION", "Temperature"); + if (heat || weakly_coupled_heat || flamelet){ + if (multicomponent) AddVolumeOutput("ENTHALPY", "Enthalpy", "SOLUTION", "Enthalpy"); + else AddVolumeOutput("TEMPERATURE", "Temperature","SOLUTION", "Temperature"); + } SetVolumeOutputFieldsScalarSolution(config); @@ -320,6 +348,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("LAMINAR_VISCOSITY", "Laminar_Viscosity", "PRIMITIVE", "Laminar viscosity"); AddVolumeOutput("HEAT_CAPACITY", "Heat_Capacity", "PRIMITIVE", "Heat capacity"); AddVolumeOutput("THERMAL_CONDUCTIVITY", "Thermal_Conductivity", "PRIMITIVE", "Thermal conductivity"); + if (heat && multicomponent) AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature"); AddVolumeOutput("SKIN_FRICTION-X", "Skin_Friction_Coefficient_x", "PRIMITIVE", "x-component of the skin friction vector"); AddVolumeOutput("SKIN_FRICTION-Y", "Skin_Friction_Coefficient_y", "PRIMITIVE", "y-component of the skin friction vector"); @@ -343,9 +372,10 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component"); if (nDim == 3) AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component"); - if (config->GetEnergy_Equation()) - AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); - + if (config->GetEnergy_Equation()){ + if (multicomponent) AddVolumeOutput("RES_ENTHALPY", "Residual_Enthalpy", "RESIDUAL", "Residual of the enthalpy"); + else AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); + } SetVolumeOutputFieldsScalarResidual(config); if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) { @@ -399,7 +429,10 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve if (nDim == 3) SetVolumeOutputValue("VELOCITY-Z", iPoint, Node_Flow->GetSolution(iPoint, 3)); - if (heat || flamelet) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetSolution(iPoint, nDim+1)); + if (heat || flamelet) { + if (multicomponent) SetVolumeOutputValue("ENTHALPY", iPoint, Node_Flow->GetSolution(iPoint, nDim+1)); + else SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetSolution(iPoint, nDim+1)); + } if (weakly_coupled_heat) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); // Radiation solver @@ -423,6 +456,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint)); SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); SetVolumeOutputValue("THERMAL_CONDUCTIVITY", iPoint, Node_Flow->GetThermalConductivity(iPoint)); + if (multicomponent && heat) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint)); } SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0)); @@ -430,8 +464,10 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("RES_VELOCITY-Y", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 2)); if (nDim == 3) SetVolumeOutputValue("RES_VELOCITY-Z", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 3)); - if (config->GetEnergy_Equation()) - SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, nDim+1)); + if (config->GetEnergy_Equation()) { + if (multicomponent) SetVolumeOutputValue("RES_ENTHALPY", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, nDim+1)); + else SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, nDim+1)); + } if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) { SetVolumeOutputValue("LIMITER_PRESSURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 0)); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 3c651269b02..979d911f3b2 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -150,6 +150,7 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry const bool compressible = config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE; const bool incompressible = config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE; const bool energy = config->GetEnergy_Equation(); + const bool multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE; const bool streamwisePeriodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); const bool species = config->GetKind_Species_Model() == SPECIES_MODEL::SPECIES_TRANSPORT; const auto nSpecies = config->GetnSpecies(); @@ -241,7 +242,11 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry sqrt(config->GetBulk_Modulus()/(flow_nodes->GetDensity(iPoint))); } Temperature = flow_nodes->GetTemperature(iPoint); - Enthalpy = flow_nodes->GetSpecificHeatCp(iPoint)*Temperature; + if (energy && multicomponent) { + Enthalpy = flow_nodes->GetEnthalpy(iPoint); + } else { + Enthalpy = flow_nodes->GetSpecificHeatCp(iPoint) * Temperature; + } TotalTemperature = Temperature + 0.5*Velocity2/flow_nodes->GetSpecificHeatCp(iPoint); TotalPressure = Pressure + 0.5*Density*Velocity2; } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 6aaf09a5183..aa6bea96a0d 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -57,6 +57,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned bool time_stepping = config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING; bool adjoint = (config->GetContinuous_Adjoint()) || (config->GetDiscrete_Adjoint()); const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED; + const bool energy_multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) && config->GetEnergy_Equation()); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); @@ -119,7 +120,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned /*--- Make sure to align the sizes with the constructor of CIncEulerVariable. ---*/ nVar = nDim + 2; - nPrimVar = nDim + 9; + nPrimVar = nDim + 10; /*--- Centered schemes only need gradients for viscous fluxes (T and v, but we need also to include P). ---*/ nPrimVarGrad = nDim + (centered ? 2 : 4); @@ -178,6 +179,12 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned Pressure_Inf = config->GetPressure_FreeStreamND(); Velocity_Inf = config->GetVelocity_FreeStreamND(); Temperature_Inf = config->GetTemperature_FreeStreamND(); + if (energy_multicomponent){ + CFluidModel *auxFluidModel = new CFluidScalar(config->GetPressure_Thermodynamic(), config); + const su2double *scalar_init = config->GetSpecies_Init(); + auxFluidModel->SetTDState_T(Temperature_Inf,scalar_init); // compute total enthalpy from temperature + Enthalpy_Inf = auxFluidModel->GetEnthalpy(); + } /*--- Initialize the secondary values for direct derivative approxiations ---*/ @@ -208,7 +215,11 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned /*--- Initialize the solution to the far-field state everywhere. ---*/ if (navier_stokes) { - nodes = new CIncNSVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config); + if (energy_multicomponent){ + nodes = new CIncNSVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); + }else{ + nodes = new CIncNSVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config); + } } else { nodes = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config); } @@ -1227,6 +1238,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const bool bounded_scalar = config->GetBounded_Scalar(); + const bool energy_multicomponent = + ((config->GetEnergy_Equation()) && (config->GetKind_FluidModel() == FLUID_MIXTURE)); /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of * variables, otherwise switch to the faster adjoint evaluation mode. ---*/ @@ -1304,6 +1317,31 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont Primitive_i[iVar] = V_i[iVar]; Primitive_j[iVar] = V_j[iVar]; } + if (energy_multicomponent) { + su2double Project_Grad_Enthalpy_i = 0.0; + su2double Project_Grad_Enthalpy_j = 0.0; + + for (iDim = 0; iDim < nDim; iDim++) { + Project_Grad_Enthalpy_i += Vector_ij[iDim] * nodes->GetAuxVarGradient(iPoint, 1, iDim); + Project_Grad_Enthalpy_j -= Vector_ij[iDim] * nodes->GetAuxVarGradient(jPoint, 1, iDim); + } + su2double lim_i = 1.0; + su2double lim_j = 1.0; + if (van_albada) { + su2double V_ij = V_j[nDim + 9] - V_i[nDim + 9]; + lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_Enthalpy_i, V_ij, EPS); + lim_j = LimiterHelpers<>::vanAlbadaFunction(-Project_Grad_Enthalpy_j, V_ij, EPS); + } else if (limiter) { + lim_i = 1.0; // nodes->GetLimiter_Primitive(iPoint, iVar); + lim_j = 1.0; // nodes->GetLimiter_Primitive(jPoint, iVar); + } + Primitive_i[nDim + 9] = V_i[nDim + 9] + lim_i * Project_Grad_Enthalpy_i; + Primitive_j[nDim + 9] = V_j[nDim + 9] + lim_j * Project_Grad_Enthalpy_j; + const su2double* scalar_i = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint); + const su2double* scalar_j = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(jPoint); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, scalar_i); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, scalar_j); + } /*--- Check for non-physical solutions after reconstruction. If found, use the cell-average value of the solution. This results in a locally @@ -1375,6 +1413,19 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont FinalizeResidualComputation(geometry, pausePreacc, counter_local, config); } +void CIncEulerSolver::ComputeConsistentExtrapolation(CFluidModel* fluidModel, unsigned short nDim, su2double* primitive, + const su2double* scalar) { + const CIncEulerVariable::CIndices prim_idx(nDim, 0); + const su2double enthalpy = primitive[prim_idx.Enthalpy()]; + su2double temperature = primitive[prim_idx.Temperature()]; + fluidModel->ComputeTempFromEnthalpy(enthalpy, &temperature, scalar); + + fluidModel->SetTDState_T(temperature, scalar); + + primitive[prim_idx.Temperature()] = temperature; + primitive[prim_idx.Density()] = fluidModel->GetDensity(); +} + void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { @@ -1396,6 +1447,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont const bool energy = config->GetEnergy_Equation(); const bool streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); const bool streamwise_periodic_temperature = config->GetStreamwise_Periodic_Temperature(); + const bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); AD::StartNoSharedReading(); @@ -1578,6 +1630,35 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); + if(multicomponent && energy){ + /*--- retrieve number of species that are solved and set maximum static array ---*/ + int n_species = config->GetnSpecies(); + static constexpr size_t MAXNVAR_SPECIES = 20UL; + /*--- Obtain fluid model for computing the enthalpy diffusion terms. ---*/ + CFluidModel* FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); + /*--- retrieve species gradient needed for multicomponent. ---*/ + CMatrixView Species_Grad_i = solver_container[SPECIES_SOL]->GetNodes()->GetGradient(iPoint); + /*--- Set thermodynamic state. ---*/ + FluidModel->SetTDState_T(nodes->GetTemperature(iPoint),solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint)); + /*--- Get enthalpy diffusion terms and its gradients(for implicit). ---*/ + su2double EnthalpyDiffusion_i[MAXNVAR_SPECIES]{0.0}; + su2double GradEnthalpyDiffusion_i[MAXNVAR_SPECIES]{0.0}; + FluidModel->GetEnthalpyDiffusivity(EnthalpyDiffusion_i); + if (implicit) FluidModel->GetGradEnthalpyDiffusivity(GradEnthalpyDiffusion_i); + /*--- Compute Enthalpy diffusion flux and its jacobian (for implicit iterations) ---*/ + su2double flux_enthalpy_diffusion = 0.0; + su2double jac_flux_enthalpy_diffusion = 0.0; + for (int i_species = 0; i_species < n_species; i_species++) { + flux_enthalpy_diffusion += EnthalpyDiffusion_i[i_species]* Species_Grad_i[i_species][1]; + if (implicit) + jac_flux_enthalpy_diffusion += GradEnthalpyDiffusion_i[i_species] * Species_Grad_i[i_species][1]; + } + + /*--- Set heat flux and jacobian (for implicit) due to enthalpy diffusion ---*/ + + numerics->SetHeatFluxDiffusion(flux_enthalpy_diffusion); + if (implicit) numerics->SetJacHeatFluxDiffusion(jac_flux_enthalpy_diffusion); + } } /*--- Compute Source term Residual ---*/ @@ -1977,12 +2058,13 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo unsigned short iDim, jDim, iVar, jVar; - su2double BetaInc2, Density, dRhodT, Temperature, oneOverCp, Cp; + su2double BetaInc2, Density, dRhodT, Temperature, oneOverCp, Cp, Enthalpy; su2double Velocity[MAXNDIM] = {0.0}; bool variable_density = (config->GetVariable_Density_Model()); bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); bool energy = config->GetEnergy_Equation(); + bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); /*--- Access the primitive variables at this node. ---*/ @@ -1991,6 +2073,7 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo Cp = nodes->GetSpecificHeatCp(iPoint); oneOverCp = 1.0/Cp; Temperature = nodes->GetTemperature(iPoint); + if (energy && multicomponent) Enthalpy = nodes->GetEnthalpy(iPoint); for (iDim = 0; iDim < nDim; iDim++) Velocity[iDim] = nodes->GetVelocity(iPoint,iDim); @@ -2000,7 +2083,11 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo law, but in the future, dRhodT should be in the fluid model. ---*/ if (variable_density) { - dRhodT = -Density/Temperature; + if (multicomponent && energy){ + dRhodT = -Density / (Cp * Temperature); + } else { + dRhodT = -Density/Temperature; + } } else { dRhodT = 0.0; } @@ -2017,8 +2104,15 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim+1][0] = Velocity[iDim]/BetaInc2; - if (energy) Preconditioner[nDim+1][0] = Cp*Temperature/BetaInc2; - else Preconditioner[nDim+1][0] = 0.0; + if (energy) { + if (multicomponent) { + Preconditioner[nDim + 1][0] = Enthalpy / BetaInc2; + } else { + Preconditioner[nDim + 1][0] = Cp * Temperature / BetaInc2; + } + } else { + Preconditioner[nDim + 1][0] = 0.0; + } for (jDim = 0; jDim < nDim; jDim++) { Preconditioner[0][jDim+1] = 0.0; @@ -2033,46 +2127,72 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim+1][nDim+1] = Velocity[iDim]*dRhodT; - if (energy) Preconditioner[nDim+1][nDim+1] = Cp*(dRhodT*Temperature + Density); - else Preconditioner[nDim+1][nDim+1] = 1.0; + if (energy) { + if (multicomponent) { + Preconditioner[nDim + 1][nDim + 1] = dRhodT * Enthalpy + Density; + } else { + Preconditioner[nDim + 1][nDim + 1] = Cp * (dRhodT * Temperature + Density); + } + } else { + Preconditioner[nDim + 1][nDim + 1] = 1.0; + } for (iVar = 0; iVar < nVar; iVar ++ ) for (jVar = 0; jVar < nVar; jVar ++ ) Preconditioner[iVar][jVar] = delta*Preconditioner[iVar][jVar]; } else { - /*--- For explicit calculations, we move the residual to the right-hand side and pre-multiply by the preconditioner inverse. Therefore, we build inv(Precon) here and multiply by the residual later in the R-K and Euler Explicit time integration schemes. ---*/ - Preconditioner[0][0] = Temperature*BetaInc2*dRhodT/Density + BetaInc2; - for (iDim = 0; iDim < nDim; iDim ++) - Preconditioner[iDim+1][0] = -1.0*Velocity[iDim]/Density; - - if (energy) Preconditioner[nDim+1][0] = -1.0*Temperature/Density; - else Preconditioner[nDim+1][0] = 0.0; + if (multicomponent && energy) { + Preconditioner[0][0] = Enthalpy * BetaInc2 * dRhodT / Density + BetaInc2; + } else { + Preconditioner[0][0] = Temperature * BetaInc2 * dRhodT / Density + BetaInc2; + } + for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][0] = -1.0 * Velocity[iDim] / Density; + if (energy) { + if (multicomponent) { + Preconditioner[nDim + 1][0] = -1.0 * Enthalpy / Density; + } else { + Preconditioner[nDim + 1][0] = -1.0 * Temperature / Density; + } + } else { + Preconditioner[nDim + 1][0] = 0.0; + } for (jDim = 0; jDim < nDim; jDim++) { - Preconditioner[0][jDim+1] = 0.0; + Preconditioner[0][jDim + 1] = 0.0; for (iDim = 0; iDim < nDim; iDim++) { - if (iDim == jDim) Preconditioner[iDim+1][jDim+1] = 1.0/Density; - else Preconditioner[iDim+1][jDim+1] = 0.0; + if (iDim == jDim) + Preconditioner[iDim + 1][jDim + 1] = 1.0 / Density; + else + Preconditioner[iDim + 1][jDim + 1] = 0.0; } - Preconditioner[nDim+1][jDim+1] = 0.0; + Preconditioner[nDim + 1][jDim + 1] = 0.0; } - Preconditioner[0][nDim+1] = -1.0*BetaInc2*dRhodT*oneOverCp/Density; - for (iDim = 0; iDim < nDim; iDim ++) - Preconditioner[iDim+1][nDim+1] = 0.0; + if (multicomponent && energy) { + Preconditioner[0][nDim + 1] = -1.0 * BetaInc2 * dRhodT / Density; + } else { + Preconditioner[0][nDim + 1] = -1.0 * BetaInc2 * dRhodT * oneOverCp / Density; + } + for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][nDim + 1] = 0.0; - if (energy) Preconditioner[nDim+1][nDim+1] = oneOverCp/Density; - else Preconditioner[nDim+1][nDim+1] = 0.0; + if (energy) { + if (multicomponent) { + Preconditioner[nDim + 1][nDim + 1] = 1.0 / Density; + } else { + Preconditioner[nDim + 1][nDim + 1] = oneOverCp / Density; + } + } else { + Preconditioner[nDim + 1][nDim + 1] = 0.0; + } } - } void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -2217,6 +2337,7 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool viscous = config->GetViscous(); + const bool energy_multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) && (config->GetEnergy_Equation())); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); @@ -2380,6 +2501,15 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[prim_idx.Pressure()] = nodes->GetPressure(iPoint); } + /*-- Enthalpy is needed for energy equation in multicomponent and reacting flows. ---*/ + if (energy_multicomponent) { + CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel(); + const su2double* scalar_inlet = config->GetInlet_SpeciesVal(config->GetMarker_All_TagBound(val_marker)); + auxFluidModel->SetTDState_T(V_inlet[prim_idx.Temperature()], + scalar_inlet); // compute total enthalpy from temperature + V_inlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy(); + } + /*--- Access density at the node. This is either constant by construction, or will be set fixed implicitly by the temperature and equation of state. ---*/ @@ -2417,7 +2547,7 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Viscous contribution, commented out because serious convergence problems ---*/ - if (!viscous) continue; + if (!viscous || energy_multicomponent) continue; /*--- Set transport properties at the inlet ---*/ @@ -2469,6 +2599,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool viscous = config->GetViscous(); + const bool energy_multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) && (config->GetEnergy_Equation())); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); su2double Normal[MAXNDIM] = {0.0}; @@ -2591,6 +2722,15 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, V_outlet[prim_idx.CpTotal()] = nodes->GetSpecificHeatCp(iPoint); + /*-- Enthalpy is needed for energy equation in multicomponent and reacting flows. ---*/ + if (energy_multicomponent) { + CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel();; + const su2double* scalar_outlet = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint); + auxFluidModel->SetTDState_T(nodes->GetTemperature(iPoint), + scalar_outlet); // compute total enthalpy from temperature + V_outlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy(); + } + /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_outlet); @@ -2615,7 +2755,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Viscous contribution, commented out because serious convergence problems ---*/ - if (!viscous) continue; + if (!viscous || energy_multicomponent) continue; /*--- Set transport properties at the outlet. ---*/ diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index e0a4fda96ef..2128f85a450 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -69,10 +69,20 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE) && (InnerIter <= config->GetLimiterIter()); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const bool wall_functions = config->GetWall_Functions(); + const bool energy_multicomponent = (config->GetEnergy_Equation()) && (config->GetKind_FluidModel() == FLUID_MIXTURE); /*--- Common preprocessing steps (implemented by CEulerSolver) ---*/ CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); + if (energy_multicomponent && muscl) { + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto i_point = 0u; i_point < nPoint; i_point++) { + solver_container[FLOW_SOL]->GetNodes()->SetAuxVar(i_point, 1, + solver_container[FLOW_SOL]->GetNodes()->GetEnthalpy(i_point)); + } + END_SU2_OMP_FOR + } /*--- Compute gradient for MUSCL reconstruction ---*/ @@ -86,6 +96,21 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container default: break; } } + if (muscl && !center && energy_multicomponent) { + /*--- Gradient computation for MUSCL reconstruction of Enthalpy for multicomponent flows. ---*/ + + switch (config->GetKind_Gradient_Method_Recon()) { + case GREEN_GAUSS: + SetAuxVar_Gradient_GG(geometry, config); + break; + case LEAST_SQUARES: + case WEIGHTED_LEAST_SQUARES: + SetAuxVar_Gradient_LS(geometry, config); + break; + default: + break; + } + } /*--- Compute gradient of the primitive variables ---*/ @@ -275,7 +300,78 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, void CIncNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config) { + const bool energy_multicomponent = ((config->GetKind_FluidModel() == FLUID_MIXTURE) && config->GetEnergy_Equation()); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + /*--- Contribution to heat flux due to enthalpy diffusion for multicomponent and reacting flows ---*/ + if (energy_multicomponent) { + CVariable* speciesNodes = solver_container[SPECIES_SOL]->GetNodes(); + /*--- Points in edge ---*/ + + auto iPoint = geometry->edges->GetNode(iEdge, 0); + auto jPoint = geometry->edges->GetNode(iEdge, 1); + + /*--- Points coordinates, and normal vector ---*/ + + const su2double* Normal = geometry->edges->GetNormal(iEdge); + const su2double* Coord_i = geometry->nodes->GetCoord(iPoint); + const su2double* Coord_j = geometry->nodes->GetCoord(jPoint); + + /*--- Obtain fluid model for computing the enthalpy diffusion terms. ---*/ + + CFluidModel* FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); + /*--- retrieve number of species that are solved and set maximum static array ---*/ + + int n_species = config->GetnSpecies(); + static constexpr size_t MAXNVAR_SPECIES = 20UL; + + /*--- Species variables, and its gradients ---*/ + const su2double* Species_i = speciesNodes->GetSolution(iPoint); + const su2double* Species_j = speciesNodes->GetSolution(jPoint); + CMatrixView Species_Grad_i = speciesNodes->GetGradient(iPoint); + CMatrixView Species_Grad_j = speciesNodes->GetGradient(jPoint); + + /*--- Compute Projected gradient for species variables ---*/ + su2double ProjGradScalarVarNoCorr[MAXNVAR_SPECIES]{0.0}; + su2double Proj_Mean_GradScalarVar[MAXNVAR_SPECIES]{0.0}; + su2double proj_vector_ij = numerics->ComputeProjectedGradient( + nDim, n_species, Normal, Coord_i, Coord_j, Species_Grad_i, Species_Grad_j, true, Species_i, Species_j, + ProjGradScalarVarNoCorr, Proj_Mean_GradScalarVar); + (void)proj_vector_ij; + + /*--- Get enthalpy diffusion terms and its gradient(for implicit) for each species at point i. ---*/ + + su2double EnthalpyDiffusion_i[MAXNVAR_SPECIES]{0.0}; + su2double GradEnthalpyDiffusion_i[MAXNVAR_SPECIES]{0.0}; + FluidModel->SetTDState_T(nodes->GetPrimitive(iPoint)[prim_idx.Temperature()], Species_i); + FluidModel->GetEnthalpyDiffusivity(EnthalpyDiffusion_i); + if (implicit) FluidModel->GetGradEnthalpyDiffusivity(GradEnthalpyDiffusion_i); + + /*--- Repeat the above computations for jPoint. ---*/ + + su2double EnthalpyDiffusion_j[MAXNVAR_SPECIES]{0.0}; + su2double GradEnthalpyDiffusion_j[MAXNVAR_SPECIES]{0.0}; + FluidModel->SetTDState_T(nodes->GetPrimitive(jPoint)[prim_idx.Temperature()], Species_j); + FluidModel->GetEnthalpyDiffusivity(EnthalpyDiffusion_j); + if (implicit) FluidModel->GetGradEnthalpyDiffusivity(GradEnthalpyDiffusion_j); + + /*--- Compute Enthalpy diffusion flux and its jacobian (for implicit iterations) ---*/ + su2double flux_enthalpy_diffusion = 0.0; + su2double jac_flux_enthalpy_diffusion = 0.0; + for (int i_species = 0; i_species < n_species; i_species++) { + flux_enthalpy_diffusion += + 0.5 * (EnthalpyDiffusion_i[i_species] + EnthalpyDiffusion_j[i_species]) * Proj_Mean_GradScalarVar[i_species]; + if (implicit) + jac_flux_enthalpy_diffusion += 0.5 * (GradEnthalpyDiffusion_i[i_species] + GradEnthalpyDiffusion_j[i_species]) * + Proj_Mean_GradScalarVar[i_species]; + } + + /*--- Set heat flux and jacobian (for implicit) due to enthalpy diffusion ---*/ + + numerics->SetHeatFluxDiffusion(flux_enthalpy_diffusion); + if (implicit) numerics->SetJacHeatFluxDiffusion(jac_flux_enthalpy_diffusion); + } Viscous_Residual_impl(iEdge, geometry, solver_container, numerics, config); } @@ -337,6 +433,7 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool energy = config->GetEnergy_Equation(); const bool py_custom = config->GetMarker_All_PyCustom(val_marker); + const bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); /*--- Variables for streamwise periodicity ---*/ const bool streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); @@ -482,7 +579,6 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con const su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij; /*--- Get thermal conductivity ---*/ - const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint); /*--- Apply a weak boundary condition for the energy equation. @@ -493,7 +589,12 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con /*--- Jacobian contribution for temperature equation. ---*/ if (implicit) { - Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity * Area / dist_ij); + if (multicomponent) { + const su2double Cp = nodes->GetSpecificHeatCp(iPoint); + Jacobian.AddVal2Diag(iPoint, nDim + 1, thermal_conductivity * Area / (dist_ij * Cp)); + } else { + Jacobian.AddVal2Diag(iPoint, nDim + 1, thermal_conductivity * Area / dist_ij); + } } break; } // switch @@ -524,6 +625,7 @@ void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **sol const su2double Temperature_Ref = config->GetTemperature_Ref(); const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool energy = config->GetEnergy_Equation(); + const bool multicomponent = (config->GetKind_FluidModel() == FLUID_MIXTURE); /*--- Identify the boundary ---*/ @@ -609,7 +711,18 @@ void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **sol /*--- Strong imposition of the temperature on the fluid zone. ---*/ LinSysRes(iPoint, nDim+1) = 0.0; - nodes->SetSolution_Old(iPoint, nDim+1, Twall); + if (multicomponent) { + /*--- Retrieve scalars at wall node. ---*/ + su2double* scalars = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint); + /*--- Retrieve fluid model. ---*/ + CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); + /*--- Set thermodynamic state given wall temperature and species composition. ---*/ + fluid_model_local->SetTDState_T(Twall, scalars); + /*--- Set enthalpy obtained from fluid model. ---*/ + nodes->SetSolution_Old(iPoint, nDim + 1, fluid_model_local->GetEnthalpy()); + } else { + nodes->SetSolution_Old(iPoint, nDim + 1, Twall); + } nodes->SetEnergy_ResTruncError_Zero(iPoint); } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 089845f6eb3..b2c6d848971 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -30,7 +30,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double temperature, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) - : CFlowVariable(npoint, ndim, nvar, ndim + 9, + : CFlowVariable(npoint, ndim, nvar, ndim + 10, ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config), indices(ndim, 0) { diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index cbcfa3b0aad..b8b7f7deb81 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -46,19 +46,32 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su AuxVar.resize(nPoint,nAuxVar) = su2double(0.0); Grad_AuxVar.resize(nPoint,nAuxVar,nDim); } + if ((config->GetKind_FluidModel() == FLUID_MIXTURE) && config->GetEnergy_Equation()) { + Energy_Multicomponent = true; + if (config->GetMUSCL_Flow()) { + nAuxVar = 2; + AuxVar.resize(nPoint, nAuxVar) = su2double(0.0); + Grad_AuxVar.resize(nPoint, nAuxVar, nDim); + } + } } bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel, const su2double *scalar) { bool physical = true; + su2double Temperature; /*--- Set the value of the pressure ---*/ SetPressure(iPoint); - /*--- Set the value of the temperature directly ---*/ - - su2double Temperature = Solution(iPoint, indices.Temperature()); + if(Energy_Multicomponent){ + su2double Enthalpy = Solution(iPoint, nDim +1); + FluidModel->ComputeTempFromEnthalpy(Enthalpy, &Temperature, scalar); + } else { + /*--- Set the value of the temperature directly ---*/ + Temperature = Solution(iPoint, indices.Temperature()); + } auto check_temp = SetTemperature(iPoint, Temperature); /*--- Use the fluid model to compute the new value of density. @@ -69,11 +82,13 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do FluidModel->SetTDState_T(Temperature, scalar); - /*--- for FLAMELET: copy the LUT temperature into the solution ---*/ - Solution(iPoint,nDim+1) = FluidModel->GetTemperature(); - /*--- for FLAMELET: update the local temperature using LUT variables ---*/ - Temperature = Solution(iPoint,indices.Temperature()); - check_temp = SetTemperature(iPoint, Temperature); + if (!Energy_Multicomponent) { + /*--- for FLAMELET: copy the LUT temperature into the solution ---*/ + Solution(iPoint, nDim + 1) = FluidModel->GetTemperature(); + /*--- for FLAMELET: update the local temperature using LUT variables ---*/ + Temperature = Solution(iPoint, indices.Temperature()); + check_temp = SetTemperature(iPoint, Temperature); + } /*--- Set the value of the density ---*/ @@ -90,7 +105,12 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do /*--- Recompute the primitive variables ---*/ - Temperature = Solution(iPoint, indices.Temperature()); + if (Energy_Multicomponent) { + su2double Enthalpy = Solution(iPoint, nDim + 1); + FluidModel->ComputeTempFromEnthalpy(Enthalpy, &Temperature, scalar); + } else { + Temperature = Solution(iPoint, indices.Temperature()); + } SetTemperature(iPoint, Temperature); FluidModel->SetTDState_T(Temperature, scalar); SetDensity(iPoint, FluidModel->GetDensity()); @@ -122,6 +142,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do SetSpecificHeatCp(iPoint, FluidModel->GetCp()); SetSpecificHeatCv(iPoint, FluidModel->GetCv()); + if (Energy_Multicomponent) SetEnthalpy(iPoint, FluidModel->GetEnthalpy()); return physical; diff --git a/SU2_PY/SU2/io/historyMap.py b/SU2_PY/SU2/io/historyMap.py index a477478359c..be17da14643 100644 --- a/SU2_PY/SU2/io/historyMap.py +++ b/SU2_PY/SU2/io/historyMap.py @@ -244,6 +244,12 @@ "HEADER": "bgs[T]", "TYPE": "RESIDUAL", }, + "BGS_ENTHALPY": { + "DESCRIPTION": "Block-Gauss-Seidel residual of the " "enthalpy", + "GROUP": "BGS_RES", + "HEADER": "bgs[h]", + "TYPE": "RESIDUAL", + }, "BGS_TKE": { "DESCRIPTION": "BGS residual of kinetic energy (SST model).", "GROUP": "BGS_RES", @@ -894,6 +900,12 @@ "HEADER": "max[T]", "TYPE": "RESIDUAL", }, + "MAX_ENTHALPY": { + "DESCRIPTION": "Maximum residual of the enthalpy", + "GROUP": "MAX_RES", + "HEADER": "max[h]", + "TYPE": "RESIDUAL", + }, "MAX_TKE": { "DESCRIPTION": "Maximum residual of kinetic energy (SST model).", "GROUP": "MAX_RES", @@ -1154,6 +1166,12 @@ "HEADER": "rms[T]", "TYPE": "RESIDUAL", }, + "RMS_ENTHALPY": { + "DESCRIPTION": "Root mean square residual of the " "enthalpy", + "GROUP": "RMS_RES", + "HEADER": "rms[h]", + "TYPE": "RESIDUAL", + }, "RMS_TKE": { "DESCRIPTION": "Root-mean square residual of kinetic energy (SST " "model).", "GROUP": "RMS_RES",