Skip to content

SST turbulence model with sustaining terms #765

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vdweide merged 9 commits into from
Aug 25, 2019
Merged

SST turbulence model with sustaining terms #765

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
b174af6
Added the SST turbulence model with sustaining terms
vdweide Aug 20, 2019
b19fadc
Fixed a bug and added a couple of regression test cases for SST_SUST
vdweide Aug 20, 2019
1ca22ae
A small change in the comments of the serial regression test
vdweide Aug 20, 2019
5512136
Merge branch 'develop' into feature_SST_SUST
vdweide Aug 20, 2019
efb7e99
Changed travis such that only I get an error message when something g…
vdweide Aug 22, 2019
0fa9f53
Reverted travis
vdweide Aug 22, 2019
f779211
Merged with develop and resolved the conflicts
vdweide Aug 23, 2019
7bbb549
Updated the config files for the SST_SUST test cases
vdweide Aug 25, 2019
9ba04e8
Merge branch 'develop' into feature_SST_SUST
vdweide Aug 25, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions Common/include/option_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -829,28 +829,30 @@ static const map<string, ENUM_LIMITER> Limiter_Map = CCreateMap<string, ENUM_LIM
*/
enum ENUM_TURB_MODEL {
NO_TURB_MODEL = 0, /*!< \brief No turbulence model. */
SA = 1, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
SA_NEG = 2, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
SST = 3, /*!< \brief Kind of Turbulence model (Menter SST). */
SA_E = 4, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards). */
SA_COMP = 5, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Compressibility Correction). */
SA_E_COMP = 6, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards with Compressibility Correction). */
SA = 1, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
SA_NEG = 2, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
SA_E = 3, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards). */
SA_COMP = 4, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Compressibility Correction). */
SA_E_COMP = 5, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards with Compressibility Correction). */
SST = 6, /*!< \brief Kind of Turbulence model (Menter SST). */
SST_SUST = 7 /*!< \brief Kind of Turbulence model (Menter SST with sustaining terms for free-stream preservation). */
};
static const map<string, ENUM_TURB_MODEL> Turb_Model_Map = CCreateMap<string, ENUM_TURB_MODEL>
("NONE", NO_TURB_MODEL)
("SA", SA)
("SA_NEG", SA_NEG)
("SST", SST)
("SA_E", SA_E)
("SA_COMP", SA_COMP)
("SA_E_COMP", SA_E_COMP);
("SA_E_COMP", SA_E_COMP)
("SST", SST)
("SST_SUST", SST_SUST);

/*!
* \brief types of transition models
*/
enum ENUM_TRANS_MODEL {
NO_TRANS_MODEL = 0, /*!< \brief No transition model. */
LM = 1, /*!< \brief Kind of transition model (LM for Spalart-Allmaras). */
LM = 1, /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */
BC = 2 /*!< \brief Kind of transition model (BAS-CAKMAKCIOGLU (BC) for Spalart-Allmaras). */
};
static const map<string, ENUM_TRANS_MODEL> Trans_Model_Map = CCreateMap<string, ENUM_TRANS_MODEL>
Expand Down
17 changes: 9 additions & 8 deletions Common/src/config_structure.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4082,7 +4082,7 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_

/* --- Throw error if UQ used for any turbulence model other that SST --- */

