Merge pull request #1127 from su2code/feature_ReynoldsStressInCNumerics

Introduce a function to compute stress tensor

Author: maxaehle

SU2_CFD/include/numerics/CNumerics.hpp

@@ -459,6 +459,86 @@ class CNumerics {
     TurbPsi_Grad_j = val_turbpsivar_grad_j;
   }
 
+  /*!
+   * \brief Compute the mean rate of strain matrix.
+   * \details The parameter primvargrad can be e.g. PrimVar_Grad_i or Mean_GradPrimVar.
+   * \param[in] nDim - 2 or 3
+   * \param[out] rateofstrain - Rate of strain matrix
+   * \param[in] velgrad - A velocity gradient matrix.
+   * \tparam TWOINDICES_1 - any type that supports the [][] interface
+   * \tparam TWOINDICES_2 - any type that supports the [][] interface
+   */
+  template<class TWOINDICES_1, class TWOINDICES_2>
+  inline static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, TWOINDICES_1& rateofstrain, const TWOINDICES_2& velgrad){
+
+    /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
+
+    if (nDim == 3){
+      rateofstrain[0][0] = velgrad[0][0];
+      rateofstrain[1][1] = velgrad[1][1];
+      rateofstrain[2][2] = velgrad[2][2];
+      rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+      rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]);
+      rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]);
+      rateofstrain[1][0] = rateofstrain[0][1];
+      rateofstrain[2][1] = rateofstrain[1][2];
+      rateofstrain[2][0] = rateofstrain[0][2];
+    }
+    else { // nDim==2
+      rateofstrain[0][0] = velgrad[0][0];
+      rateofstrain[1][1] = velgrad[1][1];
+      rateofstrain[2][2] = 0.0;
+      rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+      rateofstrain[0][2] = 0.0;
+      rateofstrain[1][2] = 0.0;
+      rateofstrain[1][0] = rateofstrain[0][1];
+      rateofstrain[2][1] = rateofstrain[1][2];
+      rateofstrain[2][0] = rateofstrain[0][2];
+    }
+  }
+
+  /*!
+   * \brief Compute the stress tensor from the velocity gradients.
+   * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result
+   * of this function by (-rho). The argument density is only used if turb_ke is not 0.
+   * To select the velocity gradient components from a primitive variable gradient PrimVar_Grad_i,
+   * write PrimVar_Grad_i+1.
+   * If <code>nDim==2</code>, we use the same formula but only only access the entries [0][0]..[1][1] of
+   * stress and velgrad. If <code>reynolds3x3</code> is true, the other non-diagonal entries of stress
+   * set to zero, and <code>stress[2][2]</code> to some value.
+   * \param[in] nDim - Dimension of the flow problem, 2 or 3
+   * \param[out] stress - Stress tensor
+   * \param[in] velgrad - A velocity gradient matrix.
+   * \param[in] viscosity - Viscosity
+   * \param[in] density - Density
+   * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor
+   * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2.
+   * \tparam TWOINDICES_1 - any type that supports the [][] interface
+   * \tparam TWOINDICES_2 - any type that supports the [][] interface
+   */
+  template<class TWOINDICES_1, class TWOINDICES_2>
+  inline static void ComputeStressTensor(unsigned short nDim, TWOINDICES_1& stress, const TWOINDICES_2& velgrad,
+                                      su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false){
+    su2double divVel = 0;
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
+      divVel += velgrad[iDim][iDim];
+    }
+    su2double pTerm = 2./3. * (divVel * viscosity + density * turb_ke);
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
+      for (unsigned short jDim = 0; jDim < nDim; jDim++){
+        stress[iDim][jDim] = viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]);
+      }
+      stress[iDim][iDim] -= pTerm;
+    }
+
+    if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix
+      stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0;
+      stress[2][2] = -pTerm;
+    }
+
+  }
+
   /*!
    * \brief Set the value of the first blending function.
    * \param[in] val_F1_i - Value of the first Menter blending function at point i.

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -190,12 +190,6 @@ class CAvgGrad_Base : public CNumerics {
    */
   void SetPerturbedRSM(su2double turb_ke, const CConfig* config);
 
-  /*!
-   * \brief Get the mean rate of strain matrix based on velocity gradients
-   * \param[in] S_ij
-   */
-  void GetMeanRateOfStrainMatrix(su2double **S_ij) const;
-
 public:
 
   /*!

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -339,12 +339,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
      */
   void SetPerturbedStrainMag(su2double turb_ke);
 
