From 46405cdd5ef8dab7ac6676052c7a38d6b4263423 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 21 Dec 2024 14:48:10 -0800 Subject: [PATCH 1/7] add thermal expansion to CConfig --- Common/include/CConfig.hpp | 11 +++++++++-- Common/src/CConfig.cpp | 14 +++++++++++--- .../numerics/elasticity/CFEAElasticity.hpp | 2 ++ .../src/numerics/elasticity/CFEAElasticity.cpp | 13 ++++++------- .../numerics/elasticity/CFEALinearElasticity.cpp | 7 +++++++ .../src/numerics/elasticity/nonlinear_models.cpp | 2 +- TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg | 16 ++++++++-------- 7 files changed, 44 insertions(+), 21 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 0c32f044288..9eb6ef74e71 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1013,10 +1013,12 @@ class CConfig { su2double *Int_Coeffs; /*!< \brief Time integration coefficients for structural method. */ unsigned short nElasticityMod, /*!< \brief Number of different values for the elasticity modulus. */ nPoissonRatio, /*!< \brief Number of different values for the Poisson ratio modulus. */ - nMaterialDensity; /*!< \brief Number of different values for the Material density. */ + nMaterialDensity, /*!< \brief Number of different values for the Material density. */ + nMaterialThermalExpansion; /*!< \brief Number of different values for thermal expansion coefficient. */ su2double *ElasticityMod, /*!< \brief Value of the elasticity moduli. */ *PoissonRatio, /*!< \brief Value of the Poisson ratios. */ - *MaterialDensity; /*!< \brief Value of the Material densities. */ + *MaterialDensity, /*!< \brief Value of the Material densities. */ + *MaterialThermalExpansion; /*!< \brief Value of the thermal expansion coefficients. */ unsigned short nElectric_Field, /*!< \brief Number of different values for the electric field in the membrane. */ nDim_Electric_Field; /*!< \brief Dimensionality of the problem. */ unsigned short nDim_RefNode; /*!< \brief Dimensionality of the vector . */ @@ -2389,6 +2391,11 @@ class CConfig { */ su2double GetMaterialDensity(unsigned short id_val) const { return MaterialDensity[id_val]; } + /*! + * \brief Get the thermal expansion coefficient. + */ + su2double GetMaterialThermalExpansion(unsigned short id_val) const { return MaterialThermalExpansion[id_val]; } + /*! * \brief Compressibility/incompressibility of the solids analysed using the structural solver. * \return Compressible or incompressible. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c5bd653b04a..90af5948d25 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2440,6 +2440,8 @@ void CConfig::SetConfig_Options() { addDoubleListOption("POISSON_RATIO", nPoissonRatio, PoissonRatio); /* DESCRIPTION: Material density */ addDoubleListOption("MATERIAL_DENSITY", nMaterialDensity, MaterialDensity); + /* DESCRIPTION: Material thermal expansion coefficient */ + addDoubleListOption("MATERIAL_THERMAL_EXPANSION_COEFF", nMaterialThermalExpansion, MaterialThermalExpansion); /* DESCRIPTION: Knowles B constant */ addDoubleOption("KNOWLES_B", Knowles_B, 1.0); /* DESCRIPTION: Knowles N constant */ @@ -4834,9 +4836,15 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854; } - if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity) { - SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, and MATERIAL_DENSITY need to have the same number " - "of entries (the number of materials).", CURRENT_FUNCTION); + if (nMaterialThermalExpansion == 0) { + nMaterialThermalExpansion = 1; + MaterialThermalExpansion = new su2double[1](); + } + + if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity || + nElasticityMod != nMaterialThermalExpansion) { + SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, MATERIAL_DENSITY, and THERMAL_EXPANSION_COEFF need " + "to have the same number of entries (the number of materials).", CURRENT_FUNCTION); } if (nElectric_Constant == 0) { diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 11470fd0fba..43d5ebd5a63 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -54,6 +54,7 @@ class CFEAElasticity : public CNumerics { su2double Nu = 0.0; /*!< \brief Aux. variable, Poisson's ratio. */ su2double Rho_s = 0.0; /*!< \brief Aux. variable, Structural density. */ su2double Rho_s_DL = 0.0; /*!< \brief Aux. variable, Structural density (for dead loads). */ + su2double Alpha = 0.0; /*!< \brief Aux. variable, thermal expansion coefficient. */ su2double Mu = 0.0; /*!< \brief Aux. variable, Lame's coeficient. */ su2double Lambda = 0.0; /*!< \brief Aux. variable, Lame's coeficient. */ @@ -63,6 +64,7 @@ class CFEAElasticity : public CNumerics { su2double *Nu_i = nullptr; /*!< \brief Poisson's ratio. */ su2double *Rho_s_i = nullptr; /*!< \brief Structural density. */ su2double *Rho_s_DL_i = nullptr; /*!< \brief Structural density (for dead loads). */ + su2double *Alpha_i = nullptr; /*!< \brief Thermal expansion coefficient. */ su2double **Ba_Mat = nullptr; /*!< \brief Matrix B for node a - Auxiliary. */ su2double **Bb_Mat = nullptr; /*!< \brief Matrix B for node b - Auxiliary. */ diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 44014591204..b91be37223e 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -44,19 +44,16 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, const auto nProp = config->GetnElasticityMat(); E_i = new su2double[nProp]; - for (iVar = 0; iVar < nProp; iVar++) - E_i[iVar] = config->GetElasticyMod(iVar); - Nu_i = new su2double[nProp]; - for (iVar = 0; iVar < nProp; iVar++) - Nu_i[iVar] = config->GetPoissonRatio(iVar); - Rho_s_i = new su2double[nProp]; // For inertial effects Rho_s_DL_i = new su2double[nProp]; // For dead loads - + Alpha_i = new su2double[nProp]; for (iVar = 0; iVar < nProp; iVar++) { + E_i[iVar] = config->GetElasticyMod(iVar); + Nu_i[iVar] = config->GetPoissonRatio(iVar); Rho_s_DL_i[iVar] = config->GetMaterialDensity(iVar); Rho_s_i[iVar] = pseudo_static ? 0.0 : config->GetMaterialDensity(iVar); + Alpha_i[iVar] = config->GetMaterialThermalExpansion(iVar); } // Initialization @@ -64,6 +61,7 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, Nu = Nu_i[0]; Rho_s = Rho_s_i[0]; Rho_s_DL = Rho_s_DL_i[0]; + Alpha = Alpha_i[0]; Compute_Lame_Parameters(); @@ -302,6 +300,7 @@ void CFEAElasticity::SetElement_Properties(const CElement *element, const CConfi Nu = Nu_i[element->Get_iProp()]; Rho_s = Rho_s_i[element->Get_iProp()]; Rho_s_DL = Rho_s_DL_i[element->Get_iProp()]; + Alpha = Alpha_i[element->Get_iProp()]; switch (config->GetDV_FEA()) { case YOUNG_MODULUS: diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index a6763bed0b8..de7ad86219f 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -99,6 +99,13 @@ void CFEALinearElasticity::Compute_Tangent_Matrix(CElement *element, const CConf for (iNode = 0; iNode < nNode; iNode++) { + su2double KAux_t_a[3] = {0.0}; + const su2double thermalStress = -E * 2e-4 / (1 - 2 * Nu); + for (iVar = 0; iVar < nDim; iVar++) { + KAux_t_a[iVar] += Weight * thermalStress * GradNi_Ref_Mat[iNode][iVar] * Jac_X; + } + element->Add_Kt_a(iNode, KAux_t_a); + if (nDim == 2) { Ba_Mat[0][0] = GradNi_Ref_Mat[iNode][0]; Ba_Mat[1][1] = GradNi_Ref_Mat[iNode][1]; diff --git a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp index d21fc228300..55d48a29010 100644 --- a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp +++ b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp @@ -112,7 +112,7 @@ void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfi for (iVar = 0; iVar < 3; iVar++) { for (jVar = 0; jVar < 3; jVar++) { su2double dij = deltaij(iVar,jVar); - Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + Lambda_J * log(J_F) * dij; + Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) - E * 2e-4 / (1 - 2 * Nu)) * dij; } } diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg index a6c9cde3c84..128958ecb19 100644 --- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg @@ -15,13 +15,13 @@ MESH_FILENAME= meshBeam_3d.su2 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 MATERIAL_DENSITY=7854 -MARKER_CLAMPED = ( left , right ) -MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0) -MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) -LINEAR_SOLVER= FGMRES -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-8 -LINEAR_SOLVER_ITER= 500 +MARKER_CLAMPED = ( left ) +MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0) +%MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-9 +LINEAR_SOLVER_ITER= 1000 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV CONV_FILENAME= history_beam @@ -29,4 +29,4 @@ VOLUME_FILENAME= beam RESTART_FILENAME= restart_beam.dat SOLUTION_FILENAME= restart_beam.dat OUTPUT_WRT_FREQ= 1 -INNER_ITER=1 +INNER_ITER=10 From 458e73224ccdee72f99d79635745197b879c4973 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 22 Dec 2024 18:00:05 -0800 Subject: [PATCH 2/7] thermal stress term, not coupled to heat solver --- Common/include/CConfig.hpp | 8 +++- Common/include/geometry/elements/CElement.hpp | 41 ++++++++++++++--- Common/src/CConfig.cpp | 2 + Common/src/geometry/elements/CElement.cpp | 2 + Common/src/geometry/elements/CQUAD4.cpp | 1 - .../numerics/elasticity/CFEAElasticity.hpp | 11 +++-- .../elasticity/CFEANonlinearElasticity.hpp | 3 +- .../numerics/elasticity/nonlinear_models.hpp | 8 ++-- .../numerics/elasticity/CFEAElasticity.cpp | 46 ++++++------------- .../elasticity/CFEALinearElasticity.cpp | 10 ++-- .../elasticity/CFEANonlinearElasticity.cpp | 6 +-- .../numerics/elasticity/nonlinear_models.cpp | 26 ++++++----- SU2_CFD/src/solvers/CFEASolver.cpp | 8 ++++ .../fea_fsi/StatBeam_3d/configBeam_3d.cfg | 8 ++-- 14 files changed, 110 insertions(+), 70 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 9eb6ef74e71..d91a6011611 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1018,7 +1018,8 @@ class CConfig { su2double *ElasticityMod, /*!< \brief Value of the elasticity moduli. */ *PoissonRatio, /*!< \brief Value of the Poisson ratios. */ *MaterialDensity, /*!< \brief Value of the Material densities. */ - *MaterialThermalExpansion; /*!< \brief Value of the thermal expansion coefficients. */ + *MaterialThermalExpansion, /*!< \brief Value of the thermal expansion coefficients. */ + MaterialReferenceTemperature; /*!< \brief Value of the reference temperature for thermal expansion. */ unsigned short nElectric_Field, /*!< \brief Number of different values for the electric field in the membrane. */ nDim_Electric_Field; /*!< \brief Dimensionality of the problem. */ unsigned short nDim_RefNode; /*!< \brief Dimensionality of the vector . */ @@ -2396,6 +2397,11 @@ class CConfig { */ su2double GetMaterialThermalExpansion(unsigned short id_val) const { return MaterialThermalExpansion[id_val]; } + /*! + * \brief Temperature at which there is no stress from thermal expansion. + */ + su2double GetMaterialReferenceTemperature() const { return MaterialReferenceTemperature; } + /*! * \brief Compressibility/incompressibility of the solids analysed using the structural solver. * \return Compressible or incompressible. diff --git a/Common/include/geometry/elements/CElement.hpp b/Common/include/geometry/elements/CElement.hpp index 535323ec0c0..6f7bd435560 100644 --- a/Common/include/geometry/elements/CElement.hpp +++ b/Common/include/geometry/elements/CElement.hpp @@ -67,6 +67,8 @@ class CElement { su2activematrix NodalExtrap; /*!< \brief Coordinates of the nodal points for Gaussian extrapolation. */ su2activematrix NodalStress; /*!< \brief Stress at the nodes. */ + su2activevector NodalTemperature; /*!< \brief Temperature at the nodes. */ + /*--- Stiffness and load matrices. ---*/ std::vector Kab; /*!< \brief Structure for the constitutive component of the tangent matrix. */ su2activematrix Mab; /*!< \brief Structure for the nodal components of the mass matrix. */ @@ -151,7 +153,7 @@ class CElement { * \param[in] iDim - Dimension. * \param[in] val_CoordRef - Value of the coordinate. */ - inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordRef) { + inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordRef) { RefCoord(iNode, iDim) = val_CoordRef; } @@ -161,10 +163,22 @@ class CElement { * \param[in] iDim - Dimension. * \param[in] val_CoordRef - Value of the coordinate. */ - inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordCurr) { + inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordCurr) { CurrentCoord(iNode, iDim) = val_CoordCurr; } + /*! + * \brief Set the value of the temperature of a node. + */ + inline void SetTemperature(unsigned short iNode, const su2double& val_Temperature) { + NodalTemperature[iNode] = val_Temperature; + } + + /*! + * \brief Set the value of the temperature of all nodes. + */ + inline void SetTemperature(const su2double& val_Temperature) { NodalTemperature = val_Temperature; } + /*! * \brief Get the value of the coordinate of the nodes in the reference configuration. * \param[in] iNode - Index of node. @@ -181,6 +195,17 @@ class CElement { */ inline su2double GetCurr_Coord(unsigned short iNode, unsigned short iDim) const { return CurrentCoord(iNode, iDim); } + /*! + * \brief Get the value of the temperature at a Gaussian integration point. + */ + inline su2double GetTemperature(unsigned short iGauss) const { + su2double Temperature = 0; + for (auto iNode = 0u; iNode < nNodes; ++iNode) { + Temperature += GetNi(iNode, iGauss) * NodalTemperature[iNode]; + } + return Temperature; + } + /*! * \brief Get the weight of the corresponding Gaussian Point. * \param[in] iGauss - index of the Gaussian point. @@ -214,7 +239,9 @@ class CElement { * \param[in] nodeB - index of Node b. * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution. */ - inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, su2double val_Mab) { Mab(nodeA, nodeB) += val_Mab; } + inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Mab) { + Mab(nodeA, nodeB) += val_Mab; + } /*! * \brief Add the value of a submatrix K relating nodes a and b, for the constitutive term. @@ -243,7 +270,7 @@ class CElement { * \param[in] nodeB - index of Node b. * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution. */ - inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, su2double val_Ks_ab) { + inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Ks_ab) { Ks_ab(nodeA, nodeB) += val_Ks_ab; } @@ -354,7 +381,7 @@ class CElement { * \param[in] iVar - Variable index. * \param[in] val_Stress - Value of the stress added. */ - inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, su2double val_Stress) { + inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, const su2double& val_Stress) { NodalStress(iNode, iVar) += val_Stress; } @@ -789,7 +816,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> { /*! * \brief Shape functions (Ni) evaluated at point Xi,Eta. */ - inline static void ShapeFunctions(su2double Xi, su2double Eta, su2double* Ni) { + inline static void ShapeFunctions(const su2double& Xi, const su2double& Eta, su2double* Ni) { Ni[0] = 0.25 * (1.0 - Xi) * (1.0 - Eta); Ni[1] = 0.25 * (1.0 + Xi) * (1.0 - Eta); Ni[2] = 0.25 * (1.0 + Xi) * (1.0 + Eta); @@ -799,7 +826,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> { /*! * \brief Shape function Jacobian (dNi) evaluated at point Xi,Eta. */ - inline static void ShapeFunctionJacobian(su2double Xi, su2double Eta, su2double dNi[][2]) { + inline static void ShapeFunctionJacobian(const su2double& Xi, const su2double& Eta, su2double dNi[][2]) { dNi[0][0] = -0.25 * (1.0 - Eta); dNi[0][1] = -0.25 * (1.0 - Xi); dNi[1][0] = 0.25 * (1.0 - Eta); diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 90af5948d25..5c220686271 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2442,6 +2442,8 @@ void CConfig::SetConfig_Options() { addDoubleListOption("MATERIAL_DENSITY", nMaterialDensity, MaterialDensity); /* DESCRIPTION: Material thermal expansion coefficient */ addDoubleListOption("MATERIAL_THERMAL_EXPANSION_COEFF", nMaterialThermalExpansion, MaterialThermalExpansion); + /* DESCRIPTION: Temperature at which there is no stress from thermal expansion */ + addDoubleOption("MATERIAL_REFERENCE_TEMPERATURE", MaterialReferenceTemperature, 288.15); /* DESCRIPTION: Knowles B constant */ addDoubleOption("KNOWLES_B", Knowles_B, 1.0); /* DESCRIPTION: Knowles N constant */ diff --git a/Common/src/geometry/elements/CElement.cpp b/Common/src/geometry/elements/CElement.cpp index c42a3707562..5acf96a3ac9 100644 --- a/Common/src/geometry/elements/CElement.cpp +++ b/Common/src/geometry/elements/CElement.cpp @@ -47,6 +47,8 @@ CElement::CElement(unsigned short ngauss, unsigned short nnodes, unsigned short NodalStress.resize(nNodes, 6) = su2double(0.0); + NodalTemperature.resize(nNodes) = su2double(0.0); + Mab.resize(nNodes, nNodes); Ks_ab.resize(nNodes, nNodes); Kab.resize(nNodes); diff --git a/Common/src/geometry/elements/CQUAD4.cpp b/Common/src/geometry/elements/CQUAD4.cpp index 07725a3af19..29d3a621281 100644 --- a/Common/src/geometry/elements/CQUAD4.cpp +++ b/Common/src/geometry/elements/CQUAD4.cpp @@ -67,7 +67,6 @@ CQUAD4::CQUAD4() : CElementWithKnownSizes() { /*--- Store the extrapolation functions (used to compute nodal stresses) ---*/ su2double ExtrapCoord[4][2], sqrt3 = 1.732050807568877; - ; ExtrapCoord[0][0] = -sqrt3; ExtrapCoord[0][1] = -sqrt3; diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index 43d5ebd5a63..bbaa805e046 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -59,6 +59,7 @@ class CFEAElasticity : public CNumerics { su2double Mu = 0.0; /*!< \brief Aux. variable, Lame's coeficient. */ su2double Lambda = 0.0; /*!< \brief Aux. variable, Lame's coeficient. */ su2double Kappa = 0.0; /*!< \brief Aux. variable, Compressibility constant. */ + su2double ThermalStressTerm = 0.0; /*!< \brief Aux. variable, Relationship between stress and delta T. */ su2double *E_i = nullptr; /*!< \brief Young's modulus of elasticity. */ su2double *Nu_i = nullptr; /*!< \brief Poisson's ratio. */ @@ -66,6 +67,8 @@ class CFEAElasticity : public CNumerics { su2double *Rho_s_DL_i = nullptr; /*!< \brief Structural density (for dead loads). */ su2double *Alpha_i = nullptr; /*!< \brief Thermal expansion coefficient. */ + su2double ReferenceTemperature = 0.0; /*!< \brief Reference temperature for thermal expansion. */ + su2double **Ba_Mat = nullptr; /*!< \brief Matrix B for node a - Auxiliary. */ su2double **Bb_Mat = nullptr; /*!< \brief Matrix B for node b - Auxiliary. */ su2double *Ni_Vec = nullptr; /*!< \brief Vector of shape functions - Auxiliary. */ @@ -74,8 +77,6 @@ class CFEAElasticity : public CNumerics { su2double **GradNi_Ref_Mat = nullptr; /*!< \brief Gradients of Ni - Auxiliary. */ su2double **GradNi_Curr_Mat = nullptr; /*!< \brief Gradients of Ni - Auxiliary. */ - su2double *FAux_Dead_Load = nullptr; /*!< \brief Auxiliar vector for the dead loads */ - 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. */ @@ -232,6 +233,8 @@ 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); } /*! @@ -240,8 +243,8 @@ class CFEAElasticity : public CNumerics { * \param[in] jVar - Index j. * \return 1 if i=j, 0 otherwise. */ - inline static su2double deltaij(unsigned short iVar, unsigned short jVar) { - return su2double(iVar==jVar); + inline static passivedouble deltaij(unsigned short iVar, unsigned short jVar) { + return static_cast(iVar == jVar); } }; diff --git a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp index 0e57f54e529..49ba78c09a2 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp @@ -141,8 +141,9 @@ class CFEANonlinearElasticity : public CFEAElasticity { * \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. */ - virtual void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) = 0; + virtual void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) = 0; /*! * \brief Update an element with Maxwell's stress. diff --git a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp b/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp index a36eee0a67d..dbd7b4db422 100644 --- a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp +++ b/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp @@ -73,7 +73,7 @@ class CFEM_NeoHookean_Comp final : public CFEANonlinearElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override; + void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override; }; @@ -124,7 +124,7 @@ class CFEM_Knowles_NearInc final : public CFEANonlinearElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override; + void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override; }; @@ -172,7 +172,7 @@ class CFEM_DielectricElastomer final : public CFEANonlinearElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override; + void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override; }; @@ -222,6 +222,6 @@ class CFEM_IdealDE final : public CFEANonlinearElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override; + 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 b91be37223e..73df569ff26 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -35,7 +35,6 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, nDim = val_nDim; nVar = val_nVar; - bool body_forces = config->GetDeadLoad(); // Body forces (dead loads). bool pseudo_static = config->GetPseudoStatic(); unsigned short iVar; @@ -62,14 +61,11 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, Rho_s = Rho_s_i[0]; Rho_s_DL = Rho_s_DL_i[0]; Alpha = Alpha_i[0]; + ReferenceTemperature = config->GetMaterialReferenceTemperature(); Compute_Lame_Parameters(); - // Auxiliary vector for body forces (dead load) - FAux_Dead_Load = nullptr; - if (body_forces) FAux_Dead_Load = new su2double [nDim]; - - plane_stress = (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRESS); + plane_stress = (nDim == 2) && (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRESS); KAux_ab = new su2double* [nDim]; for (iVar = 0; iVar < nDim; iVar++) { @@ -162,12 +158,11 @@ CFEAElasticity::~CFEAElasticity() { delete[] DV_Val; - delete [] FAux_Dead_Load; - delete [] E_i; delete [] Nu_i; delete [] Rho_s_i; delete [] Rho_s_DL_i; + delete [] Alpha_i; delete [] Ni_Vec; } @@ -242,43 +237,32 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) AD::SetPreaccIn(Rho_s_DL); element->SetPreaccIn_Coords(false); - unsigned short iGauss, nGauss; - unsigned short iNode, iDim, nNode; - - su2double Weight, Jac_X; - /* -- Gravity directionality: * -- For 2D problems, we assume the direction for gravity is -y * -- For 3D problems, we assume the direction for gravity is -z */ su2double g_force[3] = {0.0,0.0,0.0}; - - if (nDim == 2) g_force[1] = -1*STANDARD_GRAVITY; - else if (nDim == 3) g_force[2] = -1*STANDARD_GRAVITY; + g_force[nDim - 1] = -STANDARD_GRAVITY; element->ClearElement(); /*--- Restart the element to avoid adding over previous results. --*/ element->ComputeGrad_Linear(); /*--- Need to compute the gradients to obtain the Jacobian. ---*/ - nNode = element->GetnNodes(); - nGauss = element->GetnGaussPoints(); - - for (iGauss = 0; iGauss < nGauss; iGauss++) { + const auto nNode = element->GetnNodes(); + const auto nGauss = element->GetnGaussPoints(); - Weight = element->GetWeight(iGauss); - Jac_X = element->GetJ_X(iGauss); /*--- The dead load is computed in the reference configuration ---*/ + for (auto iGauss = 0u; iGauss < nGauss; iGauss++) { - /*--- Retrieve the values of the shape functions for each node ---*/ - /*--- This avoids repeated operations ---*/ - for (iNode = 0; iNode < nNode; iNode++) { - Ni_Vec[iNode] = element->GetNi(iNode,iGauss); - } + const auto Weight = element->GetWeight(iGauss); + /*--- The dead load is computed in the reference configuration ---*/ + const auto Jac_X = element->GetJ_X(iGauss); - for (iNode = 0; iNode < nNode; iNode++) { + for (auto iNode = 0; iNode < nNode; iNode++) { + const auto Ni = element->GetNi(iNode,iGauss); - for (iDim = 0; iDim < nDim; iDim++) { - FAux_Dead_Load[iDim] = Weight * Ni_Vec[iNode] * Jac_X * Rho_s_DL * g_force[iDim]; + su2double FAux_Dead_Load[3] = {0.0,0.0,0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + FAux_Dead_Load[iDim] = Weight * Ni * Jac_X * Rho_s_DL * g_force[iDim]; } - element->Add_FDL_a(iNode, FAux_Dead_Load); } diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index de7ad86219f..8e15ed0cf6b 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -97,10 +97,11 @@ void CFEALinearElasticity::Compute_Tangent_Matrix(CElement *element, const CConf } } + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); + for (iNode = 0; iNode < nNode; iNode++) { su2double KAux_t_a[3] = {0.0}; - const su2double thermalStress = -E * 2e-4 / (1 - 2 * Nu); for (iVar = 0; iVar < nDim; iVar++) { KAux_t_a[iVar] += Weight * thermalStress * GradNi_Ref_Mat[iNode][iVar] * Jac_X; } @@ -321,14 +322,17 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, } - /*--- Compute the Stress Vector as D*epsilon ---*/ + /*--- Compute the Stress Vector as D*epsilon + thermal stress ---*/ su2double Stress[DIM_STRAIN_3D] = {0.0}; + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); + for (iVar = 0; iVar < bDim; iVar++) { for (jVar = 0; jVar < bDim; jVar++) { Stress[iVar] += D_Mat[iVar][jVar]*Strain[jVar]; } + if (iVar < nDim) Stress[iVar] += thermalStress; avgStress[iVar] += Stress[iVar] / nGauss; } @@ -370,10 +374,10 @@ CFEAMeshElasticity::CFEAMeshElasticity(unsigned short val_nDim, unsigned short v unsigned long val_nElem, const CConfig *config) : CFEALinearElasticity() { DV_Val = nullptr; - FAux_Dead_Load = nullptr; Rho_s_i = nullptr; Rho_s_DL_i = nullptr; Nu_i = nullptr; + Alpha_i = nullptr; nDim = val_nDim; nVar = val_nVar; diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index 0b8ed7cb855..8f74c164dd6 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp @@ -349,7 +349,7 @@ void CFEANonlinearElasticity::Compute_Tangent_Matrix(CElement *element, const CC /*--- Compute the constitutive matrix ---*/ - Compute_Stress_Tensor(element, config); + Compute_Stress_Tensor(element, config, iGauss); // if (maxwell_stress) Add_MaxwellStress(element, config); Compute_Constitutive_Matrix(element, config); @@ -571,7 +571,7 @@ void CFEANonlinearElasticity::Compute_NodalStress_Term(CElement *element, const /*--- Compute the stress tensor ---*/ - Compute_Stress_Tensor(element, config); + Compute_Stress_Tensor(element, config, iGauss); // if (maxwell_stress) Add_MaxwellStress(element, config); for (iNode = 0; iNode < nNode; iNode++) { @@ -850,7 +850,7 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen /*--- Compute the stress tensor ---*/ - Compute_Stress_Tensor(element, config); + Compute_Stress_Tensor(element, config, iGauss); if (maxwell_stress) Add_MaxwellStress(element, config); avgStress[0] += Stress_Tensor[0][0] / nGauss; diff --git a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp index 55d48a29010..4b787db3416 100644 --- a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp +++ b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp @@ -98,7 +98,7 @@ void CFEM_NeoHookean_Comp::Compute_Constitutive_Matrix(CElement *element, const } -void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfig *config) { +void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) { unsigned short iVar,jVar; su2double Mu_J = 0.0, Lambda_J = 0.0; @@ -109,10 +109,12 @@ void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfi Lambda_J = Lambda/J_F; } + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); + for (iVar = 0; iVar < 3; iVar++) { for (jVar = 0; jVar < 3; jVar++) { su2double dij = deltaij(iVar,jVar); - Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) - E * 2e-4 / (1 - 2 * Nu)) * dij; + Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) + thermalStress) * dij; } } @@ -182,7 +184,7 @@ void CFEM_Knowles_NearInc::Compute_Constitutive_Matrix(CElement *element, const } -void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfig *config) { +void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) { /* -- Suchocki (2011) (full reference in class constructor). ---*/ @@ -199,10 +201,12 @@ void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfi Ek = Kappa * (2.0 * J_F - 1.0); Pr = Kappa * (J_F - 1.0); + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); + for (iVar = 0; iVar < 3; iVar++){ for (jVar = 0; jVar < 3; jVar++){ Stress_Tensor[iVar][jVar] = term1 * (b_Mat_Iso[iVar][jVar] - deltaij(iVar,jVar)*trbbar) + - deltaij(iVar,jVar) * Pr; + deltaij(iVar,jVar) * (Pr + thermalStress); } } @@ -234,7 +238,7 @@ void CFEM_DielectricElastomer::Compute_Constitutive_Matrix(CElement *element, co } -void CFEM_DielectricElastomer::Compute_Stress_Tensor(CElement *element, const CConfig *config) { +void CFEM_DielectricElastomer::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) { unsigned short iDim, jDim; @@ -315,12 +319,11 @@ void CFEM_IdealDE::Compute_Constitutive_Matrix(CElement *element, const CConfig } -void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *config) { +void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) { /* -- Zhao, X. and Suo, Z. (2008) (full reference in class constructor). ---*/ unsigned short iVar, jVar; - su2double dij = 0.0; /*--- Compute the isochoric deformation gradient Fbar and left Cauchy-Green tensor bbar ---*/ Compute_Isochoric_F_b(); @@ -333,13 +336,12 @@ void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *confi Pr = Kappa * (J_F - 1.0); Eg23 = 2.0 * Eg / 3.0; - // Stress tensor + const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature); for (iVar = 0; iVar < 3; iVar++){ - for (jVar = 0; jVar < 3; jVar++){ - if (iVar == jVar) dij = 1.0; - else if (iVar != jVar) dij = 0.0; - Stress_Tensor[iVar][jVar] = Eg * ( b_Mat_Iso[iVar][jVar] - dij * trbbar) + dij * Pr ; + for (jVar = 0; jVar < 3; jVar++) { + const su2double dij = deltaij(iVar, jVar); + Stress_Tensor[iVar][jVar] = Eg * ( b_Mat_Iso[iVar][jVar] - dij * trbbar) + dij * (Pr + thermalStress); } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 4cd718d6ce1..44a19e1421a 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -71,6 +71,10 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge element_container[FEA_TERM][EL_TRIA+offset] = new CTRIA1(); element_container[FEA_TERM][EL_QUAD+offset] = new CQUAD4(); + /*--- Initialize temperature ---*/ + element_container[FEA_TERM][EL_TRIA+offset]->SetTemperature(config->GetTemperature_FreeStream()); + element_container[FEA_TERM][EL_QUAD+offset]->SetTemperature(config->GetTemperature_FreeStream()); + if (de_effects) { element_container[DE_TERM][EL_TRIA+offset] = new CTRIA1(); element_container[DE_TERM][EL_QUAD+offset] = new CQUAD4(); @@ -82,6 +86,10 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge element_container[FEA_TERM][EL_PYRAM+offset] = new CPYRAM5(); element_container[FEA_TERM][EL_PRISM+offset] = new CPRISM6(); + for (const auto el : {EL_TETRA, EL_HEXA, EL_PYRAM, EL_PRISM}) { + element_container[FEA_TERM][el+offset]->SetTemperature(config->GetTemperature_FreeStream()); + } + if (de_effects) { element_container[DE_TERM][EL_TETRA+offset] = new CTETRA1(); element_container[DE_TERM][EL_HEXA +offset] = new CHEXA8 (); diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg index 128958ecb19..b08a38181d6 100644 --- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg @@ -14,13 +14,15 @@ MATERIAL_MODEL= LINEAR_ELASTIC MESH_FILENAME= meshBeam_3d.su2 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 +MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5 +MATERIAL_REFERENCE_TEMPERATURE= 288.15 +FREESTREAM_TEMPERATURE= 350 MATERIAL_DENSITY=7854 MARKER_CLAMPED = ( left ) MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0) -%MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) LINEAR_SOLVER= CONJUGATE_GRADIENT LINEAR_SOLVER_PREC= ILU -LINEAR_SOLVER_ERROR= 1E-9 +LINEAR_SOLVER_ERROR= 1E-8 LINEAR_SOLVER_ITER= 1000 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV @@ -29,4 +31,4 @@ VOLUME_FILENAME= beam RESTART_FILENAME= restart_beam.dat SOLUTION_FILENAME= restart_beam.dat OUTPUT_WRT_FREQ= 1 -INNER_ITER=10 +INNER_ITER=1 From bac0e067fe3169496b1484346be42d0d2177497a Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 22 Dec 2024 18:15:30 -0800 Subject: [PATCH 3/7] update testcase --- TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg | 5 +++-- config_template.cfg | 6 ++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg index b08a38181d6..56b57777a55 100644 --- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg @@ -18,8 +18,9 @@ MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5 MATERIAL_REFERENCE_TEMPERATURE= 288.15 FREESTREAM_TEMPERATURE= 350 MATERIAL_DENSITY=7854 -MARKER_CLAMPED = ( left ) -MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0) +MARKER_CLAMPED = ( left, right ) +MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0 ) +MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 ) LINEAR_SOLVER= CONJUGATE_GRADIENT LINEAR_SOLVER_PREC= ILU LINEAR_SOLVER_ERROR= 1E-8 diff --git a/config_template.cfg b/config_template.cfg index 8a953b7be78..2ecb8ca03b2 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1727,6 +1727,12 @@ ELASTICITY_MODULUS= 1000.0 % Poisson ratio POISSON_RATIO= 0.35 % +% Thermal expansion coefficient +MATERIAL_THERMAL_EXPANSION_COEFF= 0 +% +% Freestream temperature at which there is no stress from thermal expansion +MATERIAL_REFERENCE_TEMPERATURE= 288.15 +% % Knowles B constant KNOWLES_B= 1.0 % From f0808bd7dc0a7109b6f823ce1b854b0c18d73bad Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 23 Dec 2024 18:28:26 -0800 Subject: [PATCH 4/7] update cases --- Common/src/CConfig.cpp | 5 ++++- TestCases/hybrid_regression.py | 4 ++-- TestCases/parallel_regression.py | 4 ++-- TestCases/serial_regression.py | 4 ++-- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 5c220686271..ba2f975c699 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -914,7 +914,10 @@ void CConfig::SetPointersNull() { Inlet_Velocity = nullptr; Inflow_Mach = nullptr; Inflow_Pressure = nullptr; Outlet_Pressure = nullptr; Isothermal_Temperature = nullptr; - ElasticityMod = nullptr; PoissonRatio = nullptr; MaterialDensity = nullptr; + ElasticityMod = nullptr; + PoissonRatio = nullptr; + MaterialDensity = nullptr; + MaterialThermalExpansion = nullptr; Load_Dir = nullptr; Load_Dir_Value = nullptr; Load_Dir_Multiplier = nullptr; Disp_Dir = nullptr; Disp_Dir_Value = nullptr; Disp_Dir_Multiplier = nullptr; diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index e5b0a847e61..fb05c2513db 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -680,8 +680,8 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-2.378370, -1.585252, -2.028505, 64359.000000] - statbeam3d.test_vals_aarch64 = [-2.382650, -1.561882, -2.045083, 64433.000000] + statbeam3d.test_vals = [-2.787802, -1.721974, -2.438436, 110350] + statbeam3d.test_vals_aarch64 = [-2.787802, -1.721974, -2.438436, 110350] test_list.append(statbeam3d) # Dynamic beam, 2d diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index f19b323f761..4a86a6fcc00 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1220,8 +1220,8 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-8.396797, -8.162206, -8.156102, 64095.000000] - statbeam3d.test_vals_aarch64 = [-8.396793, -8.162255, -8.156118, 64095.0] #last 4 columns + statbeam3d.test_vals = [-6.058758, -5.750933, -5.892188, 110190] + statbeam3d.test_vals_aarch64 = [-6.058758, -5.750933, -5.892188, 110190] statbeam3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") test_list.append(statbeam3d) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 4aa5f883481..a0d531353f6 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1012,8 +1012,8 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-8.498245, -8.230816, -8.123810, 64095.000000] - statbeam3d.test_vals_aarch64 = [-8.498254, -8.230683, -8.123819, 64095.0] #last 4 columns + statbeam3d.test_vals = [-6.168640, -5.939035, -6.071159, 110190] + statbeam3d.test_vals_aarch64 = [-6.168640, -5.939035, -6.071159, 110190] #last 4 columns test_list.append(statbeam3d) # Mix elem, 3d beam, Knowles From 2cda5ec8ff45ea99466ceb59187f031cba654b41 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 24 Dec 2024 14:10:28 -0800 Subject: [PATCH 5/7] generalize body force for FEA solver and add centrifugal loads --- Common/include/CConfig.hpp | 7 ++- Common/src/CConfig.cpp | 10 ++-- SU2_CFD/include/numerics/CNumerics.hpp | 4 +- .../numerics/elasticity/CFEAElasticity.hpp | 6 +-- SU2_CFD/include/solvers/CFEASolver.hpp | 9 ++-- SU2_CFD/include/solvers/CSolver.hpp | 6 +-- .../numerics/elasticity/CFEAElasticity.cpp | 51 ++++++++++++++----- SU2_CFD/src/output/CAdjElasticityOutput.cpp | 4 +- SU2_CFD/src/solvers/CFEASolver.cpp | 26 +++------- SU2_CFD/src/variables/CFEAVariable.cpp | 2 +- TestCases/disc_adj_fea/configAD_fem.cfg | 3 +- .../flow_load_sens/configAD_fem.cfg | 1 - config_template.cfg | 5 +- 13 files changed, 73 insertions(+), 61 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index d91a6011611..b5407e71eeb 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1004,7 +1004,7 @@ class CConfig { bool ExtraOutput; /*!< \brief Check if extra output need. */ bool Wall_Functions; /*!< \brief Use wall functions with the turbulence model */ long ExtraHeatOutputZone; /*!< \brief Heat solver zone with extra screen output */ - bool DeadLoad; /*!< \brief Application of dead loads to the FE analysis */ + bool CentrifugalForce; /*!< \brief Application of centrifugal forces to the FE analysis */ bool PseudoStatic; /*!< \brief Application of dead loads to the FE analysis */ bool SteadyRestart; /*!< \brief Restart from a steady state for FSI problems. */ su2double Newmark_beta, /*!< \brief Parameter alpha for Newmark method. */ @@ -8942,10 +8942,9 @@ class CConfig { su2double GetAitkenDynMinInit(void) const { return AitkenDynMinInit; } /*! - * \brief Decide whether to apply dead loads to the model. - * \return TRUE if the dead loads are to be applied, FALSE otherwise. + * \brief Decide whether to apply centrifugal forces to the model. */ - bool GetDeadLoad(void) const { return DeadLoad; } + bool GetCentrifugalForce(void) const { return CentrifugalForce; } /*! * \brief Identifies if the mesh is matching or not (temporary, while implementing interpolation procedures). diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index ba2f975c699..757388c9cc8 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2507,11 +2507,11 @@ void CConfig::SetConfig_Options() { addEnumOption("NONLINEAR_FEM_SOLUTION_METHOD", Kind_SpaceIteScheme_FEA, Space_Ite_Map_FEA, STRUCT_SPACE_ITE::NEWTON); /* DESCRIPTION: Formulation for bidimensional elasticity solver */ addEnumOption("FORMULATION_ELASTICITY_2D", Kind_2DElasForm, ElasForm_2D, STRUCT_2DFORM::PLANE_STRAIN); - /* DESCRIPTION: Apply dead loads - * Options: NO, YES \ingroup Config */ - addBoolOption("DEAD_LOAD", DeadLoad, false); + /* DESCRIPTION: Apply centrifugal forces + * Options: NO, YES \ingroup Config */ + addBoolOption("CENTRIFUGAL_FORCE", CentrifugalForce, false); /* DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis) - * Options: NO, YES \ingroup Config */ + * Options: NO, YES \ingroup Config */ addBoolOption("PSEUDO_STATIC", PseudoStatic, false); /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */ addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25); @@ -3070,6 +3070,8 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){ newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN instead.\n\n"); else if (!option_name.compare("SPECIES_USE_STRONG_BC")) newString.append("SPECIES_USE_STRONG_BC is deprecated. Use MARKER_SPECIES_STRONG_BC= (marker1, ...) instead.\n\n"); + else if (!option_name.compare("DEAD_LOAD")) + newString.append("DEAD_LOAD is deprecated. Use GRAVITY_FORCE or BODY_FORCE instead.\n\n"); else { /*--- Find the most likely candidate for the unrecognized option, based on the length of start and end character sequences shared by candidates and the option. ---*/ diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 535335947ca..2e45916edb1 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -1483,10 +1483,10 @@ class CNumerics { inline virtual void Compute_Mass_Matrix(CElement *element_container, const CConfig* config) { } /*! - * \brief A virtual member to compute the residual component due to dead loads + * \brief A virtual member to compute the residual component due to body forces. * \param[in] element_container - Element structure for the particular element integrated. */ - inline virtual void Compute_Dead_Load(CElement *element_container, const CConfig* config) { } + inline virtual void Compute_Body_Forces(CElement *element_container, const CConfig* config) { } /*! * \brief A virtual member to compute the averaged nodal stresses diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index bbaa805e046..a1395fdf130 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -159,11 +159,11 @@ class CFEAElasticity : public CNumerics { void Compute_Mass_Matrix(CElement *element_container, const CConfig *config) final; /*! - * \brief Compute the nodal gravity loads for an element. - * \param[in,out] element_container - The element for which the dead loads are computed. + * \brief Compute the nodal inertial loads for an element. + * \param[in,out] element_container - The element for which the inertial loads are computed. * \param[in] config - Definition of the problem. */ - void Compute_Dead_Load(CElement *element_container, const CConfig *config) final; + void Compute_Body_Forces(CElement *element_container, const CConfig *config) final; /*! * \brief Build the tangent stiffness matrix of an element. diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 4d8bf74be8e..9442fa190d7 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -98,6 +98,7 @@ class CFEASolver : public CFEASolverBase { bool element_based; /*!< \brief Bool to determine if an element-based file is used. */ bool topol_filter_applied; /*!< \brief True if density filtering has been performed. */ bool initial_calc = true; /*!< \brief Becomes false after first call to Preprocessing. */ + bool body_forces = false; /*!< \brief Whether any body force is active. */ /*! * \brief The highest level in the variable hierarchy this solver can safely use, @@ -335,14 +336,14 @@ class CFEASolver : public CFEASolverBase { const CConfig *config); /*! - * \brief Compute the dead loads. + * \brief Compute the inertial loads. * \param[in] geometry - Geometrical definition of the problem. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ - void Compute_DeadLoad(CGeometry *geometry, - CNumerics **numerics, - const CConfig *config) final; + void Compute_BodyForces(CGeometry *geometry, + CNumerics **numerics, + const CConfig *config) final; /*! * \brief Clamped boundary conditions. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index b6c0d45d219..2ed78a78a9b 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3645,9 +3645,9 @@ class CSolver { * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ - inline virtual void Compute_DeadLoad(CGeometry *geometry, - CNumerics **numerics, - const CConfig *config) { } + inline virtual void Compute_BodyForces(CGeometry *geometry, + CNumerics **numerics, + const CConfig *config) { } /*! * \brief A virtual member. Set the volumetric heat source diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 73df569ff26..ca3da2b26ae 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -27,6 +27,7 @@ #include "../../../include/numerics/elasticity/CFEAElasticity.hpp" #include "../../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../../Common/include/toolboxes/geometry_toolbox.hpp" CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, @@ -226,7 +227,11 @@ void CFEAElasticity::Compute_Mass_Matrix(CElement *element, const CConfig *confi } -void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) { +void CFEAElasticity::Compute_Body_Forces(CElement *element, const CConfig *config) { + + const bool gravity = config->GetGravityForce(); + const bool body_force = config->GetBody_Force(); + const bool centrifugal = config->GetCentrifugalForce(); /*--- Initialize values for the material model considered ---*/ SetElement_Properties(element, config); @@ -237,13 +242,24 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) AD::SetPreaccIn(Rho_s_DL); element->SetPreaccIn_Coords(false); - /* -- Gravity directionality: - * -- For 2D problems, we assume the direction for gravity is -y - * -- For 3D problems, we assume the direction for gravity is -z - */ - su2double g_force[3] = {0.0,0.0,0.0}; - g_force[nDim - 1] = -STANDARD_GRAVITY; + std::array acceleration{}; + if (gravity) { + /*--- For 2D problems, we assume gravity is in the -y direction, + * and for 3D problems in the -z direction. ---*/ + acceleration[nDim - 1] = -STANDARD_GRAVITY; + } else if (body_force) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + acceleration[iDim] = config->GetBody_Force_Vector()[iDim]; + } + } + std::array center{}, omega{}; + if (centrifugal) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + center[iDim] = config->GetMotion_Origin(iDim); + omega[iDim] = config->GetRotation_Rate(iDim); + } + } element->ClearElement(); /*--- Restart the element to avoid adding over previous results. --*/ element->ComputeGrad_Linear(); /*--- Need to compute the gradients to obtain the Jacobian. ---*/ @@ -256,15 +272,26 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) /*--- The dead load is computed in the reference configuration ---*/ const auto Jac_X = element->GetJ_X(iGauss); - for (auto iNode = 0; iNode < nNode; iNode++) { + for (auto iNode = 0u; iNode < nNode; iNode++) { const auto Ni = element->GetNi(iNode,iGauss); - su2double FAux_Dead_Load[3] = {0.0,0.0,0.0}; + auto total_accel = acceleration; + if (centrifugal) { + /*--- For nonlinear this should probably use the current (deformed) + * coordinates, but it should not make a big difference. ---*/ + std::array r{}, wr{}, w2r{}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + r[iDim] = element->GetRef_Coord(iNode, iDim) - center[iDim]; + } + GeometryToolbox::CrossProduct(omega.data(), r.data(), wr.data()); + GeometryToolbox::CrossProduct(omega.data(), wr.data(), w2r.data()); + for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] += w2r[iDim]; + } + std::array force{}; for (auto iDim = 0u; iDim < nDim; iDim++) { - FAux_Dead_Load[iDim] = Weight * Ni * Jac_X * Rho_s_DL * g_force[iDim]; + force[iDim] = Weight * Ni * Jac_X * Rho_s_DL * total_accel[iDim]; } - element->Add_FDL_a(iNode, FAux_Dead_Load); - + element->Add_FDL_a(iNode, force.data()); } } diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index 7466b506d60..add084006ea 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -113,7 +113,7 @@ void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){ if (config->GetTime_Domain() && !config->GetPseudoStatic()) { AddHistoryOutput("SENS_RHO_" + iVarS, "Sens[Rho" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Material density"); } - if (config->GetDeadLoad()) { + if (config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce()) { AddHistoryOutput("SENS_RHO_DL_" + iVarS, "Sens[RhoDL" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Dead load density"); } } @@ -151,7 +151,7 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge if (config->GetTime_Domain() && !config->GetPseudoStatic()) { SetHistoryOutputValue("SENS_RHO_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar)); } - if (config->GetDeadLoad()) { + if (config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce()) { SetHistoryOutputValue("SENS_RHO_DL_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho_DL(iVar)); } } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 44a19e1421a..23bad36a8a4 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -49,7 +49,7 @@ CFEASolver::CFEASolver(LINEAR_SOLVER_MODE mesh_deform_mode) : CFEASolverBase(mes CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(geometry, config) { - bool dynamic = (config->GetTime_Domain()); + bool dynamic = config->GetTime_Domain(); config->SetDelta_UnstTimeND(config->GetDelta_UnstTime()); /*--- Test whether we consider dielectric elastomers ---*/ @@ -59,6 +59,7 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge element_based = false; topol_filter_applied = false; initial_calc = true; + body_forces = config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce(); /*--- Here is where we assign the kind of each element ---*/ @@ -121,22 +122,15 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge /*--- The length of the solution vector depends on whether the problem is static or dynamic ---*/ - unsigned short nSolVar; string text_line, filename; ifstream restart_file; - if (dynamic) nSolVar = 3 * nVar; - else nSolVar = nVar; - - auto* SolInit = new su2double[nSolVar](); - /*--- Initialize from zero everywhere ---*/ - nodes = new CFEABoundVariable(SolInit, nPoint, nDim, nVar, config); + std::array zeros{}; + nodes = new CFEABoundVariable(zeros.data(), nPoint, nDim, nVar, config); SetBaseClassPointerToNodes(); - delete [] SolInit; - /*--- Set which points are vertices and allocate boundary data. ---*/ for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) @@ -567,7 +561,6 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, const bool dynamic = config->GetTime_Domain(); const bool disc_adj_fem = (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_FEM); - const bool body_forces = config->GetDeadLoad(); const bool topology_mode = config->GetTopology_Optimization(); /* @@ -602,7 +595,7 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, * Only initialized once, at the first iteration of the first time step. */ if (body_forces && (initial_calc || disc_adj_fem)) - Compute_DeadLoad(geometry, numerics, config); + Compute_BodyForces(geometry, numerics, config); /*--- Clear the linear system solution. ---*/ SU2_OMP_PARALLEL @@ -1415,7 +1408,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } -void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { +void CFEASolver::Compute_BodyForces(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { /*--- Start OpenMP parallel region. ---*/ @@ -1465,7 +1458,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con /*--- Set the properties of the element and compute its mass matrix. ---*/ element->Set_ElProperties(element_properties[iElem]); - numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Dead_Load(element, config); + numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Body_Forces(element, config); /*--- Add contributions of this element to the mass matrix. ---*/ for (iNode = 0; iNode < nNodes; iNode++) { @@ -2113,7 +2106,6 @@ void CFEASolver::ImplicitNewmark_Iteration(const CGeometry *geometry, CNumerics const bool linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL); const bool nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE); const bool newton_raphson = (config->GetKind_SpaceIteScheme_FEA() == STRUCT_SPACE_ITE::NEWTON); - const bool body_forces = config->GetDeadLoad(); /*--- For simplicity, no incremental loading is handled with increment of 1. ---*/ const su2double loadIncr = config->GetIncrementalLoad()? loadIncrement : su2double(1.0); @@ -2301,7 +2293,6 @@ void CFEASolver::GeneralizedAlpha_Iteration(const CGeometry *geometry, CNumerics const bool linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL); const bool nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE); const bool newton_raphson = (config->GetKind_SpaceIteScheme_FEA() == STRUCT_SPACE_ITE::NEWTON); - const bool body_forces = config->GetDeadLoad(); /*--- Blend between previous and current timestep. ---*/ const su2double alpha_f = config->Get_Int_Coeffs(2); @@ -2925,9 +2916,6 @@ void CFEASolver::Compute_OFVolFrac(CGeometry *geometry, const CConfig *config) void CFEASolver::Compute_OFCompliance(CGeometry *geometry, const CConfig *config) { - /*--- Types of loads to consider ---*/ - const bool body_forces = config->GetDeadLoad(); - /*--- If the loads are being applied incrementaly ---*/ const bool incremental_load = config->GetIncrementalLoad(); diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp index 7e486af3a21..4760706a814 100644 --- a/SU2_CFD/src/variables/CFEAVariable.cpp +++ b/SU2_CFD/src/variables/CFEAVariable.cpp @@ -46,7 +46,7 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig * still work as expected for the primal solver). ---*/ const bool dynamic_analysis = config->GetTime_Domain(); - const bool body_forces = config->GetDeadLoad(); + const bool body_forces = config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce(); const bool prestretch_fem = config->GetPrestretch(); // Structure is prestretched const bool discrete_adjoint = config->GetDiscrete_Adjoint(); const bool refgeom = config->GetRefGeom(); // Reference geometry needs to be stored diff --git a/TestCases/disc_adj_fea/configAD_fem.cfg b/TestCases/disc_adj_fea/configAD_fem.cfg index 188c7e451bf..fe1d2c6c5b7 100644 --- a/TestCases/disc_adj_fea/configAD_fem.cfg +++ b/TestCases/disc_adj_fea/configAD_fem.cfg @@ -16,7 +16,7 @@ RESTART_SOL= NO ITER = 10 -OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY +OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY REFERENCE_GEOMETRY = YES REFERENCE_GEOMETRY_FILENAME = reference_geometry.dat @@ -39,7 +39,6 @@ MATERIAL_COMPRESSIBILITY= COMPRESSIBLE ELASTICITY_MODULUS=21000 POISSON_RATIO=0.4 MATERIAL_DENSITY=100 -DEAD_LOAD=NO FORMULATION_ELASTICITY_2D = PLANE_STRAIN NONLINEAR_FEM_SOLUTION_METHOD = NEWTON_RAPHSON diff --git a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg index 63edf6869cf..b5d0f15d944 100644 --- a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg +++ b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg @@ -35,7 +35,6 @@ MATERIAL_COMPRESSIBILITY= COMPRESSIBLE ELASTICITY_MODULUS=21000 POISSON_RATIO=0.4 MATERIAL_DENSITY=100 -DEAD_LOAD=NO FORMULATION_ELASTICITY_2D = PLANE_STRAIN NONLINEAR_FEM_SOLUTION_METHOD = NEWTON_RAPHSON diff --git a/config_template.cfg b/config_template.cfg index 2ecb8ca03b2..17891c34373 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -2023,10 +2023,7 @@ NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON % Formulation for bidimensional elasticity solver FORMULATION_ELASTICITY_2D= PLANE_STRAIN % -% Apply dead loads -DEAD_LOAD= NO -% -% pseudo static analysis (no density in dynamic analysis) +% Pseudo static analysis (no density in dynamic analysis) PSEUDO_STATIC= NO % % Dynamic or static structural analysis (deprecated -> use TIME_DOMAIN) From bee055e4f89e2f00f3c32d186ceb6a904a1287ad Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 24 Dec 2024 19:12:01 -0800 Subject: [PATCH 6/7] add verification case --- .../numerics/elasticity/CFEAElasticity.cpp | 2 +- .../fea_fsi/rotating_cylinder/config.cfg | 34 +++++++++++++++++++ TestCases/parallel_regression.py | 12 +++++++ 3 files changed, 47 insertions(+), 1 deletion(-) create mode 100644 TestCases/fea_fsi/rotating_cylinder/config.cfg diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index ca3da2b26ae..efe0e19a0b2 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -285,7 +285,7 @@ void CFEAElasticity::Compute_Body_Forces(CElement *element, const CConfig *confi } GeometryToolbox::CrossProduct(omega.data(), r.data(), wr.data()); GeometryToolbox::CrossProduct(omega.data(), wr.data(), w2r.data()); - for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] += w2r[iDim]; + for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] -= w2r[iDim]; } std::array force{}; for (auto iDim = 0u; iDim < nDim; iDim++) { diff --git a/TestCases/fea_fsi/rotating_cylinder/config.cfg b/TestCases/fea_fsi/rotating_cylinder/config.cfg new file mode 100644 index 00000000000..7a9d124da76 --- /dev/null +++ b/TestCases/fea_fsi/rotating_cylinder/config.cfg @@ -0,0 +1,34 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Spinning Cylinder % +% Author: Pedro Gomes % +% Institution: SU2 Foundation % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +MATH_PROBLEM= DIRECT +GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS +MATERIAL_MODEL= LINEAR_ELASTIC + +ELASTICITY_MODULUS= 2E11 +POISSON_RATIO= 0.3 +MATERIAL_DENSITY= 7854 + +MARKER_SYM= ( x_minus, per_1, per_2 ) +MARKER_PRESSURE= ( inner, 0, outer, 0, x_plus, 0 ) + +CENTRIFUGAL_FORCE= YES +ROTATION_RATE= ( 1500, 0, 0 ) +MOTION_ORIGIN= ( 0, 0, 0 ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 1000 + +SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS +OUTPUT_WRT_FREQ= 1 +INNER_ITER= 1 + +MESH_FILENAME= cylinder.su2 + diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index d60bd28e07e..b0c18f755b3 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1225,6 +1225,18 @@ def main(): statbeam3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") test_list.append(statbeam3d) + # Rotating cylinder, 3d + rotating_cylinder_fea = TestCase('rotating_cylinder_fea') + rotating_cylinder_fea.cfg_dir = "fea_fsi/rotating_cylinder" + rotating_cylinder_fea.cfg_file = "config.cfg" + 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. + rotating_cylinder_fea.test_vals = [-1.986097, -1.023250, -1.022700, 47, -8.101266, 1.6458e+08] + rotating_cylinder_fea.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") + test_list.append(rotating_cylinder_fea) + # Dynamic beam, 2d dynbeam2d = TestCase('dynbeam2d') dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d" From 28fb23ddef1d624dbbe485f122e8cab3123cfa67 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 25 Dec 2024 10:01:45 -0800 Subject: [PATCH 7/7] update values and config template --- TestCases/parallel_regression.py | 3 +- config_template.cfg | 49 +++++++++++++++----------------- 2 files changed, 24 insertions(+), 28 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index b0c18f755b3..8bc45cb6a8b 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1233,8 +1233,7 @@ def main(): # 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. - rotating_cylinder_fea.test_vals = [-1.986097, -1.023250, -1.022700, 47, -8.101266, 1.6458e+08] - rotating_cylinder_fea.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f") + rotating_cylinder_fea.test_vals = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8] test_list.append(rotating_cylinder_fea) # Dynamic beam, 2d diff --git a/config_template.cfg b/config_template.cfg index 17891c34373..49cca43332b 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1742,12 +1742,17 @@ KNOWLES_N= 1.0 % ID of the region we want to compute the sensitivities using direct differentiation FEA_ID_DIRECTDIFF= 0 % -% RESTART_STEADY_STATE= NO % % Time discretization TIME_DISCRE_FEA= NEWMARK_IMPLICIT % +% Iterative method for non-linear structural analysis +NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON +% +% Pseudo static analysis (no density in dynamic analysis) +PSEUDO_STATIC= NO +% % Parameter alpha for Newmark scheme (s) NEWMARK_BETA= 0.25 % @@ -1781,6 +1786,9 @@ MATERIAL_MODEL= LINEAR_ELASTIC % Compressibility of the material MATERIAL_COMPRESSIBILITY= COMPRESSIBLE % +% Formulation for 2-dimensional elasticity solver +FORMULATION_ELASTICITY_2D= PLANE_STRAIN +% % -------------------- Dielectric effects ------------------% % % Include DE effects @@ -1794,6 +1802,20 @@ ELECTRIC_FIELD_MOD= 20e5 % % Direction of the Electic Fields ELECTRIC_FIELD_DIR= (0.0, 1.0) +% +% ------------------------ Prestretch -----------------------% +% +% Consider a prestretch in the structural domain +PRESTRETCH= NO +% +% Filename to input for prestretching membranes +PRESTRETCH_FILENAME= prestretch_file.dat +% +% ----------------------- Body Forces -----------------------% +% +% Centrifugal forces due to ROTATION_RATE around MOTION_ORIGIN. +% GRAVITY_FORCE and BODY_FORCE can also be used with the FEA solver. +CENTRIFUGAL_FORCE= NO % -------------------- Weakly Coupled Heat ------------------% % @@ -2010,31 +2032,6 @@ STRESS_PENALTY_PARAM= (1.0, 10.0) % Preaccumulation in the AD mode. PREACC= YES -% ---------------- PRESTRETCH FOR STRUCTURES -------------------% -% Consider a prestretch in the structural domain -PRESTRETCH= NO -% -% Filename to input for prestretching membranes -PRESTRETCH_FILENAME= prestretch_file.dat -% -% Iterative method for non-linear structural analysis -NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON -% -% Formulation for bidimensional elasticity solver -FORMULATION_ELASTICITY_2D= PLANE_STRAIN -% -% Pseudo static analysis (no density in dynamic analysis) -PSEUDO_STATIC= NO -% -% Dynamic or static structural analysis (deprecated -> use TIME_DOMAIN) -DYNAMIC_ANALYSIS= NO -% -% Time Step for dynamic analysis (s) (deprecated -> use TIME_STEP) -DYN_TIMESTEP= 0.0 -% -% Total Physical Time for dual time stepping simulations (s) (deprecated -> use MAX_TIME) -DYN_TIME= 1.0 - % ---------------- MESH DEFORMATION PARAMETERS (NEW SOLVER) -------------------% % % Use the reformatted pseudo-elastic solver for grid deformation