testcase and handle non dimensional mode

Author: pcarruscag

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -101,13 +101,9 @@ class CFEASolver : public CFEASolverBase {
   bool body_forces = false;    /*!< \brief Whether any body force is active. */
 
   /*!
-   * \brief Pointer to the heat solver for coupled simulations.
-   *
-   * This member stores a pointer to the heat solver, which handles the solution of the heat equation in weakly coupled
-   * simulations. It is initialized during the preprocessing step if the configuration enables the weak coupling of
-   * heat and elasticity solvers. This solver provides temperature information to the finite element elasticity solver.
+   * \brief Pointer to the heat solver nodes to access temperature for coupled simulations.
    */
-  CSolver* heat_solver = nullptr;
+  const CVariable* heat_nodes = nullptr;
 
   /*!
    * \brief The highest level in the variable hierarchy this solver can safely use,

SU2_CFD/include/solvers/CHeatSolver.hpp

@@ -42,7 +42,6 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
   static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */
 
   const bool flow; /*!< \brief Use solver as a scalar transport equation of Temperature for the inc solver. */
-  const bool heat_equation; /*!< \brief use solver for heat conduction in solids. */
 
   su2double Global_Delta_Time = 0.0, Global_Delta_UnstTimeND = 0.0;
 

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -503,10 +503,10 @@ void CScalarSolver<VariableType>::CompleteImplicitIteration(CGeometry* geometry,
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
         /*--- Multiply the Solution var with density to get the conservative transported quantity, if necessary. ---*/
-        /* Note that for consistency with residual and jacobian calaulcations, use of current density for conservative variables 
+        /* Note that for consistency with residual and jacobian calaulcations, use of current density for conservative variables
         * of the old solution is used. see pull request https://github.com/su2code/SU2/pull/2458*/
         const su2double density = flowNodes->GetDensity(iPoint);
-        
+
         for (unsigned short iVar = 0; iVar < nVar; iVar++) {
           nodes->AddClippedSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint) * LinSysSol(iPoint, iVar),
                                     lowerlimit[iVar], upperlimit[iVar], density, density);

SU2_CFD/src/output/CHeatOutput.cpp

@@ -77,7 +77,6 @@ void CHeatOutput::LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSol
   CSolver* heat_solver = solver[HEAT_SOL];
 
   output->SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux());
-  output->SetHistoryOutputValue("MAXIMUM_HEATFLUX", heat_solver->GetTotal_MaxHeatFlux());
   output->SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature());
   output->SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0)));
   output->SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0)));
@@ -109,7 +108,6 @@ void CHeatOutput::SetHistoryOutputFieldsImpl(CConfig *config, COutput* output) {
   output->AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "Block-Gauss-Seidel residual of the temperature", HistoryFieldType::RESIDUAL);
 
   output->AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
-  output->AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
   output->AddHistoryOutput("AVG_TEMPERATURE", "AvgTemp", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT);
   output->AddHistoryOutput("CFL_NUMBER", "CFL number", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current value of the CFL number");
 }

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -566,7 +566,7 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container,
 
   /*--- Set the pointer to the heat solver so we can access temperatures. ---*/
   if (config->GetWeakly_Coupled_Heat()) {
-    heat_solver = solver_container[HEAT_SOL];
+    heat_nodes = solver_container[HEAT_SOL]->GetNodes();
   }
 
   /*
@@ -652,6 +652,7 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics,
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_exponent = config->GetSIMP_Exponent();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   /*--- Start OpenMP parallel region. ---*/
 
@@ -694,8 +695,8 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics,
             element->SetRef_Coord(iNode, iDim, val_Coord);
             element->SetCurr_Coord(iNode, iDim, val_Sol);
           }
-          if (heat_solver) {
-            element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -749,6 +750,7 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_exponent = config->GetSIMP_Exponent();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   /*--- Start OpenMP parallel region. ---*/
 
@@ -804,8 +806,8 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri
               de_elem->SetRef_Coord(iNode, iDim, val_Coord);
             }
           }