-  /*!
-   * \brief Get the mean rate of strain matrix based on velocity gradients
-   * \param[in] S_ij
-   */
-  void GetMeanRateOfStrainMatrix(su2double **S_ij);
-
 public:
   /*!
    * \brief Constructor of the class.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2030,12 +2030,11 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, div_vel, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
+  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
             RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, thermal_conductivity_tr, thermal_conductivity_ve = 0.0,
-            MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor,
-            delta[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+            MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
   su2double **Grad_PrimVar = nullptr, dTn, dTven;
 
@@ -2190,16 +2189,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       }
 
       /*--- Evaluate Tau ---*/
-
-      div_vel = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim];
-
-      for (iDim = 0; iDim < nDim; iDim++) {
-        for (jDim = 0; jDim < nDim; jDim++) {
-          Tau[iDim][jDim] = Viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) -
-                            TWO3 * Viscosity * div_vel * delta[iDim][jDim];
-        }
-      }
+      CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity);
 
       /*--- If necessary evaluate the QCR contribution to Tau ---*/
 

SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp

@@ -188,22 +188,15 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry
 
     su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow);
 
-    const su2double* const* GradVel = &flow_nodes->GetGradient_Primitive(Point_Flow)[1];
-
-    // Divergence of the velocity
-    su2double DivVel = 0.0;
-    for (auto iVar = 0u; iVar < nVar; iVar++) DivVel += GradVel[iVar][iVar];
-
+    su2double tau[3][3];
+    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity);
     for (auto iVar = 0u; iVar < nVar; iVar++) {
       for (auto jVar = 0u; jVar < nVar; jVar++) {
-        // Viscous stress
-        su2double delta_ij = (iVar == jVar);
-        su2double tau_ij = Viscosity*(GradVel[jVar][iVar] + GradVel[iVar][jVar] - TWO3*DivVel*delta_ij);
-
         // Viscous component in the tn vector --> Units of force (non-dimensional).
-        Donor_Variable[iVar] += tau_ij * Normal_Flow[jVar];
+        Donor_Variable[iVar] += tau[iVar][jVar] * Normal_Flow[jVar];
       }
     }
+
   }
 
   // Redimensionalize and take into account ramp transfer of the loads

SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp

@@ -237,7 +237,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   // rather than the standard V = [r1, ... , rn, T, Tve, ... ]
 
   unsigned short iSpecies, iVar, iDim, jDim;
-  su2double *Ds, *V, **GV, mu, ktr, kve, div_vel;
+  su2double *Ds, *V, **GV, mu, ktr, kve;
   su2double rho, T, Tve, RuSI, Ru;
   auto& Ms = fluidmodel->GetSpeciesMolarMass();
 
@@ -279,12 +279,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   //Cpve = V[RHOCVVE_INDEX]+Ru/Mass;
   //kve += Cpve*(val_eddy_viscosity/Prandtl_Turb);
 
-  /*--- Calculate the velocity divergence ---*/
-  div_vel = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++)
-    div_vel += GV[VEL_INDEX+iDim][iDim];
-
-
+  
   /*--- Pre-compute mixture quantities ---*/
   for (iDim = 0; iDim < nDim; iDim++) {
     Vector[iDim] = 0.0;
@@ -294,16 +289,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   }
 
   /*--- Compute the viscous stress tensor ---*/
-  for (iDim = 0; iDim < nDim; iDim++)
-    for (jDim = 0; jDim < nDim; jDim++)
-      tau[iDim][jDim] = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++) {
-    for (jDim = 0 ; jDim < nDim; jDim++) {
-      tau[iDim][jDim] += mu * (val_gradprimvar[VEL_INDEX+jDim][iDim] +
-          val_gradprimvar[VEL_INDEX+iDim][jDim]);
-    }
-    tau[iDim][iDim] -= TWO3*mu*div_vel;
-  }
+  ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX, mu);
 
   /*--- Populate entries in the viscous flux vector ---*/
   for (iDim = 0; iDim < nDim; iDim++) {

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

@@ -123,28 +123,22 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
                            const su2double val_laminar_viscosity,
                            const su2double val_eddy_viscosity) {
 
-  unsigned short iDim, jDim;
   const su2double Density = val_primvar[nDim+2];
-  const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity;
-
-  su2double div_vel = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++)
-    div_vel += val_gradprimvar[iDim+1][iDim];
 
