k-omega SST 2D axisymmetric source terms

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -305,8 +305,10 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   alfa_2,
   beta_1,
   beta_2,
-  sigma_omega_1,
-  sigma_omega_2,
+  sigma_k_1,
+  sigma_k_2,
+  sigma_w_1,
+  sigma_w_2,
   beta_star,
   a1;
 
@@ -320,12 +322,49 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
   bool incompressible;
   bool sustaining_terms;
+  bool axisymmetric;
 
   /*!
    * \brief A virtual member. Get strain magnitude based on perturbed reynolds stress matrix
    * \param[in] turb_ke: turbulent kinetic energy of the node
    */
   void SetPerturbedStrainMag(su2double turb_ke);
+
+  /*!
+   * \brief Add contribution due to axisymmetric formulation to 2D residual
+   */
+  inline void ResidualAxisymmetric(su2double alfa_blended, su2double zeta){
+
+    if (Coord_i[1] < EPS) return;
+
+    su2double yinv, rhov, k, w;
+    su2double sigma_k_i, sigma_w_i;
+    su2double pk_axi, pw_axi, cdk_axi, cdw_axi;
+
+    AD::SetPreaccIn(Coord_i[1]);
+
+    yinv = 1.0/Coord_i[1];
+    rhov = Density_i*V_i[2];
+    k = TurbVar_i[0];
+    w = TurbVar_i[1];
+
+    /*--- Compute blended constants ---*/
+    sigma_k_i = F1_i*sigma_k_1+(1.0-F1_i)*sigma_k_2;
+    sigma_w_i = F1_i*sigma_w_1+(1.0-F1_i)*sigma_w_2;
+
+    /*--- Production ---*/
+    pk_axi = max(0.0,2.0/3.0*rhov*k*(2.0/zeta*(yinv*V_i[2]-PrimVar_Grad_i[2][1]-PrimVar_Grad_i[1][0])-1.0));
+    pw_axi = alfa_blended*zeta/k*pk_axi;
+
+    /*--- Convection-Diffusion ---*/
+    cdk_axi = rhov*k-(Laminar_Viscosity_i+sigma_k_i*Eddy_Viscosity_i)*TurbVar_Grad_i[0][1];
+    cdw_axi = rhov*w-(Laminar_Viscosity_i+sigma_w_i*Eddy_Viscosity_i)*TurbVar_Grad_i[1][1];
+
+    /*--- Add terms to the residuals ---*/
+    Residual[0] += yinv*Volume*(pk_axi-cdk_axi);
+    Residual[1] += yinv*Volume*(pw_axi-cdw_axi);
+
+  }
 
 public:
   /*!

SU2_CFD/src/numerics/turbulent/turb_sources.cpp

@@ -760,16 +760,19 @@ CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim,
 
   incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
   sustaining_terms = (config->GetKind_Turb_Model() == SST_SUST);
+  axisymmetric = config->GetAxisymmetric();
 
   /*--- Closure constants ---*/
-  beta_star     = constants[6];
-  sigma_omega_1 = constants[2];
-  sigma_omega_2 = constants[3];
+  sigma_k_1     = constants[0];
+  sigma_k_2     = constants[1];
+  sigma_w_1     = constants[2];
+  sigma_w_2     = constants[3];
   beta_1        = constants[4];
   beta_2        = constants[5];
+  beta_star     = constants[6];
+  a1            = constants[7];
   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;
@@ -845,7 +848,6 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi
      pk = Eddy_Viscosity_i*StrainMag_i*StrainMag_i - 2.0/3.0*Density_i*TurbVar_i[0]*diverg;
    }
 
-
    pk = min(pk,20.0*beta_star*Density_i*TurbVar_i[1]*TurbVar_i[0]);
    pk = max(pk,0.0);
 
@@ -889,6 +891,10 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi
    /*--- Cross diffusion ---*/
 
    Residual[1] += (1.0 - F1_i)*CDkw_i*Volume;
+
+   /*--- Contribution due to 2D axisymmetric formulation ---*/
+
+   if (axisymmetric) ResidualAxisymmetric(alfa_blended,zeta);
 
    /*--- Implicit part ---*/
 

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -322,6 +322,8 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
 
 void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,
                                      CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
+
+  bool axisymmetric = config->GetAxisymmetric();
 
   CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
 
@@ -372,6 +374,11 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta
 
     numerics->SetCrossDiff(nodes->GetCrossDiff(iPoint),0.0);
 
+    if (axisymmetric){
+      /*--- Set y coordinate ---*/
+      numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
+    }
+
     /*--- Compute the source term ---*/
 
     auto residual = numerics->ComputeResidual(config);

