Skip to content

Commit 92a016b

Browse files
AngPassAngPass
authored
Merge branch 'develop' into develop
2 parents 14323bb + 2c0049f commit 92a016b

File tree

58 files changed

+1200
-291
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

58 files changed

+1200
-291
lines changed

‎Common/include/CConfig.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -919,6 +919,7 @@ class CConfig {
919919
Initial_BCThrust, /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
920920
Pressure_FreeStream, /*!< \brief Total pressure of the fluid. */
921921
Pressure_Thermodynamic, /*!< \brief Thermodynamic pressure of the fluid. */
922+
Standard_Ref_Temperature, /*!< \brief Standard reference temperature for multicomponent flows. */
922923
Temperature_FreeStream, /*!< \brief Total temperature of the fluid. */
923924
Temperature_ve_FreeStream; /*!< \brief Total vibrational-electronic temperature of the fluid. */
924925
unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
@@ -1967,6 +1968,12 @@ class CConfig {
19671968
*/
19681969
su2double GetPressure_Thermodynamic(void) const { return Pressure_Thermodynamic; }
19691970

1971+
/*!
1972+
* \brief Get the value of the standard reference temperature for multicomponent flows.
1973+
* \return Standard reference temperature, Non-dimensionalized if it is needed for Non-Dimensional problems.
1974+
*/
1975+
su2double GetStandard_RefTemperatureND(void) const { return Standard_Ref_Temperature / Temperature_Ref; }
1976+
19701977
/*!
19711978
* \brief Get the value of the non-dimensionalized thermodynamic pressure.
19721979
* \return Non-dimensionalized thermodynamic pressure.

‎Common/src/CConfig.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1229,6 +1229,9 @@ void CConfig::SetConfig_Options() {
12291229
addDoubleOption("GAMMA_VALUE", Gamma, 1.4);
12301230
/*!\brief THERMODYNAMIC_PRESSURE \n DESCRIPTION: Thermodynamics(operating) Pressure (101325 Pa), only for incompressible flows) \ingroup Config*/
12311231
addDoubleOption("THERMODYNAMIC_PRESSURE", Pressure_Thermodynamic, 101325.0);
1232+
/*!\brief STANDARD_REFERENCE_TEMPERATURE \n DESCRIPTION: Standard reference temperature (298.15K), only for
1233+
* multicomponent incompressible flows) \ingroup Config*/
1234+
addDoubleOption("STANDARD_REFERENCE_TEMPERATURE", Standard_Ref_Temperature, 298.15);
12321235
/*!\brief CP_VALUE \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/
12331236
addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp);
12341237
/*!\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*/

‎SU2_CFD/include/fluid/CConstantDensity.hpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,5 +56,16 @@ class CConstantDensity final : public CFluidModel {
5656
decoupled equation. Hence, we update the value.
5757
Note Cp = Cv, (gamma = 1).*/
5858
Temperature = t;
59+
Enthalpy = Cp * Temperature;
60+
}
61+
62+
/*!
63+
* \brief Set the Dimensionless State using Enthalpy.
64+
* \param[in] val_enthalpy - Enthalpy value at the point.
65+
* \param[in] val_scalars - not used here.
66+
*/
67+
void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
68+
Enthalpy = val_enthalpy;
69+
Temperature = Enthalpy / Cp;
5970
}
6071
};