if (Kind_Solver == RANS && Kind_Turb_Model != SST && using_uq){
if (Kind_Solver == RANS && Kind_Turb_Model != SST && Kind_Turb_Model != SST_SUST && using_uq){
SU2_MPI::Error("UQ capabilities only implemented for NAVIER_STOKES solver SST turbulence model", CURRENT_FUNCTION);
}

Expand Down Expand Up @@ -5047,21 +5047,22 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) {
if (Kind_Regime == INCOMPRESSIBLE) cout << "Incompressible RANS equations." << endl;
cout << "Turbulence model: ";
switch (Kind_Turb_Model) {
case SA: cout << "Spalart Allmaras" << endl; break;
case SA_NEG: cout << "Negative Spalart Allmaras" << endl; break;
case SST: cout << "Menter's SST" << endl; break;
case SA_E: cout << "Edwards Spalart Allmaras" << endl; break;
case SA: cout << "Spalart Allmaras" << endl; break;
case SA_NEG: cout << "Negative Spalart Allmaras" << endl; break;
case SA_E: cout << "Edwards Spalart Allmaras" << endl; break;
case SA_COMP: cout << "Compressibility Correction Spalart Allmaras" << endl; break;
case SA_E_COMP: cout << "Compressibility Correction Edwards Spalart Allmaras" << endl; break;
case SA_E_COMP: cout << "Compressibility Correction Edwards Spalart Allmaras" << endl; break;
case SST: cout << "Menter's SST" << endl; break;
case SST_SUST: cout << "Menter's SST with sustaining terms" << endl; break;
}
if (QCR) cout << "Using Quadratic Constitutive Relation, 2000 version (QCR2000)" << endl;
cout << "Hybrid RANS/LES: ";
switch (Kind_HybridRANSLES){
case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
case SA_DES: cout << "Detached Eddy Simulation (DES97) " << endl; break;
case SA_DES: cout << "Detached Eddy Simulation (DES97) " << endl; break;
case SA_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
}
if (using_uq){
cout << "Perturbing Reynold's Stress Matrix towards "<< eig_val_comp << " component turbulence"<< endl;
Expand Down
6 changes: 4 additions & 2 deletions Common/src/geometry_structure.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16022,8 +16022,10 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config) {

ifstream restart_file;
string filename = config->GetSolution_AdjFileName();
bool sst = config->GetKind_Turb_Model() == SST;
bool sa = (config->GetKind_Turb_Model() == SA) || (config->GetKind_Turb_Model() == SA_NEG);
bool sst = (config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST);
bool sa = (config->GetKind_Turb_Model() == SA) || (config->GetKind_Turb_Model() == SA_NEG) ||
(config->GetKind_Turb_Model() == SA_E) || (config->GetKind_Turb_Model() == SA_COMP) ||
(config->GetKind_Turb_Model() == SA_E_COMP);
bool grid_movement = config->GetGrid_Movement();
bool frozen_visc = config->GetFrozen_Visc_Disc();
unsigned short Kind_Solver = config->GetKind_Solver();
Expand Down
8 changes: 6 additions & 2 deletions SU2_CFD/include/numerics_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5074,8 +5074,11 @@ class CSourcePieceWise_TurbSST : public CNumerics {
a1;

su2double CDkw_i, CDkw_j;


su2double kAmb, omegaAmb;

bool incompressible;
bool sustaining_terms;

public:

Expand All @@ -5085,7 +5088,8 @@ class CSourcePieceWise_TurbSST : public CNumerics {
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double* constants, CConfig *config);
CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double* constants,
su2double val_kine_Inf, su2double val_omega_Inf, CConfig *config);

/*!
* \brief Destructor of the class.
Expand Down
61 changes: 38 additions & 23 deletions SU2_CFD/src/drivers/CDriver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1099,12 +1099,13 @@ void CDriver::Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolve

if (turbulent || fem_turbulent)
switch (config->GetKind_Turb_Model()) {
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E_COMP: e_comp_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SST_SUST: menter_sst = true; break;
default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
}

Expand Down Expand Up @@ -1574,12 +1575,14 @@ void CDriver::Solver_Postprocessing(CSolver ****solver, CGeometry **geometry,

if (turbulent)
switch (config->GetKind_Turb_Model()) {
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E_COMP: e_comp_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SST_SUST: menter_sst = true; break;
default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
}

/*--- Definition of the Class for the solution: solver_container[DOMAIN][MESH_LEVEL][EQUATION]. Note that euler, ns
Expand Down Expand Up @@ -1812,6 +1815,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri
numerics = new CNumerics***[config->GetnMGLevels()+1];

su2double *constants = NULL;
su2double kine_Inf = 0.0, omega_Inf = 0.0;

bool
euler, adj_euler,
Expand Down Expand Up @@ -1865,14 +1869,24 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri

if (turbulent || fem_turbulent)
switch (config->GetKind_Turb_Model()) {
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E_COMP: e_comp_spalart_allmaras = true; break;
case SST: menter_sst = true; constants = solver[MESH_0][TURB_SOL]->GetConstants(); break;
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E_COMP: e_comp_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SST_SUST: menter_sst = true; break;
default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
}

/*--- If the Menter SST model is used, store the constants of the model and determine the
free stream values of the turbulent kinetic energy and dissipation rate. ---*/

if (menter_sst) {
constants = solver[MESH_0][TURB_SOL]->GetConstants();
kine_Inf = solver[MESH_0][TURB_SOL]->GetTke_Inf();
omega_Inf = solver[MESH_0][TURB_SOL]->GetOmega_Inf();
}

/*--- Number of variables for the template ---*/

Expand Down Expand Up @@ -2269,7 +2283,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri
else if (comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config);
else if (e_comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config);
else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config);
else if (menter_sst) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, config);
else if (menter_sst) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config);
numerics[iMGlevel][TURB_SOL][SOURCE_SECOND_TERM] = new CSourceNothing(nDim, nVar_Turb, config);
}

