Fix 2D von Mises formula and thermal expansion for nonlinear material in plane stress (#2667)

* FIX: Correct 2D Von Mises formula in CFEAElasticity (#2603)

* FIX: Implement Plane Strain Von Mises formulation with Poisson ratio

* Update SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

* Update SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

* FIX: Compute and store Szz for Plane Strain in 2D linear elasticity

* DOCS: Add contributor to AUTHORS.md

* REF: Decouple physics from VonMisesStress utility

* CHORE: Revert accidental submodule updates

* CHORE: Sync opdi with develop branch

* fixes for 2D VMS and nonlinear plane stress with thermal expansion

* update regressions

* update

---------

Co-authored-by: Ayush <ayushssp74@gmail.com>

Author: pcarruscag

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

Common/src/linear_algebra/CSysSolve.cpp

@@ -626,14 +626,15 @@ unsigned long CSysSolve<ScalarType>::RFGMRES_LinSolver(const CSysVector<ScalarTy
   }
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
-  for (auto totalIter = 0ul; totalIter < MaxIter;) {
+  auto totalIter = 0ul;
+  while (totalIter < MaxIter) {
     /*--- Enforce a hard limit on total number of iterations ---*/
     auto iterLimit = min(restartIter, MaxIter - totalIter);
     auto iter = FGMRES_LinSolver(b, x, mat_vec, precond, tol, iterLimit, residual, monitoring, config);
     totalIter += iter;
-    if (residual <= tol || iter < iterLimit) return totalIter;
+    if (residual <= tol || iter < iterLimit) break;
   }
-  return 0;
+  return totalIter;
 }
 
 template <class ScalarType>

SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

@@ -80,7 +80,8 @@ class CFEAElasticity : public CNumerics {
   su2double *DV_Val = nullptr;            /*!< \brief For optimization cases, value of the design variables. */
   unsigned short n_DV = 0;                /*!< \brief For optimization cases, number of design variables. */
 
-  bool plane_stress = false;              /*!< \brief Checks if we are solving a plane stress case */
+  bool plane_stress = false;              /*!< \brief Checks if we are solving a plane stress case. */
+  bool linear = false;                    /*!< \brief Checks if we are solving a linear elasticity case. */
 
 public:
   /*!
@@ -181,28 +182,21 @@ class CFEAElasticity : public CNumerics {
 
   /*!
    * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz.
+   * \note Szz is required in 2D whereas Sxz and Syz are assumed to be 0.
    */
-  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 S1, S2; S1 = S2 = (Sxx+Syy)/2;
-      su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2));
-      S1 += tauMax;
-      S2 -= tauMax;
-
-      return sqrt(S1*S1+S2*S2-2*S1*S2);
-    }
-    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)));
+    /*--- In 2D, we only have 4 components: Sxx, Syy, Sxy, Szz. ---*/
+    const auto& Sxx = stress[0];
+    const auto& Syy = stress[1];
+    const auto& Sxy = stress[2];
+    const auto& Szz = stress[3];
+    su2double Sxz = 0, Syz = 0;
+    if (nDim == 3) {
+      Sxz = stress[4];
+      Syz = stress[5];
     }
+    return sqrt(0.5 * (pow(Sxx - Syy, 2) + pow(Syy - Szz, 2) + pow(Szz - Sxx, 2) + 6 * (Sxy * Sxy + Sxz * Sxz + Syz * Syz)));
   }
 
 protected:
@@ -233,8 +227,12 @@ class CFEAElasticity : public CNumerics {
     Mu     = E / (2.0*(1.0 + Nu));
     Lambda = Nu*E/((1.0+Nu)*(1.0-2.0*Nu));
     Kappa  = Lambda + (2/3)*Mu;
-    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php ---*/
-    ThermalStressTerm = -Alpha * E / (1 - (plane_stress ? 1 : 2) * Nu);
+    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php
+     * The stress tensor for nonlinear problems is still 3x3 and plane stress is imposed
+     * by determining the deformation that makes sigma_33 = 0, hence this denominator is
+     * only changed for linear elasticity. ---*/
+    const auto nu_mult = (linear && plane_stress) ? 1 : 2;
+    ThermalStressTerm = -Alpha * E / (1 - nu_mult * Nu);
   }
 
   /*!

SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp

@@ -134,8 +134,9 @@ class CFEANonlinearElasticity : public CFEAElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  virtual void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) = 0;
+  virtual void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) = 0;
 
   /*!
    * \brief Compute the stress tensor.

SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp

@@ -58,8 +58,9 @@ class CFEM_NeoHookean_Comp final : public CFEANonlinearElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) override;
+  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
   /*!
    * \brief Compute the constitutive matrix.
@@ -72,9 +73,9 @@ class CFEM_NeoHookean_Comp final : public CFEANonlinearElasticity {
    * \brief Compute the stress tensor.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
   void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
-
 };
 
 
@@ -109,8 +110,9 @@ class CFEM_Knowles_NearInc final : public CFEANonlinearElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) override;
+  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
   /*!
    * \brief Compute the constitutive matrix.
@@ -123,6 +125,7 @@ class CFEM_Knowles_NearInc final : public CFEANonlinearElasticity {
    * \brief Compute the stress tensor.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
   void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
@@ -157,8 +160,9 @@ class CFEM_DielectricElastomer final : public CFEANonlinearElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  inline void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) override { };
+  inline void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) override { };
 
   /*!
    * \brief Compute the constitutive matrix.
@@ -171,6 +175,7 @@ class CFEM_DielectricElastomer final : public CFEANonlinearElasticity {
    * \brief Compute the stress tensor.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
   void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
@@ -207,8 +212,9 @@ class CFEM_IdealDE final : public CFEANonlinearElasticity {
    * \brief Compute the plane stress term.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config) override;
+  void Compute_Plane_Stress_Term(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
   /*!
    * \brief Compute the constitutive matrix.
@@ -221,6 +227,7 @@ class CFEM_IdealDE final : public CFEANonlinearElasticity {
    * \brief Compute the stress tensor.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
   void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 

SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp

@@ -64,9 +64,10 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   Alpha = Alpha_i[0];
   ReferenceTemperature = config->GetMaterialReferenceTemperature();
 
-  Compute_Lame_Parameters();
-
   plane_stress = (nDim == 2) && (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRESS);
+  linear = config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL;
+
+  Compute_Lame_Parameters();
 
   KAux_ab = new su2double* [nDim];
   for (iVar = 0; iVar < nDim; iVar++) {

SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp

@@ -206,26 +206,18 @@ void CFEALinearElasticity::Compute_Constitutive_Matrix(CElement *element_contain
 
   /*--- Compute the D Matrix (for plane stress and 2-D)---*/
 
-
   if (nDim == 2) {
     if (plane_stress) {
-      /*--- We enable plane stress cases ---*/
-
-      D_Mat[0][0] = E/(1-Nu*Nu);        D_Mat[0][1] = (E*Nu)/(1-Nu*Nu);  D_Mat[0][2] = 0.0;
-      D_Mat[1][0] = (E*Nu)/(1-Nu*Nu);   D_Mat[1][1] = E/(1-Nu*Nu);       D_Mat[1][2] = 0.0;
-      D_Mat[2][0] = 0.0;                D_Mat[2][1] = 0.0;               D_Mat[2][2] = ((1-Nu)*E)/(2*(1-Nu*Nu));
-    }
-    else {
-      /*--- Assuming plane strain as a general case ---*/
-
+      su2double D = E / (1 - Nu * Nu);
+      D_Mat[0][0] = D;        D_Mat[0][1] = Nu * D;  D_Mat[0][2] = 0.0;
+      D_Mat[1][0] = Nu * D;   D_Mat[1][1] = D;       D_Mat[1][2] = 0.0;
+      D_Mat[2][0] = 0.0;      D_Mat[2][1] = 0.0;     D_Mat[2][2] = (1 - Nu) * D / 2;
+    } else {
       D_Mat[0][0] = Lambda + 2.0*Mu;  D_Mat[0][1] = Lambda;           D_Mat[0][2] = 0.0;
       D_Mat[1][0] = Lambda;           D_Mat[1][1] = Lambda + 2.0*Mu;  D_Mat[1][2] = 0.0;
       D_Mat[2][0] = 0.0;              D_Mat[2][1] = 0.0;              D_Mat[2][2] = Mu;
     }
-
-  }
-  else {
-
+  } else {
     su2double Lbda_2Mu = Lambda + 2.0*Mu;
 
     D_Mat[0][0] = Lbda_2Mu;  D_Mat[0][1] = Lambda;    D_Mat[0][2] = Lambda;    D_Mat[0][3] = 0.0;  D_Mat[0][4] = 0.0;  D_Mat[0][5] = 0.0;
@@ -234,7 +226,6 @@ void CFEALinearElasticity::Compute_Constitutive_Matrix(CElement *element_contain
     D_Mat[3][0] = 0.0;       D_Mat[3][1] = 0.0;       D_Mat[3][2] = 0.0;       D_Mat[3][3] = Mu;   D_Mat[3][4] = 0.0;  D_Mat[3][5] = 0.0;
     D_Mat[4][0] = 0.0;       D_Mat[4][1] = 0.0;       D_Mat[4][2] = 0.0;       D_Mat[4][3] = 0.0;  D_Mat[4][4] = Mu;   D_Mat[4][5] = 0.0;
     D_Mat[5][0] = 0.0;       D_Mat[5][1] = 0.0;       D_Mat[5][2] = 0.0;       D_Mat[5][3] = 0.0;  D_Mat[5][4] = 0.0;  D_Mat[5][5] = Mu;
-
   }
 
 }
@@ -336,30 +327,27 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element,
       avgStress[iVar] += Stress[iVar] / nGauss;
     }
 