-  /* --- If UQ methodology is used, calculate tau using the perturbed reynolds stress tensor --- */
+  /* --- If UQ methodology is used, use the perturbed Reynolds stress tensor
+   * for the turbulent part of tau. Otherwise both the laminar and turbulent
+   * parts of tau can be computed with the total viscosity. --- */
 
   if (using_uq){
-    for (iDim = 0 ; iDim < nDim; iDim++)
-      for (jDim = 0 ; jDim < nDim; jDim++)
-        tau[iDim][jDim] = val_laminar_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] )
-        - TWO3*val_laminar_viscosity*div_vel*delta[iDim][jDim] - Density * MeanPerturbedRSM[iDim][jDim];
-
+    ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity);  // laminar part
+    // add turbulent part which was perturbed
+    for (unsigned short iDim = 0 ; iDim < nDim; iDim++)
+      for (unsigned short jDim = 0 ; jDim < nDim; jDim++)
+        tau[iDim][jDim] += (-Density) * MeanPerturbedRSM[iDim][jDim];
   } else {
-
-    for (iDim = 0 ; iDim < nDim; iDim++)
-      for (jDim = 0 ; jDim < nDim; jDim++)
-        tau[iDim][jDim] = total_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] )
-                        - TWO3*total_viscosity*div_vel*delta[iDim][jDim];
+    // compute both parts in one step
+    const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity;
+    ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, 0.0); // TODO why ignore turb_ke?
   }
 }
 
@@ -217,67 +211,14 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal,
       tau[iDim][jDim] = tau[iDim][jDim]*(val_tau_wall/WallShearStress);
 }
 
-void CAvgGrad_Base::GetMeanRateOfStrainMatrix(su2double **S_ij) const
-{
-  /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
-
-  if (nDim == 3){
-    S_ij[0][0] = Mean_GradPrimVar[1][0];
-    S_ij[1][1] = Mean_GradPrimVar[2][1];
-    S_ij[2][2] = Mean_GradPrimVar[3][2];
-    S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]);
-    S_ij[0][2] = 0.5 * (Mean_GradPrimVar[1][2] + Mean_GradPrimVar[3][0]);
-    S_ij[1][2] = 0.5 * (Mean_GradPrimVar[2][2] + Mean_GradPrimVar[3][1]);
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-  }
-  else {
-    S_ij[0][0] = Mean_GradPrimVar[1][0];
-    S_ij[1][1] = Mean_GradPrimVar[2][1];
-    S_ij[2][2] = 0.0;
-    S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]);
-    S_ij[0][2] = 0.0;
-    S_ij[1][2] = 0.0;
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-
-  }
-}
-
 void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){
-
-  unsigned short iDim, jDim;
-  su2double **S_ij = new su2double* [3];
-  su2double muT = Mean_Eddy_Viscosity;
-  su2double divVel = 0;
-  su2double density;
-  su2double TWO3 = 2.0/3.0;
-  density = Mean_PrimVar[nDim+2];
-
-  for (iDim = 0; iDim < 3; iDim++){
-    S_ij[iDim] = new su2double [3];
-  }
-
-  GetMeanRateOfStrainMatrix(S_ij);
-
-  /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */
-
-  for (iDim = 0; iDim < 3; iDim++){
-    divVel += S_ij[iDim][iDim];
-  }
-
-  for (iDim = 0; iDim < 3; iDim++){
-    for (jDim = 0; jDim < 3; jDim++){
-      MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim]
-      - muT / density * (2 * S_ij[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]);
+  su2double meandensity = Mean_PrimVar[nDim+2];
+  ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar+1, Mean_Eddy_Viscosity, meandensity, turb_ke, true);
+  for(unsigned short iDim=0; iDim<3; iDim++){
+    for(unsigned short jDim=0; jDim<3; jDim++){
+      MeanReynoldsStress[iDim][jDim] /= (-meandensity);
     }
   }
-
-  for (iDim = 0; iDim < 3; iDim++)
-    delete [] S_ij[iDim];
-  delete [] S_ij;
 }
 
 void CAvgGrad_Base::SetPerturbedRSM(su2double turb_ke, const CConfig* config){

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -136,7 +136,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim
 CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) {
 
   su2double yinv, Velocity_i[3];
-  unsigned short iDim, jDim, iVar, jVar;
+  unsigned short iDim, iVar, jVar;
 
   if (Coord_i[1] > EPS) {
 
@@ -197,21 +197,12 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo
       Eddy_Viscosity_i       = V_i[nDim+5];
       Thermal_Conductivity_i = V_i[nDim+6];
 
-      su2double total_viscosity, div_vel;
+      su2double total_viscosity;
 
       total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i);
 
       /*--- The full stress tensor is needed for variable density ---*/
-
-      div_vel = 0.0;
-      for (iDim = 0 ; iDim < nDim; iDim++)
-        div_vel += PrimVar_Grad_i[iDim+1][iDim];
-
-      for (iDim = 0 ; iDim < nDim; iDim++)
-        for (jDim = 0 ; jDim < nDim; jDim++)
-          tau[iDim][jDim] = (total_viscosity*(PrimVar_Grad_i[jDim+1][iDim] +
-                                              PrimVar_Grad_i[iDim+1][jDim] )
-                             -TWO3*total_viscosity*div_vel*delta[iDim][jDim]);
+      ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity);
 
       /*--- Viscous terms. ---*/