Merge pull request #2538 from JanRotti/develop

Camberline Hicks-Henne Parameterization as design feature

Author: pcarruscag

Common/include/geometry/CPhysicalGeometry.hpp

@@ -620,15 +620,15 @@ class CPhysicalGeometry final : public CGeometry {
 
   /*!
    * \brief Compute the maximum thickness of an airfoil.
-   * \return Maximum thickness at a particular seccion.
+   * \return Maximum thickness at a particular section.
    */
   su2double Compute_MaxThickness(su2double* Plane_P0, su2double* Plane_Normal, CConfig* config,
                                  vector<su2double>& Xcoord_Airfoil, vector<su2double>& Ycoord_Airfoil,
                                  vector<su2double>& Zcoord_Airfoil) override;
 
   /*!
    * \brief Compute the twist of an airfoil.
-   * \return Twist at a particular seccion.
+   * \return Twist at a particular section.
    */
   su2double Compute_Twist(su2double* Plane_P0, su2double* Plane_Normal, vector<su2double>& Xcoord_Airfoil,
                           vector<su2double>& Ycoord_Airfoil, vector<su2double>& Zcoord_Airfoil) override;
@@ -708,7 +708,7 @@ class CPhysicalGeometry final : public CGeometry {
 
   /*!
    * \brief Compute the dihedral of a wing.
-   * \return Dihedral at a particular seccion.
+   * \return Dihedral at a particular section.
    */
   su2double Compute_Dihedral(su2double* LeadingEdge_im1, su2double* TrailingEdge_im1, su2double* LeadingEdge_i,
                              su2double* TrailingEdge_i) override;

Common/include/grid_movement/CSurfaceMovement.hpp

@@ -95,6 +95,13 @@ class CSurfaceMovement : public CGridMovement {
    */
   void SetCST(CGeometry* boundary, CConfig* config, unsigned short iDV, bool ResetDef);
 
+  /*!
+   * \brief Set a Hicks-Henne deformation bump function on the camberline of an airfoil.
+   * \param[in] boundary - Geometry of the boundary.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetHicksHenneCamber(CGeometry* boundary, CConfig* config);
+
   /*!
    * \brief Set a NACA 4 digits airfoil family for airfoil deformation.
    * \param[in] boundary - Geometry of the boundary.

Common/include/option_structure.hpp

@@ -2287,6 +2287,7 @@ enum ENUM_PARAM {
   CST = 34,                   /*!< \brief CST method with Kulfan parameters for airfoil deformation. */
   SURFACE_BUMP = 35,          /*!< \brief Surfacebump function for flat surfaces deformation. */
   SURFACE_FILE = 36,          /*!< \brief Nodal coordinates for surface set using a file (external parameterization). */
+  HICKS_HENNE_CAMBER = 37,    /*!< \brief Hicks-Henne bump function for camber line deformation. */
   DV_EFIELD = 40,             /*!< \brief Electric field in deformable membranes. */
   DV_YOUNG = 41,
   DV_POISSON = 42,
@@ -2304,6 +2305,7 @@ static const MapType<std::string, ENUM_PARAM> Param_Map = {
   MakePair("FFD_CAMBER_2D", FFD_CAMBER_2D)
   MakePair("FFD_THICKNESS_2D", FFD_THICKNESS_2D)
   MakePair("HICKS_HENNE", HICKS_HENNE)
+  MakePair("HICKS_HENNE_CAMBER", HICKS_HENNE_CAMBER)
   MakePair("SURFACE_BUMP", SURFACE_BUMP)
   MakePair("ANGLE_OF_ATTACK", ANGLE_OF_ATTACK)
   MakePair("NACA_4DIGITS", NACA_4DIGITS)

Common/include/option_structure.inl

@@ -26,6 +26,7 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "option_structure.hpp"
 #include "parallelization/mpi_structure.hpp"
 using namespace std;
 
@@ -570,6 +571,9 @@ class COptionDVParam : public COptionBase {
         case HICKS_HENNE:
           nParamDV[iDV] = 2;
           break;
+        case HICKS_HENNE_CAMBER:
+          nParamDV[iDV] = 1;
+          break;
         case SURFACE_BUMP:
           nParamDV[iDV] = 3;
           break;

Common/src/CConfig.cpp

@@ -2372,6 +2372,7 @@ void CConfig::SetConfig_Options() {
    - FFD_CAMBER_2D ( FFDBox ID, i_Ind )
    - FFD_THICKNESS_2D ( FFDBox ID, i_Ind )
    - HICKS_HENNE ( Lower Surface (0)/Upper Surface (1)/Only one Surface (2), x_Loc )
+   - HICKS_HENNE_CAMBER ( x_Loc )
    - SURFACE_BUMP ( x_start, x_end, x_Loc )
    - CST ( Lower Surface (0)/Upper Surface (1), Kulfan parameter number, Total number of Kulfan parameters for surface )
    - NACA_4DIGITS ( 1st digit, 2nd digit, 3rd and 4th digit )
@@ -6695,6 +6696,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           case FFD_CAMBER_2D:         cout << "FFD 2D (camber) <-> "; break;
           case FFD_THICKNESS_2D:      cout << "FFD 2D (thickness) <-> "; break;
           case HICKS_HENNE:           cout << "Hicks Henne <-> " ; break;
+          case HICKS_HENNE_CAMBER:    cout << "Hicks Henne (camber) <-> " ; break;
           case SURFACE_BUMP:          cout << "Surface bump <-> " ; break;
           case ANGLE_OF_ATTACK:       cout << "Angle of attack <-> " ; break;
           case CST:                   cout << "Kulfan parameter number (CST) <-> " ; break;
@@ -6730,7 +6732,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         if ((Design_Variable[iDV] == NO_DEFORMATION) ||
             (Design_Variable[iDV] == FFD_SETTING) ||
             (Design_Variable[iDV] == SCALE) ) nParamDV = 0;
-        if (Design_Variable[iDV] == ANGLE_OF_ATTACK) nParamDV = 1;
+        if ((Design_Variable[iDV] == ANGLE_OF_ATTACK) ||
+            (Design_Variable[iDV] == HICKS_HENNE_CAMBER)) nParamDV = 1;
         if ((Design_Variable[iDV] == FFD_CAMBER_2D) ||
             (Design_Variable[iDV] == FFD_THICKNESS_2D) ||
             (Design_Variable[iDV] == HICKS_HENNE) ||

Common/src/grid_movement/CSurfaceMovement.cpp

@@ -584,7 +584,7 @@ vector<vector<su2double> > CSurfaceMovement::SetSurface_Deformation(CGeometry* g
   else if ((config->GetDesign_Variable(0) == ROTATION) || (config->GetDesign_Variable(0) == TRANSLATION) ||
            (config->GetDesign_Variable(0) == SCALE) || (config->GetDesign_Variable(0) == HICKS_HENNE) ||
            (config->GetDesign_Variable(0) == SURFACE_BUMP) || (config->GetDesign_Variable(0) == ANGLE_OF_ATTACK) ||
-           (config->GetDesign_Variable(0) == CST)) {
+           (config->GetDesign_Variable(0) == CST) || (config->GetDesign_Variable(0) == HICKS_HENNE_CAMBER)) {
     /*--- Apply rotation, displacement and stretching design variables (this
      should be done before the bump function design variables) ---*/
 
@@ -640,6 +640,11 @@ vector<vector<su2double> > CSurfaceMovement::SetSurface_Deformation(CGeometry* g
       }
     }
 
+    /*--- HICKS_HENNE_CAMBER design variable ---*/
+
+    if (config->GetDesign_Variable(0) == HICKS_HENNE_CAMBER) {
+      SetHicksHenneCamber(geometry, config);
+    }
   }
 
   /*--- NACA_4Digits design variable ---*/
@@ -3696,6 +3701,130 @@ void CSurfaceMovement::SetExternal_Deformation(CGeometry* geometry, CConfig* con
   }
 }
 
+void CSurfaceMovement::SetHicksHenneCamber(CGeometry* boundary, CConfig* config) {
+  unsigned long iVertex;
+  unsigned short iMarker, nDV_Camber = 0;
+  su2double VarCoord[3] = {0.0, 0.0, 0.0}, VarCoordTrans[3] = {0.0, 0.0, 0.0}, *CoordTrans, *NormalTrans, ek, fk,
+            Coord[3] = {0.0, 0.0, 0.0}, TPCoord[2] = {0.0, 0.0}, LPCoord[2] = {0.0, 0.0}, USTPCoord[2] = {0.0, 0.0},
+            LSTPCoord[2] = {0.0, 0.0}, Distance, Chord, AoA, ValCos, ValSin;
+  vector<su2double> positions, values, X_Coord_upper, Y_Coord_upper, X_Coord_lower, Y_Coord_lower;
+
+  // --- Check if the type of design variables is only HICKS_HENNE_CAMBER ---// // TODO: Extend to CAMBER + THICKNESS
+  for (unsigned short iDV = 0; iDV < config->GetnDV(); iDV++)
+    if (config->GetDesign_Variable(iDV) != HICKS_HENNE_CAMBER) {
+      SU2_MPI::Error("'HicksHenneCamber' can not be combined with other design variables.", CURRENT_FUNCTION);
+    } else {
+      nDV_Camber++;
+    }
+
+  positions.resize(nDV_Camber);
+  values.resize(nDV_Camber);
+
+  /*--- Collect the values of the Hicks-Henne camber design variables ---*/
+  for (unsigned short iDV = 0, counter = 0; iDV < config->GetnDV(); iDV++) {
+    if (config->GetDesign_Variable(iDV) == HICKS_HENNE_CAMBER) {
+      positions[counter] =
+          config->GetParamDV(iDV, 0);             /*--- Position of the camber point as a fraction of the chord ---*/
+      values[counter] = config->GetDV_Value(iDV); /*--- Value of the deformation ---*/
+      counter++;
+    }
+  }
+
+  /*--- Compute the angle of attack, Leading-edge (LP) and Trailing-edge (TP) point.
+        We do this for upper and lower surface separately to detect blunt trailing edges ---*/
+
+  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    if (config->GetMarker_All_DV(iMarker) == YES) {
+      CoordTrans = boundary->vertex[iMarker][0]->GetCoord();
+      LSTPCoord[0] = CoordTrans[0];
+      LSTPCoord[1] = CoordTrans[1];
+      USTPCoord[0] = CoordTrans[0];
+      USTPCoord[1] = CoordTrans[1];
+      for (iVertex = 1; iVertex < boundary->nVertex[iMarker]; iVertex++) {
+        CoordTrans = boundary->vertex[iMarker][iVertex]->GetCoord();
+        NormalTrans = boundary->vertex[iMarker][0]->GetNormal();
+        if ((CoordTrans[0] > TPCoord[0]) && (abs(NormalTrans[1]) > abs(NormalTrans[0]))) {
+          if (NormalTrans[1] >= 0.0) {
+            USTPCoord[0] = CoordTrans[0];
+            USTPCoord[1] = CoordTrans[1];
+          } else {
+            LSTPCoord[0] = CoordTrans[0];
+            LSTPCoord[1] = CoordTrans[1];
+          }
+        }
+      }
+    }
+  }
+
+  /*--- Correct the TP coordinates if the trailing edge is blunt ---*/
+  if ((fabs(USTPCoord[0] - LSTPCoord[0]) > EPS) || (fabs(USTPCoord[1] - LSTPCoord[1]) > EPS)) {
+    TPCoord[0] = (USTPCoord[0] + LSTPCoord[0]) / 2.0;
+    TPCoord[1] = (USTPCoord[1] + LSTPCoord[1]) / 2.0;
+  } else {
+    TPCoord[0] = USTPCoord[0];
+    TPCoord[1] = USTPCoord[1];
+  }
+
+  Chord = 0.0;
+  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    if (config->GetMarker_All_DV(iMarker) == YES) {
+      for (iVertex = 0; iVertex < boundary->nVertex[iMarker]; iVertex++) {
+        CoordTrans = boundary->vertex[iMarker][iVertex]->GetCoord();
+        Distance = GeometryToolbox::Distance(2, CoordTrans, TPCoord);
+        if (Chord < Distance) {
+          Chord = Distance;
+          LPCoord[0] = CoordTrans[0];
+          LPCoord[1] = CoordTrans[1];
+        }
+      }
+    }
+  }
+
+  AoA = atan((LPCoord[1] - TPCoord[1]) / (TPCoord[0] - LPCoord[0])) * 180 / PI_NUMBER;
+
+  /*--- Apply deformation based on the camberline Hicks-Henne deformation and the existing airfoil thickness
+   * distribution ---*/
+  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    for (iVertex = 0; iVertex < boundary->nVertex[iMarker]; iVertex++) {
+      if (config->GetMarker_All_DV(iMarker) == YES) {
+        CoordTrans = boundary->vertex[iMarker][iVertex]->GetCoord();
+        NormalTrans = boundary->vertex[iMarker][iVertex]->GetNormal();
+
+        ValCos = cos(AoA * PI_NUMBER / 180.0);
+        ValSin = sin(AoA * PI_NUMBER / 180.0);
+
+        /*--- Transform the coordinates ---*/
+        Coord[0] = (CoordTrans[0] - LPCoord[0]) * ValCos - (CoordTrans[1] - LPCoord[1]) * ValSin;
+        Coord[0] = max(0.0, Coord[0]);  // Coord x should be always positive
+        Coord[1] = (CoordTrans[1] - LPCoord[1]) * ValCos + (CoordTrans[0] - LPCoord[0]) * ValSin;
+
+        /*--- Special case: surface point is part of trailing edge ---*/
+        if (Coord[0] > 0.99) {
+          continue;
+        }
+
+        /*--- Compute the Hicks-Henne bumps ---*/
+        fk = 0.0;
+        for (unsigned short iDV = 0; iDV < nDV_Camber; iDV++) {
+          ek = log10(0.5) / log10(positions[iDV]);
+          if (Coord[0] > 10 * EPS) fk += pow(sin(PI_NUMBER * pow(Coord[0], ek)), 3) * values[iDV];
+        }
+        VarCoord[1] = fk;
+      }
+
+      /*--- Apply the transformation to the coordinate variation ---*/
+
+      ValCos = cos(-AoA * PI_NUMBER / 180.0);
+      ValSin = sin(-AoA * PI_NUMBER / 180.0);
+
+      VarCoordTrans[0] = VarCoord[0] * ValCos - VarCoord[1] * ValSin;
+      VarCoordTrans[1] = VarCoord[1] * ValCos + VarCoord[0] * ValSin;
+
+      boundary->vertex[iMarker][iVertex]->SetVarCoord(VarCoordTrans);
+    }
+  }
+}
+
 void CSurfaceMovement::SetNACA_4Digits(CGeometry* boundary, CConfig* config) {
   unsigned long iVertex;
   unsigned short iMarker;

SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp

@@ -563,6 +563,10 @@ void CDiscAdjDeformationDriver::SetProjection_FD(CGeometry* geometry, CConfig* c
       surface_movement->SetRotation(geometry, config, iDV, true);
     }
 
+    else if (config->GetDesign_Variable(iDV) == HICKS_HENNE_CAMBER) {
+      surface_movement->SetHicksHenneCamber(geometry, config);
+    }
+
     /*--- NACA_4Digits design variable. ---*/
 
     else if (config->GetDesign_Variable(iDV) == NACA_4DIGITS) {

SU2_GEO/src/SU2_GEO.cpp

@@ -934,6 +934,17 @@ int main(int argc, char* argv[]) {
         surface_movement->SetRotation(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true);
       }
 
+      /*--- HICKS_HENNE_CAMBER design variable ---*/
+
+      else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == HICKS_HENNE_CAMBER) {
+        if (rank == MASTER_NODE) {
+          cout << endl << "Design variable number " << iDV << "." << endl;
+          cout << "Perform 2D deformation of the surface." << endl;
+        }
+        MoveSurface = true;
+        surface_movement->SetHicksHenneCamber(geometry_container[ZONE_0], config_container[ZONE_0]);
+      }
+
       /*--- NACA_4Digits design variable ---*/
 
       else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == NACA_4DIGITS) {

SU2_PY/SU2/io/tools.py

@@ -600,6 +600,7 @@ def get_dvMap():
         34: "CST",
         35: "SURFACE_BUMP",
         36: "SURFACE_FILE",
+        37: "HICKS_HENNE_CAMBER",
         40: "DV_EFIELD",
         41: "DV_YOUNG",
         42: "DV_POISSON",
@@ -722,6 +723,9 @@ def get_gradFileFormat(grad_type, plot_format, kindID, special_cases=[]):
     elif kindID == "HICKS_HENNE":
         header.append(r',"Up/Down","Loc_Max"')
         write_format.append(r", %s, %s")
+    if kindID == "HICKS_HENNE_CAMBER":
+        header.append(r',"Loc_Max"')
+        write_format.append(r", %s")
     elif kindID == "SURFACE_BUMP":
         header.append(r',"Loc_Start","Loc_End","Loc_Max"')
         write_format.append(r", %s, %s, %s")