‎SU2_CFD/include/fluid/CFluidFlamelet.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,13 @@ class CFluidFlamelet final : public CFluidModel {
118118
*/
119119
void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override;
120120

121+
/*!
122+
* \brief Set the thermodynamic state.
123+
* \param[in] val_enthalpy - enthalpy
124+
* \param[in] val_scalars - pointer to species mass fractions
125+
*/
126+
void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
127+
121128
/*!
122129
* \brief Evaluate data-set for flamelet simulations.
123130
* \param[in] input_scalar - controlling variables used to interpolate manifold.

‎SU2_CFD/include/fluid/CFluidModel.hpp‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ class CLookUpTable;
4747
class CFluidModel {
4848
protected:
4949
su2double StaticEnergy{0.0}; /*!< \brief Internal Energy. */
50+
su2double Enthalpy{0.0}; /*!<*\brief Enthalpy. */
5051
su2double Entropy{0.0}; /*!< \brief Entropy. */
5152
su2double Density{0.0}; /*!< \brief Density. */
5253
su2double Pressure{0.0}; /*!< \brief Pressure. */
@@ -113,6 +114,11 @@ class CFluidModel {
113114
*/
114115
su2double GetStaticEnergy() const { return StaticEnergy; }
115116

117+
/*!
118+
* \brief Get fluid enthalpy.
119+
*/
120+
su2double GetEnthalpy() const { return Enthalpy; }
121+
116122
/*!
117123
* \brief Get fluid density.
118124
*/
@@ -186,6 +192,30 @@ class CFluidModel {
186192
return mass_diffusivity;
187193
}
188194

195+
/*!
196+
* \brief Get the enthalpy diffusivity terms for all species being solved.
197+
*
198+
* This function computes and retrieves the enthalpy diffusion terms required in the energy equation
199+
* for multicomponent flows.
200+
*
201+
* \param[in,out] enthalpy_diffusions - Array containing the enthalpy diffusion terms for all
202+
* species to be solved. The size of \p enthalpy_diffusions must be at least (n_species_mixture - 1),
203+
* corresponding to the number of species transport equations in the system.
204+
*/
205+
virtual void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions = nullptr) const {}
206+
207+
/*!
208+
* \brief Get the gradient of enthalpy diffusivity terms for all species being solved.
209+
*
210+
* This function computes and retrieves the gradient of the enthalpy diffusion terms with respect to temperature.
211+
* These terms are required for implicit computations when solving the energy equation for multicomponent flows.
212+
*
213+
* \param[in,out] grad_enthalpy_diffusions - Array containing the gradient of enthalpy diffusion terms for all
214+
* species to be solved. The size of \p grad_enthalpy_diffusions must be at least (n_species_mixture - 1),
215+
* corresponding to the number of species transport equations in the system.
216+
*/
217+
virtual void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions = nullptr) const {}
218+
189219
/*!
190220
* \brief Get fluid pressure partial derivative.
191221
*/
@@ -339,6 +369,13 @@ class CFluidModel {
339369
*/
340370
virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {}
341371

372+
/*!
373+
* \brief Virtual member.
374+
* \param[in] val_enthalpy - Enthalpy value at the point.
375+
* \param[in] val_scalars - Scalar mass fractions.
376+
*/
377+
virtual void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) {}
378+
342379
/*!
343380
* \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
344381
*/

‎SU2_CFD/include/fluid/CFluidScalar.hpp‎

