FIX: Correct 2D Von Mises Stress Calculation in CFEAElasticity (#2603)

AUTHORS.md

@@ -56,6 +56,7 @@ Aniket C. Aranake
 Antonio Rubino
 Arne Bachmann
 Arne Voß
+Ayush Kumar
 Beckett Y. Zhou
 Benjamin S. Kirk
 Brendan Tracey

SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

@@ -183,16 +183,20 @@ class CFEAElasticity : public CNumerics {
    * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz.
    */
   template<class T>
-  static su2double VonMisesStress(unsigned short nDim, const T& stress) {
+  static su2double VonMisesStress(unsigned short nDim, const T& stress, su2double Nu, bool isPlaneStrain) {
     if (nDim == 2) {
       su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2];
 
-      su2double S1, S2; S1 = S2 = (Sxx+Syy)/2;
-      su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2));
-      S1 += tauMax;
-      S2 -= tauMax;
+      /*--- In Plane Strain, Szz is not zero. It is determined by Poisson's ratio. ---*/
+      su2double Szz = 0.0;
+      if (isPlaneStrain) {
+        Szz = Nu * (Sxx + Syy);
+      }
 
-      return sqrt(S1*S1+S2*S2-2*S1*S2);
+      /*--- General 3D Von Mises formula reduced to 2D components + Szz ---*/
+      /*--- Sigma_vm = sqrt( 0.5 * [ (Sxx-Syy)^2 + (Syy-Szz)^2 + (Szz-Sxx)^2 + 6*Sxy^2 ] ) ---*/
+
+      return sqrt(0.5 * (pow(Sxx - Syy, 2) + pow(Syy - Szz, 2) + pow(Szz - Sxx, 2) + 6.0 * pow(Sxy, 2)));
     }
     else {
       su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3];
@@ -201,7 +205,7 @@ class CFEAElasticity : public CNumerics {
       return sqrt(0.5*(pow(Sxx - Syy, 2) +
                        pow(Syy - Szz, 2) +
                        pow(Szz - Sxx, 2) +
-                       6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz)));
+                       6.0*(pow(Sxy, 2) + pow(Syz, 2) + pow(Sxz, 2))));
     }
   }
 

SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp

@@ -242,6 +242,9 @@ void CFEALinearElasticity::Compute_Constitutive_Matrix(CElement *element_contain
 
 su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) {
 
+  su2double Nu = config->GetPoissonRatio(0);
+  bool isPlaneStrain = (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRAIN);
+
   unsigned short iVar, jVar;
   unsigned short iGauss, nGauss;
   unsigned short iNode, nNode;
@@ -341,9 +344,14 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element,
       su2double Ni_Extrap = element->GetNi_Extrap(iNode, iGauss);
 
       if (nDim == 2) {
-        for(iVar = 0; iVar < 3; ++iVar)
-          element->Add_NodalStress(iNode, iVar, Stress[iVar] * Ni_Extrap);
+      for(iVar = 0; iVar < 3; ++iVar)
+        element->Add_NodalStress(iNode, iVar, Stress[iVar] * Ni_Extrap);
+
+      if (isPlaneStrain) {
+        su2double Szz = Nu * (Stress[0] + Stress[1]);
+        element->Add_NodalStress(iNode, 3, Szz * Ni_Extrap);
       }
+    }
       else {
         /*--- If nDim is 3 and we compute it this way, the 3rd component is the Szz,
          *    while in the output it is the 4th component for practical reasons. ---*/
@@ -359,7 +367,7 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element,
   }
 
   if (nDim == 3) std::swap(avgStress[2], avgStress[3]);
-  auto elStress = VonMisesStress(nDim, avgStress);
+  auto elStress = VonMisesStress(nDim, avgStress, Nu, isPlaneStrain);
 
   /*--- We only differentiate w.r.t. an avg VM stress for the element as
    * considering all nodal stresses would use too much memory. ---*/

SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp

@@ -746,6 +746,9 @@ void CFEANonlinearElasticity::Assign_cijkl_D_Mat() {
 
 su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) {
 
+  su2double Nu = config->GetPoissonRatio(0);
+  bool isPlaneStrain = (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRAIN);
+
   unsigned short iVar, jVar, kVar;
   unsigned short iGauss, nGauss;
   unsigned short iDim, iNode, nNode;
@@ -893,7 +896,7 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
 
   }
 
-  auto elStress = VonMisesStress(nDim, avgStress);
+  auto elStress = VonMisesStress(nDim, avgStress, Nu, isPlaneStrain);
 
   /*--- We only differentiate w.r.t. an avg VM stress for the element as
    * considering all nodal stresses would use too much memory. ---*/

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -1286,10 +1286,17 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
     /*--- Compute the von Misses stress at each point, and the maximum for the domain. ---*/
     su2double maxVonMises = 0.0;
 
+    /*--- Pre-fetch configuration for 2D Plane Strain check ---*/
+    bool isPlaneStrain = (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRAIN);
+
+    /*--- Note: For multi-zone/material problems, Nu should vary per point.
+     * Here we use the first material's Poisson ratio as the reference. ---*/
+    su2double Nu = config->GetPoissonRatio(0);
+
     SU2_OMP_FOR_(schedule(static,omp_chunk_size) SU2_NOWAIT)
     for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
 
-      const auto vms = CFEAElasticity::VonMisesStress(nDim, nodes->GetStress_FEM(iPoint));
+      const auto vms = CFEAElasticity::VonMisesStress(nDim, nodes->GetStress_FEM(iPoint), Nu, isPlaneStrain);
 
       nodes->SetVonMises_Stress(iPoint, vms);