add thermal expansion to CConfig

Author: pcarruscag

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.

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

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

SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp

@@ -44,26 +44,24 @@ 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
   E   = E_i[0];
   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:

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

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

TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg

@@ -15,18 +15,18 @@ 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
 VOLUME_FILENAME= beam
 RESTART_FILENAME= restart_beam.dat
 SOLUTION_FILENAME= restart_beam.dat
 OUTPUT_WRT_FREQ= 1
-INNER_ITER=1
+INNER_ITER=10