diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index d5d118aa1de4..30207e6301f7 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -326,6 +326,7 @@ class CConfig { su2double *Engine_Area; /*!< \brief Specified engine area for nacelle boundaries. */ su2double *Outlet_Pressure; /*!< \brief Specified back pressures (static) for outlet boundaries. */ su2double *Isothermal_Temperature; /*!< \brief Specified isothermal wall temperatures (static). */ + mutable std::map> Marker_Custom_Temperature; /*!< \brief Spatially varying isothermal wall temperatures. */ su2double *HeatTransfer_Coeff; /*!< \brief Specified heat transfer coefficients. */ su2double *HeatTransfer_WallTemp; /*!< \brief Specified temperatures at infinity alongside heat transfer coefficients. */ su2double *Heat_Flux; /*!< \brief Specified wall heat fluxes. */ @@ -911,6 +912,9 @@ class CConfig { 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. */ array KtPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */ + array NASA_CpLowPolyCoefficients{{0.0}}; /*!< \brief Definition of the NASA temperature polynomial coefficients for specific heat Cp (low temperature range). */ + array NASA_CpHighPolyCoefficients{{0.0}}; /*!< \brief Definition of the NASA temperature polynomial coefficients for specific heat Cp (high temperature range). */ + su2double NASA_TransitionTemperature; /*!< \brief Transition temperature for the NASA polynomial gas model. */ su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ su2double Energy_FreeStream, /*!< \brief Free-stream total energy of the fluid. */ ModVel_FreeStream, /*!< \brief Magnitude of the free-stream velocity of the fluid. */ @@ -1124,6 +1128,9 @@ class CConfig { array cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ + array nasa_cp_low_polycoeffs{{0.0}}; /*!< \brief Array for NASA high temperature polynomial coefficients. */ + array nasa_cp_high_polycoeffs{{0.0}}; /*!< \brief Array for NASA high temperature polynomial coefficients. */ + su2double nasa_transition_temperature; /*!< \brief Transition temperature for the NASA polynomial gas model. */ bool Body_Force; /*!< \brief Flag to know if a body force is included in the formulation. */ struct CMUSCLRampParam { @@ -7426,6 +7433,22 @@ class CConfig { */ su2double GetIsothermal_Temperature(const string& val_index) const; + /*! + * \brief Get the custom wall temperature (static) at an isothermal boundary node. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Vertex index on the marker. + * \return The custom wall temperature (returns -1.0 if not set). + */ + su2double GetCustom_Temperature(unsigned short iMarker, unsigned long iVertex) const; + + /*! + * \brief Set the custom wall temperature (static) at an isothermal boundary node. + * \param[in] iMarker - Marker index. + * \param[in] iVertex - Vertex index on the marker. + * \param[in] val - The custom wall temperature. + */ + void SetCustom_Temperature(unsigned short iMarker, unsigned long iVertex, su2double val); + /*! * \brief Get the wall heat flux on a constant heat flux boundary. * \param[in] val_index - Index corresponding to the constant heat flux boundary. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1361a12072da..acf81cb9f949 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -202,6 +202,7 @@ constexpr passivedouble COLORING_EFF_THRESH = 0.875; /*!< \brief Below this val Cp(T) = b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4. By default, all coeffs are set to zero and will be properly non-dim. in the solver. ---*/ constexpr int N_POLY_COEFFS = 5; /*!< \brief Number of coefficients in temperature polynomial fits for fluid models. */ +constexpr int N_NASA_POLY_COEFFS = 7; /*!< \brief Number of coefficients in NASA temperature polynomial fits. */ /*! * \brief Boolean answers @@ -551,6 +552,7 @@ enum ENUM_FLUIDMODEL { COOLPROP = 10, /*!< \brief Thermodynamics library. */ FLUID_FLAMELET = 11, /*!< \brief lookup table (LUT) method for premixed flamelets. */ DATADRIVEN_FLUID = 12, /*!< \brief multi-layer perceptron driven fluid model. */ + NASA_GAS = 13, /*!< \brief NASA polynomial ideal gas model. */ }; static const MapType FluidModel_Map = { MakePair("STANDARD_AIR", STANDARD_AIR) @@ -566,6 +568,7 @@ static const MapType FluidModel_Map = { MakePair("COOLPROP", COOLPROP) MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID) MakePair("FLUID_FLAMELET", FLUID_FLAMELET) + MakePair("NASA_GAS", NASA_GAS) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 79fe6d6aab1e..d60aa0813b50 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1334,6 +1334,12 @@ void CConfig::SetConfig_Options() { addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, false, mu_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data()); + /* DESCRIPTION: Definition of the NASA temperature polynomial coefficients for specific heat Cp (low temperature range). */ + addDoubleArrayOption("NASA_CP_LOW_COEFFS", N_NASA_POLY_COEFFS, false, nasa_cp_low_polycoeffs.data()); + /* DESCRIPTION: Definition of the NASA temperature polynomial coefficients for specific heat Cp (high temperature range). */ + addDoubleArrayOption("NASA_CP_HIGH_COEFFS", N_NASA_POLY_COEFFS, false, nasa_cp_high_polycoeffs.data()); + /* DESCRIPTION: Transition temperature for the NASA polynomial gas model. */ + addDoubleOption("NASA_T_TRANS", nasa_transition_temperature, 1000.0); /*!\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); @@ -9608,6 +9614,21 @@ su2double CConfig::GetIsothermal_Temperature(const string& val_marker) const { return Isothermal_Temperature[0]; } +su2double CConfig::GetCustom_Temperature(unsigned short iMarker, unsigned long iVertex) const { + auto it_marker = Marker_Custom_Temperature.find(iMarker); + if (it_marker != Marker_Custom_Temperature.end()) { + auto it_vertex = it_marker->second.find(iVertex); + if (it_vertex != it_marker->second.end()) { + return it_vertex->second; + } + } + return -1.0; +} + +void CConfig::SetCustom_Temperature(unsigned short iMarker, unsigned long iVertex, su2double val) { + Marker_Custom_Temperature[iMarker][iVertex] = val; +} + su2double CConfig::GetWall_HeatFlux(const string& val_marker) const { for (unsigned short iMarker_HeatFlux = 0; iMarker_HeatFlux < nMarker_HeatFlux; iMarker_HeatFlux++) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index b01f560fa281..252566b1274e 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -555,6 +555,30 @@ class CDriver : public CDriverBase { */ void SetMarkerTranslationRate(unsigned short iMarker, passivedouble vel_x, passivedouble vel_y, passivedouble vel_z); + /*! + * \brief Set a custom temperature at a boundary vertex. + * \param[in] iMarker - Index of the boundary marker. + * \param[in] iVertex - Index of the boundary vertex. + * \param[in] Temperature - Temperature value. + */ + void SetMarkerCustomTemperature(unsigned short iMarker, unsigned long iVertex, su2double Temperature); + + /*! + * \brief Get the local speed of sound at a boundary vertex. + * \param[in] iMarker - Index of the boundary marker. + * \param[in] iVertex - Index of the boundary vertex. + * \return Local speed of sound. + */ + passivedouble GetMarkerLocalSpeedOfSound(unsigned short iMarker, unsigned long iVertex); + + /*! + * \brief Get the Mach number at a boundary vertex. + * \param[in] iMarker - Index of the boundary marker. + * \param[in] iVertex - Index of the boundary vertex. + * \return Mach number. + */ + passivedouble GetMarkerMachNumber(unsigned short iMarker, unsigned long iVertex); + /*! * \brief Get the Freestream Density for nondimensionalization * \return Freestream Density diff --git a/SU2_CFD/include/fluid/CNasaGas.hpp b/SU2_CFD/include/fluid/CNasaGas.hpp new file mode 100644 index 000000000000..fd74bbdd2d88 --- /dev/null +++ b/SU2_CFD/include/fluid/CNasaGas.hpp @@ -0,0 +1,64 @@ +/*! + * \file CNasaGas.hpp + * \brief Defines the NASA polynomial gas model. + * \author Om Giri + * \version 8.4.0 "Harrier" + */ + +#pragma once + +#include "CIdealGas.hpp" +#include + +/*! + * \class CNasaGas + * \brief Child class for defining the NASA polynomial gas model. + */ +class CNasaGas : public CIdealGas { + protected: + array CpLowPolyCoefficientsND{{0.0}}; /*!< \brief NASA low temp polynomial coefficients. */ + array CpHighPolyCoefficientsND{{0.0}}; /*!< \brief NASA high temp polynomial coefficients. */ + su2double TransitionTemperatureND{0.0}; /*!< \brief NASA transition temperature. */ + + public: + /*! + * \brief Constructor of the class. + */ + CNasaGas(su2double gamma, su2double R, bool CompEntropy = true); + + /*! + * \brief Set specific heat Cp model. + */ + void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override; + + /*! + * \brief Set the Dimensionless State using Density and Internal Energy + */ + void SetTDState_rhoe(su2double rho, su2double e) override; + + /*! + * \brief Set the Dimensionless State using Pressure and Temperature + */ + void SetTDState_PT(su2double P, su2double T) override; + + /*! + * \brief Set the Dimensionless State using Pressure and Density + */ + void SetTDState_Prho(su2double P, su2double rho) override; + + /*! + * \brief Set the Dimensionless State using Density and Temperature + */ + void SetTDState_rhoT(su2double rho, su2double T) override; + + /*! + * \brief Set the Dimensionless State using Pressure and Entropy + */ + void SetTDState_Ps(su2double P, su2double s) override; + + protected: + /*! + * \brief Calculate Cp, Enthalpy and Entropy from Temperature. + */ + void ComputeProperties_T(su2double T); +}; diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp index 992a9adb491a..3a16e5844756 100644 --- a/SU2_CFD/src/fluid/CFluidModel.cpp +++ b/SU2_CFD/src/fluid/CFluidModel.cpp @@ -45,6 +45,7 @@ #include "../../include/fluid/CCoolPropViscosity.hpp" #include "../../include/fluid/CConstantLewisDiffusivity.hpp" #include "../../include/fluid/CCoolPropConductivity.hpp" +#include "../../include/fluid/CNasaGas.hpp" unique_ptr CFluidModel::MakeLaminarViscosityModel(const CConfig* config, unsigned short iSpecies) { switch (config->GetKind_ViscosityModel()) { diff --git a/SU2_CFD/src/fluid/CNasaGas.cpp b/SU2_CFD/src/fluid/CNasaGas.cpp new file mode 100644 index 000000000000..9e3da6ef6554 --- /dev/null +++ b/SU2_CFD/src/fluid/CNasaGas.cpp @@ -0,0 +1,146 @@ +/*! + * \file CNasaGas.cpp + * \brief Source of the NASA polynomial gas model. + * \author Om Giri + * \version 8.4.0 "Harrier" + */ + +#include "../../include/fluid/CNasaGas.hpp" + +CNasaGas::CNasaGas(su2double gamma, su2double R, bool CompEntropy) : CIdealGas(gamma, R, CompEntropy) { + // Constructor inherited from CIdealGas +} + +void CNasaGas::SetCpModel(const CConfig* config, su2double val_Temperature_Ref) { + for (int i = 0; i < 7; ++i) { + CpLowPolyCoefficientsND[i] = config->GetNASA_CpLowPolyCoeff(i); + CpHighPolyCoefficientsND[i] = config->GetNASA_CpHighPolyCoeff(i); + } + TransitionTemperatureND = config->GetNASA_TransitionTemperature() / val_Temperature_Ref; +} + +void CNasaGas::ComputeProperties_T(su2double T) { + Temperature = T; + const array* c; + + if (T < TransitionTemperatureND) { + c = &CpLowPolyCoefficientsND; + } else { + c = &CpHighPolyCoefficientsND; + } + + /*--- NASA polynomials: + Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4 + H/RT = -a1*T^-2 + a2*log(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T + S/R = -a1*T^-2/2 - a2*T^-1 + a3*log(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9 + + Note: SU2 usually uses a 7-coefficient form: + Cp/R = a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4 + H/RT = a0 + a1*T/2 + a2*T^2/3 + a3*T^3/4 + a4*T^4/5 + a5/T + S/R = a0*log(T) + a1*T + a2*T^2/2 + a3*T^3/3 + a4*T^4/4 + a6 + + Wait, standard NASA 7-coefficient format is: + Cp/R = a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4 + H/RT = a0 + a1*T/2 + a2*T^2/3 + a3*T^3/4 + a4*T^4/5 + a5/T + S/R = a0*lnT + a1*T + a2*T^2/2 + a3*T^3/3 + a4*T^4/4 + a6 + + I will use this standard 7-coefficient format. + ---*/ + + Cp = ((*c)[0] + (*c)[1]*T + (*c)[2]*T*T + (*c)[3]*T*T*T + (*c)[4]*T*T*T*T) * Gas_Constant; + Enthalpy = ((*c)[0] + (*c)[1]*T/2.0 + (*c)[2]*T*T/3.0 + (*c)[3]*T*T*T/4.0 + (*c)[4]*T*T*T*T/5.0 + (*c)[5]/T) * T * Gas_Constant; + + if (ComputeEntropy) { + Entropy = ((*c)[0]*log(T) + (*c)[1]*T + (*c)[2]*T*T/2.0 + (*c)[3]*T*T*T/3.0 + (*c)[4]*T*T*T*T/4.0 + (*c)[6]) * Gas_Constant; + // Note: This is partial entropy (temperature dependent part). Pressure part added in SetTDState_rhoe if needed. + // However, SU2 CIdealGas includes log(1/rho) term. + } + + // Cv = Cp - R + Cv = Cp - Gas_Constant; + // Gamma = Cp / Cv + Gamma = Cp / Cv; + Gamma_Minus_One = Gamma - 1.0; +} + +void CNasaGas::SetTDState_rhoe(su2double rho, su2double e) { + Density = rho; + StaticEnergy = e; + + /*--- Solve for T using Newton-Raphson: e(T) = h(T) - R*T ---*/ + su2double T_iter = Temperature; + if (T_iter <= 0.0) T_iter = 1.0; // Initial guess + + for (int i = 0; i < 10; ++i) { + ComputeProperties_T(T_iter); + su2double e_iter = Enthalpy - Gas_Constant * T_iter; + su2double de_dT = Cv; // de/dT = Cv for ideal gas even with variable Cp + su2double deltaT = (e - e_iter) / de_dT; + T_iter += deltaT; + if (abs(deltaT) < 1e-8 * T_iter) break; + } + + Temperature = T_iter; + ComputeProperties_T(Temperature); + + Pressure = Density * Gas_Constant * Temperature; + SoundSpeed2 = Gamma * Pressure / Density; + + /*--- Derivatives ---*/ + dPde_rho = Density * Gas_Constant / Cv; + dPdrho_e = Gas_Constant * Temperature - StaticEnergy * dPde_rho; // From P = rho*R*T and e = e(T) + + dTde_rho = 1.0 / Cv; + dTdrho_e = 0.0; + + if (ComputeEntropy) { + // Entropy was S(T). Now add pressure/density part: S = S(T) - R*ln(rho*R) + R*ln(R_ref?) + // CIdealGas uses: Entropy = (1.0 / Gamma_Minus_One * log(Temperature) + log(1.0 / Density)) * Gas_Constant; + // For NASA, S(T) already has the log(T) part. + Entropy += log(1.0 / Density) * Gas_Constant; + } +} + +void CNasaGas::SetTDState_PT(su2double P, su2double T) { + Pressure = P; + Temperature = T; + ComputeProperties_T(T); + Density = P / (Gas_Constant * T); + StaticEnergy = Enthalpy - Gas_Constant * T; + SetTDState_rhoe(Density, StaticEnergy); +} + +void CNasaGas::SetTDState_Prho(su2double P, su2double rho) { + Pressure = P; + Density = rho; + Temperature = P / (Density * Gas_Constant); + ComputeProperties_T(Temperature); + StaticEnergy = Enthalpy - Gas_Constant * Temperature; + SetTDState_rhoe(Density, StaticEnergy); +} + +void CNasaGas::SetTDState_rhoT(su2double rho, su2double T) { + Density = rho; + Temperature = T; + ComputeProperties_T(T); + Pressure = Density * Gas_Constant * T; + StaticEnergy = Enthalpy - Gas_Constant * T; + SetTDState_rhoe(Density, StaticEnergy); +} + +void CNasaGas::SetTDState_Ps(su2double P, su2double s) { + /*--- This would require another nested iteration to find T from P and s. + For now, we can use an approximate T or a simple Newton solver. ---*/ + su2double T_iter = Temperature; + if (T_iter <= 0.0) T_iter = 1.0; + + for (int i = 0; i < 10; ++i) { + ComputeProperties_T(T_iter); + su2double s_iter = Entropy + log(Gas_Constant * T_iter / P) * Gas_Constant; + su2double ds_dT = Cp / T_iter; + su2double deltaT = (s - s_iter) / ds_dT; + T_iter += deltaT; + if (abs(deltaT) < 1e-8 * T_iter) break; + } + SetTDState_PT(P, T_iter); +} diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3b40822a34f4..4083a71c8c0b 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -13,7 +13,8 @@ su2_cfd_src += files(['fluid/CFluidModel.cpp', 'fluid/CNEMOGas.cpp', 'fluid/CMutationTCLib.cpp', 'fluid/CSU2TCLib.cpp', - 'fluid/CDataDrivenFluid.cpp']) + 'fluid/CDataDrivenFluid.cpp', + 'fluid/CNasaGas.cpp']) su2_cfd_src += files(['output/COutputFactory.cpp', 'output/CAdjElasticityOutput.cpp', diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 54a30d1b3a2b..d9d1fdaa5804 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -166,6 +166,8 @@ void CFlowCompOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("MIN_DELTA_TIME", "Min DT", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current minimum local time step"); AddHistoryOutput("MAX_DELTA_TIME", "Max DT", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current maximum local time step"); + AddHistoryOutput("AVG_SOUND_SPEED", "Avg[a]", ScreenOutputFormat::FIXED, "PRIMITIVE", "Average local speed of sound.", HistoryFieldType::DEFAULT); + AddHistoryOutput("MIN_CFL", "Min CFL", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current minimum of the local CFL numbers"); AddHistoryOutput("MAX_CFL", "Max CFL", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current maximum of the local CFL numbers"); AddHistoryOutput("AVG_CFL", "Avg CFL", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current average of the local CFL numbers"); @@ -238,6 +240,7 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("PRESSURE", "Pressure", "PRIMITIVE", "Pressure"); AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature"); AddVolumeOutput("MACH", "Mach", "PRIMITIVE", "Mach number"); + AddVolumeOutput("LOCAL_SPEED_OF_SOUND", "Local_Speed_of_Sound", "PRIMITIVE", "Local speed of sound"); AddVolumeOutput("PRESSURE_COEFF", "Pressure_Coefficient", "PRIMITIVE", "Pressure coefficient"); AddVolumeOutput("VELOCITY-X", "Velocity_x", "PRIMITIVE", "x-component of the velocity vector"); AddVolumeOutput("VELOCITY-Y", "Velocity_y", "PRIMITIVE", "y-component of the velocity vector"); @@ -329,6 +332,8 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv SetVolumeOutputValue("ENERGY", iPoint, Node_Flow->GetSolution(iPoint, 3)); } + SetVolumeOutputValue("LOCAL_SPEED_OF_SOUND", iPoint, Node_Flow->GetSoundSpeed(iPoint)); + if (gridMovement){ SetVolumeOutputValue("GRID_VELOCITY-X", iPoint, Node_Geo->GetGridVel(iPoint)[0]); SetVolumeOutputValue("GRID_VELOCITY-Y", iPoint, Node_Geo->GetGridVel(iPoint)[1]); @@ -437,6 +442,26 @@ void CFlowCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSol SetHistoryOutputValue("MIN_DELTA_TIME", flow_solver->GetMin_Delta_Time()); SetHistoryOutputValue("MAX_DELTA_TIME", flow_solver->GetMax_Delta_Time()); + /*--- Compute average local speed of sound for screen output ---*/ + su2double avg_a = 0.0; + unsigned long nPoint = geometry->GetnPointDomain(); + const auto* Node_Flow = flow_solver->GetNodes(); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + avg_a += Node_Flow->GetSoundSpeed(iPoint); + } + su2double total_a = 0.0; + unsigned long total_points = 0; +#ifdef HAVE_MPI + unsigned long nPoint_local = nPoint; + SU2_MPI::Allreduce(&avg_a, &total_a, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&nPoint_local, &total_points, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); +#else + total_a = avg_a; + total_points = nPoint; +#endif + if (total_points > 0) total_a /= (su2double)total_points; + SetHistoryOutputValue("AVG_SOUND_SPEED", total_a); + SetHistoryOutputValue("MIN_CFL", flow_solver->GetMin_CFL_Local()); SetHistoryOutputValue("MAX_CFL", flow_solver->GetMax_CFL_Local()); SetHistoryOutputValue("AVG_CFL", flow_solver->GetAvg_CFL_Local()); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 446db73dc13f..a9f939530422 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -250,8 +250,19 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry Mach = sqrt(Velocity2)/SoundSpeed; Temperature = Pressure / (Gas_Constant * Density); Enthalpy = flow_nodes->GetEnthalpy(iPoint); - TotalTemperature = Temperature * (1.0 + Mach * Mach * 0.5 * (Gamma - 1.0)); - TotalPressure = Pressure * pow( 1.0 + Mach * Mach * 0.5 * (Gamma - 1.0), Gamma / (Gamma - 1.0)); + if (config->GetKind_FluidModel() == NASA_GAS) { + TotalTemperature = Temperature; + TotalPressure = Pressure; + // For variable Cp, we'd need an iterative or approximate approach for stagnation properties. + // For now, use a simplified approach or leave as static if complex. + // A common approximation for small Mach: Tt = T + V^2/(2*Cp) + TotalTemperature = Temperature + 0.5 * Velocity2 / flow_nodes->GetSpecificHeatCp(iPoint); + // Pt approximation: Pt = P * (Tt/T)^(Cp/R) + TotalPressure = Pressure * pow(TotalTemperature / Temperature, flow_nodes->GetSpecificHeatCp(iPoint) / Gas_Constant); + } else { + TotalTemperature = Temperature * (1.0 + Mach * Mach * 0.5 * (Gamma - 1.0)); + TotalPressure = Pressure * pow( 1.0 + Mach * Mach * 0.5 * (Gamma - 1.0), Gamma / (Gamma - 1.0)); + } } /*--- Compute the mass Surface_MassFlow ---*/ diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 8c5f939305da..3a05e028164b 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -28,6 +28,7 @@ #include "../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../include/drivers/CDriver.hpp" #include "../include/drivers/CSinglezoneDriver.hpp" +#include void CDriver::PreprocessPythonInterface(CConfig** config, CGeometry**** geometry, CSolver***** solver) { int rank = MASTER_NODE; @@ -179,3 +180,23 @@ void CDriver::SetMarkerTranslationRate(unsigned short iMarker, passivedouble vel config_container[selected_zone]->SetMarkerTranslationRate(iMarker, 2, vel_z); } +void CDriver::SetMarkerCustomTemperature(unsigned short iMarker, unsigned long iVertex, su2double Temperature) { + geometry_container[selected_zone][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, Temperature); +} + +passivedouble CDriver::GetMarkerLocalSpeedOfSound(unsigned short iMarker, unsigned long iVertex) { + const auto* geometry = geometry_container[selected_zone][INST_0][MESH_0]; + const auto* flow_solver = solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]; + unsigned long iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + return SU2_TYPE::GetValue(flow_solver->GetNodes()->GetSoundSpeed(iPoint)); +} + +passivedouble CDriver::GetMarkerMachNumber(unsigned short iMarker, unsigned long iVertex) { + const auto* geometry = geometry_container[selected_zone][INST_0][MESH_0]; + const auto* flow_solver = solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]; + unsigned long iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + su2double V2 = flow_solver->GetNodes()->GetVelocity2(iPoint); + su2double a = flow_solver->GetNodes()->GetSoundSpeed(iPoint); + return SU2_TYPE::GetValue(sqrt(V2) / a); +} + diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 3f8d642e9c6e..302e7f32db10 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -881,6 +881,13 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes auxFluidModel = new CDataDrivenFluid(config); + break; + + case NASA_GAS: + + auxFluidModel = new CNasaGas(Gamma, config->GetGas_Constant()); + auxFluidModel->SetCpModel(config, 1.0); + break; case COOLPROP: @@ -1133,6 +1140,11 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes FluidModel[thread] = new CDataDrivenFluid(config, false); break; + case NASA_GAS: + FluidModel[thread] = new CNasaGas(Gamma, Gas_ConstantND); + GetFluidModel()->SetCpModel(config, config->GetTemperature_Ref()); + break; + case COOLPROP: FluidModel[thread] = new CCoolProp(config->GetFluid_Name()); break; @@ -1311,6 +1323,9 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes case PR_GAS: ModelTable << "PR_GAS"; break; + case NASA_GAS: + ModelTable << "NASA_GAS"; + break; case COOLPROP: ModelTable << "CoolProp library"; break; diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 3f7694a0cbee..afab2eee2766 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -33,6 +33,7 @@ #include "../../include/fluid/CPengRobinson.hpp" #include "../../include/fluid/CCoolProp.hpp" #include "../../include/fluid/CDataDrivenFluid.hpp" +#include "../../include/fluid/CNasaGas.hpp" enum { SIZE_ARR_NORM = 8 @@ -914,6 +915,21 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, } break; + + case NASA_GAS: + FluidModel = new CNasaGas(Gamma, config->GetGas_Constant(), config->GetCompute_Entropy()); + FluidModel->SetCpModel(config, 1.0); + if (free_stream_temp) { + FluidModel->SetTDState_PT(Pressure_FreeStream, Temperature_FreeStream); + Density_FreeStream = FluidModel->GetDensity(); + config->SetDensity_FreeStream(Density_FreeStream); + } + else { + FluidModel->SetTDState_Prho(Pressure_FreeStream, Density_FreeStream ); + Temperature_FreeStream = FluidModel->GetTemperature(); + config->SetTemperature_FreeStream(Temperature_FreeStream); + } + break; } Mach2Vel_FreeStream = FluidModel->GetSoundSpeed(); @@ -1121,6 +1137,12 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, FluidModel = new CDataDrivenFluid(config); FluidModel->SetEnergy_Prho(Pressure_FreeStreamND, Density_FreeStreamND); break; + + case NASA_GAS: + FluidModel = new CNasaGas(Gamma, Gas_ConstantND, config->GetCompute_Entropy()); + FluidModel->SetCpModel(config, config->GetTemperature_Ref()); + FluidModel->SetEnergy_Prho(Pressure_FreeStreamND, Density_FreeStreamND); + break; } Energy_FreeStreamND = FluidModel->GetStaticEnergy() + 0.5*ModVel_FreeStreamND*ModVel_FreeStreamND; @@ -1291,6 +1313,9 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, case PR_GAS: ModelTable << "PR_GAS"; break; + case NASA_GAS: + ModelTable << "NASA_GAS"; + break; case COOLPROP: ModelTable << "CoolProp library"; break; diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index f7b5d71f358d..eea8e104a4a9 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -33,7 +33,8 @@ unsigned long EulerNPrimVarGrad(const CConfig *config, unsigned long ndim) { if (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) return ndim + 1; const bool ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || - config->GetKind_FluidModel() == IDEAL_GAS; + config->GetKind_FluidModel() == IDEAL_GAS || + config->GetKind_FluidModel() == NASA_GAS; const bool low_mach = config->Low_Mach_Correction(); if (ideal_gas && !low_mach && (config->GetKind_Upwind_Flow() == UPWIND::ROE || config->GetKind_Upwind_Flow() == UPWIND::MSW)) {