Bas-Cakmakcioglu (BC) Transition Model

Common/include/option_structure.hpp

@@ -641,11 +641,13 @@ static const map<string, ENUM_TURB_MODEL> Turb_Model_Map = CCreateMap<string, EN
  */
 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 (LM for 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>
 ("NONE", NO_TRANS_MODEL)
-("LM", LM);
+("LM", LM)
+("BC", BC); //BAS-CAKMAKCIOGLU
 
 /*!
  * \brief type of time integration schemes

SU2_CFD/include/numerics_structure.hpp

@@ -4202,6 +4202,7 @@ class CSourcePieceWise_TurbSA : public CNumerics {
   su2double dr, dg, dfw;
   bool incompressible;
   bool rotating_frame;
+  bool transition;
   su2double intermittency;
   su2double Production, Destruction, CrossProduction;
 

SU2_CFD/src/numerics_direct_turbulent.cpp

@@ -506,6 +506,7 @@ CSourcePieceWise_TurbSA::CSourcePieceWise_TurbSA(unsigned short val_nDim, unsign
 
   incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
   rotating_frame = config->GetRotating_Frame();
+  transition = (config->GetKind_Trans_Model() == BC);
 
   /*--- Spalart-Allmaras closure constants ---*/
 
@@ -526,7 +527,7 @@ CSourcePieceWise_TurbSA::CSourcePieceWise_TurbSA(unsigned short val_nDim, unsign
 CSourcePieceWise_TurbSA::~CSourcePieceWise_TurbSA(void) { }
 
 void CSourcePieceWise_TurbSA::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) {
-  
+	
 //  AD::StartPreacc();
 //  AD::SetPreaccIn(V_i, nDim+6);
 //  AD::SetPreaccIn(Vorticity_i, nDim);
@@ -535,6 +536,10 @@ void CSourcePieceWise_TurbSA::ComputeResidual(su2double *val_residual, su2double
 //  AD::SetPreaccIn(TurbVar_Grad_i[0], nDim);
 //  AD::SetPreaccIn(Volume); AD::SetPreaccIn(dist_i);
 
+//  BC Transition Model variables
+  su2double vmag, rey, re_theta, re_theta_t, re_v;
+  su2double tu , nu_cr, nu_t, nu_BC, chi_1, chi_2, gamma_BC, term1, term2, term_exponential;
+
   if (incompressible) {
     Density_i = V_i[nDim+1];
     Laminar_Viscosity_i = V_i[nDim+3];
@@ -550,6 +555,18 @@ void CSourcePieceWise_TurbSA::ComputeResidual(su2double *val_residual, su2double
   CrossProduction = 0.0;
   val_Jacobian_i[0][0] = 0.0;
 
+  gamma_BC = 0.0;
+  vmag = 0.0;
+  tu   = config->GetTurbulenceIntensity_FreeStream();
+  rey  = config->GetReynolds();
+
+  if (nDim==2) {
+    vmag = sqrt(V_i[1]*V_i[1]+V_i[2]*V_i[2]);
+  }
+  else if (nDim==3) {
+    vmag = sqrt(V_i[1]*V_i[1]+V_i[2]*V_i[2]+V_i[3]*V_i[3]);
+  }
+
   /*--- Evaluate Omega ---*/
 
   Omega = sqrt(Vorticity_i[0]*Vorticity_i[0] + Vorticity_i[1]*Vorticity_i[1] + Vorticity_i[2]*Vorticity_i[2]);
@@ -576,14 +593,36 @@ void CSourcePieceWise_TurbSA::ComputeResidual(su2double *val_residual, su2double
     Shat = S + TurbVar_i[0]*fv2*inv_k2_d2;
     Shat = max(Shat, 1.0e-10);
     inv_Shat = 1.0/Shat;
-
-    /*--- Production term ---*/;
 
 //    Original SA model
 //    Production = cb1*(1.0-ft2)*Shat*TurbVar_i[0]*Volume;
 
-    Production = cb1*Shat*TurbVar_i[0]*Volume;
+    if (transition) {
+
+//    BC model constants	  
+      chi_1 = 0.002;
+      chi_2 = 5.0;
+
+      nu_t = (TurbVar_i[0]*fv1); //S-A variable
+      nu_cr = chi_2/rey;
+      nu_BC = (nu_t)/(vmag*dist_i);
 
+      re_v   = ((Density_i*pow(dist_i,2.))/(Laminar_Viscosity_i))*Omega;
+      re_theta = re_v/2.193;
+      re_theta_t = (803.73 * pow((tu + 0.6067),-1.027)); //MENTER correlation
+      //re_theta_t = 163.0 + exp(6.91-tu); //ABU-GHANNAM & SHAW correlation
+
+      term1 = sqrt(max(re_theta-re_theta_t,0.)/(chi_1*re_theta_t));
+      term2 = sqrt(max(nu_BC-nu_cr,0.)/(nu_cr));
+      term_exponential = (term1 + term2);
+      gamma_BC = 1.0 - exp(-term_exponential);
+
+      Production = gamma_BC*cb1*Shat*TurbVar_i[0]*Volume;
+    }
+    else {
+      Production = cb1*Shat*TurbVar_i[0]*Volume;
+    }
+
     /*--- Destruction term ---*/
 
     r = min(TurbVar_i[0]*inv_Shat*inv_k2_d2,10.0);
@@ -613,7 +652,13 @@ void CSourcePieceWise_TurbSA::ComputeResidual(su2double *val_residual, su2double
     dfv2 = -(1/nu-Ji_2*dfv1)/pow(1.+Ji*fv1,2.);
     if ( Shat <= 1.0e-10 ) dShat = 0.0;
     else dShat = (fv2+TurbVar_i[0]*dfv2)*inv_k2_d2;
-    val_Jacobian_i[0][0] += cb1*(TurbVar_i[0]*dShat+Shat)*Volume;
+
+    if (transition) {
+        val_Jacobian_i[0][0] += gamma_BC*cb1*(TurbVar_i[0]*dShat+Shat)*Volume;
+    }
+    else {
+        val_Jacobian_i[0][0] += cb1*(TurbVar_i[0]*dShat+Shat)*Volume;
+    }
 
     /*--- Implicit part, destruction term ---*/
 

SU2_CFD/src/solver_direct_turbulent.cpp

@@ -1255,13 +1255,19 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
   Viscosity_Inf = config->GetViscosity_FreeStreamND();
 
   /*--- Factor_nu_Inf in [3.0, 5.0] ---*/
-  
+
   Factor_nu_Inf = config->GetNuFactor_FreeStream();
   nu_tilde_Inf  = Factor_nu_Inf*Viscosity_Inf/Density_Inf;
+  if (config->GetKind_Trans_Model() == BC) {
+    nu_tilde_Inf  = 0.005*Factor_nu_Inf*Viscosity_Inf/Density_Inf;
+  }
 
   /*--- Factor_nu_Engine ---*/
   Factor_nu_Engine = config->GetNuFactor_Engine();
   nu_tilde_Engine  = Factor_nu_Engine*Viscosity_Inf/Density_Inf;
+  if (config->GetKind_Trans_Model() == BC) {
+    nu_tilde_Engine  = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf;
+  }
 
   /*--- Factor_nu_ActDisk ---*/
   Factor_nu_ActDisk = config->GetNuFactor_Engine();

TestCases/parallel_regression.py

@@ -234,6 +234,21 @@ def main():
     inc_turb_naca0012.tol       = 0.00001
     test_list.append(inc_turb_naca0012)
 
+    ############################
+    ###      Transition      ###
+    ############################
+
+    # Schubauer-Klebanoff Natural Transition Case
+    schubauer_klebanoff_transition              = TestCase('Schubauer_Klebanoff')
+    schubauer_klebanoff_transition.cfg_dir      = "transition/Schubauer_Klebanoff"
+    schubauer_klebanoff_transition.cfg_file     = "transitional_BC_model_ConfigFile.cfg"
+    schubauer_klebanoff_transition.test_iter    = 250
+    schubauer_klebanoff_transition.test_vals    = [-9.467876, -15.630835, 0.000520, 0.001062] #last 4 columns
+    schubauer_klebanoff_transition.su2_exec     = "parallel_computation.py -f"
+    schubauer_klebanoff_transition.timeout      = 1600
+    schubauer_klebanoff_transition.tol          = 0.00001
+    test_list.append(schubauer_klebanoff_transition)
+
     #####################################
     ### Cont. adj. compressible Euler ###
     #####################################

TestCases/serial_regression.py

@@ -233,6 +233,21 @@ def main():
     inc_turb_naca0012.tol       = 0.00001
     test_list.append(inc_turb_naca0012)
 
+    #########################
+    ###    Transition     ###
+    #########################
+
+    # Schubauer-Klebanoff Natural Transition
+    schubauer_klebanoff_transition              = TestCase('Schubauer_Klebanoff')
+    schubauer_klebanoff_transition.cfg_dir      = "transition/Schubauer_Klebanoff"
+    schubauer_klebanoff_transition.cfg_file     = "transitional_BC_model_ConfigFile.cfg"
+    schubauer_klebanoff_transition.test_iter    = 250
+    schubauer_klebanoff_transition.test_vals    = [-9.474452, -15.630834, 0.000520, 0.001062] #last 4 columns
+    schubauer_klebanoff_transition.su2_exec     = "SU2_CFD"
+    schubauer_klebanoff_transition.timeout      = 1600
+    schubauer_klebanoff_transition.tol          = 0.00001
+    test_list.append(schubauer_klebanoff_transition)
+
     #####################################
     ### Cont. adj. compressible Euler ###
     #####################################