Lines changed: 41 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -40,10 +40,12 @@
4040
class CFluidScalar final : public CFluidModel {
4141
private:
4242
const int n_species_mixture; /*!< \brief Number of species in mixture. */
43-
su2double Gas_Constant; /*!< \brief Specific gas constant. */
43+
su2double Gas_Constant; /*!< \brief Specific gas constant. */
4444
const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
45+
const su2double Ref_Temperature; /*!< \brief Standard Reference temperature, usually set to 298.15 K. */
4546
const su2double GasConstant_Ref; /*!< \brief Gas constant reference needed for Nondimensional problems. */
46-
const su2double Prandtl_Number; /*!< \brief Prandlt number.*/
47+
const su2double Prandtl_Turb_Number; /*!< \brief Prandlt turbulent number.*/
48+
const su2double Schmidt_Turb_Number; /*!< \brief Schmidt turbulent number.*/
4749

4850
const bool wilke;
4951
const bool davidson;
@@ -91,6 +93,11 @@ class CFluidScalar final : public CFluidModel {
9193
*/
9294
su2double ComputeMeanSpecificHeatCp(const su2double* val_scalars);
9395

96+
/*!
97+
* \brief Compute Enthalpy given the temperature and scalars.
98+
*/
99+
su2double ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars);
100+
94101
/*!
95102
* \brief Compute gas constant for mixture.
96103
*/
@@ -130,16 +137,47 @@ class CFluidScalar final : public CFluidModel {
130137
/*!
131138
* \brief Get fluid thermal conductivity.
132139
*/
133-
inline su2double GetThermalConductivity() override { return Kt + Mu_Turb * Cp / Prandtl_Number; }
140+
inline su2double GetThermalConductivity() override { return Kt + Mu_Turb * Cp / Prandtl_Turb_Number; }
134141

135142
/*!
136143
* \brief Get fluid mass diffusivity.
137144
*/
138145
inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[ivar]; }
139146

147+
/*!
148+
* \brief Get the enthalpy diffusivity terms for all species being solved.
149+
*
150+
* This function computes and retrieves the enthalpy diffusion terms required in the energy equation
151+
* for multicomponent flows.
152+
*
153+
* \param[in,out] enthalpy_diffusions - Array containing the enthalpy diffusion terms for all
154+
* species to be solved. The size of \p enthalpy_diffusions must be at least (n_species_mixture - 1),
155+
* corresponding to the number of species transport equations in the system.
156+
*/
157+
void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const override;
158+
159+
/*!
160+
* \brief Get the gradient of enthalpy diffusivity terms for all species being solved.
161+
*
162+
* This function computes and retrieves the gradient of the enthalpy diffusion terms with respect to temperature.
163+
* These terms are required for implicit computations when solving the energy equation for multicomponent flows.
164+
*
165+
* \param[in,out] grad_enthalpy_diffusions - Array containing the gradient of enthalpy diffusion terms for all
166+
* species to be solved. The size of \p grad_enthalpy_diffusions must be at least (n_species_mixture - 1),
167+
* corresponding to the number of species transport equations in the system.
168+
*/
169+
void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) const override;
170+
140171
/*!
141172
* \brief Set the Dimensionless State using Temperature.
142173
* \param[in] t - Temperature value at the point.
143174
*/
144175
void SetTDState_T(su2double val_temperature, const su2double* val_scalars) override;
176+
177+
/*!
178+
* \brief Virtual member.
179+
* \param[in] val_enthalpy - Enthalpy value at the point.
180+
* \param[in] val_scalars - Scalar mass fractions.
181+
*/
182+
void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
145183
};

‎SU2_CFD/include/fluid/CIncIdealGas.hpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,17 @@ class CIncIdealGas final : public CFluidModel {
5858
/*--- The EoS only depends upon temperature. ---*/
5959
Temperature = t;
6060
Density = Pressure / (Temperature * Gas_Constant);
61+
Enthalpy = Cp * Temperature;
62+
}
63+
64+
/*!
65+
* \brief Set the Dimensionless State using Enthalpy.
66+
* \param[in] val_enthalpy - Enthalpy value at the point.
67+
*/
68+
void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
69+
Enthalpy = val_enthalpy;
70+
Temperature = Enthalpy / Cp;
71+
Density = Pressure / (Temperature * Gas_Constant);
6172
}
6273

6374
private:

‎SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp‎

