diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 0c32f044288..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. */
@@ -1013,10 +1013,13 @@ 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. */
+ 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 . */
@@ -2389,6 +2392,16 @@ 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 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.
@@ -8929,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/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 c5bd653b04a..757388c9cc8 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;
@@ -2440,6 +2443,10 @@ 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: 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 */
@@ -2500,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);
@@ -3063,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. ---*/
@@ -4834,9 +4843,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/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/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 11470fd0fba..a1395fdf130 100644
--- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
+++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
@@ -54,15 +54,20 @@ 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. */
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. */
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 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. */
@@ -72,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. */
@@ -156,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.
@@ -230,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);
}
/*!
@@ -238,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/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 44014591204..efe0e19a0b2 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,
@@ -35,7 +36,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;
@@ -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,14 +61,12 @@ 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];
+ 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++) {
@@ -164,12 +159,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;
}
@@ -233,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);
@@ -244,45 +242,56 @@ 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;
+ 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. ---*/
- 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 = 0u; 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];
+ 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];
}
-
- element->Add_FDL_a(iNode, FAux_Dead_Load);
-
+ std::array force{};
+ for (auto iDim = 0u; iDim < nDim; iDim++) {
+ force[iDim] = Weight * Ni * Jac_X * Rho_s_DL * total_accel[iDim];
+ }
+ element->Add_FDL_a(iNode, force.data());
}
}
@@ -302,6 +311,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..8e15ed0cf6b 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
@@ -97,8 +97,16 @@ 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};
+ 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];
@@ -314,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;
}
@@ -363,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 d21fc228300..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) * 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/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 4cd718d6ce1..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 ---*/
@@ -71,6 +72,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 +87,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 ();
@@ -113,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++)
@@ -559,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();
/*
@@ -594,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
@@ -1407,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. ---*/
@@ -1457,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++) {
@@ -2105,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);
@@ -2293,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);
@@ -2917,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/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
index a6c9cde3c84..56b57777a55 100644
--- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
+++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
@@ -14,14 +14,17 @@ 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 , 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
+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
-LINEAR_SOLVER_ITER= 500
+LINEAR_SOLVER_ITER= 1000
MESH_FORMAT= SU2
TABULAR_FORMAT= CSV
CONV_FILENAME= history_beam
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/hybrid_regression.py b/TestCases/hybrid_regression.py
index 80c709f8d67..da3af4f5c9d 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 3ff60005cee..8bc45cb6a8b 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1220,11 +1220,22 @@ 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)
+ # 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 = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8]
+ test_list.append(rotating_cylinder_fea)
+
# Dynamic beam, 2d
dynbeam2d = TestCase('dynbeam2d')
dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d"
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/TestCases/serial_regression.py b/TestCases/serial_regression.py
index d37d6562534..b4b9cd473a9 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
diff --git a/config_template.cfg b/config_template.cfg
index 8a953b7be78..49cca43332b 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
%
@@ -1736,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
%
@@ -1775,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
@@ -1788,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 ------------------%
%
@@ -2004,34 +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
-%
-% Apply dead loads
-DEAD_LOAD= NO
-%
-% 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