Expand Down Expand Up @@ -2680,13 +2694,14 @@ void CDriver::Numerics_Postprocessing(CNumerics *****numerics,

if (turbulent || fem_turbulent)
switch (config->GetKind_Turb_Model()) {
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SST: menter_sst = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA: spalart_allmaras = true; break;
case SA_NEG: neg_spalart_allmaras = true; break;
case SA_COMP: comp_spalart_allmaras = true; break;
case SA_E: e_spalart_allmaras = true; break;
case SA_E_COMP: e_comp_spalart_allmaras = true; break;

case SST: menter_sst = true; break;
case SST_SUST: menter_sst = true; break;
default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
}

/*--- Solver definition for the template problem ---*/
Expand Down
31 changes: 27 additions & 4 deletions SU2_CFD/src/numerics_direct_turbulent.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1166,9 +1166,11 @@ void CAvgGrad_TurbSST::FinishResidualCalc(su2double *val_residual, su2double **J
}

CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double *constants,
CConfig *config) : CNumerics(val_nDim, val_nVar, config) {
su2double val_kine_Inf, su2double val_omega_Inf, CConfig *config)
: CNumerics(val_nDim, val_nVar, config) {

incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
sustaining_terms = (config->GetKind_Turb_Model() == SST_SUST);

/*--- Closure constants ---*/
beta_star = constants[6];
Expand All @@ -1179,6 +1181,10 @@ CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, unsi
alfa_1 = constants[8];
alfa_2 = constants[9];
a1 = constants[7];

/*--- Set the ambient values of k and omega to the free stream values. ---*/
kAmb = val_kine_Inf;
omegaAmb = val_omega_Inf;
}

CSourcePieceWise_TurbSST::~CSourcePieceWise_TurbSST(void) { }
Expand Down Expand Up @@ -1256,10 +1262,27 @@ void CSourcePieceWise_TurbSST::ComputeResidual(su2double *val_residual, su2doubl
else {
pw = StrainMag_i*StrainMag_i - 2.0/3.0*zeta*diverg;
}
pw = max(pw,0.0);
pw = alfa_blended*Density_i*max(pw,0.0);

/*--- Sustaining terms, if desired. Note that if the production terms are
larger equal than the sustaining terms, the original formulation is
obtained again. This is in contrast to the version in literature
where the sustaining terms are simply added. This latter approach could
lead to problems for very big values of the free-stream turbulence
intensity. ---*/

if ( sustaining_terms ) {
const su2double sust_k = beta_star*Density_i*kAmb*omegaAmb;

@pcarruscag pcarruscag Aug 22, 2019

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes are so slight I missed them while scrolling to get to the meat of the feature.

Why not simply add a config option for USE_SST_SUSTAINING_TERMS instead of through a new turb solver type? (which forces you to change lots of booleans in the code)
Or maybe use the turb solver option but in the config post processing set a different boolean to indicate sustaining_terms and reset KindTurbSolver to SST.

const su2double sust_w = beta_blended*Density_i*omegaAmb*omegaAmb;

pk = max(pk, sust_k);
pw = max(pw, sust_w);
}

/*--- Add the production terms to the residuals. ---*/

val_residual[0] += pk*Volume;
val_residual[1] += alfa_blended*Density_i*pw*Volume;
val_residual[1] += pw*Volume;

/*--- Dissipation ---*/

Expand Down
Loading