Implementing the gust encounter for 3D meshes and for free-flying, elastic aircraft

Common/include/option_structure.hpp

@@ -779,11 +779,13 @@ static const MapType<std::string, ENUM_GUST_TYPE> Gust_Type_Map = {
  */
 enum ENUM_GUST_DIR {
   X_DIR = 0,  /*!< \brief Gust direction-X. */
-  Y_DIR = 1   /*!< \brief Gust direction-Y. */
+  Y_DIR = 1,  /*!< \brief Gust direction-Y. */
+  Z_DIR = 2   /*!< \brief Gust direction-Z. */
 };
 static const MapType<std::string, ENUM_GUST_DIR> Gust_Dir_Map = {
   MakePair("X_DIR", X_DIR)
   MakePair("Y_DIR", Y_DIR)
+  MakePair("Z_DIR", Z_DIR)
 };
 
 // If you add to ENUM_CENTERED, you must also add the option to ENUM_CONVECTIVE

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -302,13 +302,7 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
   // derivatives needed for the source term are calculated when applicable. If the gust derivatives are zero the source
   // term is also zero. The source term itself is implemented in the class CSourceWindGust
 
-  if (rank == MASTER_NODE) cout << endl << "Running simulation with a Wind Gust." << endl;
-  unsigned short iDim, nDim = geometry[MESH_0]->GetnDim();  // We assume nDim = 2
-  if (nDim != 2) {
-    if (rank == MASTER_NODE) {
-      cout << endl << "WARNING - Wind Gust capability is only verified for 2 dimensional simulations." << endl;
-    }
-  }
+  unsigned short iDim, nDim = geometry[MESH_0]->GetnDim();
 
   /*--- Gust Parameters from config ---*/
   unsigned short Gust_Type = config->GetGust_Type();
@@ -324,7 +318,7 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
   unsigned long iPoint;
   unsigned short iMGlevel, nMGlevel = config->GetnMGLevels();
 
-  su2double x, y, x_gust, dgust_dx, dgust_dy, dgust_dt;
+  su2double x, y, x_gust, dgust_dx, dgust_dy, dgust_dz, dgust_dt;
   su2double *Gust, *GridVel, *NewGridVel, *GustDer;
 
   su2double Physical_dt = config->GetDelta_UnstTime();
@@ -335,16 +329,17 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
 
   su2double Uinf = solver[MESH_0][FLOW_SOL]->GetVelocity_Inf(0);  // Assumption gust moves at infinity velocity
 
-  Gust = new su2double[nDim];
-  NewGridVel = new su2double[nDim];
-  for (iDim = 0; iDim < nDim; iDim++) {
-    Gust[iDim] = 0.0;
-    NewGridVel[iDim] = 0.0;
-  }
-
-  GustDer = new su2double[3];
-  for (unsigned short i = 0; i < 3; i++) {
-    GustDer[i] = 0.0;
+  Gust       = new su2double[nDim]();
+  NewGridVel = new su2double[nDim]();
+  GustDer    = new su2double[nDim+1]();
+
+  // Print some information to check that we are doing the right thing. Not sure how to convert the index back to a string...
+  if (rank == MASTER_NODE) {
+	  cout << endl << "Setting up a wind gust type " << Gust_Type << " with amplitude of " << gust_amp << " in direction " << GustDir << endl;
+	  cout << " U_inf      = " << Uinf << endl;
+	  cout << " Physical_t = " << Physical_t << endl;
+	  su2double loc_x = (xbegin + L + Uinf * (Physical_t - tbegin));
+	  cout << " Location_x = " << loc_x << endl;
   }
 
   // Vortex variables
@@ -366,7 +361,7 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
 
     for (iPoint = 0; iPoint < geometry[iMGlevel]->GetnPoint(); iPoint++) {
       /*--- Reset the Grid Velocity to zero if there is no grid movement ---*/
-      if (Kind_Grid_Movement == GUST && !(config->GetFSI_Simulation())) {
+      if (Kind_Grid_Movement == GUST && !(config->GetFSI_Simulation()) && !(config->GetDeform_Mesh())) {
         for (iDim = 0; iDim < nDim; iDim++) geometry[iMGlevel]->nodes->SetGridVel(iPoint, iDim, 0.0);
       }
 
@@ -377,6 +372,7 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
       }
       dgust_dx = 0.0;
       dgust_dy = 0.0;
+      dgust_dz = 0.0;
       dgust_dt = 0.0;
 
       /*--- Begin applying the gust ---*/
@@ -413,7 +409,7 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
           case ONE_M_COSINE:
             // Check if we are in the region where the gust is active
             if (x_gust > 0 && x_gust < n) {
-              Gust[GustDir] = gust_amp * (1 - cos(2 * PI_NUMBER * x_gust));
+              Gust[GustDir] = gust_amp * 0.5 * (1 - cos(2 * PI_NUMBER * x_gust));
 
               // Gust derivatives
               // dgust_dx = gust_amp*2*PI_NUMBER*(sin(2*PI_NUMBER*x_gust))/L;
@@ -454,11 +450,18 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
       }
 
       /*--- Set the Wind Gust, Wind Gust Derivatives and the Grid Velocities ---*/
-
-      GustDer[0] = dgust_dx;
-      GustDer[1] = dgust_dy;
-      GustDer[2] = dgust_dt;
-
+      if (nDim == 2) {
+        GustDer[0] = dgust_dx;
+        GustDer[1] = dgust_dy;
+        GustDer[2] = dgust_dt;
+      }
+      else {
+        GustDer[0] = dgust_dx;
+        GustDer[1] = dgust_dy;
+        GustDer[2] = dgust_dz;
+        GustDer[3] = dgust_dt;
+      }
+      // I think we don't need to set any source terms because they depend on the derivatives, which are zero in all cases from above.
       solver[iMGlevel][FLOW_SOL]->GetNodes()->SetWindGust(iPoint, Gust);
       solver[iMGlevel][FLOW_SOL]->GetNodes()->SetWindGustDer(iPoint, GustDer);
 

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -783,8 +783,13 @@ CNumerics::ResidualType<> CSourceWindGust::ComputeResidual(const CConfig* config
     residual[2] = smy*Volume;
     //residual[3] = smz*Volume;
     residual[3] = se*Volume;
-  } else {
-    SU2_MPI::Error("You should only be in the gust source term in two dimensions", CURRENT_FUNCTION);
+  }
+  else {
+    residual[0] = 0.0;
+    residual[1] = 0.0;
+    residual[2] = 0.0;
+    residual[3] = 0.0;
+    residual[4] = 0.0;
   }
 
   /*--- For now the source term Jacobian is just set to zero ---*/

TestCases/gust/cosine_gust_zdir.cfg

@@ -0,0 +1,111 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Cosine gust in z-direction of a 3D mesh                    %
+% Author: Arne Voß                                                             %
+% Institution: DLR                                                             %
+% Date: 25.05.2023                                                             %
+% File Version 7.5.1 "Blackbird"                                               %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= EULER
+MATH_PROBLEM= DIRECT
+%
+RESTART_SOL= YES
+RESTART_ITER= 72
+%
+% ------------------------------- SOLVER CONTROL ------------------------------%
+%
+INNER_ITER= 30
+CONV_RESIDUAL_MINVAL= -6
+CONV_STARTITER= 0
+%
+% ------------------------- UNSTEADY SIMULATION -------------------------------%
+%
+TIME_DOMAIN=YES
+TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+TIME_STEP= 0.001
+TIME_ITER= 80
+%
+% --------------------------- GUST SIMULATION ---------------------------------%
+%
+% The gust simulation requires the GRID_MOVEMENT flag to be set to YES.
+% and the GRID_MOVEMENT_KIND to be set to GUST or any of the other options.
+% Apply a wind gust (NO, YES)
+GRID_MOVEMENT=GUST
+WIND_GUST= YES
+GUST_TYPE= ONE_M_COSINE
+GUST_DIR= Z_DIR
+GUST_WAVELENGTH= 5.0
+GUST_PERIODS= 1.0
+% Gust amplitude corresponds to ~2.0 deg AoA
+GUST_AMPL= 2.37
+%
+GUST_BEGIN_TIME= 0.0
+% Gust is placed 5m ahead of the wing
+GUST_BEGIN_LOC=-10.0
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.2
+AOA= 0.0
+SIDESLIP_ANGLE= 0.0
+FREESTREAM_OPTION= DENSITY_FS
+FREESTREAM_PRESSURE= 101325.0
+FREESTREAM_DENSITY= 1.225
+FREESTREAM_TEMPERATURE= 288.15
+%
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 3.0
+SEMI_SPAN= 1.5
+REF_DIMENSIONALIZATION= DIMENSIONAL
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_EULER= ( wing, tip)
+MARKER_FAR= ( far_away )
+
+% ------------------------ SURFACES IDENTIFICATION ----------------------------%
+%
+MARKER_PLOTTING = ( wing, tip )
+MARKER_MONITORING = ( wing, tip )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+CFL_NUMBER= 1000.0
+CFL_ADAPT= YES
+%
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+MGLEVEL= 3
+MGCYCLE= W_CYCLE
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= JST
+LINEAR_SOLVER_ITER= 20
+%
+% ------------------------- SCREEN/HISTORY VOLUME OUTPUT --------------------------%
+%
+SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, LIFT, MOMENT_X, MOMENT_Y, MOMENT_Z)
+HISTORY_OUTPUT= (ITER, RMS_RES, AERO_COEFF, CAUCHY, WALL_TIME)
+%
+% ------------------------- INPUT/OUTPUT FILE INFORMATION --------------------------%
+%
+MESH_FILENAME= mesh_rectangular_wing.su2
+MESH_FORMAT= SU2
+SOLUTION_FILENAME= restart_gust.dat
+RESTART_FILENAME= restart_gust.dat
+SURFACE_FILENAME= surface_gust
+VOLUME_FILENAME= volume_gust
+CONV_FILENAME= history_gust
+TABULAR_FORMAT= CSV
+OUTPUT_FILES= (RESTART, RESTART_ASCII)
+%

TestCases/gust/gust_with_mesh_deformation.cfg

@@ -0,0 +1,122 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Cosine gust combined with mesh deformation                 %
+% Author: Arne Voß                                                             %
+% Institution: DLR                                                             %
+% Date: 25.05.2023                                                             %
+% File Version 7.5.1 "Blackbird"                                               %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= EULER
+MATH_PROBLEM= DIRECT
+% ------------------------------- SOLVER CONTROL ------------------------------%
+%
+INNER_ITER= 100
+CONV_RESIDUAL_MINVAL= -6
+CONV_STARTITER= 0
+%
+% ------------------------- UNSTEADY SIMULATION -------------------------------%
+%
+TIME_DOMAIN=YES
+TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+TIME_STEP= 0.001
+TIME_ITER= 7
+%
+% --------------------------- GUST SIMULATION ---------------------------------%
+%
+% The gust simulation requires the GRID_MOVEMENT flag to be set to YES.
+% and the GRID_MOVEMENT_KIND to be set to GUST or any of the other options.
+% Apply a wind gust (NO, YES)
+GRID_MOVEMENT=GUST
+WIND_GUST= YES
+GUST_TYPE= ONE_M_COSINE
+GUST_DIR= Y_DIR
+GUST_WAVELENGTH= 5.0
+GUST_PERIODS= 1.0
+GUST_AMPL= 1.0
+GUST_BEGIN_TIME= 0.0
+GUST_BEGIN_LOC=-5.0
+%
+% --------------------------- MESH DEFORMATION--------------------------------%
+% Type of dynamic surface movement (NONE, DEFORMING, MOVING_WALL,
+% AEROELASTIC, AEROELASTIC_RIGID_MOTION EXTERNAL, EXTERNAL_ROTATION)
+SURFACE_MOVEMENT= DEFORMING
+%
+% Moving wall boundary marker(s) (NONE = no marker, ignored for RIGID_MOTION)
+MARKER_MOVING= ( airfoil )
+%
+% Plunging angular freq. (rad/s) in x, y, & z directions
+SURFACE_PLUNGING_OMEGA= 0.0 125.6 0.0
+%
+% Plunging amplitude (m or ft) in x, y, & z directions
+SURFACE_PLUNGING_AMPL= 0.0 0.0001 0.0
+%
+% Move Motion Origin for marker moving (1 or 0)
+MOVE_MOTION_ORIGIN = 0
+%
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.2
+AOA= 0.0
+SIDESLIP_ANGLE= 0.0
+FREESTREAM_OPTION= DENSITY_FS
+FREESTREAM_PRESSURE= 101325.0
+FREESTREAM_DENSITY= 1.225
+FREESTREAM_TEMPERATURE= 288.15
+%
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+SEMI_SPAN= 1.0
+REF_DIMENSIONALIZATION= DIMENSIONAL
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_EULER= ( airfoil )
+MARKER_FAR= ( farfield )
+%
+% ------------------------ SURFACES IDENTIFICATION ----------------------------%
+%
+MARKER_PLOTTING = ( airfoil )
+MARKER_MONITORING = ( airfoil )
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+CFL_NUMBER= 100.0
+CFL_ADAPT= YES
+%
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+MGLEVEL= 3
+MGCYCLE= W_CYCLE
+%
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= JST
+LINEAR_SOLVER_ITER= 20
+%
+% ------------------------- SCREEN/HISTORY VOLUME OUTPUT --------------------------%
+%
+SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, LIFT, MOMENT_Z)
+HISTORY_OUTPUT= (ITER, RMS_RES, AERO_COEFF, CAUCHY, WALL_TIME)
+%
+% ------------------------- INPUT/OUTPUT FILE INFORMATION --------------------------%
+%
+MESH_FILENAME= mesh_NACA0012_inv.su2
+MESH_FORMAT= SU2
+SOLUTION_FILENAME= restart_gust.dat
+RESTART_FILENAME= restart_gust.dat
+SURFACE_FILENAME= surface_gust
+VOLUME_FILENAME= volume_gust
+CONV_FILENAME= history_gust
+TABULAR_FORMAT= CSV
+OUTPUT_FILES= (NONE)
+%

TestCases/hybrid_regression.py

@@ -470,6 +470,24 @@ def main():
     sine_gust.test_vals = [-1.977520, 3.481804, -0.012402, -0.007454]
     sine_gust.unsteady  = True
     test_list.append(sine_gust)
+
+    # Cosine gust in z-direction
+    cosine_gust           = TestCase('cosine_gust_zdir')
+    cosine_gust.cfg_dir   = "gust"
+    cosine_gust.cfg_file  = "cosine_gust_zdir.cfg"
+    cosine_gust.test_iter = 79
+    cosine_gust.test_vals = [-2.418813, 0.004650, -0.001878, -0.000637, -0.000271]
+    cosine_gust.unsteady  = True
+    test_list.append(cosine_gust)
+
+    # Gust with mesh deformation
+    gust_mesh_defo           = TestCase('gust_with_mesh_deformation')
+    gust_mesh_defo.cfg_dir   = "gust"
+    gust_mesh_defo.cfg_file  = "gust_with_mesh_deformation.cfg"
+    gust_mesh_defo.test_iter = 6
+    gust_mesh_defo.test_vals = [-1.844778, 0.000846, -0.000408]
+    gust_mesh_defo.unsteady  = True
+    test_list.append(gust_mesh_defo)
 
     # Aeroelastic
     aeroelastic           = TestCase('aeroelastic')