Unsteady incompressible flow simulation with grid movement

Common/src/config_structure.cpp

@@ -4299,10 +4299,10 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
     }
   }
 
-  /*--- Grid motion is not yet supported with the incompressible solver. ---*/
+  /*--- Rotating frame is not yet supported with the incompressible solver. ---*/
 
-  if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (GetGrid_Movement())) {
-    SU2_MPI::Error("Support for grid movement not yet implemented for incompressible flows.", CURRENT_FUNCTION);
+  if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (Kind_GridMovement == ROTATING_FRAME)) {
+    SU2_MPI::Error("Support for rotating frame simulation not yet implemented for incompressible flows.", CURRENT_FUNCTION);
   }
 
   /*--- Assert that there are two markers being analyzed if the

SU2_CFD/include/numerics_structure.hpp

@@ -5263,6 +5263,42 @@ class CSourceIncBodyForce : public CNumerics {
 
 };
 
+/*!
+ * \class CSourceIncRotatingFrame_Flow
+ * \brief Class for a rotating frame source term.
+ * \ingroup SourceDiscr
+ */
+class CSourceIncRotatingFrame_Flow : public CNumerics {
+
+private:
+  su2double Omega[3];  /*!< \brief Angular velocity */
+  bool implicit; /*!< \brief Implicit calculation. */
+
+public:
+
+  /*!
+   * \brief Constructor of the class.
+   * \param[in] val_nDim - Number of dimensions of the problem.
+   * \param[in] val_nVar - Number of variables of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  CSourceIncRotatingFrame_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config);
+
+  /*!
+   * \brief Destructor of the class.
+   */
+  ~CSourceIncRotatingFrame_Flow(void);
+
+  /*!
+   * \brief Residual of the rotational frame source term.
+   * \param[out] val_residual - Pointer to the total residual.
+   * \param[out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation).
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config);
+
+};
+
 /*!
  * \class CSourceBoussinesq
  * \brief Class for the source term integration of the Boussinesq approximation for incompressible flow.

SU2_CFD/include/solver_structure.hpp

@@ -567,6 +567,7 @@ class CSolver {
   /*!
    * \brief Load the geometries at the previous time states n and nM1.
    * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
    */
   virtual void Restart_OldGeometry(CGeometry *geometry, CConfig *config);
 

SU2_CFD/src/drivers/CDriver.cpp

@@ -2205,7 +2205,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol
       else if (incompressible && (config->GetKind_DensityModel() == BOUSSINESQ))
         numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceBoussinesq(nDim, nVar_Flow, config);
       else if (config->GetRotating_Frame() == YES)
-        numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config);
+        if (incompressible) numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncRotatingFrame_Flow(nDim, nVar_Flow, config);
+        else numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config);
       else if (config->GetAxisymmetric() == YES)
         if (incompressible) numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config);
       else numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config);
@@ -3831,7 +3832,7 @@ void CDriver::StartSolver(){
 
     /*--- Perform a dynamic mesh update if required. ---*/
 
-      if (!fem_solver) {
+      if (!fem_solver && !(config_container[ZONE_0]->GetGrid_Movement() && config_container[ZONE_0]->GetDiscrete_Adjoint())) {
         DynamicMeshUpdate(ExtIter);
       }
 

SU2_CFD/src/numerics_direct_mean_inc.cpp

@@ -86,8 +86,15 @@ CUpwFDSInc_Flow::~CUpwFDSInc_Flow(void) {
 
 void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) {
 
+  su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0};
+  su2double ProjGridVel = 0.0;
+
   AD::StartPreacc();
   AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); AD::SetPreaccIn(Normal, nDim);
+  if (dynamic_grid) {
+    AD::SetPreaccIn(GridVel_i, nDim);
+    AD::SetPreaccIn(GridVel_j, nDim);
+  }
 
   /*--- Face area (norm or the normal vector) ---*/
 
@@ -120,6 +127,16 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J
     ProjVelocity       += MeanVelocity[iDim]*Normal[iDim];
   }
 
+  /*--- Projected velocity adjustment due to mesh motion ---*/
+
+  if (dynamic_grid) { 
+    ProjGridVel = 0.0; 
+    for (iDim = 0; iDim < nDim; iDim++) { 
+      ProjGridVel   += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; 
+    } 
+    ProjVelocity   -= ProjGridVel; 
+  }
+
   /*--- Mean variables at points iPoint and jPoint ---*/
 
   MeanDensity     = 0.5*(DensityInc_i  + DensityInc_j);
