thermal stress term, not coupled to heat solver

Author: pcarruscag

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.

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

@@ -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 */

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

@@ -59,13 +59,16 @@ 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. */
   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. */
   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<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;
 
 };

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

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;