From 925c734780a0e1d65cec3fd71f21e49d8c1f44bf Mon Sep 17 00:00:00 2001 From: Ayush Date: Sun, 28 Dec 2025 17:17:39 +0000 Subject: [PATCH 01/12] 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 02/12] 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 03/12] 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 04/12] 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 05/12] 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 06/12] 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 07/12] 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 08/12] 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 09/12] 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 From 4bfee6d249f4f68a8c6c0dc607bf7c34c4880c84 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 31 Dec 2025 11:57:20 -0800 Subject: [PATCH 10/12] fixes for 2D VMS and nonlinear plane stress with thermal expansion --- Common/src/linear_algebra/CSysSolve.cpp | 7 +- .../numerics/elasticity/CFEAElasticity.hpp | 45 +++++------ .../elasticity/CFEANonlinearElasticity.hpp | 3 +- .../numerics/elasticity/nonlinear_models.hpp | 17 +++-- .../numerics/elasticity/CFEAElasticity.cpp | 5 +- .../elasticity/CFEALinearElasticity.cpp | 74 +++++-------------- .../elasticity/CFEANonlinearElasticity.cpp | 38 ++-------- .../numerics/elasticity/nonlinear_models.cpp | 47 ++++++------ SU2_CFD/src/output/CElasticityOutput.cpp | 29 ++++---- SU2_CFD/src/solvers/CFEASolver.cpp | 30 +------- SU2_CFD/src/variables/CFEAVariable.cpp | 5 +- .../VonMissesVerif/linear_plane_strain_2d.cfg | 45 +++++++++++ .../nonlinear_plane_stress_2d.cfg | 54 ++++++++++++++ TestCases/parallel_regression.py | 19 ++++- 14 files changed, 222 insertions(+), 196 deletions(-) create mode 100644 TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg create mode 100644 TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index c13df2d36a52..a444556f8f1a 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -626,14 +626,15 @@ unsigned long CSysSolve::RFGMRES_LinSolver(const CSysVector diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 5da75ff06d6a..2a7ea4e36381 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/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,33 +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 static su2double VonMisesStress(unsigned short nDim, const T& stress) { - su2double VM = 0.0; - - 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! - - VM = sqrt( 0.5 * ( (Sxx-Syy)*(Sxx-Syy) + (Syy-Szz)*(Syy-Szz) + (Szz-Sxx)*(Szz-Sxx) + 6.0*Sxy*Sxy ) ); + /*--- 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]; } - else { - /*--- 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; + return sqrt(0.5 * (pow(Sxx - Syy, 2) + pow(Syy - Szz, 2) + pow(Szz - Sxx, 2) + 6 * (Sxy * Sxy + Sxz * Sxz + Syz * Syz))); } protected: @@ -238,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); } /*! diff --git a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp index d3fc1d796679..51f8801a9779 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp +++ b/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. diff --git a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp b/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp index 97a332a88831..9c82a99c3623 100644 --- a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp +++ b/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; diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 2b82f0ca285f..6c9bf0a57934 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/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++) { diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index 27f433120458..5635a9329ee1 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/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; - } } @@ -242,9 +233,6 @@ 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; @@ -339,53 +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) + 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); - - 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. ---*/ - 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); } } } + /*--- See note regarding output order above. ---*/ if (nDim == 3) std::swap(avgStress[2], avgStress[3]); - /*--- 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); + 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. ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index a18f28383c2f..f7df2516c4cf 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/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 @@ -746,9 +746,6 @@ 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; @@ -824,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; } } @@ -859,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; } @@ -886,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); } @@ -896,26 +893,7 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen } - /*--- 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); + 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. ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp index f7d99ddb058b..9623e0a9867d 100644 --- a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp +++ b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp @@ -34,38 +34,33 @@ CFEM_NeoHookean_Comp::CFEM_NeoHookean_Comp(unsigned short val_nDim, CFEANonlinearElasticity(val_nDim, val_nVar, config) { } -void CFEM_NeoHookean_Comp::Compute_Plane_Stress_Term(CElement *element, const CConfig *config) { +void CFEM_NeoHookean_Comp::Compute_Plane_Stress_Term(CElement *element, const CConfig *config, unsigned short iGauss) { - su2double j_red = 1.0; - su2double fx = 0.0, fpx = 1.0; - su2double xkm1 = 1.0, xk = 1.0; - su2double cte = 0.0; - - unsigned short iNR, nNR; - su2double NRTOL; + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); - // Maximum number of iterations and tolerance (relative) - nNR = 10; - NRTOL = 1E-25; + // Tolerate for NR method. + const su2double NRTOL = 1E-14; - // j_red: reduced jacobian, for the 2x2 submatrix of F - j_red = F_Mat[0][0] * F_Mat[1][1] - F_Mat[1][0] * F_Mat[0][1]; + // j_red: reduced jacobian, for the 2x2 submatrix of F. + const su2double j_red = F_Mat[0][0] * F_Mat[1][1] - F_Mat[1][0] * F_Mat[0][1]; // cte: constant term in the NR method - cte = Lambda*log(j_red) - Mu; + const su2double cte = Lambda*log(j_red) - Mu + thermalStress; - // f(f33) = mu*f33^2 + lambda*ln(f33) + (lambda*ln(j_red)-mu) = 0 + // f(f33) = mu*f33^2 + lambda*ln(f33) + (lambda*ln(j_red)-mu+thermalStress) = 0 // f'(f33) = 2*mu*f33 + lambda/f33 - for (iNR = 0; iNR < nNR; iNR++) { - fx = Mu*pow(xk,2.0) + Lambda*log(xk) + cte; - fpx = 2*Mu*xk + (Lambda / xk); - xkm1 = xk - fx / fpx; - if (((xkm1 - xk) / xk) < NRTOL) break; - xk = xkm1; + // Initialize as plane strain. + su2double xk = 1.0; + for (auto iNR = 0; iNR < 20; ++iNR) { + const su2double fx = Mu * pow(xk, 2) + Lambda * log(xk) + cte; + const su2double fpx = 2 * Mu * xk + Lambda / xk; + f33 = xk - fx / fpx; + if (fabs(f33 - xk) < NRTOL) break; + xk = f33; + } + if (fabs(f33 - xk) >= NRTOL) { + std::cout << "WARNING: Plane Stress term did not converge (residual = " << fabs(f33 - xk) << ")\n"; } - - f33 = xkm1; - } void CFEM_NeoHookean_Comp::Compute_Constitutive_Matrix(CElement *element, const CConfig *config) { @@ -140,7 +135,7 @@ CFEM_Knowles_NearInc::CFEM_Knowles_NearInc(unsigned short val_nDim, unsigned sho } -void CFEM_Knowles_NearInc::Compute_Plane_Stress_Term(CElement *element, const CConfig *config) { +void CFEM_Knowles_NearInc::Compute_Plane_Stress_Term(CElement *element, const CConfig *config, unsigned short iGauss) { SU2_MPI::Error("This material model cannot (yet) be used for plane stress.",CURRENT_FUNCTION); @@ -286,7 +281,7 @@ CFEM_IdealDE::CFEM_IdealDE(unsigned short val_nDim, unsigned short val_nVar, } -void CFEM_IdealDE::Compute_Plane_Stress_Term(CElement *element, const CConfig *config) { +void CFEM_IdealDE::Compute_Plane_Stress_Term(CElement *element, const CConfig *config, unsigned short iGauss) { SU2_MPI::Error("This material model cannot (yet) be used for plane stress.", CURRENT_FUNCTION); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index f40b6b9cc21e..98e5837479b3 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -232,9 +232,9 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo SetVolumeOutputValue("STRESS-XX", iPoint, Node_Struc->GetStress_FEM(iPoint)[0]); SetVolumeOutputValue("STRESS-YY", iPoint, Node_Struc->GetStress_FEM(iPoint)[1]); + SetVolumeOutputValue("STRESS-ZZ", iPoint, Node_Struc->GetStress_FEM(iPoint)[3]); SetVolumeOutputValue("STRESS-XY", iPoint, Node_Struc->GetStress_FEM(iPoint)[2]); - if (nDim == 3){ - SetVolumeOutputValue("STRESS-ZZ", iPoint, Node_Struc->GetStress_FEM(iPoint)[3]); + if (nDim == 3) { SetVolumeOutputValue("STRESS-XZ", iPoint, Node_Struc->GetStress_FEM(iPoint)[4]); SetVolumeOutputValue("STRESS-YZ", iPoint, Node_Struc->GetStress_FEM(iPoint)[5]); } @@ -260,17 +260,17 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector"); if (nDim == 3) AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector"); - AddVolumeOutput("DISPLACEMENT-X", "Displacement_x", "SOLUTION", "x-component of the displacement vector"); - AddVolumeOutput("DISPLACEMENT-Y", "Displacement_y", "SOLUTION", "y-component of the displacement vector"); + AddVolumeOutput("DISPLACEMENT-X", "Displacement_x", "SOLUTION", "x-component of the displacement vector"); + AddVolumeOutput("DISPLACEMENT-Y", "Displacement_y", "SOLUTION", "y-component of the displacement vector"); if (nDim == 3) AddVolumeOutput("DISPLACEMENT-Z", "Displacement_z", "SOLUTION", "z-component of the displacement vector"); if (dynamic) { - AddVolumeOutput("VELOCITY-X", "Velocity_x", "SOLUTION", "x-component of the velocity vector"); - AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); + AddVolumeOutput("VELOCITY-X", "Velocity_x", "SOLUTION", "x-component of the velocity vector"); + AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector"); - AddVolumeOutput("ACCELERATION-X", "Acceleration_x", "SOLUTION", "x-component of the acceleration vector"); - AddVolumeOutput("ACCELERATION-Y", "Acceleration_y", "SOLUTION", "y-component of the acceleration vector"); + AddVolumeOutput("ACCELERATION-X", "Acceleration_x", "SOLUTION", "x-component of the acceleration vector"); + AddVolumeOutput("ACCELERATION-Y", "Acceleration_y", "SOLUTION", "y-component of the acceleration vector"); if (nDim == 3) AddVolumeOutput("ACCELERATION-Z", "Acceleration_z", "SOLUTION", "z-component of the acceleration vector"); } @@ -278,14 +278,13 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); } - AddVolumeOutput("STRESS-XX", "Sxx", "STRESS", "x-component of the normal stress vector"); - AddVolumeOutput("STRESS-YY", "Syy", "STRESS", "y-component of the normal stress vector"); - AddVolumeOutput("STRESS-XY", "Sxy", "STRESS", "xy shear stress component"); - + AddVolumeOutput("STRESS-XX", "Sxx", "STRESS", "x-component of the normal stress vector"); + AddVolumeOutput("STRESS-YY", "Syy", "STRESS", "y-component of the normal stress vector"); + AddVolumeOutput("STRESS-ZZ", "Szz", "STRESS", "z-component of the normal stress vector"); + AddVolumeOutput("STRESS-XY", "Sxy", "STRESS", "xy shear stress component"); if (nDim == 3) { - AddVolumeOutput("STRESS-ZZ", "Szz", "STRESS", "z-component of the normal stress vector"); - AddVolumeOutput("STRESS-XZ", "Sxz", "STRESS", "xz shear stress component"); - AddVolumeOutput("STRESS-YZ", "Syz", "STRESS", "yz shear stress component"); + AddVolumeOutput("STRESS-XZ", "Sxz", "STRESS", "xz shear stress component"); + AddVolumeOutput("STRESS-YZ", "Syz", "STRESS", "yz shear stress component"); } AddVolumeOutput("VON_MISES_STRESS", "Von_Mises_Stress", "STRESS", "von-Mises stress"); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index c8ecd94985e6..8b31c874bdce 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1163,7 +1163,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const su2double stress_scale = 1.0 / stressParam[0]; const su2double ks_mult = stressParam[1]; - const unsigned short nStress = (nDim == 2) ? 3 : 6; + const unsigned short nStress = 2 * nDim; su2double StressPenalty = 0.0; su2double MaxVonMises_Stress = 0.0; @@ -1286,36 +1286,10 @@ 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++) { - /*--- 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); + const auto vms = CFEAElasticity::VonMisesStress(nDim, nodes->GetStress_FEM(iPoint)); nodes->SetVonMises_Stress(iPoint, vms); diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp index d921fd1dc0e4..f00e8f85f4e2 100644 --- a/SU2_CFD/src/variables/CFEAVariable.cpp +++ b/SU2_CFD/src/variables/CFEAVariable.cpp @@ -54,9 +54,8 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig const bool fsi_analysis = config->GetFSI_Simulation() || multizone; VonMises_Stress.resize(nPoint) = su2double(0.0); - - if (nDim==2) Stress.resize(nPoint,3); - else Stress.resize(nPoint,6); + /*--- For completeness we also output sigma_zz in 2D. ---*/ + Stress.resize(nPoint, 2 * nDim); /*--- Initialization of variables ---*/ for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) diff --git a/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg b/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg new file mode 100644 index 000000000000..f9f2ea27b60b --- /dev/null +++ b/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg @@ -0,0 +1,45 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Cantilever 2D beam in plane strain. % +% File Version 8.3.0 "Harrier" % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS +MATERIAL_MODEL= LINEAR_ELASTIC +FORMULATION_ELASTICITY_2D= PLANE_STRAIN + +ELASTICITY_MODULUS= 3E7 +POISSON_RATIO= 0.3 +MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5 +MATERIAL_REFERENCE_TEMPERATURE= 300 +FREESTREAM_TEMPERATURE= 350 +MATERIAL_DENSITY= 7854 + +MARKER_CLAMPED= ( x_minus ) +MARKER_PRESSURE= ( y_minus, 0, y_plus, 0 ) +MARKER_LOAD= ( x_plus, 1, 1000, 0, -1, 0 ) +% For the 3D version: +% MARKER_SYM= ( z_minus, z_plus ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 500 + +% RECTANGLE for 2D and BOX for 3D. +MESH_FORMAT= RECTANGLE +MESH_BOX_LENGTH= ( 1, 0.05, 0.01 ) +MESH_BOX_SIZE= ( 129, 33, 2 ) + +% The Von-Misses stress was checked against a 3D version of this case. +% 2D and 3D must match exactly. +SCREEN_OUTPUT= INNER_ITER, RMS_RES, VMS, LINSOL + +TABULAR_FORMAT= CSV +CONV_FILENAME= history +VOLUME_FILENAME= beam +RESTART_FILENAME= restart +SOLUTION_FILENAME= restart +OUTPUT_WRT_FREQ= 1 +ITER= 1 diff --git a/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg b/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg new file mode 100644 index 000000000000..453fcf847f8d --- /dev/null +++ b/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg @@ -0,0 +1,54 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Cantilever 2D beam in plane stress. % +% File Version 8.3.0 "Harrier" % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +MATERIAL_MODEL= NEO_HOOKEAN +FORMULATION_ELASTICITY_2D= PLANE_STRESS + +ELASTICITY_MODULUS= 3E7 +POISSON_RATIO= 0.3 +MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5 +MATERIAL_REFERENCE_TEMPERATURE= 300 +FREESTREAM_TEMPERATURE= 350 +MATERIAL_DENSITY= 7854 + +MARKER_CLAMPED= ( x_minus ) +MARKER_PRESSURE= ( y_minus, 0, y_plus, 0 ) +% For the 3D version: +% MARKER_PRESSURE= ( y_minus, 0, y_plus, 0, z_minus, 0, z_plus, 0 ) +MARKER_LOAD= ( x_plus, 1, 1000, 0, -1, 0 ) + +LINEAR_SOLVER= FGCRODR +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-4 +LINEAR_SOLVER_ITER= 100 +LINEAR_SOLVER_RESTART_FREQUENCY= 100 +LINEAR_SOLVER_RESTART_DEFLATION= 20 + +% RECTANGLE for 2D and BOX for 3D. +MESH_FORMAT= RECTANGLE +MESH_BOX_LENGTH= ( 1, 0.05, 0.005 ) +MESH_BOX_SIZE= ( 129, 33, 2 ) + +% The Von-Misses stress was checked against a 3D version of this case. +% Note that plane stress conditions are difficult to replicate in 3D +% because linear systems become stiff for very thin elements. +% However we can see the 3D result approach the 2D result as we decrease +% the box size in the z direction. +SCREEN_OUTPUT= INNER_ITER, RMS_RES, VMS, LINSOL + +CONV_FIELD= REL_RMS_RTOL +CONV_STARTITER= 0 +CONV_RESIDUAL_MINVAL= -6 +INNER_ITER= 30 + +TABULAR_FORMAT= CSV +CONV_FILENAME= history +VOLUME_FILENAME= beam +RESTART_FILENAME= restart +SOLUTION_FILENAME= restart +OUTPUT_WRT_FREQ= 1 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 8c9279df366a..f4b5df1cacaf 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1262,6 +1262,23 @@ def main(): rotating_cylinder_fea.test_vals_aarch64 = [-6.861939, -6.835539, -6.895498, 22, -8.313847, 1.6502e+08] test_list.append(rotating_cylinder_fea) + # 2D beam in plain strain with thermal expansion. This tests fixes to the 2D Von Misses stress calculation. + linear_plane_strain = TestCase('linear_plane_strain') + linear_plane_strain.cfg_dir = "fea_fsi/VonMissesVerif" + linear_plane_strain.cfg_file = "linear_plane_strain_2d.cfg" + linear_plane_strain.test_iter = 0 + linear_plane_strain.test_vals = [-6.038376, -5.972195, 0, 1.2014e+05, 125, -8.009476] + test_list.append(linear_plane_strain) + + # 2D beam in plain stress with thermal expansion. This tests fixes to the 2D Von Misses stress calculation, + # and to the plane stress formulation with thermal expansion for nonlinear materials. + nonlinear_plane_stress = TestCase('nonlinear_plane_stress') + nonlinear_plane_stress.cfg_dir = "fea_fsi/VonMissesVerif" + nonlinear_plane_stress.cfg_file = "nonlinear_plane_stress_2d.cfg" + nonlinear_plane_stress.test_iter = 19 + nonlinear_plane_stress.test_vals = [-7.444178, -3.365607, -13.998049, 1.6248e+05, 43, -4.095522] + test_list.append(nonlinear_plane_stress) + # Dynamic beam, 2d dynbeam2d = TestCase('dynbeam2d') dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d" @@ -1662,7 +1679,7 @@ def main(): species3_primitiveVenturi_inletFile.test_vals = [-5.537438, -4.503863, -4.553632, -5.400874, -0.945967, -5.818774, -5.945211, 5.000000, -0.544749, 5.000000, -2.599435, 5.000000, -0.596360] test_list.append(species3_primitiveVenturi_inletFile) - # 3 species (2 eq) primitive venturi mixing with new flux and value boundary conditions + # 3 species (2 eq) primitive venturi mixing with new flux and value boundary conditions species3_primitiveVenturi_fluxvalue = TestCase('species3_primitiveVenturi_fluxvalue') species3_primitiveVenturi_fluxvalue.cfg_dir = "species_transport/venturi_primitive_3species" species3_primitiveVenturi_fluxvalue.cfg_file = "species3_primitiveVenturi_flux_value.cfg" From 6a88d3ca80c2e5e9ee24b08d54dd6058f1a0bb20 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 31 Dec 2025 14:06:44 -0800 Subject: [PATCH 11/12] update regressions --- .../VonMissesVerif/linear_plane_strain_2d.cfg | 2 +- .../nonlinear_plane_stress_2d.cfg | 2 +- TestCases/hybrid_regression.py | 6 ++--- TestCases/parallel_regression.py | 23 +++++++++---------- TestCases/parallel_regression_AD.py | 4 ++-- TestCases/py_wrapper/custom_load_fea/run.py | 2 +- TestCases/serial_regression.py | 10 ++++---- 7 files changed, 24 insertions(+), 25 deletions(-) diff --git a/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg b/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg index f9f2ea27b60b..399c03dc5f29 100644 --- a/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg +++ b/TestCases/fea_fsi/VonMissesVerif/linear_plane_strain_2d.cfg @@ -32,7 +32,7 @@ MESH_FORMAT= RECTANGLE MESH_BOX_LENGTH= ( 1, 0.05, 0.01 ) MESH_BOX_SIZE= ( 129, 33, 2 ) -% The Von-Misses stress was checked against a 3D version of this case. +% The von Mises stress was checked against a 3D version of this case. % 2D and 3D must match exactly. SCREEN_OUTPUT= INNER_ITER, RMS_RES, VMS, LINSOL diff --git a/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg b/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg index 453fcf847f8d..5d91de45440a 100644 --- a/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg +++ b/TestCases/fea_fsi/VonMissesVerif/nonlinear_plane_stress_2d.cfg @@ -34,7 +34,7 @@ MESH_FORMAT= RECTANGLE MESH_BOX_LENGTH= ( 1, 0.05, 0.005 ) MESH_BOX_SIZE= ( 129, 33, 2 ) -% The Von-Misses stress was checked against a 3D version of this case. +% The von Mises stress was checked against a 3D version of this case. % Note that plane stress conditions are difficult to replicate in 3D % because linear systems become stiff for very thin elements. % However we can see the 3D result approach the 2D result as we decrease diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 726c2f491438..3b13cfb22f4a 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -688,7 +688,7 @@ def main(): dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d" dynbeam2d.cfg_file = "configBeam_2d.cfg" dynbeam2d.test_iter = 6 - dynbeam2d.test_vals = [-3.240016, 2.895057, -0.353147, 66127.000000] + dynbeam2d.test_vals = [-3.240013, 2.895060, -0.353141, 76220] dynbeam2d.unsteady = True test_list.append(dynbeam2d) @@ -697,7 +697,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4.000000, 0.000000, -3.726028, -4.277769] + fsi2d.test_vals = [4, 0, -3.726029, -4.277531] fsi2d.multizone= True fsi2d.unsteady = True fsi2d.enabled_with_tsan = False @@ -708,7 +708,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.330725, -4.152983, 0.000000, 102.000000] + dyn_fsi.test_vals = [-4.330725, -4.152808, 0, 102] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index f4b5df1cacaf..6067df66f371 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1257,26 +1257,26 @@ def main(): rotating_cylinder_fea.test_iter = 0 # For a thin disk with the inner and outer radius of this geometry, from # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4, - # the maximum stress is 165.6MPa, we get a Von Misses stress very close to that. + # the maximum stress is 165.6MPa, we get a von Mises stress very close to that. rotating_cylinder_fea.test_vals = [-6.886145, -6.917148, -6.959634, 23, -8.369804, 1.6502e+08] rotating_cylinder_fea.test_vals_aarch64 = [-6.861939, -6.835539, -6.895498, 22, -8.313847, 1.6502e+08] test_list.append(rotating_cylinder_fea) - # 2D beam in plain strain with thermal expansion. This tests fixes to the 2D Von Misses stress calculation. + # 2D beam in plain strain with thermal expansion. This tests fixes to the 2D Von Mises stress calculation. linear_plane_strain = TestCase('linear_plane_strain') linear_plane_strain.cfg_dir = "fea_fsi/VonMissesVerif" linear_plane_strain.cfg_file = "linear_plane_strain_2d.cfg" linear_plane_strain.test_iter = 0 - linear_plane_strain.test_vals = [-6.038376, -5.972195, 0, 1.2014e+05, 125, -8.009476] + linear_plane_strain.test_vals = [-6.157686, -6.051735, 0, 120140, 325, -8.105017] test_list.append(linear_plane_strain) - # 2D beam in plain stress with thermal expansion. This tests fixes to the 2D Von Misses stress calculation, + # 2D beam in plain stress with thermal expansion. This tests fixes to the 2D von Mises stress calculation, # and to the plane stress formulation with thermal expansion for nonlinear materials. nonlinear_plane_stress = TestCase('nonlinear_plane_stress') nonlinear_plane_stress.cfg_dir = "fea_fsi/VonMissesVerif" nonlinear_plane_stress.cfg_file = "nonlinear_plane_stress_2d.cfg" nonlinear_plane_stress.test_iter = 19 - nonlinear_plane_stress.test_vals = [-7.444178, -3.365607, -13.998049, 1.6248e+05, 43, -4.095522] + nonlinear_plane_stress.test_vals = [7.433102, -3.355151, -13.983258, 162480, 43, -4.071408] test_list.append(nonlinear_plane_stress) # Dynamic beam, 2d @@ -1284,7 +1284,7 @@ def main(): dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d" dynbeam2d.cfg_file = "configBeam_2d.cfg" dynbeam2d.test_iter = 6 - dynbeam2d.test_vals = [-3.240015, 2.895057, -0.353146, 66127.000000] + dynbeam2d.test_vals = [-3.240012, 2.895060, -0.353140, 76220] dynbeam2d.unsteady = True test_list.append(dynbeam2d) @@ -1293,7 +1293,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4.000000, 0.000000, -3.726013, -4.277768] + fsi2d.test_vals = [4, 0, -3.726013, -4.277529] fsi2d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") fsi2d.multizone= True fsi2d.unsteady = True @@ -1313,7 +1313,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.330741, -4.153001, 0.000000, 97.000000] + dyn_fsi.test_vals = [-4.330741, -4.152826, 0, 97] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -1443,8 +1443,7 @@ def main(): pywrapper_custom_fea_load.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_custom_fea_load.cfg_file = "config.cfg" pywrapper_custom_fea_load.test_iter = 13 - pywrapper_custom_fea_load.test_vals = [-7.263559, -4.946814, -14.165142, 34.000000, -6.380144, 320.580000] - pywrapper_custom_fea_load.test_vals_aarch64 = [-7.263558, -4.946814, -14.165142, 35.000000, -6.802790, 320.580000] + pywrapper_custom_fea_load.test_vals = [-7.262040, -4.945686, -14.163208, 34, -6.009642, 362.23] pywrapper_custom_fea_load.command = TestCase.Command("mpirun -np 2", "python", "run.py") test_list.append(pywrapper_custom_fea_load) @@ -1453,7 +1452,7 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.726013, -4.277768] + pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.726013, -4.277529] pywrapper_fsi2d.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--nZone 2 --fsi True --parallel -f") pywrapper_fsi2d.unsteady = True pywrapper_fsi2d.multizone = True @@ -1464,7 +1463,7 @@ def main(): pywrapper_unsteadyFSI.cfg_dir = "py_wrapper/dyn_fsi" pywrapper_unsteadyFSI.cfg_file = "config.cfg" pywrapper_unsteadyFSI.test_iter = 4 - pywrapper_unsteadyFSI.test_vals = [0.000000, 31.000000, 5.000000, 58.000000, -1.756677, -2.828286, -7.638545, -6.863959, 0.000156] + pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 58, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156] pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") pywrapper_unsteadyFSI.unsteady = True pywrapper_unsteadyFSI.multizone = True diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index f800dbdc7697..0aada2937bc9 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -276,7 +276,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-4.772623, 0.917211, -3.863369, 0.295450, 3.841200] + discadj_fsi2.test_vals = [-4.772676, 0.917733, -3.863369, 0.295450, 3.839800] discadj_fsi2.test_vals_aarch64 = [-4.773008, 0.916312, -3.863369, 0.295450, 3.841200] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) @@ -472,7 +472,7 @@ def main(): pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" pywrapper_Unst_FEA_AD.test_iter = 100 - pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.319877, 0.320149, -0.184491, -0.184509] + pywrapper_Unst_FEA_AD.test_vals = [0.262441, 0.262444, 0.324361, 0.324640, -0.188221, -0.188238] pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") pywrapper_Unst_FEA_AD.timeout = 1600 pywrapper_Unst_FEA_AD.tol = 0.00001 diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py index f5389e70fdb5..948226ab6aca 100755 --- a/TestCases/py_wrapper/custom_load_fea/run.py +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -83,7 +83,7 @@ def main(): if NodeFound: print(f"Vertical displacement of tip: {Disp}") # Test the value against expected. - assert abs(Disp / 0.095439 - 1) < 1e-5, "Test FAILED" + assert abs(Disp / 0.095439 - 1) < 1e-5, f"Test FAILED, computed displacement = {Disp}" # Finalize the solver and exit cleanly. SU2Driver.Finalize() diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 2498aa5e4922..be32fc412224 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1037,7 +1037,7 @@ def main(): dynbeam2d.cfg_file = "configBeam_2d.cfg" dynbeam2d.unsteady = True dynbeam2d.test_iter = 6 - dynbeam2d.test_vals = [-3.240015, 2.895057, -0.353146, 66127.000000] + dynbeam2d.test_vals = [-3.240012, 2.895060, -0.353140, 76220] test_list.append(dynbeam2d) # # FSI, 2d @@ -1045,7 +1045,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4.000000, 0.000000, -3.726014, -4.277767] + fsi2d.test_vals = [4, 0, -3.726015, -4.277528] fsi2d.multizone = True fsi2d.unsteady = True test_list.append(fsi2d) @@ -1073,7 +1073,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.330728, -4.152995, 0.000000, 86.000000] + dyn_fsi.test_vals = [-4.330728, -4.152820, 0, 85] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -1576,10 +1576,10 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4.000000, 0.000000, -3.726014, -4.277767] + pywrapper_fsi2d.test_vals = [4, 0, -3.726015, -4.277528] pywrapper_fsi2d.command = TestCase.Command(exec = "SU2_CFD.py", param = "--nZone 2 --fsi True -f") pywrapper_fsi2d.unsteady = True - pywrapper_fsi2d.multizone = True + pywrapper_fsi2d.multizone = True pywrapper_fsi2d.timeout = 1600 pywrapper_fsi2d.tol = 0.00001 pywrapper_fsi2d.enabled_with_asan = False From a37a8e739f0a382c960c12a8723afa47ad82663e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 31 Dec 2025 19:26:43 -0800 Subject: [PATCH 12/12] update --- TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref | 16 ++++++++-------- TestCases/parallel_regression.py | 2 +- TestCases/py_wrapper/custom_load_fea/run.py | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref index 0dc991272850..60e4d2cf956e 100644 --- a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref +++ b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref @@ -1,9 +1,9 @@ INDEX GRAD -0 -4.570869246220687e-04 -1 -2.401466765893676e-04 -2 -9.134698429099325e-05 -3 -1.628087033585022e-05 -4 -1.741052439794770e-05 -5 -9.462489834101617e-05 -6 -2.452466275277050e-04 -7 -4.635483662879789e-04 +0 -4.570880257851463e-04 +1 -2.401474161436175e-04 +2 -9.134742016285272e-05 +3 -1.628103644960537e-05 +4 -1.741046823504659e-05 +5 -9.462481396563234e-05 +6 -2.452469693841206e-04 +7 -4.635498187330369e-04 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6067df66f371..dfb3229a092e 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1276,7 +1276,7 @@ def main(): nonlinear_plane_stress.cfg_dir = "fea_fsi/VonMissesVerif" nonlinear_plane_stress.cfg_file = "nonlinear_plane_stress_2d.cfg" nonlinear_plane_stress.test_iter = 19 - nonlinear_plane_stress.test_vals = [7.433102, -3.355151, -13.983258, 162480, 43, -4.071408] + nonlinear_plane_stress.test_vals = [-7.433102, -3.355151, -13.983258, 162480, 43, -4.071408] test_list.append(nonlinear_plane_stress) # Dynamic beam, 2d diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py index 948226ab6aca..65c089a79025 100755 --- a/TestCases/py_wrapper/custom_load_fea/run.py +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -83,7 +83,7 @@ def main(): if NodeFound: print(f"Vertical displacement of tip: {Disp}") # Test the value against expected. - assert abs(Disp / 0.095439 - 1) < 1e-5, f"Test FAILED, computed displacement = {Disp}" + assert abs(Disp / 0.0954409442 - 1) < 1e-5, f"Test FAILED, computed displacement = {Disp}" # Finalize the solver and exit cleanly. SU2Driver.Finalize()