-
Notifications
You must be signed in to change notification settings - Fork 843
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 Simplified LM Transition model #1901
base: develop
Are you sure you want to change the base?
Conversation
Technically, this model has only one equation, the one for intermittency. Most of the functions are the same as the original LM model. Should I consider the Simplified model (SLM) as an option for the LM model to avoid duplicates? This is how I started, but I am open to discussions and suggestions. |
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.
If the changes to the solver don't become too intrusive this approach sounds good to me
- Added correlations for Simplified LM. - There is a bug on the computation of grad(n*U)*n
I have a question: for each point in the mesh I am trying to compute the dot product between the velocity vector and the normal to the wall of the nearest point on the wall. How do I access such information? I've found that in the CPoint class I have the ClosestWall_Elem variable which stores the index of the closest element on a wall. However, when I try to assess the information with a number of cores greater than 2, it crashes. Moreover, to recover the normal of the element I perform a mean of the normals on the nodes of that element. Is there a structure that has the normals saved for each element of the primal grid? The part that I am referring to is from line 208 in CTransLMSolver.cpp . |
switch (options.Correlation_SLM) { | ||
case TURB_TRANS_CORRELATION_SLM::MENTER_SLM: { | ||
|
||
/*-- Thwaites parameter ---*/ | ||
su2double lambda_theta_local = 7.57e-3 * du_ds * wall_dist * wall_dist * Density / Laminar_Viscosity + 0.0128; | ||
lambda_theta_local = min(max(lambda_theta_local, -1.0), 1.0); | ||
|
||
/*-- Function to sensitize the transition onset to the streamwise pressure gradient ---*/ | ||
su2double FPG = 0.0; | ||
const su2double C_PG1 = 14.68; | ||
const su2double C_PG1_lim = 1.5; | ||
const su2double C_PG2 = -7.34; | ||
const su2double C_PG2_lim = 3.0; | ||
const su2double C_PG3 = 0.0; | ||
if (lambda_theta_local >= 0.0) { | ||
FPG = min(1+ C_PG1 * lambda_theta_local, C_PG1_lim); | ||
} else { | ||
const su2double FirstTerm = C_PG2 * lambda_theta_local; | ||
const su2double SecondTerm = C_PG3 * min(lambda_theta_local + 0.0681, 0.0); | ||
FPG = min(1 + FirstTerm + SecondTerm, C_PG2_lim); | ||
} | ||
|
||
FPG = max(FPG, 0.0); | ||
|
||
const su2double C_TU1 = 100.0; | ||
const su2double C_TU2 = 1000.0; | ||
const su2double C_TU3 = 1.0; | ||
rethetac = C_TU1 + C_TU2 * exp(-C_TU3 * Tu_L * FPG); | ||
|
||
break; | ||
} case TURB_TRANS_CORRELATION_SLM::CODER_SLM: { | ||
|
||
/*-- Local pressure gradient parameter ---*/ | ||
const su2double H_c = max(min(wall_dist * VorticityMag / VelocityMag, 1.1542), 0.3823); | ||
|
||
/*-- Thwaites parameter ---*/ | ||
su2double lambda_theta_local = 0.0; | ||
const su2double H_c_delta = 0.587743 - H_c; | ||
if ( H_c >= 0.587743 ) { | ||
const su2double FirstTerm = 0.1919 * pow(H_c_delta, 3.0); | ||
const su2double SecondTerm = 0.4182 * pow(H_c_delta, 2.0); | ||
const su2double ThirdTerm = 0.2959 * H_c_delta; | ||
lambda_theta_local = FirstTerm + SecondTerm + ThirdTerm; | ||
} else { | ||
const su2double FirstTerm = 4.7596 * pow(H_c_delta, 3.0); | ||
const su2double SecondTerm = -0.3837 * pow(H_c_delta, 2.0); | ||
const su2double ThirdTerm = 0.3575 * H_c_delta; | ||
lambda_theta_local = FirstTerm + SecondTerm + ThirdTerm; | ||
} | ||
|
||
/*-- Function to sensitize the transition onset to the streamwise pressure gradient ---*/ | ||
su2double FPG = 0.0; | ||
if (lambda_theta_local <= 0.0) { | ||
const su2double FirstTerm = -12.986 * lambda_theta_local; | ||
const su2double SecondTerm = -123.66 * pow(lambda_theta_local, 2.0); | ||
const su2double ThirdTerm = -405.689 * pow(lambda_theta_local, 3.0); | ||
FPG = 1 - (FirstTerm + SecondTerm + ThirdTerm) * exp(-pow(Tu_L/1.5,1.5)); | ||
} else { | ||
FPG = 1 + 0.275 * (1 - exp(-35.0 * lambda_theta_local)) * exp(-Tu_L/0.5); | ||
} | ||
|
||
// This is not reported in the paper | ||
//FPG = max(FPG, 0.0); | ||
|
||
const su2double C_TU1 = 100.0; | ||
const su2double C_TU2 = 1000.0; | ||
const su2double C_TU3 = 1.0; | ||
rethetac = C_TU1 + C_TU2 * exp(-C_TU3 * Tu_L * FPG); | ||
|
||
break; | ||
} case TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM: { | ||
|
||
/*-- Local pressure gradient parameter ---*/ | ||
const su2double H_c = max(min(wall_dist * VorticityMag / VelocityMag, 1.1542), 0.3823); | ||
|
||
/*-- H_32 Shape factor --*/ | ||
const su2double H_32 = 1.515095 + 0.2041 * pow((1.1542 - H_c), 2.0956); | ||
|
||
rethetac = exp(127.94 * pow((H_32-1.515095), 2.0) + 6.774224); | ||
|
||
break; | ||
} | ||
case TURB_TRANS_CORRELATION_SLM::DEFAULT: | ||
SU2_MPI::Error("Transition correlation for Simplified LM model is set to DEFAULT but no default value has ben set in the code.", | ||
CURRENT_FUNCTION); | ||
break; | ||
} |
Check notice
Code scanning / CodeQL
Long switch case Note
CODER_SLM (43 lines)
@@ -25,6 +25,7 @@ | |||
*/ | |||
|
|||
#pragma once | |||
//#include <cmath> |
Check notice
Code scanning / CodeQL
Commented-out code Note
- Modified LM_OPTIONS to include cross-flow effects: from LM2015 to CROSSFLOW
See what is done at the bottom of CGeometry.cpp in CGeometry::ComputeWallDistance. |
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.
CodeQL found more than 10 potential problems in the proposed changes. Check the Files changed tab for more details.
* \brief Get the index of the closest wall element. | ||
* \param[in] iPoint - Index of the point. | ||
*/ | ||
inline unsigned long GetClosestWall_Elem(unsigned long iPoint) {return ClosestWall_Elem(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.
Can you add the param[out] for these 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.
Sure, I'll do it now.
Common/include/option_structure.hpp
Outdated
|
||
if (NFoundCorrelations > 1) { | ||
SU2_MPI::Error("Two correlations selected for LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); | ||
} | ||
if (NFoundCorrelations_SLM > 1) { | ||
SU2_MPI::Error("Two correlations selected for Simplified model into LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); |
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.
SU2_MPI::Error("Two correlations selected for Simplified model into LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); | |
SU2_MPI::Error("Two correlations selected for simplified LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); |
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.
I'd stick with my change since the options are for the simplified model. They are not simplified options.
@@ -3471,9 +3471,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i | |||
if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { | |||
lmParsedOptions = ParseLMOptions(LM_Options, nLM_Options, rank, Kind_Turb_Model); | |||
|
|||
/*--- Check if problem is 2D and LM2015 has been selected ---*/ | |||
if (lmParsedOptions.LM2015 && val_nDim == 2) { | |||
SU2_MPI::Error("LM2015 is available only for 3D problems", CURRENT_FUNCTION); |
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 LM2015 gone? or is crossflow the same?
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.
I have changed the option to CROSSFLOW, since I will use it also for the Simplified model.
// This is not reported in the paper | ||
//FPG = max(FPG, 0.0); | ||
|
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.
// This is not reported in the paper | |
//FPG = max(FPG, 0.0); |
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.
I have left it there since I do not know if I have to keep it or not.
There is the geometry toolbox for dot product and normal: |
The problem is more related to the finding of the wall-normal for a point within the volume mesh, not to the computations that it will be involved in. |
The solution I suggested didn’t work? |
I still have to check if the implementation is correct but with more than two cores the code breaks. However, I found out that the wall-normal of a volume point can be computed as the normalized gradient of the wall-distance. Does this sound correct to you? However, there is a problem: I am using the aux variables to compute these gradients, but to compute dot(n, U) I first need n, thus I cannot compute them simultaneously. Since these computations are performed in the Preprocessing of the solvers, I was thinking to compute the normal within the FLOW_SOL preprocessing and the dot(n, U) in the TRANS_SOL preprocessing since the flow solver comes before the trans solver. Is this right? |
It's not correct. You need to follow the pattern from CGeometry::ComputeWallDistance |
- Added variables only for debug
I managed to implement the computation of grad(n*U)*n, but at the moment it is located into CTransLMSolver::PreProcessing. It seems to work with a structured mesh on a flat plate. Currently, I am testing with a 2D profile too. However, being into the PreProcessing of the transition solver, the normals are computed at each iteration, thus it is not computationally efficient if non-deforming meshes are used. I am looking into where to put it such that it is computed just if the mesh is updated (and at the first iteration of course). Plus, I have added a whole lot of variables to the output, but they will be removed in the final version. They are just used as debug. |
Do what we do with wall roughness and do it in the same place (computewalldistance and store in CPoint) |
Ok, I think I managed to add it to the ComputeWallDistance function. I tested it on a 2D cylinder and it works just fine. However, when I use it on a profile with a sharp trailing edge it gives me strange results on the region behind the TE. I guess it is intrinsic within the computation of the wall distances. |
My guess is that the closest element for those regions is somewhere on the trailing edge and the normals there change very quickly. |
Yeah, I do not think that it is solvable. Plus, I do not think that it will affect the solution either. |
Any updates @rois1995? Is this almost ready? |
I am sorry for leaving this hanging. I am currently back on working on it. These are the results on the Eppler E367 case with the new SLM model implemented and with three different correlation functions. Using the L2Roe scheme Using the Roe scheme I will update the branch to the last develop commit and proceed from there. EDIT: I think this one is too far behind with respect to the develop. Should I create another branch, make the changes there and then remove this one? |
@rois1995 I'd just merge it with develop and have a look at the conflicts in vscode. It's probably just some lines that have moved around. |
The problem is that I am behind more than 1100 commits. I have already created another branch starting from the current develop and I am running the test cases from the paper in the description. I may be able to merge the new branch into this one and then I'll be on par with the develop, right? |
if (config->GetViscous_Wall(iMarker)) { | ||
WallNormal_container[iZone][iMarker].resize(geometry->GetnElem_Bound(iMarker), 3); | ||
|
||
// cout << "geometry->nVertex[iMarker] = " << geometry->nVertex[iMarker] << endl; |
Check notice
Code scanning / CodeQL
Commented-out code Note
if (dist_i > 1e-10) { | ||
su2double Tu_L = 1.0; | ||
if (TurbFamily == TURB_FAMILY::KW) Tu_L = min(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / (ScalarVar_i[1]*dist_i), 100.0); | ||
// if (TurbFamily == TURB_FAMILY::KW) Tu_L = min(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / (Velocity_Mag), 100.0); |
Check notice
Code scanning / CodeQL
Commented-out code Note
// for (int j = 0; j < nDim; j++) | ||
// du_ds = du_ds + Velocity[i]*Velocity[j]*PrimVar_Grad_i[i][j]; | ||
|
||
// du_ds = du_ds/(Velocity_Mag*Velocity_Mag); |
Check notice
Code scanning / CodeQL
Commented-out code Note
|
||
// du_ds = du_ds/(Velocity_Mag*Velocity_Mag); | ||
|
||
// const su2double lambda_theta = -7.57e-3 * du_ds * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0128; |
Check notice
Code scanning / CodeQL
Commented-out code Note
// const su2double lambda_theta = du_ds * dist_i * dist_i * Density_i / Laminar_Viscosity_i; | ||
// duds_Here = du_ds; |
Check notice
Code scanning / CodeQL
Commented-out code Note
duds_Here = AuxVar; | ||
lambda_theta_Here = lambda_theta; | ||
|
||
// const su2double lambda_theta = 7.57e-3 * AuxVar * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0128; |
Check notice
Code scanning / CodeQL
Commented-out code Note
// cout << Re_t << " " << Corr_Rec << endl; | ||
// if (geometry->nodes->GetDomain(iPoint)) cout << "Point is on boundary" << endl; |
Check notice
Code scanning / CodeQL
Commented-out code Note
if(TurbFamily == TURB_FAMILY::SA) | ||
f_wake = 1.0; | ||
|
||
//cout << "StrainMag = " << StrainMag << " rho = " << rho << " dist = " << dist << " Re_v = " << Re_v << " Corr_Rec = " << Corr_Rec << endl; |
Check notice
Code scanning / CodeQL
Commented-out code Note
const su2double var2 = 1.0 - pow(var1,2.0); | ||
const su2double f_theta = min(max(f_wake*exp(-pow(dist/delta, 4)), var2), 1.0); | ||
su2double Intermittency_Sep = 2.0*max(0.0, Re_v/(3.235*Corr_Rec)-1.0)*f_reattach; | ||
//if (Intermittency_Sep>1.0) cout << "StrainMag = " << StrainMag << " rho = " << rho << " dist = " << dist << " Re_v = " << Re_v << " Corr_Rec = " << Corr_Rec << " Intermittency: " << Intermittency_Sep << " f_reattach = " << f_reattach << endl; |
Check notice
Code scanning / CodeQL
Commented-out code Note
Proposed Changes
Hi everyone! I am starting to implement the Simplified version of the Langtry-Menter transition model following the work in https://doi.org/10.2514/6.2012-672. More updates will follow.
Related Work
This follows the implementation of the LM model as #1751.