TestCases/axisymmetric_rans/air_nozzle/air_nozzle.cfg

@@ -0,0 +1,214 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Axisymmetric supersonic converging-diverging air nozzle    %
+% Author: Florian Dittmann                                                     %
+% Date: 2021.12.02                                                             %
+% File Version 7.10 "Blackbird"                                                %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+% Physical governing equations (EULER, NAVIER_STOKES,
+%                               FEM_EULER, FEM_NAVIER_STOKES, FEM_RANS, FEM_LES,
+%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
+%                               POISSON_EQUATION)
+SOLVER= RANS
+%
+% Specify turbulence model (NONE, SA, SA_NEG, SST, SA_E, SA_COMP, SA_E_COMP)
+KIND_TURB_MODEL= SST
+%
+% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT)
+MATH_PROBLEM= DIRECT
+%
+% Restart solution (NO, YES)
+RESTART_SOL= YES
+%
+% System of measurements (SI, US)
+% International system of units (SI): ( meters, kilograms, Kelvins,
+%                                       Newtons = kg m/s^2, Pascals = N/m^2,
+%                                       Density = kg/m^3, Speed = m/s,
+%                                       Equiv. Area = m^2 )
+% United States customary units (US): ( inches, slug, Rankines, lbf = slug ft/s^2,
+%                                       psf = lbf/ft^2, Density = slug/ft^3,
+%                                       Speed = ft/s, Equiv. Area = ft^2 )
+SYSTEM_MEASUREMENTS= SI
+%
+AXISYMMETRIC= YES
+%
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+% Mach number (non-dimensional, based on the free-stream values)
+MACH_NUMBER= 1E-9
+%
+% Angle of attack (degrees, only for compressible flows)
+AOA= 0.0
+%
+% Side-slip angle (degrees, only for compressible flows)
+SIDESLIP_ANGLE= 0.0
+%
+% Init option to choose between Reynolds (default) or thermodynamics quantities
+% for initializing the solution (REYNOLDS, TD_CONDITIONS)
+INIT_OPTION= TD_CONDITIONS
+%
+% Free-stream option to choose between density and temperature (default) for
+% initializing the solution (TEMPERATURE_FS, DENSITY_FS)
+FREESTREAM_OPTION= TEMPERATURE_FS
+%
+% Free-stream pressure (101325.0 N/m^2, 2116.216 psf by default)
+FREESTREAM_PRESSURE= 1400000
+%
+% Free-stream temperature (288.15 K, 518.67 R by default)
+FREESTREAM_TEMPERATURE= 373.15
+%
+% Compressible flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE,
+%                              FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE)
+REF_DIMENSIONALIZATION= DIMENSIONAL
+
+% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------%
+%
+% Fluid model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS,
+%              CONSTANT_DENSITY, INC_IDEAL_GAS, INC_IDEAL_GAS_POLY)
+FLUID_MODEL= STANDARD_AIR
+
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY, POLYNOMIAL_VISCOSITY).
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+%
+% Molecular Viscosity that would be constant (1.716E-5 by default)
+MU_CONSTANT= 1.716E-5 
+
+% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------%
+%
+% Laminar Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL,
+% POLYNOMIAL_CONDUCTIVITY).
+CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
+%
+% Laminar Prandtl number (0.72 (air), only for CONSTANT_PRANDTL)
+PRANDTL_LAM= 0.72
+%
+% Turbulent Prandtl number (0.9 (air) by default)
+PRANDTL_TURB= 0.90
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Navier-Stokes (no-slip), constant heat flux wall  marker(s) (NONE = no marker)
+% Format: ( marker name, constant heat flux (J/m^2), ... )
+MARKER_HEATFLUX= ( WALL, 0.0 )
+%
+% Symmetry boundary marker(s) (NONE = no marker)
+MARKER_SYM= ( SYMMETRY )
+%
+% Riemann boundary marker(s) (NONE = no marker)
+% Format: (marker, data kind flag, list of data)
+MARKER_RIEMANN= ( INFLOW, TOTAL_CONDITIONS_PT, 1400000.0, 373.15, 1.0, 0.0, 0.0, OUTFLOW, STATIC_PRESSURE, 100000.0, 0.0, 0.0, 0.0, 0.0 )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
+NUM_METHOD_GRAD= GREEN_GAUSS
+%
+% CFL number (initial value for the adaptive CFL number)
+CFL_NUMBER= 1000.0
+%
+% Adaptive CFL number (NO, YES)
+CFL_ADAPT= NO
+%
+% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
+%                                        CFL max value )
+CFL_ADAPT_PARAM= ( 0.1, 2.0, 10.0, 1000.0 )
+%
+% Maximum Delta Time in local time stepping simulations
+MAX_DELTA_TIME= 1E6
+
+% ----------- SLOPE LIMITER AND DISSIPATION SENSOR DEFINITION -----------------%
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_FLOW= YES
+%
+% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG,
+%                BARTH_JESPERSEN, VAN_ALBADA_EDGE)
+SLOPE_LIMITER_FLOW= NONE
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI,
+%                                                      SMOOTHER_ILU, SMOOTHER_LUSGS,
+%                                                      SMOOTHER_LINELET)
+LINEAR_SOLVER= FGMRES
+%
+% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI)
+LINEAR_SOLVER_PREC= ILU
+%
+% Linael solver ILU preconditioner fill-in level (0 by default)
+LINEAR_SOLVER_ILU_FILL_IN= 0
+%
+% Minimum error of the linear solver for implicit formulations
+LINEAR_SOLVER_ERROR= 0.01
+%
+% Max number of iterations of the linear solver for the implicit formulation
+LINEAR_SOLVER_ITER= 10
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+% Multi-grid levels (0 = no multi-grid)
+MGLEVEL= 0
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, AUSMPLUSUP, AUSMPLUSUP2, HLLC,
+%                              TURKEL_PREC, MSW, FDS)
+CONV_NUM_METHOD_FLOW= ROE
+%
+% Entropy fix coefficient (0.0 implies no entropy fixing, 1.0 implies scalar
+%                          artificial dissipation)
+ENTROPY_FIX_COEFF= 0.1
+%
+% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+%
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+% Convective numerical method (SCALAR_UPWIND)
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+%
+% Time discretization (EULER_IMPLICIT)
+TIME_DISCRE_TURB= EULER_IMPLICIT
+%
+% Reduction factor of the CFL coefficient in the turbulence problem
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Number of total iterations
+ITER= 1000
+%
+% Convergence criteria (CAUCHY, RESIDUAL)
+CONV_CRITERIA= RESIDUAL
+%
+% Min value of the residual (log10 of the residual)
+CONV_RESIDUAL_MINVAL= -12
+%
+% Start convergence criteria at iteration number
+CONV_STARTITER= 10
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+% Mesh input file
+MESH_FILENAME= nozzle.su2
+%
+% Mesh input file format (SU2, CGNS)
+MESH_FORMAT= SU2
+%
+% Restart flow input file
+SOLUTION_FILENAME= solution_flow.dat
+%
+% Writing solution file frequency
+OUTPUT_WRT_FREQ= 1000
+%
+% Screen output
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION)