@@ -214,6 +231,38 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J
     }
   }
 
+  /*--- Corrections due to grid motion ---*/
+  if (dynamic_grid) {
+
+    /*--- Recompute conservative variables ---*/
+
+    U_i[0] = DensityInc_i; U_j[0] = DensityInc_j;
+    for (iDim = 0; iDim < nDim; iDim++) {
+      U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim];
+    }
+    U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j;
+
+    ProjVelocity = 0.0;
+    for (iDim = 0; iDim < nDim; iDim++)
+      ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
+
+    /*--- Residual contributions ---*/
+    for (iVar = 0; iVar < nVar; iVar++) {
+      val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
+
+      /*--- Jacobian contributions ---*/
+      /*--- Implicit terms ---*/
+      if (implicit) {
+        for (iDim = 0; iDim < nDim; iDim++){
+          val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
+          val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+        }
+        val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i;
+        val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j;
+      }
+    }
+  }
+
   if (!energy) {
     val_residual[nDim+1] = 0.0;
     if (implicit) {
@@ -277,6 +326,9 @@ CCentJSTInc_Flow::~CCentJSTInc_Flow(void) {
 
 void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) {
 
+  su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0};
+  su2double ProjGridVel = 0.0;
+
   /*--- Primitive variables at point i and j ---*/
 
   Pressure_i    = V_i[0];             Pressure_j    = V_j[0];
@@ -336,7 +388,39 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
       }
     }
   }
-
+
+  /*--- Corrections due to grid motion ---*/
+  if (dynamic_grid) {
+
+    /*--- Recompute conservative variables ---*/
+
+    U_i[0] = DensityInc_i; U_j[0] = DensityInc_j;
+    for (iDim = 0; iDim < nDim; iDim++) {
+      U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim];
+    }
+    U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j;
+
+    su2double ProjVelocity = 0.0;
+    for (iDim = 0; iDim < nDim; iDim++)
+      ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
+
+    /*--- Residual contributions ---*/
+    for (iVar = 0; iVar < nVar; iVar++) {
+      val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
+
+      /*--- Jacobian contributions ---*/
+      /*--- Implicit terms ---*/
+      if (implicit) {
+        for (iDim = 0; iDim < nDim; iDim++){
+          val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
+          val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+        }
+        val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i;
+        val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j;
+      }
+    }
+  }
+
   /*--- Computes differences between Laplacians and conservative variables ---*/
 
   for (iVar = 0; iVar < nVar; iVar++) {
@@ -351,6 +435,17 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
   /*--- Compute the local spectral radius of the preconditioned system
    and the stretching factor. ---*/
 
+  /*--- Projected velocity adjustment due to mesh motion ---*/
+
+  if (dynamic_grid) {
+    ProjGridVel = 0.0;
+    for (iDim = 0; iDim < nDim; iDim++) {
+      ProjGridVel   += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
+    }
+    ProjVelocity_i -= ProjGridVel;
+    ProjVelocity_j -= ProjGridVel;
+  }
+
   SoundSpeed_i = sqrt(BetaInc2_i*Area*Area);
   SoundSpeed_j = sqrt(BetaInc2_j*Area*Area);
 
@@ -396,7 +491,6 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
       }
     }
   }
-
 }
 
 CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) {
@@ -442,6 +536,9 @@ CCentLaxInc_Flow::~CCentLaxInc_Flow(void) {
 
 void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) {
 
+  su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0};
+  su2double ProjGridVel = 0.0, ProjVelocity = 0.0;
+
   /*--- Primitive variables at point i and j ---*/
 
   Pressure_i    = V_i[0];             Pressure_j    = V_j[0];
@@ -503,7 +600,38 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
       }
     }
   }
-
+
+  /*--- Corrections due to grid motion ---*/
+  if (dynamic_grid) {
+
+    /*--- Recompute conservative variables ---*/
+
+    U_i[0] = DensityInc_i; U_j[0] = DensityInc_j;
+    for (iDim = 0; iDim < nDim; iDim++) {
+      U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim];
+    }
+    U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j;
+
+    for (iDim = 0; iDim < nDim; iDim++)
+      ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
+
+    /*--- Residual contributions ---*/
+    for (iVar = 0; iVar < nVar; iVar++) {
+      val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
+
+      /*--- Jacobian contributions ---*/
+      /*--- Implicit terms ---*/
+      if (implicit) {
+        for (iDim = 0; iDim < nDim; iDim++){
+          val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
+          val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
+        }
+        val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i;
+        val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j;
+      }
+    }
+  }
+
   /*--- Computes differences btw. conservative variables ---*/
 
   for (iVar = 0; iVar < nVar; iVar++)