-          if (heat_solver) {
-            fea_elem->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            fea_elem->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -890,6 +892,7 @@ void CFEASolver::Compute_MassMatrix(const CGeometry *geometry, CNumerics **numer
 
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   /*--- Never record this method as the mass matrix is passive (but the mass residual is not). ---*/
   const bool wasActive = AD::BeginPassive();
@@ -930,8 +933,8 @@ void CFEASolver::Compute_MassMatrix(const CGeometry *geometry, CNumerics **numer
             su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
-          if (heat_solver) {
-            element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -980,6 +983,7 @@ void CFEASolver::Compute_MassRes(const CGeometry *geometry, CNumerics **numerics
 
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   /*--- Clear vector before calculation. ---*/
   TimeRes.SetValZero();
@@ -1014,8 +1018,8 @@ void CFEASolver::Compute_MassRes(const CGeometry *geometry, CNumerics **numerics
           su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
           element->SetRef_Coord(iNode, iDim, val_Coord);
         }
-        if (heat_solver) {
-          element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+        if (heat_nodes) {
+          element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
         }
       }
 
@@ -1062,6 +1066,7 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_exponent = config->GetSIMP_Exponent();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   /*--- Start OpenMP parallel region. ---*/
 
@@ -1111,8 +1116,8 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric
             element->SetCurr_Coord(iNode, iDim, val_Sol);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
-          if (heat_solver) {
-            element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -1158,6 +1163,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_exponent = config->GetSIMP_Exponent();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
+  const su2double t_ref = config->GetTemperature_Ref();
 
   const auto stressParam = config->GetStressPenaltyParam();
   const su2double stress_scale = 1.0 / stressParam[0];
@@ -1230,8 +1236,8 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
             element->SetCurr_Coord(iNode, iDim, val_Sol);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
-          if (heat_solver) {
-            element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -1434,6 +1440,8 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
 
 void CFEASolver::Compute_BodyForces(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
+  const su2double t_ref = config->GetTemperature_Ref();
+
   /*--- Start OpenMP parallel region. ---*/
 
   SU2_OMP_PARALLEL
@@ -1473,8 +1481,8 @@ void CFEASolver::Compute_BodyForces(CGeometry *geometry, CNumerics **numerics, c
             su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
             element->SetRef_Coord(iNode, iDim, val_Coord);
           }
-          if (heat_solver) {
-            element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+          if (heat_nodes) {
+            element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
           }
         }
 
@@ -3017,6 +3025,7 @@ void CFEASolver::Stiffness_Penalty(CGeometry *geometry, CNumerics **numerics, CC
     PenaltyValue = 0.0;
     return;
   }
+  const su2double t_ref = config->GetTemperature_Ref();
 
   su2double weightedValue = 0.0;
   su2double weightedValue_reduce = 0.0;
@@ -3048,8 +3057,8 @@ void CFEASolver::Stiffness_Penalty(CGeometry *geometry, CNumerics **numerics, CC
         su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim);
         element->SetRef_Coord(iNode, iDim, val_Coord);
       }
-      if (heat_solver) {
-        element->SetTemperature(iNode, heat_solver->GetNodes()->GetSolution(indexNode[iNode], 0));
+      if (heat_nodes) {
+        element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref);
       }
     }
 

SU2_CFD/src/solvers/CHeatSolver.cpp

@@ -35,7 +35,7 @@ template class CScalarSolver<CHeatVariable>;
 
 CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh)
   : CScalarSolver<CHeatVariable>(geometry, config, false),
-    flow(config->GetFluidProblem()), heat_equation(config->GetHeatProblem()) {
+    flow(config->GetFluidProblem()) {
 
   /*--- Dimension of the problem --> temperature is the only conservative variable ---*/
 
@@ -101,14 +101,13 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iM
   /*--- Set the reference values for heat fluxes. If the heat solver runs stand-alone,
    *    thermal conductivity is read directly from config file ---*/
 
-  if (heat_equation) {
+  if (!flow) {
     su2double rho_cp = config->GetMaterialDensity(0)*config->GetSpecific_Heat_Cp();
     config->SetThermalDiffusivity(config->GetThermal_Conductivity_Constant() / rho_cp);
 
     /*--- Fluxes are computed via thermal diffusivity (not conductivity), so we have to divide by rho*cp ---*/
     config->SetHeat_Flux_Ref(rho_cp*Temperature_Ref);
-  }
-  else if (flow) {
+  } else {
     config->SetHeat_Flux_Ref(config->GetViscosity_Ref()*config->GetSpecific_Heat_Cp());
   }
 
@@ -551,8 +550,7 @@ void CHeatSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solv
       }
     }
     END_SU2_OMP_FOR
-  }
-  else if (heat_equation) {
+  } else {
     SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
     for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) {
 

TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg

@@ -0,0 +1,48 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SU2 configuration file                                                       %
+% Case description: 3D beam with thermal expansion                             %
+% File Version 8.1.0 "Harrier"                                                 %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+SOLVER= ELASTICITY
+MATH_PROBLEM= DIRECT
+GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS
+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
+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= CONJUGATE_GRADIENT
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-8
+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
+
+% Coupling with heat solver.
+WEAKLY_COUPLED_HEAT_EQUATION= YES
+FREESTREAM_TEMPERATURE= 300
+SPECIFIC_HEAT_CP= 460
+THERMAL_CONDUCTIVITY_CONSTANT= 45
+% NOTE: These markers a duplicates of "left" and "right" to allow specifying
+% boundary conditions for both solvers. This is work in progress.
+MARKER_ISOTHERMAL= ( left_heat, 400, right_heat, 300 )
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+CFL_NUMBER= 1e8
+
+MARKER_MONITORING= ( left_heat )
+SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS, TOTAL_HEATFLUX
+

TestCases/parallel_regression.py

@@ -1233,6 +1233,15 @@ def main():
     statbeam3d.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
     test_list.append(statbeam3d)
 
+    # Static beam, 3d with coupled temperature
+    thermal_beam_3d = TestCase('thermal_beam_3d')
+    thermal_beam_3d.cfg_dir = "fea_fsi/ThermalBeam_3d"
+    thermal_beam_3d.cfg_file = "configBeam_3d.cfg"
+    thermal_beam_3d.test_iter = 0
+    thermal_beam_3d.test_vals = [-6.123939, -5.890462, -5.999947, -8.091358, 249, -8.27818, 81, -8.40525, 1.3562e+05, 1.4465e+02]
+    thermal_beam_3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
+    test_list.append(thermal_beam_3d)
+
     # Rotating cylinder, 3d
     rotating_cylinder_fea           = TestCase('rotating_cylinder_fea')
     rotating_cylinder_fea.cfg_dir   = "fea_fsi/rotating_cylinder"