-    for (iNode = 0; iNode < nNode; iNode++) {
+    if (nDim == 2 && config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRAIN) {
+      Stress[3] = Nu * (Stress[0] + Stress[1]) + (1 - 2 * Nu) * thermalStress;
+      avgStress[3] += Stress[3] / nGauss;
+    }
 
-      su2double Ni_Extrap = element->GetNi_Extrap(iNode, iGauss);
+    /*--- 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. ---*/
+    if (nDim == 3) std::swap(Stress[2], Stress[3]);
 
-      if (nDim == 2) {
-        for(iVar = 0; iVar < 3; ++iVar)
-          element->Add_NodalStress(iNode, iVar, Stress[iVar] * 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. ---*/
-        element->Add_NodalStress(iNode, 0, Stress[0] * Ni_Extrap);
-        element->Add_NodalStress(iNode, 1, Stress[1] * Ni_Extrap);
-        element->Add_NodalStress(iNode, 2, Stress[3] * Ni_Extrap);
-        element->Add_NodalStress(iNode, 3, Stress[2] * Ni_Extrap);
-        element->Add_NodalStress(iNode, 4, Stress[4] * Ni_Extrap);
-        element->Add_NodalStress(iNode, 5, Stress[5] * Ni_Extrap);
+    for (iNode = 0; iNode < nNode; iNode++) {
+      su2double Ni_Extrap = element->GetNi_Extrap(iNode, iGauss);
+      for (iVar = 0; iVar < DIM_STRAIN_3D; ++iVar) {
+        element->Add_NodalStress(iNode, iVar, Stress[iVar] * Ni_Extrap);
       }
     }
 
   }
 
+  /*--- See note regarding output order above. ---*/
   if (nDim == 3) std::swap(avgStress[2], avgStress[3]);
-  auto elStress = VonMisesStress(nDim, avgStress);
+  auto elStress = CFEAElasticity::VonMisesStress(nDim, avgStress);
 
   /*--- 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

@@ -320,7 +320,7 @@ void CFEANonlinearElasticity::Compute_Tangent_Matrix(CElement *element, const CC
     if (nDim == 2) {
       if (plane_stress) {
         // Compute the value of the term 33 for the deformation gradient
-        Compute_Plane_Stress_Term(element, config);
+        Compute_Plane_Stress_Term(element, config, iGauss);
         F_Mat[2][2] = f33;
       }
       else { // plane strain
@@ -542,7 +542,7 @@ void CFEANonlinearElasticity::Compute_NodalStress_Term(CElement *element, const
     if (nDim == 2) {
       if (plane_stress) {
         // Compute the value of the term 33 for the deformation gradient
-        Compute_Plane_Stress_Term(element, config);
+        Compute_Plane_Stress_Term(element, config, iGauss);
         F_Mat[2][2] = f33;
       }
       else { // plane strain
@@ -821,10 +821,10 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
     if (nDim == 2) {
       if (plane_stress) {
         // Compute the value of the term 33 for the deformation gradient
-        Compute_Plane_Stress_Term(element, config);
+        Compute_Plane_Stress_Term(element, config, iGauss);
         F_Mat[2][2] = f33;
-      }
-      else { // plane strain
+      } else {
+        // Plane strain
         F_Mat[2][2] = 1.0;
       }
     }
@@ -856,8 +856,8 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
     avgStress[0] += Stress_Tensor[0][0] / nGauss;
     avgStress[1] += Stress_Tensor[1][1] / nGauss;
     avgStress[2] += Stress_Tensor[0][1] / nGauss;
+    avgStress[3] += Stress_Tensor[2][2] / nGauss;
     if (nDim == 3) {
-      avgStress[3] += Stress_Tensor[2][2] / nGauss;
       avgStress[4] += Stress_Tensor[0][2] / nGauss;
       avgStress[5] += Stress_Tensor[1][2] / nGauss;
     }
@@ -883,8 +883,8 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
       element->Add_NodalStress(iNode, 0, Stress_Tensor[0][0] * Ni_Extrap);
       element->Add_NodalStress(iNode, 1, Stress_Tensor[1][1] * Ni_Extrap);
       element->Add_NodalStress(iNode, 2, Stress_Tensor[0][1] * Ni_Extrap);
+      element->Add_NodalStress(iNode, 3, Stress_Tensor[2][2] * Ni_Extrap);
       if (nDim == 3) {
-        element->Add_NodalStress(iNode, 3, Stress_Tensor[2][2] * Ni_Extrap);
         element->Add_NodalStress(iNode, 4, Stress_Tensor[0][2] * Ni_Extrap);
         element->Add_NodalStress(iNode, 5, Stress_Tensor[1][2] * Ni_Extrap);
       }
@@ -893,7 +893,7 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
 
   }
 
-  auto elStress = VonMisesStress(nDim, avgStress);
+  auto elStress = CFEAElasticity::VonMisesStress(nDim, avgStress);
 
   /*--- We only differentiate w.r.t. an avg VM stress for the element as
    * considering all nodal stresses would use too much memory. ---*/