@@ -519,6 +647,17 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
   SoundSpeed_i = sqrt(BetaInc2_i*Area*Area);
   SoundSpeed_j = sqrt(BetaInc2_j*Area*Area);
 
+  /*--- Projected velocity adjustment due to mesh motion ---*/
+
+  if (dynamic_grid) {
+    ProjGridVel = 0.0;
+    for (iDim = 0; iDim < nDim; iDim++) {
+      ProjGridVel   += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim];
+    }
+    ProjVelocity_i -= ProjGridVel;
+    ProjVelocity_j -= ProjGridVel;
+  }
+
   Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i;
   Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j;
 
@@ -558,7 +697,6 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_
       }
     }
   }
-
 }
 
 CAvgGradInc_Flow::CAvgGradInc_Flow(unsigned short val_nDim,
@@ -852,6 +990,72 @@ void CSourceIncBodyForce::ComputeResidual(su2double *val_residual, CConfig *conf
 
 }
 
+CSourceIncRotatingFrame_Flow::CSourceIncRotatingFrame_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) {
+
+  implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
+
+  Gamma = config->GetGamma();
+  Gamma_Minus_One = Gamma - 1.0;
+
+  /*--- Retrieve the angular velocity vector from config. ---*/
+  for (unsigned short iDim = 0; iDim < 3; iDim++)
+    Omega[iDim] = config->GetRotation_Rate(iDim)/config->GetOmega_Ref();
+
+}
+
+CSourceIncRotatingFrame_Flow::~CSourceIncRotatingFrame_Flow(void) { }
+
+void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) {
+
+  unsigned short iDim, iVar, jVar;
+  su2double Momentum[3] = {0,0,0},
+            Velocity_i[3] = {0,0,0};
+
+  /*--- Primitive variables plus momentum at the node (point i) ---*/
+
+  DensityInc_i  = V_i[nDim+2];
+
+  for (iDim = 0; iDim < nDim; iDim++) {
+    Velocity_i[iDim] = V_i[iDim+1];
+    Momentum[iDim]   = DensityInc_i*Velocity_i[iDim];
+  }
+
+  /*--- Calculate rotating frame source term residual as ( Omega X Rho-U ) ---*/
+
+  if (nDim == 2) {
+    val_residual[0] = 0.0;
+    val_residual[1] = (Omega[1]*Momentum[2] - Omega[2]*Momentum[1])*Volume;
+    val_residual[2] = (Omega[2]*Momentum[0] - Omega[0]*Momentum[2])*Volume;
+    val_residual[3] = 0.0;
+  } else { 
+    val_residual[0] = 0.0;
+    val_residual[1] = (Omega[1]*Momentum[2] - Omega[2]*Momentum[1])*Volume;
+    val_residual[2] = (Omega[2]*Momentum[0] - Omega[0]*Momentum[2])*Volume;
+    val_residual[3] = (Omega[0]*Momentum[1] - Omega[1]*Momentum[0])*Volume;
+    val_residual[4] = 0.0;
+  } 
+
+  /*--- Calculate the source term Jacobian ---*/
+
+  if (implicit) {
+    for (iVar = 0; iVar < nVar; iVar++)
+      for (jVar = 0; jVar < nVar; jVar++)
+        val_Jacobian_i[iVar][jVar] = 0.0;
+    if (nDim == 2) {
+      val_Jacobian_i[1][2] = -DensityInc_i*Omega[2]*Volume;
+      val_Jacobian_i[2][1] =  DensityInc_i*Omega[2]*Volume;
+    } else {
+      val_Jacobian_i[1][2] = -DensityInc_i*Omega[2]*Volume;
+      val_Jacobian_i[1][3] =  DensityInc_i*Omega[1]*Volume;
+      val_Jacobian_i[2][1] =  DensityInc_i*Omega[2]*Volume;
+      val_Jacobian_i[2][3] = -DensityInc_i*Omega[0]*Volume;
+      val_Jacobian_i[3][1] = -DensityInc_i*Omega[1]*Volume;
+      val_Jacobian_i[3][2] =  DensityInc_i*Omega[0]*Volume;
+    }
+  }
+
+}
+
 CSourceBoussinesq::CSourceBoussinesq(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) {
 
   /*--- Store the pointer to the constant body force vector. ---*/