Lines changed: 60 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,9 +58,15 @@ class CIncIdealGasPolynomial final : public CFluidModel {
5858
* \param[in] config - configuration container for the problem.
5959
*/
6060
void SetCpModel(const CConfig* config) override {
61+
const su2double t_ref = config->GetStandard_RefTemperatureND();
62+
Enthalpy_Ref = 0.0;
63+
su2double t_i = 1.0;
6164
for (int i = 0; i < N; ++i) {
65+
t_i *= t_ref;
6266
coeffs_[i] = config->GetCp_PolyCoeffND(i);
67+
Enthalpy_Ref += coeffs_[i] * t_i / (i + 1);
6368
}
69+
Temperature_Min = config->GetTemperatureLimits(0);
6470
}
6571

6672
/*!
@@ -72,18 +78,70 @@ class CIncIdealGasPolynomial final : public CFluidModel {
7278
Temperature = t;
7379
Density = Pressure / (Temperature * Gas_Constant);
7480

75-
/* Evaluate the new Cp from the coefficients and temperature. */
81+
/* Evaluate the new Cp and enthalpy from the coefficients and temperature. */
7682
Cp = coeffs_[0];
83+
Enthalpy = coeffs_[0] * t - Enthalpy_Ref;
7784
su2double t_i = 1.0;
7885
for (int i = 1; i < N; ++i) {
7986
t_i *= t;
8087
Cp += coeffs_[i] * t_i;
88+
Enthalpy += coeffs_[i] * t_i * t / (i + 1);
89+
}
90+
Cv = Cp / Gamma;
91+
}
92+
93+
/*!
94+
* \brief Set the Dimensionless State using enthalpy.
95+
* \param[in] val_enthalpy - Enthalpy value at the point.
96+
*/
97+
void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
98+
Enthalpy = val_enthalpy;
99+
/*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/
100+
const su2double toll = 1e-5;
101+
su2double temp_iter = 300.0;
102+
su2double Cp_iter = 0.0;
103+
su2double delta_temp_iter = 1e10;
104+
su2double delta_enthalpy_iter;
105+
const int counter_limit = 20;
106+
107+
int counter = 0;
108+
109+
/*--- Computing temperature given enthalpy using Newton-Raphson. ---*/
110+
while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) {
111+
/* Evaluate the new Cp and enthalpy from the coefficients and temperature. */
112+
Cp_iter = coeffs_[0];
113+
su2double Enthalpy_iter = coeffs_[0] * temp_iter - Enthalpy_Ref;
114+
su2double t_i = 1.0;
115+
for (int i = 1; i < N; ++i) {
116+
t_i *= temp_iter;
117+
Cp_iter += coeffs_[i] * t_i;
118+
Enthalpy_iter += coeffs_[i] * t_i * temp_iter / (i + 1);
119+
}
120+
121+
delta_enthalpy_iter = Enthalpy - Enthalpy_iter;
122+
123+
delta_temp_iter = delta_enthalpy_iter / Cp_iter;
124+
125+
temp_iter += delta_temp_iter;
126+
if (temp_iter < Temperature_Min) {
127+
cout << "Warning: Negative temperature has been found during Newton-Raphson" << endl;
128+
temp_iter = Temperature_Min;
129+
break;
130+
}
131+
}
132+
Temperature = temp_iter;
133+
Cp = Cp_iter;
134+
if (counter == counter_limit) {
135+
cout << "Warning Newton-Raphson exceed number of max iteration in temperature computation" << endl;
81136
}
137+
Density = Pressure / (Temperature * Gas_Constant);
82138
Cv = Cp / Gamma;
83139
}
84140

85141
private:
86142
su2double Gas_Constant{0.0}; /*!< \brief Specific Gas Constant. */
87143
su2double Gamma{0.0}; /*!< \brief Ratio of specific heats. */
88-
array<su2double, N> coeffs_; /*!< \brief Polynomial coefficients for conductivity as a function of temperature. */
144+
array<su2double, N> coeffs_; /*!< \brief Polynomial coefficients for heat capacity as a function of temperature. */
145+
su2double Enthalpy_Ref; /*!< \brief Enthalpy computed at the reference temperature. */
146+
su2double Temperature_Min; /*!< \brief Minimum temperature value allowed in Newton-Raphson iterations. */
89147
};

‎SU2_CFD/include/numerics/CNumerics.hpp‎

