Add thermal expansion effects to FEA solver

Common/include/CConfig.hpp

@@ -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.

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<su2activematrix> 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);

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 */
@@ -4834,9 +4841,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) {

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);

Common/src/geometry/elements/CQUAD4.cpp

@@ -67,7 +67,6 @@ CQUAD4::CQUAD4() : CElementWithKnownSizes<NGAUSS, NNODE, NDIM>() {
   /*--- Store the extrapolation functions (used to compute nodal stresses) ---*/
 
   su2double ExtrapCoord[4][2], sqrt3 = 1.732050807568877;
-  ;
 
   ExtrapCoord[0][0] = -sqrt3;
   ExtrapCoord[0][1] = -sqrt3;

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. */
 
@@ -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<passivedouble>(iVar == jVar);
   }
 
 };

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.

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;
 
 };