diff --git a/AUTHORS.md b/AUTHORS.md index d82ccb599fe..f45e0812596 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -50,6 +50,7 @@ Alessandro Gastaldi Aman uz zaman Baig Amit Sachdeva Ana Lourenco +Andrea Rausa Andrew Burkett Andrew Wendorff Angelo Passariello diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index ae07bef8455..a9b84958313 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -411,6 +411,22 @@ class CFVMFlowSolverBase : public CSolver { void LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo, su2double* RestartSolution = nullptr, unsigned short nVar_Restart = 0); + /*! + * \brief Compute the reference values and the non-dimensional free-stream values that are common + * to the compressible solvers, and store them in the config. + * \note The reference pressure, density and temperature follow from the option REF_DIMENSIONALIZATION. + * All other reference values are derived from these three, with a reference length of one. + * \param[in,out] config - Definition of the particular problem. + * \param[in] Pressure_FreeStream - Dimensional free-stream pressure. + * \param[in] Density_FreeStream - Dimensional free-stream density. + * \param[in] Temperature_FreeStream - Dimensional free-stream temperature. + * \param[in] Viscosity_FreeStream - Dimensional free-stream laminar viscosity. + * \param[in] ModVel_FreeStream - Dimensional magnitude of the free-stream velocity. + */ + void SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream, su2double Density_FreeStream, + su2double Temperature_FreeStream, su2double Viscosity_FreeStream, + su2double ModVel_FreeStream) const; + /*! * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c628..865eb0d57d2 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -850,6 +850,112 @@ void CFVMFlowSolverBase::UpdateCustomBoundaryConditions( } } +template +void CFVMFlowSolverBase::SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream, + su2double Density_FreeStream, + su2double Temperature_FreeStream, + su2double Viscosity_FreeStream, + su2double ModVel_FreeStream) const { + + const su2double Mach = config->GetMach(); + + /*--- Reference pressure, density and temperature. ---*/ + + su2double Pressure_Ref = 0.0, Density_Ref = 0.0, Temperature_Ref = 0.0; + + if (config->GetRef_NonDim() == DIMENSIONAL) { + Pressure_Ref = 1.0; + Density_Ref = 1.0; + Temperature_Ref = 1.0; + } + else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { + Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 + } + else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { + Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 + } + else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { + Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) + Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 + Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 + } + config->SetPressure_Ref(Pressure_Ref); + config->SetDensity_Ref(Density_Ref); + config->SetTemperature_Ref(Temperature_Ref); + + /*--- Reference values that are derived from the three above. ---*/ + + const su2double Length_Ref = 1.0; + const su2double Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); + const su2double Time_Ref = Length_Ref/Velocity_Ref; + const su2double Omega_Ref = Velocity_Ref/Length_Ref; + const su2double Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; + const su2double Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; + const su2double Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); + const su2double Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; + const su2double Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; + const su2double Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); + + config->SetLength_Ref(Length_Ref); + config->SetVelocity_Ref(Velocity_Ref); + config->SetTime_Ref(Time_Ref); + config->SetOmega_Ref(Omega_Ref); + config->SetForce_Ref(Force_Ref); + config->SetHeat_Flux_Ref(Heat_Flux_Ref); + config->SetGas_Constant_Ref(Gas_Constant_Ref); + config->SetViscosity_Ref(Viscosity_Ref); + config->SetConductivity_Ref(Conductivity_Ref); + config->SetFroude(Froude); + + /*--- Divide by reference values, to compute the non-dimensional free-stream values. ---*/ + + const su2double Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); + config->SetPressure_FreeStreamND(Pressure_FreeStreamND); + + const su2double Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); + config->SetDensity_FreeStreamND(Density_FreeStreamND); + + su2double Velocity_FreeStreamND[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; + config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); + } + + const su2double Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); + config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); + + const su2double Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; + config->SetGas_ConstantND(Gas_ConstantND); + + const su2double ModVel_FreeStreamND = GeometryToolbox::Norm(nDim, Velocity_FreeStreamND); + config->SetModVel_FreeStreamND(ModVel_FreeStreamND); + + const su2double Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; + config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); + + /*--- Free-stream turbulence quantities, from the turbulence intensity and the viscosity ratio. ---*/ + + const su2double TurbIntensity = config->GetTurbulenceIntensity_FreeStream(); + + const su2double Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*TurbIntensity*TurbIntensity); + config->SetTke_FreeStream(Tke_FreeStream); + + const su2double Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*TurbIntensity*TurbIntensity); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + const su2double Omega_FreeStream = + Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); + config->SetOmega_FreeStream(Omega_FreeStream); + + const su2double Omega_FreeStreamND = + Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + config->SetOmega_FreeStreamND(Omega_FreeStreamND); +} + template void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo, su2double* SolutionRestart, diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index f1955c46a1a..cdd8acac0a4 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -60,6 +60,48 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseSetdPdU (nodes_i->GetdPdU(iPoint), nodes_j->GetdPdU(jPoint)); + numerics->SetdTdU (nodes_i->GetdTdU(iPoint), nodes_j->GetdTdU(jPoint)); + numerics->SetdTvedU(nodes_i->GetdTvedU(iPoint), nodes_j->GetdTvedU(jPoint)); + numerics->SetEve (nodes_i->GetEve(iPoint), nodes_j->GetEve(jPoint)); + numerics->SetCvve (nodes_i->GetCvve(iPoint), nodes_j->GetCvve(jPoint)); + } + + /*! + * \brief Pass the transport properties of two points of the solver to a viscous numerics object. + * \param[in,out] numerics - Numerics object that receives the data. + * \param[in] iPoint - Index of point i. + * \param[in] jPoint - Index of point j. + */ + void SetNumericsTransportState(CNumerics *numerics, unsigned long iPoint, unsigned long jPoint); + + /*! + * \brief Compute the viscous flux between a boundary point and a prescribed boundary state, + * and add it to the residual and to the Jacobian. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] visc_numerics - Description of the numerical method for viscous terms. + * \param[in] config - Definition of the particular problem. + * \param[in] iPoint - Index of the boundary point. + * \param[in] Point_Normal - Index of the closest interior neighbor of the boundary point. + * \param[in] Normal - Normal vector of the boundary face, the same that is passed to the convective numerics. + * \param[in] nodes_bound - Variable container that holds the boundary state. + * \param[in] iPoint_bound - Index of the boundary state in its container. + */ + void AddViscousBoundaryResidual(CGeometry *geometry, CNumerics *visc_numerics, const CConfig *config, + unsigned long iPoint, unsigned long Point_Normal, const su2double *Normal, + CNEMOEulerVariable *nodes_bound, unsigned long iPoint_bound); + /*! * \brief Set the maximum value of the eigenvalue. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CNEMONSSolver.hpp b/SU2_CFD/include/solvers/CNEMONSSolver.hpp index 0e4fe68a5b2..75bc17906d4 100644 --- a/SU2_CFD/include/solvers/CNEMONSSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMONSSolver.hpp @@ -72,6 +72,29 @@ class CNEMONSSolver final : public CNEMOEulerSolver { */ void SetTau_Wall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config); + /*! + * \brief Impose the velocity of a wall point as a strong (Dirichlet) boundary condition. + * The velocity is fixed and the momentum residual and truncation error are set to zero. + * \param[in] iPoint - Index of the wall point. + * \param[in] velocity - Velocity to impose (zero for no-slip, slip velocity otherwise). + */ + void SetWallVelocity_Strong(unsigned long iPoint, const su2double *velocity); + + /*! + * \brief Compute the residual of the two energy equations that weakly imposes the wall temperatures. + * \param[in] iPoint - Index of the wall point. + * \param[in] jPoint - Index of the closest normal neighbor of the wall point. + * \param[in] Twall - Target translational-rotational temperature at the wall. + * \param[in] Tvewall - Target vibrational-electronic temperature at the wall. + * \param[in] pcontrol - "Proportional control" constant of the temperature difference term. + * \param[in] area - Area of the dual-grid face on the wall. + * \param[in] dist - Distance between the wall point and its normal neighbor. + * \param[out] residual - Residual, only the two energy entries are set. + */ + void SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double Twall, su2double Tvewall, + su2double pcontrol, su2double area, su2double dist, + su2double *residual) const; + public: /*! diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e5565..055c8785674 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -807,19 +807,11 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes SU2_ZONE_SCOPED su2double Temperature_FreeStream = 0.0, Mach2Vel_FreeStream = 0.0, ModVel_FreeStream = 0.0, - Energy_FreeStream = 0.0, ModVel_FreeStreamND = 0.0, Velocity_Reynolds = 0.0, - Omega_FreeStream = 0.0, Omega_FreeStreamND = 0.0, Viscosity_FreeStream = 0.0, + Energy_FreeStream = 0.0, Velocity_Reynolds = 0.0, Viscosity_FreeStream = 0.0, Thermal_Conductivity_FreeStream = 0.0, SpecificHeat_Cp_FreeStream = 0.0, Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, Tke_FreeStream = 0.0, Re_ThetaT_FreeStream = 0.0, - Length_Ref = 0.0, Density_Ref = 0.0, Pressure_Ref = 0.0, Velocity_Ref = 0.0, - Temperature_Ref = 0.0, Time_Ref = 0.0, Omega_Ref = 0.0, Force_Ref = 0.0, - Gas_Constant_Ref = 0.0, Viscosity_Ref = 0.0, Conductivity_Ref = 0.0, Energy_Ref= 0.0, - Froude = 0.0, Pressure_FreeStreamND = 0.0, Density_FreeStreamND = 0.0, - Temperature_FreeStreamND = 0.0, Gas_ConstantND = 0.0, - Velocity_FreeStreamND[3] = {0.0, 0.0, 0.0}, Viscosity_FreeStreamND = 0.0, - Thermal_Conductivity_FreeStreamND = 0.0, SpecificHeat_Cp_FreeStreamND = 0.0, - Tke_FreeStreamND = 0.0, Energy_FreeStreamND = 0.0, - Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, TgammaR = 0.0, Heat_Flux_Ref = 0.0; + Energy_Ref= 0.0, Thermal_Conductivity_FreeStreamND = 0.0, SpecificHeat_Cp_FreeStreamND = 0.0, + Energy_FreeStreamND = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, TgammaR = 0.0; unsigned short iDim; @@ -1023,80 +1015,20 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); - /*--- Compute non dimensional quantities. By definition, - Lref is one because we have converted the grid to meters. ---*/ + /*--- Compute the reference values and the non-dimensional free-stream values, + and store them in the config. By definition, Lref is one because we have + converted the grid to meters. ---*/ - if (config->GetRef_NonDim() == DIMENSIONAL) { - Pressure_Ref = 1.0; - Density_Ref = 1.0; - Temperature_Ref = 1.0; - } - else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { - Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { - Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { - Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - } - config->SetPressure_Ref(Pressure_Ref); - config->SetDensity_Ref(Density_Ref); - config->SetTemperature_Ref(Temperature_Ref); - - Length_Ref = 1.0; config->SetLength_Ref(Length_Ref); - Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); config->SetVelocity_Ref(Velocity_Ref); - Time_Ref = Length_Ref/Velocity_Ref; config->SetTime_Ref(Time_Ref); - Omega_Ref = Velocity_Ref/Length_Ref; config->SetOmega_Ref(Omega_Ref); - Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; config->SetForce_Ref(Force_Ref); - Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; config->SetHeat_Flux_Ref(Heat_Flux_Ref); - Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); config->SetGas_Constant_Ref(Gas_Constant_Ref); - Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; config->SetViscosity_Ref(Viscosity_Ref); - Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; config->SetConductivity_Ref(Conductivity_Ref); - Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); config->SetFroude(Froude); - - /*--- Divide by reference values, to compute the non-dimensional free-stream values ---*/ - - Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); config->SetPressure_FreeStreamND(Pressure_FreeStreamND); - Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); config->SetDensity_FreeStreamND(Density_FreeStreamND); - - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); - } - - Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); + SetCompressibleNonDimValues(config, Pressure_FreeStream, Density_FreeStream, Temperature_FreeStream, + Viscosity_FreeStream, ModVel_FreeStream); - Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; config->SetGas_ConstantND(Gas_ConstantND); + /*--- Non-dimensional free-stream values that are specific to this solver ---*/ - - ModVel_FreeStreamND = 0.0; - for (iDim = 0; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; - ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); config->SetModVel_FreeStreamND(ModVel_FreeStreamND); - - Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); - Thermal_Conductivity_FreeStreamND = Thermal_Conductivity_FreeStream / Conductivity_Ref; + Thermal_Conductivity_FreeStreamND = Thermal_Conductivity_FreeStream / config->GetThermal_Conductivity_Ref(); config->SetThermalConductivity_FreeStreamND(Thermal_Conductivity_FreeStreamND); - SpecificHeat_Cp_FreeStreamND = SpecificHeat_Cp_FreeStream / Gas_Constant_Ref; + SpecificHeat_Cp_FreeStreamND = SpecificHeat_Cp_FreeStream / config->GetGas_Constant_Ref(); config->SetSpecificHeatCp_FreeStreamND(SpecificHeat_Cp_FreeStreamND); - Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); - - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStreamND(Omega_FreeStreamND); - if (config->GetTurbulenceIntensity_FreeStream() *100 <= 1.3) { if (config->GetTurbulenceIntensity_FreeStream() *100 >=0.027) { Re_ThetaT_FreeStream = (1173.51-589.428*config->GetTurbulenceIntensity_FreeStream() *100+0.2196/ @@ -1111,9 +1043,20 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes } config->SetReThetaT_FreeStream(Re_ThetaT_FreeStream); - const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); + const su2double MassDiffusivityND = + config->GetDiffusivity_Constant() / (config->GetVelocity_Ref() * config->GetLength_Ref()); config->SetDiffusivity_ConstantND(MassDiffusivityND); + /*--- Local copies of values stored in the config ---*/ + + const su2double Time_Ref = config->GetTime_Ref(); + const su2double Froude = config->GetFroude(); + const su2double Gas_ConstantND = config->GetGas_ConstantND(); + const su2double Pressure_FreeStreamND = config->GetPressure_FreeStreamND(); + const su2double Density_FreeStreamND = config->GetDensity_FreeStreamND(); + const su2double ModVel_FreeStreamND = config->GetModVel_FreeStreamND(); + const su2double Tke_FreeStreamND = config->GetTke_FreeStreamND(); + /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ /*--- Auxilary (dimensional) FluidModel no longer needed. ---*/ diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 2a09a7a1f42..449ea5cf46b 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -439,12 +439,8 @@ void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_c numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint)); /*--- Pass supplementary information to CNumerics ---*/ - numerics->SetdPdU ( nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU ( nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU( nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve ( nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve ( nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); - numerics->SetGamma ( nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); + numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); /*--- Set the largest convective eigenvalue ---*/ numerics->SetLambda(nodes->GetLambda(iPoint), nodes->GetLambda(jPoint)); @@ -520,12 +516,8 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con numerics->SetPrimitive (V_i, V_j); numerics->SetConservative(U_i, U_j); - numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); - numerics->SetGamma (nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); + numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(jPoint)); } else { @@ -991,21 +983,11 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short su2double Temperature_FreeStream = 0.0, Temperature_ve_FreeStream = 0.0, Mach2Vel_FreeStream = 0.0, - ModVel_FreeStream = 0.0, Energy_FreeStream = 0.0, ModVel_FreeStreamND = 0.0, - Velocity_Reynolds = 0.0, Omega_FreeStream = 0.0, Omega_FreeStreamND = 0.0, + ModVel_FreeStream = 0.0, Energy_FreeStream = 0.0, Velocity_Reynolds = 0.0, Viscosity_FreeStream = 0.0, Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, - Tke_FreeStream = 0.0, Length_Ref = 0.0, Density_Ref = 0.0, - Pressure_Ref = 0.0, Velocity_Ref = 0.0, Temperature_Ref = 0.0, - Temperature_ve_Ref = 0.0, Time_Ref = 0.0, Omega_Ref = 0.0, - Force_Ref = 0.0, Gas_Constant_Ref = 0.0, Viscosity_Ref = 0.0, - Conductivity_Ref = 0.0, Energy_Ref = 0.0, Pressure_FreeStreamND = 0.0, - Energy_FreeStreamND = 0.0, Temperature_FreeStreamND = 0.0, Temperature_ve_FreeStreamND = 0.0, - Gas_ConstantND = 0.0, Viscosity_FreeStreamND = 0.0, sqvel = 0.0, - Tke_FreeStreamND = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, - soundspeed = 0.0, GasConstant_Inf = 0.0, Froude = 0.0, - Density_FreeStreamND = 0.0, Heat_Flux_Ref = 0.0; - - su2double Velocity_FreeStreamND[3] = {0.0, 0.0, 0.0}; + Tke_FreeStream = 0.0, Energy_Ref = 0.0, Energy_FreeStreamND = 0.0, + sqvel = 0.0, Total_UnstTimeND = 0.0, Delta_UnstTimeND = 0.0, + soundspeed = 0.0, GasConstant_Inf = 0.0; /*--- Local variables ---*/ const su2double Alpha= config->GetAoA()*PI_NUMBER/180.0; @@ -1129,79 +1111,27 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short /*-- Compute the freestream energy. ---*/ if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); - /*--- Compute non dimensional quantities. By definition, - Lref is one because we have converted the grid to meters. ---*/ + /*--- Compute the reference values and the non-dimensional free-stream values, + and store them in the config. By definition, Lref is one because we have + converted the grid to meters. ---*/ - if (config->GetRef_NonDim() == DIMENSIONAL) { - Pressure_Ref = 1.0; - Density_Ref = 1.0; - Temperature_Ref = 1.0; - Temperature_ve_Ref = 1.0; - } - else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) { - Pressure_Ref = Pressure_FreeStream; // Pressure_FreeStream = 1.0 - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temperature_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temperature_ve_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) { - Pressure_Ref = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temp_ve_FreeStream = 1.0 - } - else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) { - Pressure_Ref = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2) - Density_Ref = Density_FreeStream; // Density_FreeStream = 1.0 - Temperature_Ref = Temperature_FreeStream; // Temp_FreeStream = 1.0 - Temperature_ve_Ref = Temperature_ve_FreeStream; // Temp_ve_FreeStream = 1.0 - } - config->SetPressure_Ref(Pressure_Ref); - config->SetDensity_Ref(Density_Ref); - config->SetTemperature_Ref(Temperature_Ref); - config->SetTemperature_ve_Ref(Temperature_ve_Ref); - - Length_Ref = 1.0; config->SetLength_Ref(Length_Ref); - Velocity_Ref = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref()); config->SetVelocity_Ref(Velocity_Ref); - Time_Ref = Length_Ref/Velocity_Ref; config->SetTime_Ref(Time_Ref); - Omega_Ref = Velocity_Ref/Length_Ref; config->SetOmega_Ref(Omega_Ref); - Force_Ref = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; config->SetForce_Ref(Force_Ref); - Heat_Flux_Ref = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref; config->SetHeat_Flux_Ref(Heat_Flux_Ref); - Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref(); config->SetGas_Constant_Ref(Gas_Constant_Ref); - Viscosity_Ref = config->GetDensity_Ref()*Velocity_Ref*Length_Ref; config->SetViscosity_Ref(Viscosity_Ref); - Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref; config->SetConductivity_Ref(Conductivity_Ref); - Froude = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref); config->SetFroude(Froude); - - /*--- Divide by reference values, to compute the non-dimensional free-stream values ---*/ - - Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref(); config->SetPressure_FreeStreamND(Pressure_FreeStreamND); - Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref(); config->SetDensity_FreeStreamND(Density_FreeStreamND); - - for (auto iDim = 0ul; iDim < nDim; iDim++) { - Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref; config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim); - } + SetCompressibleNonDimValues(config, Pressure_FreeStream, Density_FreeStream, Temperature_FreeStream, + Viscosity_FreeStream, ModVel_FreeStream); - Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref(); config->SetTemperature_FreeStreamND(Temperature_FreeStreamND); - Temperature_ve_FreeStreamND = Temperature_ve_FreeStream/config->GetTemperature_ve_Ref(); config->SetTemperature_ve_FreeStreamND(Temperature_ve_FreeStreamND); - Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref; config->SetGas_ConstantND(Gas_ConstantND); + /*--- Reference value and non-dimensional free-stream value of the vib.-el. temperature ---*/ - ModVel_FreeStreamND = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) ModVel_FreeStreamND += Velocity_FreeStreamND[iDim]*Velocity_FreeStreamND[iDim]; - ModVel_FreeStreamND = sqrt(ModVel_FreeStreamND); config->SetModVel_FreeStreamND(ModVel_FreeStreamND); - - Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); - - Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); + const su2double Temperature_ve_Ref = (config->GetRef_NonDim() == DIMENSIONAL) ? 1.0 : Temperature_ve_FreeStream; + config->SetTemperature_ve_Ref(Temperature_ve_Ref); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); + const su2double Temperature_ve_FreeStreamND = Temperature_ve_FreeStream/config->GetTemperature_ve_Ref(); + config->SetTemperature_ve_FreeStreamND(Temperature_ve_FreeStreamND); - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); + /*--- Local copies of values stored in the config ---*/ - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStreamND(Omega_FreeStreamND); + const su2double Time_Ref = config->GetTime_Ref(); + const su2double Froude = config->GetFroude(); + const su2double ModVel_FreeStreamND = config->GetModVel_FreeStreamND(); + const su2double Tke_FreeStreamND = config->GetTke_FreeStreamND(); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ @@ -1405,6 +1335,67 @@ void CNEMOEulerSolver::SetReferenceValues(const CConfig& config) { } +void CNEMOEulerSolver::SetNumericsTransportState(CNumerics *numerics, unsigned long iPoint, unsigned long jPoint) { + + /*--- Species diffusion coefficients ---*/ + numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), + nodes->GetDiffusionCoeff(jPoint)); + + /*--- Laminar viscosity ---*/ + numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), + nodes->GetLaminarViscosity(jPoint)); + + /*--- Eddy viscosity ---*/ + numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), + nodes->GetEddyViscosity(jPoint)); + + /*--- Thermal conductivity ---*/ + numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), + nodes->GetThermalConductivity(jPoint)); + + /*--- Vib-el. thermal conductivity ---*/ + numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), + nodes->GetThermalConductivity_ve(jPoint)); +} + +void CNEMOEulerSolver::AddViscousBoundaryResidual(CGeometry *geometry, CNumerics *visc_numerics, const CConfig *config, + unsigned long iPoint, unsigned long Point_Normal, + const su2double *Normal, CNEMOEulerVariable *nodes_bound, + unsigned long iPoint_bound) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + /*--- Set the normal vector and the coordinates. The second point is the + reflection of the closest interior node with respect to the boundary point. ---*/ + su2double Coord_Reflected[MAXNDIM]; + GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), + geometry->nodes->GetCoord(iPoint), Coord_Reflected); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected); + visc_numerics->SetNormal(Normal); + + /*--- Conservative and primitive variables, and gradient ---*/ + visc_numerics->SetConservative(nodes->GetSolution(iPoint), + nodes_bound->GetSolution(iPoint_bound)); + visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), + nodes_bound->GetPrimitive(iPoint_bound)); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes_bound->GetGradient_Primitive(iPoint_bound)); + + /*--- Pass supplementary information to CNumerics ---*/ + SetNumericsThermoState(visc_numerics, nodes, iPoint, nodes_bound, iPoint_bound); + + /*--- Transport properties, the values of the boundary point are used for both states ---*/ + SetNumericsTransportState(visc_numerics, iPoint, iPoint); + + /*--- Compute and update residual ---*/ + auto residual = visc_numerics->ComputeResidual(config); + LinSysRes.SubtractBlock(iPoint, residual); + + /*--- Jacobian contribution for implicit integration ---*/ + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); +} + void CNEMOEulerSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_ZONE_SCOPED @@ -1539,12 +1530,8 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, node_infty->GetPrimitive(0)); /*--- Pass supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU (nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - conv_numerics->SetdTdU (nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_infty->GetdTvedU(0)); - conv_numerics->SetEve (nodes->GetEve(iPoint), node_infty->GetEve(0)); - conv_numerics->SetCvve (nodes->GetCvve(iPoint), node_infty->GetCvve(0)); - conv_numerics->SetGamma (nodes->GetGamma(iPoint), node_infty->GetGamma(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); + conv_numerics->SetGamma(nodes->GetGamma(iPoint), node_infty->GetGamma(0)); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -1561,62 +1548,8 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - const bool viscous = config->GetViscous(); - if (viscous) { - su2double Coord_Reflected[MAXNDIM]; - GeometryToolbox::PointPointReflect( - nDim, geometry->nodes->GetCoord(Point_Normal), - geometry->nodes->GetCoord(iPoint), Coord_Reflected); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - Coord_Reflected); - visc_numerics->SetNormal(Normal); - - /*--- Primitive variables, and gradient ---*/ - visc_numerics->SetConservative(nodes->GetSolution(iPoint), - node_infty->GetSolution(0)); - visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), - node_infty->GetPrimitive(0)); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - node_infty->GetGradient_Primitive(0)); - - /*--- Pass supplementary information to CNumerics ---*/ - visc_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - visc_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - visc_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - node_infty->GetdTvedU(0)); - visc_numerics->SetEve(nodes->GetEve(iPoint), node_infty->GetEve(0)); - visc_numerics->SetCvve(nodes->GetCvve(iPoint), node_infty->GetCvve(0)); - - /*--- Species diffusion coefficients ---*/ - visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(iPoint)); - - /*--- Laminar viscosity ---*/ - visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(iPoint)); - - /*--- Eddy viscosity ---*/ - visc_numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(iPoint)); - - /*--- Thermal conductivity ---*/ - visc_numerics->SetThermalConductivity( - nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(iPoint)); - - /*--- Vib-el. thermal conductivity ---*/ - visc_numerics->SetThermalConductivity_ve( - nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(iPoint)); - - /*--- Compute and update residual ---*/ - auto residual = visc_numerics->ComputeResidual(config); - - LinSysRes.SubtractBlock(iPoint, residual); - if (implicit) { - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); - } - } + if (config->GetViscous()) + AddViscousBoundaryResidual(geometry, visc_numerics, config, iPoint, Point_Normal, Normal, node_infty, 0); } } } @@ -1624,283 +1557,8 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { SU2_ZONE_SCOPED + /*--- The subsonic inlet is not implemented for a two-temperature mixture. ---*/ SU2_MPI::Error("BC_INLET: Not operational in NEMO.", CURRENT_FUNCTION); - - unsigned short RHO_INDEX, nSpecies; - - su2double T_Total, P_Total, Velocity[3], Velocity2, H_Total, Temperature, Riemann, - Pressure, Density, Energy, Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag, - alpha, aa, bb, cc, dd, Area, UnitNormal[3] = {0.0}; - - su2double Flow_Dir[MAXNDIM] = {}; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - bool dynamic_grid = config->GetGrid_Movement(); - su2double Two_Gamma_M1 = 2.0/Gamma_Minus_One; - su2double Gas_Constant = config->GetGas_ConstantND(); - INLET_TYPE Kind_Inlet = config->GetKind_Inlet(); - - auto *U_domain = new su2double[nVar]; auto *U_inlet = new su2double[nVar]; - auto *V_domain = new su2double[nPrimVar]; auto *V_inlet = new su2double[nPrimVar]; - auto *Normal = new su2double[nDim]; - - nSpecies = config->GetnSpecies(); - auto *Spec_Density = new su2double[nSpecies]; - for (auto iSpecies = 0ul; iSpeciesGetRhoIndex(); - - /*--- Loop over all the vertices on this boundary marker ---*/ - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (auto iDim = 0ul; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); - - Area = GeometryToolbox::Norm(nDim, Normal); - - for (auto iDim = 0ul; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; - - /*--- Retrieve solution at this boundary node ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) U_domain[iVar] = nodes->GetSolution(iPoint, iVar); - for (auto iVar = 0ul; iVar < nPrimVar; iVar++) V_domain[iVar] = nodes->GetPrimitive(iPoint,iVar); - - const su2double* dir = Inlet_FlowDir[val_marker][iVertex]; - const su2double mag = GeometryToolbox::Norm(nDim, dir); - if (config->GetInletUseNormal()) { - for (auto iDim = 0u; iDim < nDim; iDim++) - Flow_Dir[iDim] = -UnitNormal[iDim]; - } else { - for (auto iDim = 0u; iDim < nDim; iDim++) - Flow_Dir[iDim] = dir[iDim]/mag; - } - - /*--- Build the fictitious intlet state based on characteristics ---*/ - - /*--- Subsonic inflow: there is one outgoing characteristic (u-c), - therefore we can specify all but one state variable at the inlet. - The outgoing Riemann invariant provides the final piece of info. - Adapted from an original implementation in the Stanford University - multi-block (SUmb) solver in the routine bcSubsonicInflow.f90 - written by Edwin van der Weide, last modified 04-20-2009. ---*/ - - switch (Kind_Inlet) { - - /*--- Total properties have been specified at the inlet. ---*/ - case INLET_TYPE::TOTAL_CONDITIONS: - - /*--- Retrieve the specified total conditions for this inlet. ---*/ - P_Total = Inlet_Ptotal[val_marker][iVertex]; - T_Total = Inlet_Ttotal[val_marker][iVertex]; - - /*--- Non-dim. the inputs if necessary. ---*/ - P_Total /= config->GetPressure_Ref(); - T_Total /= config->GetTemperature_Ref(); - - /*--- Store primitives and set some variables for clarity. ---*/ - Density = V_domain[RHO_INDEX]; - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = U_domain[nSpecies+iDim]/Density; - - Velocity2 = GeometryToolbox::SquaredNorm(nDim, Velocity); - Energy = U_domain[nVar-2]/Density; - Pressure = Gamma_Minus_One*Density*(Energy-0.5*Velocity2); - H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; - SoundSpeed2 = Gamma*Pressure/Density; - - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ - Riemann = 2.0*sqrt(SoundSpeed2)/Gamma_Minus_One; - for (auto iDim = 0ul; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; - - /*--- Total speed of sound ---*/ - SoundSpeed_Total2 = Gamma_Minus_One*(H_Total - (Energy + Pressure/Density)+0.5*Velocity2) + SoundSpeed2; - - /*--- Dot product of normal and flow direction. This should - be negative due to outward facing boundary normal convention. ---*/ - alpha = GeometryToolbox::DotProduct(nDim, UnitNormal, Flow_Dir); - - /*--- Coefficients in the quadratic equation for the velocity ---*/ - aa = 1.0 + 0.5*Gamma_Minus_One*alpha*alpha; - bb = -1.0*Gamma_Minus_One*alpha*Riemann; - cc = 0.5*Gamma_Minus_One*Riemann*Riemann - -2.0*SoundSpeed_Total2/Gamma_Minus_One; - - /*--- Solve quadratic equation for velocity magnitude. Value must - be positive, so the choice of root is clear. ---*/ - dd = bb*bb - 4.0*aa*cc; - dd = sqrt(max(0.0,dd)); - Vel_Mag = (-bb + dd)/(2.0*aa); - Vel_Mag = max(0.0,Vel_Mag); - Velocity2 = Vel_Mag*Vel_Mag; - - /*--- Compute speed of sound from total speed of sound eqn. ---*/ - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - - /*--- Mach squared (cut between 0-1), use to adapt velocity ---*/ - Mach2 = Velocity2/SoundSpeed2; - Mach2 = min(1.0,Mach2); - Velocity2 = Mach2*SoundSpeed2; - Vel_Mag = sqrt(Velocity2); - SoundSpeed2 = SoundSpeed_Total2 - 0.5*Gamma_Minus_One*Velocity2; - - /*--- Compute new velocity vector at the inlet ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = Vel_Mag*Flow_Dir[iDim]; - - /*--- Static temperature from the speed of sound relation ---*/ - Temperature = SoundSpeed2/(Gamma*Gas_Constant); - //NEED TVE AS WELL - - /*--- Static pressure using isentropic relation at a point ---*/ - Pressure = P_Total*pow((Temperature/T_Total),Gamma/Gamma_Minus_One); - - /*--- Density at the inlet from the gas law ---*/ - Density = Pressure/(Gas_Constant*Temperature); - //NEED SPECIES DENSITIES - - /*--- Using pressure, density, & velocity, compute the energy ---*/ - Energy = Pressure/(Density*Gamma_Minus_One)+0.5*Velocity2; - //NEED EVE AS WELL - - /*--- Conservative variables, using the derived quantities ---*/ - for (auto iSpecies = 0ul; iSpeciesGetDensity_Ref(); - Vel_Mag /= config->GetVelocity_Ref(); - - /*--- Get primitives from current inlet state. ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) - Velocity[iDim] = nodes->GetVelocity(iPoint, iDim); - Pressure = nodes->GetPressure(iPoint); - SoundSpeed2 = Gamma*Pressure/U_domain[0]; - - /*--- Compute the acoustic Riemann invariant that is extrapolated - from the domain interior. ---*/ - Riemann = Two_Gamma_M1*sqrt(SoundSpeed2); - for (auto iDim = 0ul; iDim < nDim; iDim++) - Riemann += Velocity[iDim]*UnitNormal[iDim]; - - /*--- Speed of sound squared for fictitious inlet state ---*/ - SoundSpeed2 = Riemann; - for (auto iDim = 0ul; iDim < nDim; iDim++) - SoundSpeed2 -= Vel_Mag*Flow_Dir[iDim]*UnitNormal[iDim]; - - SoundSpeed2 = max(0.0,0.5*Gamma_Minus_One*SoundSpeed2); - SoundSpeed2 = SoundSpeed2*SoundSpeed2; - - /*--- Pressure for the fictitious inlet state ---*/ - Pressure = SoundSpeed2*Density/Gamma; - - /*--- Energy for the fictitious inlet state ---*/ - Energy = Pressure/(Density*Gamma_Minus_One)+0.5*Vel_Mag*Vel_Mag; - - /*--- Conservative variables, using the derived quantities ---*/ - U_inlet[0] = Density; - for (auto iDim = 0ul; iDim < nDim; iDim++) - U_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]*Density; - U_inlet[nDim+1] = Energy*Density; - - /*--- Primitive variables, using the derived quantities ---*/ - V_inlet[0] = Pressure / ( Gas_Constant * Density); - for (auto iDim = 0ul; iDim < nDim; iDim++) - V_inlet[iDim+1] = Vel_Mag*Flow_Dir[iDim]; - V_inlet[nDim+1] = Pressure; - V_inlet[nDim+2] = Density; - - break; - - default: - SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION); - break; - } - - /*--- Set various quantities in the solver class ---*/ - conv_numerics->SetConservative(U_domain, U_inlet); - - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - - /*--- Compute the residual using an upwind scheme ---*/ - auto residual = conv_numerics->ComputeResidual(config); - - LinSysRes.AddBlock(iPoint, residual); - - /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - -// /*--- Viscous contribution ---*/ -// if (viscous) { - -// /*--- Set the normal vector and the coordinates ---*/ -// visc_numerics->SetNormal(Normal); -// su2double Coord_Reflected[MAXNDIM]; -// GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), -// geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// visc_numerics->SetCoord(geometry->nodes->GetCoord(), Coord_Reflected); - -// /*--- Primitive variables, and gradient ---*/ -// visc_numerics->SetPrimitive(V_domain, V_inlet); -// visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(), nodes->GetGradient_Primitive()); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), nodes->GetLaminarViscosity()); - -// /*--- Compute and update residual ---*/ -// visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); -// LinSysRes.SubtractBlock(iPoint, Residual); - -// /*--- Jacobian contribution for implicit integration ---*/ -// if (implicit) -// Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); -// } - } - } - - /*--- Free locally allocated memory ---*/ - delete [] U_domain; - delete [] U_inlet; - delete [] V_domain; - delete [] V_inlet; - delete [] Normal; - delete [] Spec_Density; } void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, @@ -1917,10 +1575,10 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); su2double *U_domain; - auto *U_outlet = new su2double[nVar]; + su2double U_outlet[MAXNVAR] = {0.0}; su2double *V_domain, *V_outlet; - auto *Normal = new su2double[nDim]; - auto *Ys = new su2double[nSpecies]; + su2double Normal[MAXNDIM] = {0.0}; + su2double Ys[MAXNVAR] = {0.0}; const unsigned short T_INDEX = nodes->GetTIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); @@ -2079,12 +1737,8 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); /*--- Passing supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU (nodes->GetdPdU(iPoint), node_infty->GetdPdU(0)); - conv_numerics->SetdTdU (nodes->GetdTdU(iPoint), node_infty->GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_infty->GetdTvedU(0)); - conv_numerics->SetEve (nodes->GetEve(iPoint), node_infty->GetEve(0)); - conv_numerics->SetCvve (nodes->GetCvve(iPoint), node_infty->GetCvve(0)); - conv_numerics->SetGamma (nodes->GetGamma(iPoint), node_infty->GetGamma(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, node_infty, 0); + conv_numerics->SetGamma(nodes->GetGamma(iPoint), node_infty->GetGamma(0)); if(config->GetKind_Upwind_Flow() == UPWIND::AUSMPLUSM) conv_numerics->SetSensor (nodes->GetSensor(iPoint), nodes->GetSensor(iPoint)); @@ -2095,64 +1749,8 @@ void CNEMOEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container /*--- Jacobian contribution for implicit integration ---*/ if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - - /*--- Viscous contribution ---*/ -// if (viscous) { - -// /*--- Set the normal vector and the coordinates ---*/ -// visc_numerics->SetNormal(Normal); -// su2double Coord_Reflected[MAXNDIM]; -// GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), -// geometry->nodes->GetCoord(iPoint), Coord_Reflected); -// visc_numerics->SetCoord(geometry->nodes->GetCoord(), Coord_Reflected); - -// /*--- Primitive variables, and gradient ---*/ -// visc_numerics->SetPrimitive(V_domain, V_outlet); -// visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(), nodes->GetGradient_Primitive()); - -// /*--- Conservative variables, and gradient ---*/ -// visc_numerics->SetConservative(U_domain, U_outlet); - -// /*--- Pass supplementary information to CNumerics ---*/ -// visc_numerics->SetdPdU(nodes->GetdPdU(), node_infty->GetdPdU()); -// visc_numerics->SetdTdU(nodes->GetdTdU(), node_infty->GetdTdU()); -// visc_numerics->SetdTvedU(nodes->GetdTvedU(), node_infty->GetdTvedU()); - -// /*--- Species diffusion coefficients ---*/ -// visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(), -// node_infty->GetDiffusionCoeff() ); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), -// node_infty->GetLaminarViscosity() ); - -// /*--- Thermal conductivity ---*/ -// visc_numerics->SetThermalConductivity(nodes->GetThermalConductivity(), -// node_infty->GetThermalConductivity()); - -// /*--- Vib-el. thermal conductivity ---*/ -// visc_numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(), -// node_infty->GetThermalConductivity_ve() ); - -// /*--- Laminar viscosity ---*/ -// visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(), nodes->GetLaminarViscosity()); - -// /*--- Compute and update residual ---*/ -// visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); -// LinSysRes.SubtractBlock(iPoint, Residual); - -// /*--- Jacobian contribution for implicit integration ---*/ -// if (implicit) -// Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); -// } } } - - /*--- Free locally allocated memory ---*/ - delete [] U_outlet; - delete [] Normal; - delete [] Ys; - } void CNEMOEulerSolver::BC_Supersonic_Inlet( @@ -2222,11 +1820,7 @@ void CNEMOEulerSolver::BC_Supersonic_Inlet( node_inlet.GetPrimitive(0)); /*--- Pass supplementary info to CNumerics ---*/ - conv_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_inlet.GetdPdU(0)); - conv_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_inlet.GetdTdU(0)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), node_inlet.GetdTvedU(0)); - conv_numerics->SetEve(nodes->GetEve(iPoint), node_inlet.GetEve(0)); - conv_numerics->SetCvve(nodes->GetCvve(iPoint), node_inlet.GetCvve(0)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, &node_inlet, 0); const bool dynamic_grid = config->GetGrid_Movement(); if (dynamic_grid) @@ -2243,60 +1837,8 @@ void CNEMOEulerSolver::BC_Supersonic_Inlet( Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); /*--- Viscous contribution ---*/ - const bool viscous = config->GetViscous(); - if (viscous) { - - /*--- Set the normal vector and the coordinates ---*/ - su2double Coord_Reflected[MAXNDIM]; - GeometryToolbox::PointPointReflect( - nDim, geometry->nodes->GetCoord(Point_Normal), - geometry->nodes->GetCoord(iPoint), Coord_Reflected); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - Coord_Reflected); - visc_numerics->SetNormal(Normal); - - /*--- Primitive variables, and gradient ---*/ - visc_numerics->SetConservative(nodes->GetSolution(iPoint), - node_inlet.GetSolution(0)); - visc_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), - node_inlet.GetPrimitive(0)); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - node_inlet.GetGradient_Primitive(0)); - - /*--- Pass supplementary information to CNumerics ---*/ - visc_numerics->SetdPdU(nodes->GetdPdU(iPoint), node_inlet.GetdPdU(0)); - visc_numerics->SetdTdU(nodes->GetdTdU(iPoint), node_inlet.GetdTdU(0)); - visc_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - node_inlet.GetdTvedU(0)); - visc_numerics->SetEve(nodes->GetEve(iPoint), node_inlet.GetEve(0)); - visc_numerics->SetCvve(nodes->GetCvve(iPoint), node_inlet.GetCvve(0)); - - /*--- Species diffusion coefficients ---*/ - visc_numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(iPoint)); - /*--- Laminar viscosity ---*/ - visc_numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(iPoint)); - /*--- Eddy viscosity ---*/ - visc_numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(iPoint)); - /*--- Thermal conductivity ---*/ - visc_numerics->SetThermalConductivity( - nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(iPoint)); - /*--- Vib-el. thermal conductivity ---*/ - visc_numerics->SetThermalConductivity_ve( - nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(iPoint)); - - /*--- Compute and update residual ---*/ - auto residual = visc_numerics->ComputeResidual(config); - LinSysRes.SubtractBlock(iPoint, residual); - - /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); - } + if (config->GetViscous()) + AddViscousBoundaryResidual(geometry, visc_numerics, config, iPoint, Point_Normal, Normal, &node_inlet, 0); } } @@ -2336,12 +1878,7 @@ void CNEMOEulerSolver::BC_Supersonic_Outlet( nodes->GetSolution(iPoint)); /*--- Pass supplementary information to CNumerics ---*/ - conv_numerics->SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint)); - conv_numerics->SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint)); - conv_numerics->SetdTvedU(nodes->GetdTvedU(iPoint), - nodes->GetdTvedU(iPoint)); - conv_numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint)); - conv_numerics->SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(iPoint)); + SetNumericsThermoState(conv_numerics, nodes, iPoint, nodes, iPoint); conv_numerics->SetGamma(nodes->GetGamma(iPoint), nodes->GetGamma(iPoint)); const bool dynamic_grid = config->GetGrid_Movement(); diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index e29934e8803..a229d707b64 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -117,18 +117,16 @@ unsigned long CNEMONSSolver::SetPrimitive_Variables(CSolver **solver_container,C unsigned long nonPhysicalPoints = 0; const TURB_MODEL turb_model = config->GetKind_Turb_Model(); - //const bool tkeNeeded = (turb_model == TURB_MODEL::SST); SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint ++) { /*--- Retrieve the value of the kinetic energy (if needed). ---*/ - su2double eddy_visc = 0.0; //su2double turb_ke = 0.0; + su2double eddy_visc = 0.0; if (turb_model != TURB_MODEL::NONE && solver_container[TURB_SOL] != nullptr) { eddy_visc = solver_container[TURB_SOL]->GetNodes()->GetmuT(iPoint); - //if (tkeNeeded) turb_ke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); nodes->SetEddyViscosity(iPoint, eddy_visc); } @@ -201,31 +199,10 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, nodes->GetGradient_Primitive(jPoint) ); /*--- Pass supplementary information to CNumerics ---*/ - numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); - numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); - numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); - numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(jPoint)); - numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); + SetNumericsThermoState(numerics, nodes, iPoint, nodes, jPoint); - /*--- Species diffusion coefficients ---*/ - numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), - nodes->GetDiffusionCoeff(jPoint) ); - - /*--- Laminar viscosity ---*/ - numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), - nodes->GetLaminarViscosity(jPoint) ); - - /*--- Eddy viscosity ---*/ - numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), - nodes->GetEddyViscosity(jPoint) ); - - /*--- Thermal conductivity ---*/ - numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), - nodes->GetThermalConductivity(jPoint)); - - /*--- Vib-el. thermal conductivity ---*/ - numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), - nodes->GetThermalConductivity_ve(jPoint) ); + /*--- Transport properties ---*/ + SetNumericsTransportState(numerics, iPoint, jPoint); /*--- Compute and update residual ---*/ auto residual = numerics->ComputeResidual(config); @@ -244,6 +221,39 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, } //iEdge } +void CNEMONSSolver::SetWallVelocity_Strong(unsigned long iPoint, const su2double *velocity) { + + /*--- Store the wall velocity in the old solution, it is not changed by the update. ---*/ + nodes->SetVelocity_Old(iPoint, velocity); + + /*--- Remove any contribution to the momentum residual and truncation error at this node. ---*/ + for (auto iDim = 0ul; iDim < nDim; iDim++) + LinSysRes(iPoint, nSpecies+iDim) = 0.0; + nodes->SetVel_ResTruncError_Zero(iPoint); +} + +void CNEMONSSolver::SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double Twall, + su2double Tvewall, su2double pcontrol, su2double area, + su2double dist, su2double *residual) const { + + /*--- Temperatures at the wall point (i) and at its closest normal neighbor (j) ---*/ + const su2double Ti = nodes->GetTemperature(iPoint); + const su2double Tj = nodes->GetTemperature(jPoint); + const su2double Tvei = nodes->GetTemperature_ve(iPoint); + const su2double Tvej = nodes->GetTemperature_ve(jPoint); + + /*--- Thermal conductivities at the wall point ---*/ + const su2double ktr = nodes->GetThermalConductivity(iPoint); + const su2double kve = nodes->GetThermalConductivity_ve(iPoint); + + /*--- Conductive heat flux from the normal neighbor, plus a term that is proportional + to the difference between the target wall temperature and the current one. The + second term drives the solution at the wall towards the target temperature. ---*/ + residual[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + + (ktr*(Twall-Ti) + kve*(Tvewall-Tvei))*pcontrol)*area/dist; + residual[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Tvewall-Tvei)*pcontrol)*area/dist; +} + void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -306,12 +316,8 @@ void CNEMONSSolver::BC_HeatFluxNonCatalytic_Wall(CGeometry *geometry, /*--- Impose the value of the velocity as a strong boundary condition (Dirichlet). Fix the velocity and remove any contribution to the residual at this node. ---*/ - su2double zero[MAXNDIM] = {0.0}; - nodes->SetVelocity_Old(iPoint, zero); - - for (auto iDim = 0ul; iDim < nDim; iDim++) - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); + const su2double zero[MAXNDIM] = {0.0}; + SetWallVelocity_Strong(iPoint, zero); /*--- Apply viscous residual to the linear system ---*/ LinSysRes.SubtractBlock(iPoint, Res_Visc); @@ -353,164 +359,8 @@ void CNEMONSSolver::BC_HeatFluxCatalytic_Wall(CGeometry *geometry, unsigned short val_marker) { SU2_ZONE_SCOPED + /*--- The catalytic contribution to a heat-flux wall is not implemented. ---*/ SU2_MPI::Error("BC_HEATFLUX with catalytic wall: Not operational in NEMO.", CURRENT_FUNCTION); - //TODO: SCALE WITH EDDY VISC - /*--- Local variables ---*/ - unsigned long iPoint; - su2double rho, Ys; - su2double *Normal, Area; - su2double *Ds, *dYdn, SdYdn; - su2double **GradY; - - /*--- Assign booleans ---*/ - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool catalytic = false; - - /*--- Set "Proportional control" coefficient ---*/ - su2double pcontrol = 0.6; - - /*--- Identify the boundary by string name ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Get the specified wall heat flux from config ---*/ - su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag) / config->GetHeat_Flux_Ref(); - if (config->GetIntegrated_HeatFlux()) { - Wall_HeatFlux /= geometry->GetSurfaceArea(config, val_marker); - } - - /*--- Get the locations of the primitive variables ---*/ - const unsigned short T_INDEX = nodes->GetTIndex(); - const unsigned short TVE_INDEX = nodes->GetTveIndex(); - const unsigned short RHOS_INDEX = nodes->GetRhosIndex(); - const unsigned short RHO_INDEX = nodes->GetRhoIndex(); - - /*--- Allocate arrays ---*/ - dYdn = new su2double[nSpecies]; - GradY = new su2double*[nSpecies]; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - GradY[iSpecies] = new su2double[nDim]; - - /*--- Loop over all of the vertices on this boundary marker ---*/ - for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Compute dual-grid area and boundary normal ---*/ - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Initialize the convective & viscous residuals to zero ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) { - Res_Visc[iVar] = 0.0; - Res_Sour[iVar] = 0.0; - } - - /*--- Assign wall velocity to "Vector" array ---*/ - for (auto iDim = 0ul; iDim < nDim; iDim++) Vector[iDim] = 0.0; - - /*--- Set the residual, truncation error, and velocity value ---*/ - nodes->SetVelocity_Old(iPoint,Vector); - for (auto iDim = 0ul; iDim < nDim; iDim++) { - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVal_ResTruncError_Zero(iPoint,nSpecies+iDim); - } - - /*--- Get temperature gradient information ---*/ - const auto V = nodes->GetPrimitive(iPoint); - const auto GradV = nodes->GetGradient_Primitive(iPoint); - su2double dTdn = 0.0; - su2double dTvedn = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) { - dTdn += GradV[T_INDEX][iDim]*Normal[iDim]; - dTvedn += GradV[TVE_INDEX][iDim]*Normal[iDim]; - } - - if (catalytic) { - SU2_MPI::Error("NEED TO IMPLEMENT CATALYTIC BOUNDARIES IN HEATFLUX!!!",CURRENT_FUNCTION); - } - else { - - /*--- Rename for convenience ---*/ - rho = V[RHO_INDEX]; - Ds = nodes->GetDiffusionCoeff(iPoint); - - /*--- Calculate normal derivative of mass fraction ---*/ - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) { - Ys = V[RHOS_INDEX+iSpecies]/rho; - dYdn[iSpecies] = 0.0; - for (auto iDim = 0ul; iDim < nDim; iDim++) - dYdn[iSpecies] += 1.0/rho * (GradV[RHOS_INDEX+iSpecies][iDim] - - Ys*GradV[RHO_INDEX][iDim])*Normal[iDim]; - } - - /*--- Calculate supplementary quantities ---*/ - SdYdn = 0.0; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - SdYdn += rho*Ds[iSpecies]*dYdn[iSpecies]; - - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) { - Ys = V[RHOS_INDEX+iSpecies]/rho; - //eves = nodes->CalcEve(config, V[TVE_INDEX], iSpecies); - //hs = nodes->CalcHs(config, V[T_INDEX], eves, iSpecies); - // Res_Visc[iSpecies] = -rho*Ds[iSpecies]*dYdn[iSpecies] + Ys*SdYdn; - // Res_Visc[nSpecies+nDim] += Res_Visc[iSpecies]*hs; - // Res_Visc[nSpecies+nDim+1] += Res_Visc[iSpecies]*eves; - } - } - - /*--- Get node thermal conductivity ---*/ - su2double ktr = nodes->GetThermalConductivity(iPoint); - su2double kve = nodes->GetThermalConductivity_ve(iPoint); - - /*--- Set the residual on the boundary with the specified heat flux ---*/ - // Note: Contributions from qtr and qve are used for proportional control - // to drive the solution toward the specified heatflux more quickly. - Res_Visc[nSpecies+nDim] += pcontrol*(ktr*dTdn+kve*dTvedn) + - Wall_HeatFlux*Area; - Res_Visc[nSpecies+nDim+1] += pcontrol*(kve*dTvedn) + - Wall_HeatFlux*Area; - - /*--- Viscous contribution to the residual at the wall ---*/ - LinSysRes.SubtractBlock(iPoint, Res_Visc); - - // /*--- Apply the non-catalytic wall boundary ---*/ - // // Note: We are re-calculating the chemistry residual and adding it to - // // the linear system to eliminate the contribution from the solution - // // (convention is to subtract sources) - // sour_numerics->SetConservative(nodes->GetSolution(), - // nodes->GetSolution() ); - // sour_numerics->SetPrimitive (nodes->GetPrimVar() , - // nodes->GetPrimVar() ); - // sour_numerics->SetdPdU (nodes->GetdPdU() , - // nodes->GetdPdU() ); - // sour_numerics->SetdTdU (nodes->GetdTdU() , - // nodes->GetdTdU() ); - // sour_numerics->SetdTvedU (nodes->GetdTvedU() , - // nodes->GetdTvedU() ); - // sour_numerics->SetVolume (geometry->nodes->GetVolume()); - // sour_numerics->ComputeChemistry(Res_Sour, Jacobian_i, config); - // LinSysRes.AddBlock(iPoint, Res_Sour); - // if (implicit) - // Jacobian.AddBlock2Diag(iPoint, iPoint, Jacobian_i); - - /*--- Only change velocity-rows of the Jacobian (includes 1 in the diagonal)/ - Note that we need to add a contribution for moving walls to the Jacobian. ---*/ - if (implicit) { - /*--- Enforce the no-slip boundary condition in a strong way ---*/ - for (auto iVar = nSpecies; iVar < nSpecies+nDim; iVar++) { - Jacobian.DeleteValsRowi(iPoint, iVar); - } - } - - } - } - - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - delete [] GradY[iSpecies]; - delete [] GradY; - delete [] dYdn; } void CNEMONSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, @@ -538,7 +388,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool ionization = config->GetIonization(); - su2double UnitNormal[MAXNDIM] = {0.0}; if (ionization) { SU2_MPI::Error("NEED TO TAKE A CLOSER LOOK AT THE JACOBIAN W/ IONIZATION",CURRENT_FUNCTION); @@ -553,13 +402,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Retrieve the specified wall temperature ---*/ const su2double Twall = config->GetIsothermal_Temperature(Marker_Tag); - su2double **Jacobian_i = nullptr; - if (implicit) { - Jacobian_i = new su2double* [nVar]; - for (auto iVar = 0ul; iVar < nVar; iVar++) - Jacobian_i[iVar] = new su2double [nVar] (); - } - /*--- Loop over boundary points to calculate energy flux ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -571,6 +413,7 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Compute dual-grid area and boundary normal ---*/ const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); const su2double Area = GeometryToolbox::Norm(nDim, Normal); + su2double UnitNormal[MAXNDIM] = {0.0}; for (auto iDim = 0ul; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim]/Area; @@ -585,39 +428,23 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, /*--- Store the corrected velocity at the wall which will be zero (v = 0), unless there is grid motion (v = u_wall)---*/ - su2double zero[MAXNDIM] = {0.0}; - nodes->SetVelocity_Old(iPoint, zero); + const su2double zero[MAXNDIM] = {0.0}; + SetWallVelocity_Strong(iPoint, zero); /*--- Initialize viscous residual to zero ---*/ for (auto iVar = 0ul; iVar < nVar; iVar ++) {Res_Visc[iVar] = 0.0;} - for (auto iDim = 0ul; iDim < nDim; iDim++) - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); - - /*--- Calculate the gradient of temperature ---*/ - const su2double Ti = nodes->GetTemperature(iPoint); - const su2double Tj = nodes->GetTemperature(Point_Normal); - const su2double Tvei = nodes->GetTemperature_ve(iPoint); - const su2double Tvej = nodes->GetTemperature_ve(Point_Normal); - /*--- Rename variables for convenience ---*/ const su2double ktr = nodes->GetThermalConductivity(iPoint); const su2double kve = nodes->GetThermalConductivity_ve(iPoint); - /*--- Apply to the linear system ---*/ - Res_Visc[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + - (ktr*(Twall-Ti) + kve*(Twall-Tvei))*C)*Area/dist_ij; - Res_Visc[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Twall-Tvei) *C)*Area/dist_ij; + /*--- Weak enforcement of the wall temperature on both energy equations ---*/ + SetWallEnergyFlux_Weak(iPoint, Point_Normal, Twall, Twall, C, Area, dist_ij, Res_Visc); /*--- Calculate Jacobian for implicit time stepping ---*/ + su2double Jacobian_i[MAXNVAR][MAXNVAR] = {{0.0}}; if (implicit) { - /*--- Initialize Jacobian to zero ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) - for (auto jVar = 0ul; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - /*--- Add contributions to the Jacobian from the weak enforcement of the energy equations. ---*/ const auto dTdU = nodes->GetdTdU(iPoint); const auto dTvedU = nodes->GetdTvedU(iPoint); @@ -645,12 +472,6 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry, } } END_SU2_OMP_FOR - - if (Jacobian_i) - for (auto iVar = 0ul; iVar < nVar; iVar++) - delete [] Jacobian_i[iVar]; - delete [] Jacobian_i; - } void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, @@ -665,8 +486,8 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, BC_IsothermalNonCatalytic_Wall(geometry, solver_container, conv_numerics, sour_numerics, config, val_marker); - /*--- Local variables ---*/ - su2double **GradY, **dVdU; + /*--- Transformation matrix from conservative to primitive variables ---*/ + su2double dVdU[MAXNVAR][MAXNVAR] = {{0.0}}; /*--- Assign booleans ---*/ const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -683,14 +504,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, const unsigned short T_INDEX = nodes->GetTIndex(); const unsigned short TVE_INDEX = nodes->GetTveIndex(); - /*--- Allocate arrays ---*/ - GradY = new su2double*[nSpecies]; - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - GradY[iSpecies] = new su2double[nDim]; - dVdU = new su2double*[nVar]; - for (auto iVar = 0ul; iVar < nVar; iVar++) - dVdU[iVar] = new su2double[nVar]; - /*--- Loop over all of the vertices on this boundary marker ---*/ for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -840,13 +653,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry, LinSysRes.SubtractBlock(iPoint, Res_Visc); } } - - for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) - delete [] GradY[iSpecies]; - delete [] GradY; - for (auto iVar = 0ul; iVar < nVar; iVar++) - delete [] dVdU[iVar]; - delete [] dVdU; } void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, @@ -912,11 +718,8 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Calculate Pressure ---*/ const su2double Pi = nodes->GetPressure(iPoint); - /*--- Calculate the gradient of temperature ---*/ + /*--- Temperature at the wall point ---*/ const su2double Ti = nodes->GetTemperature(iPoint); - const su2double Tj = nodes->GetTemperature(jPoint); - const su2double Tvei = nodes->GetTemperature_ve(iPoint); - const su2double Tvej = nodes->GetTemperature_ve(jPoint); /*--- Rename variables for convenience ---*/ const su2double ktr = nodes->GetThermalConductivity(iPoint); @@ -984,18 +787,13 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, Vector[iDim] = (1.0-alpha_V)*nodes->GetVelocity(iPoint,iDim)+(alpha_V)*Vector[iDim]; } - nodes->SetVelocity_Old(iPoint,Vector); + /*--- Impose the slip velocity as a strong boundary condition ---*/ + SetWallVelocity_Strong(iPoint, Vector); - for (auto iDim = 0ul; iDim < nDim; iDim++) { - LinSysRes(iPoint, nSpecies+iDim) = 0.0; - nodes->SetVal_ResTruncError_Zero(iPoint,nSpecies+iDim); - } + /*--- Weak enforcement of the slip temperatures on both energy equations ---*/ + SetWallEnergyFlux_Weak(iPoint, jPoint, Tslip, Tslip_ve, C, Area, dij, Res_Visc); /*--- Apply to the linear system ---*/ - Res_Visc[nSpecies+nDim] = ((ktr*(Ti-Tj) + kve*(Tvei-Tvej)) + - (ktr*(Tslip-Ti) + kve*(Tslip_ve-Tvei))*C)*Area/dij; - Res_Visc[nSpecies+nDim+1] = (kve*(Tvei-Tvej) + kve*(Tslip_ve-Tvei)*C)*Area/dij; - LinSysRes.SubtractBlock(iPoint, Res_Visc); } END_SU2_OMP_FOR