-
Notifications
You must be signed in to change notification settings - Fork 842
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.
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
[WIP] Implementation of DDES and SAS formulations for SST #2150
base: develop
Are you sure you want to change the base?
Conversation
…ons based on k and w - Added computation of k and w starting from SA variables - Added first implementation of SAS (Scale Adaptive Simulation) model
- Changed previous SAS model to SAS_SIMPLE
- Fixed DataType for Periodic Comms
…lified IDDES (SIDDES)
switch(kind_hybridRANSLES){ | ||
case SST_DDES: { | ||
|
||
const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); | ||
const su2double C_d1 = 20.0; | ||
const su2double C_d2 = 3.0; | ||
|
||
const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2)); | ||
|
||
const su2double l_LES = C_DES * h_max; | ||
|
||
DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES); | ||
|
||
nodes->SetDebug_Quantities(config, iPoint, f_d, l_RANS, l_LES, r_d); | ||
|
||
break; | ||
} | ||
case SST_IDDES: { | ||
|
||
// Constants | ||
const su2double C_w = 0.15; | ||
const su2double C_dt1 = 20.0; | ||
const su2double C_dt2 = 3.0; | ||
const su2double C_l = 5.0; | ||
const su2double C_t = 1.87; | ||
|
||
const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; | ||
const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); | ||
const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); | ||
const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); | ||
const su2double ftilda_d = max(1.0 - f_dt, f_b); | ||
|
||
const su2double r_dl = lamVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); | ||
const su2double f_l = tanh(pow(C_l*C_l*r_dl, 10.0)); | ||
const su2double f_t = tanh(pow(C_t*C_t*r_dt, 3.0)); | ||
const su2double f_e2 = 1.0 - max(f_t, f_l); | ||
const su2double f_e1 = alpha >= 0.0 ? 2.0 * exp(-11.09*alpha*alpha) : 2.0 * exp(-9.0*alpha*alpha); | ||
const su2double f_e = f_e2 * max((f_e1 - 1.0), 0.0); | ||
|
||
|
||
const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); | ||
const su2double l_LES = C_DES * Delta; | ||
|
||
nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dl, r_dt); | ||
|
||
DES_lengthScale = ftilda_d *(1.0+f_e)*l_RANS + (1.0 - ftilda_d) * l_LES; | ||
|
||
break; | ||
} | ||
case SST_SIDDES: { | ||
|
||
// Constants | ||
const su2double C_w = 0.15; | ||
const su2double C_dt1 = 20.0; | ||
const su2double C_dt2 = 3.0; | ||
const su2double C_l = 5.0; | ||
const su2double C_t = 1.87; | ||
|
||
const su2double alpha = 0.25 - sqrt(wallDist2) / h_max; | ||
const su2double f_b = min(2.0 * exp(-9.0 * alpha*alpha), 1.0); | ||
const su2double r_dt = eddyVisc / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10); | ||
const su2double f_dt = 1 - tanh(pow(C_dt1 * r_dt, C_dt2)); | ||
const su2double ftilda_d = max(1.0 - f_dt, f_b); | ||
|
||
const su2double Delta = min(C_w * max(sqrt(wallDist2), h_max), h_max); | ||
const su2double l_LES = C_DES * Delta; | ||
|
||
nodes->SetDebug_Quantities(config, iPoint, ftilda_d, l_RANS, l_LES, r_dt); | ||
|
||
DES_lengthScale = ftilda_d*l_RANS + (1.0 - ftilda_d) * l_LES; | ||
|
||
break; | ||
} | ||
} |
Check notice
Code scanning / CodeQL
Long switch case Note
SST_IDDES (31 lines)
@@ -1060,3 +1075,112 @@ | |||
} | |||
|
|||
} | |||
|
|||
|
|||
void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CConfig *config){ |
Check warning
Code scanning / CodeQL
Poorly documented large function Warning
Common/include/option_structure.hpp
Outdated
SAS_SIMPLE, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ | ||
SAS_COMPLICATED, /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are these the names in the literature?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, these are just placeholders. There are no names in literature, that is why I named them according to the complexity of the model
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe something after the authors since the models come from different papers.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sure, why not
Common/include/option_structure.hpp
Outdated
SST_OPTIONS sasModel = SST_OPTIONS::SAS_SIMPLE; /*!< \brief Enum SST base model. */ | ||
bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ | ||
bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ | ||
bool sas = false; /*!< \brief Bool for using Scale Adaptive Simulations. */ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Probably the type is enough and "none" indicates no sas.
SST_OPTIONS sasModel = SST_OPTIONS::SAS_SIMPLE; /*!< \brief Enum SST base model. */ | |
bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ | |
bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ | |
bool sas = false; /*!< \brief Bool for using Scale Adaptive Simulations. */ | |
SST_OPTIONS sasModel = SST_OPTIONS::NONE; /*!< \brief Enum SST base model. */ | |
bool sust = false; /*!< \brief Bool for SST model with sustaining terms. */ | |
bool uq = false; /*!< \brief Bool for using uncertainty quantification. */ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sure, I'll modify it.
Common/src/CConfig.cpp
Outdated
case SST_DDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; | ||
case SST_IDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; | ||
case SST_SIDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The description is slightly different right?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I planned to modify the details on the final implementation. I'll do it now
@@ -190,6 +190,10 @@ class CNumerics { | |||
|
|||
bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ | |||
|
|||
su2double lengthScale_i, lengthScale_j; | |||
su2double FTrans; /*!< \brief SAS function */ | |||
su2double VelLapl_X, VelLapl_Y, VelLapl_Z; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Use a matrix type for the laplacian in CVariable and here you can have a pointer instead.
Put these variables close to other turbulence variables please.
su2double VelLaplMag = VelLapl_X*VelLapl_X + VelLapl_Y*VelLapl_Y; | ||
if (nDim == 3) VelLaplMag += VelLapl_Z*VelLapl_Z; | ||
const su2double L_vK_1 = KolmConst * StrainMag_i / sqrt(VelLaplMag); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is the velocity laplacian only used for the source term? Or will you have changes to the convection or diffusion fluxes?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is just for the source term
/*! | ||
* \brief Get the value of the turbulence kinetic energy. | ||
* \return the value of the turbulence kinetic energy. | ||
*/ | ||
inline su2double GetSSTVariables_k(unsigned long iPoint) const { return k(iPoint); } | ||
|
||
/*! | ||
* \brief Get the value of the turbulence frequency Omega. | ||
* \return the value of the turbulence frequency Omega. | ||
*/ | ||
inline su2double GetSSTVariables_omega(unsigned long iPoint) const { return Omega(iPoint); } | ||
|
||
/*! | ||
* \brief Set the value of the SST variables computed with SA solution. | ||
* \param[in] val_k | ||
* \param[in] val_Omega | ||
*/ | ||
void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) { | ||
k(iPoint) = val_k; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are these changes related to the main theme? or were just debugging? If they can be separated to another PR let's do that please.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These are for computing the grid size for scale resolving simulations from an SA solution. Probably it has to be revised a little bit, since I found out that the formulation for the TKE is different if the QCR mod is used, thus I think I'll remove them for now.
externals/codi
Outdated
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You need to undo this change to the codi version.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sure, I'll do it
if ( config->GetKind_HybridRANSLES() == SST_DDES){ | ||
SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); | ||
SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); | ||
SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); | ||
SetVolumeOutputValue("R_D", iPoint, Node_Turb->Get_r_d(iPoint)); | ||
} else if ( config->GetKind_HybridRANSLES() == SST_IDDES){ | ||
SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); | ||
SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); | ||
SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); | ||
SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint)); | ||
SetVolumeOutputValue("R_DL", iPoint, Node_Turb->Get_r_dl(iPoint)); | ||
} else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){ | ||
SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint)); | ||
SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint)); | ||
SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint)); | ||
SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint)); | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Cleanup the repetition here, most statement are shared across the 3 cases
} | ||
const su2double mut = Node_Flow->GetEddyViscosity(iPoint); | ||
const su2double mu = Node_Flow->GetLaminarViscosity(iPoint); | ||
const su2double LESIQ = 1.0/(1.0+0.05*pow((mut+mu)/mu, 0.53)); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are the constant here also used somewhere else? It would be good to avoid repeating hard-coded constants
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
They are used just here
inline virtual su2double GetVelLapl_X(unsigned long iPoint) const { return 0.0; } | ||
/*! | ||
* \brief Get the value of the value of FTrans. | ||
*/ | ||
inline virtual su2double GetVelLapl_Y(unsigned long iPoint) const { return 0.0; } | ||
/*! | ||
* \brief Get the value of the value of FTrans. | ||
*/ | ||
inline virtual su2double GetVelLapl_Z(unsigned long iPoint) const { return 0.0; } | ||
|
||
/*! | ||
* \brief Set the value of the value of FTrans. | ||
*/ | ||
inline virtual void AddVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} | ||
/*! | ||
* \brief Set the value of the value of FTrans. | ||
*/ | ||
inline virtual void AddVelLapl_Z(unsigned long iPoint, su2double val_VelLapl_Z) {} | ||
/*! | ||
* \brief Set the value of the value of FTrans. | ||
*/ | ||
inline virtual void SetVelLapl(unsigned long iPoint, su2double val_VelLapl_X, su2double val_VelLapl_Y) {} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Use a pointer for velocity and then you only need one of each function instead of handling xyz separately.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, I'll try using the matrix type
- Changed names of SAS models - Added descriptions of functions - Clean up of output functions - Removed SA functions or SST variables
@@ -136,6 +136,9 @@ class CNumerics { | |||
const su2double | |||
*TurbPsi_i, /*!< \brief Vector of adjoint turbulent variables at point i. */ | |||
*TurbPsi_j; /*!< \brief Vector of adjoint turbulent variables at point j. */ | |||
su2double lengthScale_i, lengthScale_j; /*!< \brief length scale for SST */ | |||
su2double FTrans; /*!< \brief SAS function */ | |||
su2double *VelLapl; /*!< \brief Laplacian of the velocity */ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
su2double *VelLapl; /*!< \brief Laplacian of the velocity */ | |
const su2double *VelLapl; /*!< \brief Laplacian of the velocity */ |
/*! | ||
* \brief Get the value of the value of FTrans. | ||
*/ | ||
inline void SetVelLapl(su2double* val_VelLapl) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline void SetVelLapl(su2double* val_VelLapl) { | |
inline void SetVelLapl(const su2double* val_VelLapl) { |
AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); | ||
break; | ||
|
||
case TURB_FAMILY::KW: | ||
AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); | ||
AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); | ||
AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "SOLUTION", "desired grid size for Scale Resolving Simulations"); | ||
if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_TRAVIS) AddVolumeOutput("FTRANS", "FTrans", "SOLUTION", "value of FTrans for SAS simulation"); | ||
if (config->GetSSTParsedOptions().sasModel == SST_OPTIONS::SAS_BABU){ | ||
AddVolumeOutput("VEL-LAPLACIAN_X", "Vel Laplacian x", "SOLUTION", "value of laplacian of x-velocity for SAS simulation"); | ||
AddVolumeOutput("VEL-LAPLACIAN_Y", "Vel Laplacian y", "SOLUTION", "value of laplacian of y-velocity for SAS simulation"); | ||
if (nDim == 3) AddVolumeOutput("VEL-LAPLACIAN_Z", "Vel Laplacian z", "SOLUTION", "value of laplacian of z-velocity for SAS simulation"); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These outputs need to go somewhere else or they will break the restarts.
Solution is exclusively for solution variables, these are derived, put them where we have similar things (e.g. Q criterion)
/*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ | ||
if (boundary_i && !boundary_j) continue; | ||
|
||
/*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Check if comments need updates please
Hi @rois1995 any progress? Would be very interesting to add this to develop, looking forward to using it actually. |
Hi @bigfooted, unfortunately I do not have enough time/computational power to perform a proper validation. I only tried something with a NACA0021 at 17 deg and 60 deg and they were promising, at least for the DDES implementation. The SAS ones are not that good. I'll upload some post-processing as soon as I have the time to do so. |
I find your work really interesting. I've been studying the internal flow field in compressors and have had good results using SU2's SA_EDDES for calculating the cantilevered stator with a tip clearance. If you need help with code verification, I'd be glad to assist. |
Hi @Linnnnnn23, every help on the validation/verification is gladly accepted! Let me know if you need anything by my side. |
|
Hi @rois1995 , I'd like to help with validation/verification if I can. Let me know if you are interested. You can contact me at [email protected]. |
Hi @rois1995, I believe there is an error or a typo in the Babu's paper in the Q_SAS source term and in turn in your implementation, which may be the reason for the poor performance in SAS results. This term changes units depending on the output of the max function, which doesn't make any sense. I believe the terms in the max function should be divided by omega^2 and k^2 instead of omega and k, respectively. Which is exactly how it is done in https://resolver.tudelft.nl/uuid:5d23e2a6-5675-450d-bf3d-1dd40d736cae I will try some of the benchmark cases in these papers when I have the time. Let me know what you think. |
hi @BerkeCan97, thank you for your correction and for helping with the verification part. Indeed, that term seems really wrong, thus I expect that your correction would improve the results. However, the link that you have provided for the corrected version is broken, can you fix it please? Regarding the test-cases, we were trying to run the NACA0021 at 60deg AoA, but with relatively poor results. Maybe we should start with the flatplate and the mixing layer proposed originally in the implementation of the DDES formulation coupled to the SA model in https://www.researchgate.net/profile/Cleber-Spode-2/publication/318143234_Hybrid_RANSLES_Calculations_in_SU2/links/59c07bfba6fdcca8e56fe3c3/Hybrid-RANS-LES-Calculations-in-SU2.pdf (or any other test-case in the paper you proposed, but I am not able to read that). |
Fixed it! I think we can start with flat plate and maybe BFS. It might take a week or so for me to get anything done though. |
Sorry but it still says "Error 404 page not found". Maybe I cannot access it. Can you share the title please? A colleague of mine is working on the verification too on the flatplate, we will soon share the results and we will build on top of these ones. |
Evaluation of the SST-SAS Model, L. Davidson I can also start working on BFS if you think you got the flatplate covered or we both can work on flatplate to cross check each other. What do you think? |
I think I got the flatplate covered, it should be easy (if the equations have been implemented correctly). The backward facing step and the diffuser I think are the most natural continuation as test cases, what do you think about it? |
I'll start working on BFS then. I'll let you know if anything comes up. |
I am currently working with @rois1995 on the flat plate case, taking as a reference result the skin friction coefficient computed using the coarse grid and the SA_DDES model (∆=∆max) from figure 3a in Hybrid RANS-LES Calculations in SU2. I ran the simulations using the already available SA_DDES model, to be compared with the reference result, and the models under development for SST (SST_DDES,SST_IDDES,SST_SIDDES). The reference data were extracted using a plot digitalizer, so they may be not so accurate, particularly in the zoom region. I attach also the density residuals for all the simulations and the TKE residuals for the simulations conducted with SST. The options SST_IDDES and SST_SIDDES do not allow a complete convergence and this issue is associated to a region close to the inlet and the symmetry boundary conditions. Although seemingly the implementation of the models does not have a direct dependence on the boundary conditions, this issue is absent for SST_DDES. As evidence of this undesirable and unexpected behavior, I attach the plot of the vertical velocity computed along a line parallel to the symmetry boundary (y=1e-6 , x between -0.33 and -0.05, with the flat plate beginning at x=0). I am currently working for understanding its cause. |
Proposed Changes
Hi everyone,
I have been working on the implementation of the DDES formulations for the SST model, following the work in "Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model" (DOI:10.1007/s10494-011-9378-4). The implementation is easy, whereas the validation may take some time. I am currently working on the backward facing step and I will work on the flatplate too. The choice of the validation test cases is key since the computational cost is quite large, thus if you have any suggestions that may speed up this part then feel free to write them down.
Implemented versions:
I am also working on the Scale Adaptive Simulations (SAS) implementations for SST, following the work in "ON URANS SOLUTIONS WITH LES-LIKE BEHAVIOUR", Travis et al., and "[Evaluation of scale-adaptive simulation for transonic cavity flows", Babu et al., (https://www.inderscienceonline.com/doi/abs/10.1504/IJESMS.2016.075510). The last one, especially, is quite tricky due to the computation of the velocity laplacian. I tried computing it in Paraview as the divergence of the gradient field and it seems quite similar, but it is the first time that I've touched that section of the code, thus some work might still be needed.
Implemented versions:
PR Checklist
Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.
pre-commit run --all
to format old commits.