TestCases/hybrid_regression.py

@@ -210,6 +210,18 @@ def main():
     propeller.test_vals = [-3.389576, -8.409529, 0.000048, 0.056329]
     test_list.append(propeller)
 
+    #######################################
+    ### Axisymmetric Compressible RANS  ###
+    #######################################
+
+    # Axisymmetric air nozzle (transonic)
+    axi_rans_air_nozzle           = TestCase('axi_rans_air_nozzle')
+    axi_rans_air_nozzle.cfg_dir   = "axisymmetric_rans/air_nozzle"
+    axi_rans_air_nozzle.cfg_file  = "air_nozzle.cfg"
+    axi_rans_air_nozzle.test_iter = 10
+    axi_rans_air_nozzle.test_vals = [-12.094937, -6.622043, -8.814412, -2.393288]
+    test_list.append(axi_rans_air_nozzle)
+
     #################################
     ## Compressible RANS Restart  ###
     #################################

TestCases/parallel_regression.py

@@ -347,6 +347,21 @@ def main():
     propeller.tol       = 0.00001
     test_list.append(propeller)
 
+    #######################################
+    ### Axisymmetric Compressible RANS  ###
+    #######################################
+
+    # Axisymmetric air nozzle (transonic)
+    axi_rans_air_nozzle           = TestCase('axi_rans_air_nozzle')
+    axi_rans_air_nozzle.cfg_dir   = "axisymmetric_rans/air_nozzle"
+    axi_rans_air_nozzle.cfg_file  = "air_nozzle.cfg"
+    axi_rans_air_nozzle.test_iter = 10
+    axi_rans_air_nozzle.test_vals = [ -12.096569, -6.625843, -8.807541, -2.393279]
+    axi_rans_air_nozzle.su2_exec  = "mpirun -n 2 SU2_CFD"
+    axi_rans_air_nozzle.timeout   = 1600
+    axi_rans_air_nozzle.tol       = 0.0001
+    test_list.append(axi_rans_air_nozzle)
+
     #################################
     ## Compressible RANS Restart  ###
     #################################