From 925c734780a0e1d65cec3fd71f21e49d8c1f44bf Mon Sep 17 00:00:00 2001 From: Ayush Date: Sun, 28 Dec 2025 17:17:39 +0000 Subject: [PATCH 1/9] FIX: Correct 2D Von Mises formula in CFEAElasticity (#2603) --- SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 0eac022032eb..7e62b63e07db 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -192,7 +192,7 @@ class CFEAElasticity : public CNumerics { S1 += tauMax; S2 -= tauMax; - return sqrt(S1*S1+S2*S2-2*S1*S2); + return sqrt(S1*S1 + S2*S2 - S1*S2); } else { su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3]; From d406e77c06e9418c5be8ab92360a4d35545e7aee Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 10:26:29 +0000 Subject: [PATCH 2/9] FIX: Implement Plane Strain Von Mises formulation with Poisson ratio --- .../numerics/elasticity/CFEAElasticity.hpp | 19 ++++++++++++------- SU2_CFD/src/solvers/CFEASolver.cpp | 9 ++++++++- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 7e62b63e07db..ce60d5284461 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -181,18 +181,23 @@ class CFEAElasticity : public CNumerics { /*! * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz. + * \note Uses default arguments to maintain compatibility with legacy calls. */ template - static su2double VonMisesStress(unsigned short nDim, const T& stress) { + static su2double VonMisesStress(unsigned short nDim, const T& stress, su2double Nu = 0.0, bool isPlaneStrain = false) { 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 - 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 +206,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)))); } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index da425eff6995..c3a04b7c8856 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/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); From d5c4ea5d226484c7d7d6d354f6179cc75f2718b3 Mon Sep 17 00:00:00 2001 From: Ayush Kumar Date: Mon, 29 Dec 2025 21:42:17 +0530 Subject: [PATCH 3/9] Update SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index ce60d5284461..9367637293c1 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -181,7 +181,6 @@ class CFEAElasticity : public CNumerics { /*! * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz. - * \note Uses default arguments to maintain compatibility with legacy calls. */ template static su2double VonMisesStress(unsigned short nDim, const T& stress, su2double Nu = 0.0, bool isPlaneStrain = false) { From 040dcae22659c3f79ca20622f71cb046fb677b7f Mon Sep 17 00:00:00 2001 From: Ayush Kumar Date: Mon, 29 Dec 2025 21:42:30 +0530 Subject: [PATCH 4/9] Update SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 9367637293c1..d1c9986a9291 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -183,7 +183,7 @@ class CFEAElasticity : public CNumerics { * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz. */ template - static su2double VonMisesStress(unsigned short nDim, const T& stress, su2double Nu = 0.0, bool isPlaneStrain = false) { + 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]; From a86abc967e9d5b3dece8158ad0b84f0bca5703ad Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 16:38:20 +0000 Subject: [PATCH 5/9] FIX: Compute and store Szz for Plane Strain in 2D linear elasticity --- .../numerics/elasticity/CFEALinearElasticity.cpp | 14 +++++++++++--- .../elasticity/CFEANonlinearElasticity.cpp | 5 ++++- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index 381f7fc368da..fe03cbe98467 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/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. ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index 835af075d0e1..7e8d3a9d6537 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/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. ---*/ From 0a4df17160b6d351262d7224f636b7bb029c45ad Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 16:48:13 +0000 Subject: [PATCH 6/9] DOCS: Add contributor to AUTHORS.md --- AUTHORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS.md b/AUTHORS.md index 7f8afb021de6..3d7e7eb8416b 100644 --- a/AUTHORS.md +++ b/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 From f94ea51e12b310ebba5edc7dbcbfd4053bb3b270 Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 19:19:26 +0000 Subject: [PATCH 7/9] REF: Decouple physics from VonMisesStress utility --- .../numerics/elasticity/CFEAElasticity.hpp | 43 ++++++++++--------- .../elasticity/CFEALinearElasticity.cpp | 20 ++++++++- .../elasticity/CFEANonlinearElasticity.cpp | 21 ++++++++- SU2_CFD/src/solvers/CFEASolver.cpp | 21 ++++++++- externals/opdi | 2 +- subprojects/MLPCpp | 2 +- 6 files changed, 83 insertions(+), 26 deletions(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index d1c9986a9291..5da75ff06d6a 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -182,31 +182,32 @@ class CFEAElasticity : public CNumerics { /*! * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz. */ - template - static su2double VonMisesStress(unsigned short nDim, const T& stress, su2double Nu, bool isPlaneStrain) { + template + static su2double VonMisesStress(unsigned short nDim, const T& stress) { + su2double VM = 0.0; + if (nDim == 2) { - su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2]; - - /*--- In Plane Strain, Szz is not zero. It is determined by Poisson's ratio. ---*/ - su2double Szz = 0.0; - if (isPlaneStrain) { - Szz = Nu * (Sxx + Syy); - } - - /*--- 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))); + /*--- 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! + + 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*(pow(Sxy, 2) + pow(Syz, 2) + pow(Sxz, 2)))); + /*--- 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: diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index fe03cbe98467..27f433120458 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -367,7 +367,25 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, } if (nDim == 3) std::swap(avgStress[2], avgStress[3]); - auto elStress = VonMisesStress(nDim, avgStress, Nu, isPlaneStrain); + /*--- 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. ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index 7e8d3a9d6537..a18f28383c2f 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp @@ -896,7 +896,26 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen } - auto elStress = VonMisesStress(nDim, avgStress, Nu, isPlaneStrain); + /*--- 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. ---*/ diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index c3a04b7c8856..c8ecd94985e6 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1296,7 +1296,26 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, 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), Nu, isPlaneStrain); + /*--- 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); diff --git a/externals/opdi b/externals/opdi index 294807b0111c..a5e2ac47035b 160000 --- a/externals/opdi +++ b/externals/opdi @@ -1 +1 @@ -Subproject commit 294807b0111ce241cda97db62f80cdd5012d9381 +Subproject commit a5e2ac47035b6b3663f60d5f80b7a9fe62084867 diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index e19ca0cafb28..4936b4500ddc 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit e19ca0cafb28c4b7ba5b8cffef42883259b00dc0 +Subproject commit 4936b4500ddc6bd637f58ac15e8818a2c3fdd288 From f3d296b2f9c6dd8407032f0eb5c721b2481ce34f Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 19:29:07 +0000 Subject: [PATCH 8/9] CHORE: Revert accidental submodule updates --- subprojects/MLPCpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp index 4936b4500ddc..e19ca0cafb28 160000 --- a/subprojects/MLPCpp +++ b/subprojects/MLPCpp @@ -1 +1 @@ -Subproject commit 4936b4500ddc6bd637f58ac15e8818a2c3fdd288 +Subproject commit e19ca0cafb28c4b7ba5b8cffef42883259b00dc0 From 97251c8ad53646107b5e8fe6a79401271f20b2e1 Mon Sep 17 00:00:00 2001 From: Ayush Date: Mon, 29 Dec 2025 19:43:58 +0000 Subject: [PATCH 9/9] CHORE: Sync opdi with develop branch --- externals/opdi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/externals/opdi b/externals/opdi index a5e2ac47035b..294807b0111c 160000 --- a/externals/opdi +++ b/externals/opdi @@ -1 +1 @@ -Subproject commit a5e2ac47035b6b3663f60d5f80b7a9fe62084867 +Subproject commit 294807b0111ce241cda97db62f80cdd5012d9381