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

@@ -182,27 +182,32 @@ class CFEAElasticity : public CNumerics {
   /*!
    * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz.
    */
-  template<class T>
+  template <class T>
   static su2double VonMisesStress(unsigned short nDim, const T& stress) {
-    if (nDim == 2) {
-      su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2];
+    su2double VM = 0.0;
 
-      su2double S1, S2; S1 = S2 = (Sxx+Syy)/2;
-      su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2));
-      S1 += tauMax;
-      S2 -= tauMax;
+    if (nDim == 2) {
+      /*--- In 2D, we expect 4 components: Sxx, Syy, Sxy, Szz ---*/
+      su2double Sxx = stress[0];
+      su2double Syy = stress[1];
+      su2double Sxy = stress[2];
+      su2double Szz = stress[3]; // <--- Critical: Input must have size 4!
 
-      return sqrt(S1*S1+S2*S2-2*S1*S2);
+      VM = sqrt( 0.5 * ( (Sxx-Syy)*(Sxx-Syy) + (Syy-Szz)*(Syy-Szz) + (Szz-Sxx)*(Szz-Sxx) + 6.0*Sxy*Sxy ) );
     }
     else {
-      su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3];
-      su2double Sxy = stress[2], Sxz = stress[4], Syz = stress[5];
-
-      return sqrt(0.5*(pow(Sxx - Syy, 2) +
-                       pow(Syy - Szz, 2) +
-                       pow(Szz - Sxx, 2) +
-                       6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz)));
+      /*--- 3D: Sxx, Syy, Szz, Sxy, Syz, Szx ---*/
+      su2double Sxx = stress[0];
+      su2double Syy = stress[1];
+      su2double Szz = stress[2];
+      su2double Sxy = stress[3];
+      su2double Syz = stress[4];
+      su2double Szx = stress[5];
+
+      VM = sqrt( 0.5 * ( (Sxx-Syy)*(Sxx-Syy) + (Syy-Szz)*(Syy-Szz) + (Szz-Sxx)*(Szz-Sxx) + 6.0*(Sxy*Sxy + Syz*Syz + Szx*Szx) ) );
     }
+
+    return VM;
   }
 
 protected:

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,25 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element,
   }
 
   if (nDim == 3) std::swap(avgStress[2], avgStress[3]);
-  auto elStress = VonMisesStress(nDim, avgStress);
+  /*--- Pack Average Stress Vector ---*/
+  su2double avgStress_VM[6] = {0.0};
+
+  if (nDim == 2) {
+      avgStress_VM[0] = avgStress[0];
+      avgStress_VM[1] = avgStress[1];
+      avgStress_VM[2] = avgStress[2];
+
+      if (isPlaneStrain) {
+          avgStress_VM[3] = Nu * (avgStress[0] + avgStress[1]);
+      } else {
+          avgStress_VM[3] = 0.0;
+      }
+  } 
+  else {
+      for (unsigned short k = 0; k < 6; k++) avgStress_VM[k] = avgStress[k];
+  }
+
+  auto elStress = CFEAElasticity::VonMisesStress(nDim, avgStress_VM);
 
   /*--- 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,26 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
 
   }
 
-  auto elStress = VonMisesStress(nDim, avgStress);
+  /*--- Pack Average Stress Vector (Handle Szz locally) ---*/
+  su2double avgStress_VM[6] = {0.0};
+
+  if (nDim == 2) {
+      avgStress_VM[0] = avgStress[0];
+      avgStress_VM[1] = avgStress[1];
+      avgStress_VM[2] = avgStress[2];
+
+      if (isPlaneStrain) {
+          // Nu is available in this scope (it was passed to the old function)
+          avgStress_VM[3] = Nu * (avgStress[0] + avgStress[1]);
+      } else {
+          avgStress_VM[3] = 0.0;
+      }
+  } 
+  else {
+      for (unsigned short k = 0; k < 6; k++) avgStress_VM[k] = avgStress[k];
+  }
+
+  auto elStress = CFEAElasticity::VonMisesStress(nDim, avgStress_VM);
 
   /*--- 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,36 @@ 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));
+      /*--- Pack Stress Vector (Handle Szz logic locally) ---*/
+      const su2double* rawStress = nodes->GetStress_FEM(iPoint);
+      su2double Stress_VM[6] = {0.0};
+
+      if (nDim == 2) {
+          Stress_VM[0] = rawStress[0]; // Sxx
+          Stress_VM[1] = rawStress[1]; // Syy
+          Stress_VM[2] = rawStress[2]; // Sxy
+
+          if (isPlaneStrain) {
+              Stress_VM[3] = Nu * (rawStress[0] + rawStress[1]);
+          } else {
+              Stress_VM[3] = 0.0;
+          }
+      } 
+      else {
+          for (unsigned short k = 0; k < 6; k++) Stress_VM[k] = rawStress[k];
+      }
+
+      const auto vms = CFEAElasticity::VonMisesStress(nDim, Stress_VM);
 
       nodes->SetVonMises_Stress(iPoint, vms);