Lines changed: 26 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,10 @@ class CNumerics {
130130
const su2double
131131
*ScalarVar_i, /*!< \brief Vector of scalar variables at point i. */
132132
*ScalarVar_j; /*!< \brief Vector of scalar variables at point j. */
133+
su2double
134+
HeatFluxDiffusion; /*!< \brief Heat flux due to enthalpy diffusion for multicomponent. */
135+
su2double
136+
JacHeatFluxDiffusion; /*!< \brief Heat flux jacobian due to enthalpy diffusion for multicomponent. */
133137
const su2double
134138
*TransVar_i, /*!< \brief Vector of turbulent variables at point i. */
135139
*TransVar_j; /*!< \brief Vector of turbulent variables at point j. */
@@ -196,6 +200,8 @@ class CNumerics {
196200

197201
bool nemo; /*!< \brief Flag for NEMO problems */
198202

203+
bool energy_multicomponent = false; /*!< \brief Flag for multicomponent and reacting flow */
204+
199205
bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */
200206

201207
public:
@@ -799,6 +805,20 @@ class CNumerics {
799805
Diffusion_Coeff_j = val_diffusioncoeff_j;
800806
}
801807

808+
/*!
809+
* \brief Set the heat flux due to enthalpy diffusion
810+
* \param[in] val_heatfluxdiffusion - Value of the heat flux due to enthalpy diffusion.
811+
*/
812+
inline void SetHeatFluxDiffusion(su2double val_heatfluxdiffusion) { HeatFluxDiffusion = val_heatfluxdiffusion; }
813+
814+
/*!
815+
* \brief Set Jacobian of the heat flux due to enthalpy diffusion
816+
* \param[in] val_jacheatfluxdiffusion - Value of the heat flux jacobian due to enthalpy diffusion.
817+
*/
818+
inline void SetJacHeatFluxDiffusion(su2double val_jacheatfluxdiffusion) {
819+
JacHeatFluxDiffusion = val_jacheatfluxdiffusion;
820+
}
821+
802822
/*!
803823
* \brief Set the laminar viscosity.
804824
* \param[in] val_laminar_viscosity_i - Value of the laminar viscosity at point i.
@@ -1117,8 +1137,7 @@ class CNumerics {
11171137
* \param[in] val_density - Value of the density.
11181138
* \param[in] val_velocity - Pointer to the velocity.
11191139
* \param[in] val_betainc2 - Value of the artificial compresibility factor.
1120-
* \param[in] val_cp - Value of the specific heat at constant pressure.
1121-
* \param[in] val_temperature - Value of the temperature.
1140+
* \param[in] val_enthalpy - Value of the enthalpy.
11221141
* \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
11231142
* \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
11241143
* \param[in] val_scale - Scale of the projection.
@@ -1127,8 +1146,7 @@ class CNumerics {
11271146
void GetInviscidIncProjJac(const su2double *val_density,
11281147
const su2double *val_velocity,
11291148
const su2double *val_betainc2,
1130-
const su2double *val_cp,
1131-
const su2double *val_temperature,
1149+
const su2double *val_enthalpy,
11321150
const su2double *val_dRhodT,
11331151
const su2double *val_normal,
11341152
su2double val_scale,
@@ -1139,17 +1157,15 @@ class CNumerics {
11391157
* \param[in] val_density - Value of the density.
11401158
* \param[in] val_velocity - Pointer to the velocity.
11411159
* \param[in] val_betainc2 - Value of the artificial compresibility factor.
1142-
* \param[in] val_cp - Value of the specific heat at constant pressure.
1143-
* \param[in] val_temperature - Value of the temperature.
1144-
* \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
1160+
* \param[in] val_enthalpy - Value of the enthalpy.
1161+
* \param[in] val_dRhodh - Value of the derivative of density w.r.t. enthalpy.
11451162
* \param[out] val_Precon - Pointer to the preconditioning matrix.
11461163
*/
11471164
void GetPreconditioner(const su2double *val_density,
11481165
const su2double *val_velocity,
11491166
const su2double *val_betainc2,
1150-
const su2double *val_cp,
1151-
const su2double *val_temperature,
1152-
const su2double *val_drhodt,
1167+
const su2double *val_enthalpy,
1168+
const su2double *val_drhodh,
11531169
su2double **val_Precon) const;
11541170

11551171
/*!

0 commit comments

Comments
 (0)