diff --git a/SU2_CFD/include/output/CElasticityOutput.hpp b/SU2_CFD/include/output/CElasticityOutput.hpp index ca65cab7a9c2..528858b68d75 100644 --- a/SU2_CFD/include/output/CElasticityOutput.hpp +++ b/SU2_CFD/include/output/CElasticityOutput.hpp @@ -40,6 +40,7 @@ class CElasticityOutput final: public COutput { unsigned short nVar_FEM; //!< Number of FEM variables bool linear_analysis, //!< Boolean indicating a linear analysis nonlinear_analysis, //!< Boolean indicating a nonlinear analysis + coupled_heat, //!< Boolean indicating a thermoelastic analysis dynamic; //!< Boolean indicating a dynamic analysis public: @@ -82,6 +83,11 @@ class CElasticityOutput final: public COutput { * \param[in] config - Definition of the particular problem. * \return if the residuals should be initialized. */ - bool SetInitResiduals(const CConfig *config) override ; + bool SetInitResiduals(const CConfig *config) override; + /*! + * \brief LoadSurfaceData + */ + void LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, + unsigned short iMarker, unsigned long iVertex) override; }; diff --git a/SU2_CFD/include/output/CHeatOutput.hpp b/SU2_CFD/include/output/CHeatOutput.hpp index 04babb1b7c33..ca7bc61df25e 100644 --- a/SU2_CFD/include/output/CHeatOutput.hpp +++ b/SU2_CFD/include/output/CHeatOutput.hpp @@ -50,12 +50,22 @@ class CHeatOutput final: public CFVMOutput { */ void SetHistoryOutputFields(CConfig *config) override; + /*! + * \brief Set the available history output fields in another output instance. + */ + static void SetHistoryOutputFieldsImpl(CConfig *config, COutput* output); + /*! * \brief Load the history output field values * \param[in] config - Definition of the particular problem. */ void LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) override; + /*! + * \brief Set the history output field values in another output instance. + */ + static void LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output); + /*! * \brief Set the available volume output fields * \param[in] config - Definition of the particular problem. diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 14b0bbec9b6a..fd9cd7b4f929 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -55,6 +55,7 @@ class CSolver; class CFileWriter; class CParallelDataSorter; class CConfig; +class CHeatOutput; using namespace std; @@ -65,6 +66,7 @@ using namespace std; */ class COutput { protected: + friend class CHeatOutput; /*----------------------------- General ----------------------------*/ diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 9442fa190d7a..ada30a36fa37 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -100,6 +100,11 @@ class CFEASolver : public CFEASolverBase { bool initial_calc = true; /*!< \brief Becomes false after first call to Preprocessing. */ bool body_forces = false; /*!< \brief Whether any body force is active. */ + /*! + * \brief Pointer to the heat solver nodes to access temperature for coupled simulations. + */ + const CVariable* heat_nodes = nullptr; + /*! * \brief The highest level in the variable hierarchy this solver can safely use, * CVariable is the common denominator between the FEA and Mesh deformation variables. diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp index 61a750ebc422..1d1de687a13b 100644 --- a/SU2_CFD/include/solvers/CHeatSolver.hpp +++ b/SU2_CFD/include/solvers/CHeatSolver.hpp @@ -42,7 +42,6 @@ class CHeatSolver final : public CScalarSolver { 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; diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index bf3d585ab7b9..7eae805673b7 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -503,10 +503,10 @@ void CScalarSolver::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); diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 047650701075..9bd268932fc8 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1527,7 +1527,9 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver case MAIN_SOLVER::FEM_ELASTICITY: case MAIN_SOLVER::DISC_ADJ_FEM: - fem = true; break; + fem = true; + heat = config->GetWeakly_Coupled_Heat(); + break; case MAIN_SOLVER::ADJ_EULER: adj_euler = euler = compressible = true; break; diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index 2aa19d27de9d..d0cf9159d617 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -48,6 +48,13 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL]; CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; + /*--- Add heat solver integration step. ---*/ + if (config[val_iZone]->GetWeakly_Coupled_Heat()) { + config[val_iZone]->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS); + integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, + RUNTIME_HEAT_SYS, val_iZone, val_iInst); + } + /*--- FEA equations ---*/ config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 1af62426cce0..0327cf36a1d3 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -24,18 +24,18 @@ * You should have received a copy of the GNU Lesser General Public * License along with SU2. If not, see . */ - - #include "../../include/output/CElasticityOutput.hpp" +#include "../../include/output/CHeatOutput.hpp" #include "../../../Common/include/geometry/CGeometry.hpp" #include "../../include/solvers/CSolver.hpp" CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COutput(config, nDim, false) { - linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL); - nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE); - dynamic = (config->GetTime_Domain()); + linear_analysis = config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL; + nonlinear_analysis = config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE; + coupled_heat = config->GetWeakly_Coupled_Heat(); + dynamic = config->GetTime_Domain(); /*--- Initialize number of variables ---*/ if (linear_analysis) nVar_FEM = nDim; @@ -48,20 +48,21 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu } /*--- Default fields for screen output ---*/ - if (nRequestedScreenFields == 0){ + if (nRequestedScreenFields == 0) { if (dynamic) requestedScreenFields.emplace_back("TIME_ITER"); if (multiZone) requestedScreenFields.emplace_back("OUTER_ITER"); requestedScreenFields.emplace_back("INNER_ITER"); - if(linear_analysis){ + if (linear_analysis) { requestedScreenFields.emplace_back("RMS_DISP_X"); requestedScreenFields.emplace_back("RMS_DISP_Y"); requestedScreenFields.emplace_back("RMS_DISP_Z"); } - if(nonlinear_analysis){ + if (nonlinear_analysis) { requestedScreenFields.emplace_back("RMS_UTOL"); requestedScreenFields.emplace_back("RMS_RTOL"); requestedScreenFields.emplace_back("RMS_ETOL"); } + if (coupled_heat) requestedScreenFields.emplace_back("RMS_TEMPERATURE"); requestedScreenFields.emplace_back("VMS"); nRequestedScreenFields = requestedScreenFields.size(); } @@ -71,11 +72,8 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("STRESS"); - if (dynamic) { - requestedVolumeFields.emplace_back("VELOCITY"); - requestedVolumeFields.emplace_back("ACCELERATION"); - } if (config->GetTopology_Optimization()) requestedVolumeFields.emplace_back("TOPOLOGY"); + if (coupled_heat) requestedVolumeFields.emplace_back("PRIMITIVE"); nRequestedVolumeFields = requestedVolumeFields.size(); } @@ -106,6 +104,7 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { CSolver* fea_solver = solver[FEA_SOL]; + CSolver* heat_solver = solver[HEAT_SOL]; /*--- Residuals: ---*/ /*--- Linear analysis: RMS of the displacements in the nDim coordinates ---*/ @@ -147,6 +146,13 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("TOPOL_DISCRETENESS", fea_solver->GetTotal_OFDiscreteness()); } + /*--- Add heat solver data if available. ---*/ + if (coupled_heat) { + CHeatOutput::LoadHistoryDataImpl(config, geometry, solver, this); + SetHistoryOutputValue("LINSOL_ITER_HEAT", heat_solver->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_RESIDUAL_HEAT", log10(heat_solver->GetResLinSolver())); + } + ComputeSimpleCustomOutputs(config); /*--- Keep this as last, since it uses the history values that were set. ---*/ @@ -189,6 +195,11 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config) { } AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); + if (coupled_heat) { + CHeatOutput::SetHistoryOutputFieldsImpl(config, this); + AddHistoryOutput("LINSOL_ITER_HEAT", "LinSolIterHeat", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); + AddHistoryOutput("LINSOL_RESIDUAL_HEAT", "LinSolResHeat", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); + } } void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){ @@ -214,6 +225,10 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo SetVolumeOutputValue("ACCELERATION-Y", iPoint, Node_Struc->GetSolution_Accel(iPoint, 1)); if (nDim == 3) SetVolumeOutputValue("ACCELERATION-Z", iPoint, Node_Struc->GetSolution_Accel(iPoint, 2)); } + if (coupled_heat) { + CVariable* Node_Heat = solver[HEAT_SOL]->GetNodes(); + SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); + } SetVolumeOutputValue("STRESS-XX", iPoint, Node_Struc->GetStress_FEM(iPoint)[0]); SetVolumeOutputValue("STRESS-YY", iPoint, Node_Struc->GetStress_FEM(iPoint)[1]); @@ -228,6 +243,14 @@ void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSo if (config->GetTopology_Optimization()) { SetVolumeOutputValue("TOPOL_DENSITY", iPoint, Node_Struc->GetAuxVar(iPoint)); } + + CSolver* heat_solver = solver[HEAT_SOL]; + if (heat_solver) { + const auto Node_Heat = heat_solver->GetNodes(); + SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Heat->GetSolution(iPoint, 0)); + SetVolumeOutputValue("RES_TEMPERATURE", iPoint, heat_solver->LinSysRes(iPoint, 0)); + } + } void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ @@ -251,6 +274,10 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ if (nDim == 3) AddVolumeOutput("ACCELERATION-Z", "Acceleration_z", "SOLUTION", "z-component of the acceleration vector"); } + if (coupled_heat) { + AddVolumeOutput("TEMPERATURE", "Temperature", "SOLUTION", "Temperature"); + } + AddVolumeOutput("STRESS-XX", "Sxx", "STRESS", "x-component of the normal stress vector"); AddVolumeOutput("STRESS-YY", "Syy", "STRESS", "y-component of the normal stress vector"); AddVolumeOutput("STRESS-XY", "Sxy", "STRESS", "xy shear stress component"); @@ -266,6 +293,12 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){ if (config->GetTopology_Optimization()) { AddVolumeOutput("TOPOL_DENSITY", "Topology_Density", "TOPOLOGY", "filtered topology density"); } + + if (coupled_heat) { + AddVolumeOutput("HEAT_FLUX", "Heat_Flux", "PRIMITIVE", "Heatflux"); + AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); + } + } bool CElasticityOutput::SetInitResiduals(const CConfig *config){ @@ -273,3 +306,12 @@ bool CElasticityOutput::SetInitResiduals(const CConfig *config){ return (config->GetTime_Domain() == NO && (curInnerIter == 0)); } + +void CElasticityOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, + unsigned short iMarker, unsigned long iVertex) { + if (!coupled_heat || !config->GetViscous_Wall(iMarker)) return; + + /* Heat flux value at each surface grid node. */ + SetVolumeOutputValue("HEAT_FLUX", iPoint, solver[HEAT_SOL]->GetHeatFlux(iMarker, iVertex)); + +} diff --git a/SU2_CFD/src/output/CHeatOutput.cpp b/SU2_CFD/src/output/CHeatOutput.cpp index f7acc239e9b1..50c5c1714ece 100644 --- a/SU2_CFD/src/output/CHeatOutput.cpp +++ b/SU2_CFD/src/output/CHeatOutput.cpp @@ -72,21 +72,28 @@ CHeatOutput::CHeatOutput(CConfig *config, unsigned short nDim) : CFVMOutput(conf } +void CHeatOutput::LoadHistoryDataImpl(CConfig *config, CGeometry *geometry, CSolver **solver, COutput* output) { + + CSolver* heat_solver = solver[HEAT_SOL]; + + output->SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux()); + 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))); + if (config->GetMultizone_Problem()) { + output->SetHistoryOutputValue("BGS_TEMPERATURE", log10(heat_solver->GetRes_BGS(0))); + } + output->SetHistoryOutputValue("CFL_NUMBER", config->GetCFL(MESH_0)); +} + void CHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { CSolver* heat_solver = solver[HEAT_SOL]; - SetHistoryOutputValue("TOTAL_HEATFLUX", heat_solver->GetTotal_HeatFlux()); - SetHistoryOutputValue("MAXIMUM_HEATFLUX", heat_solver->GetTotal_MaxHeatFlux()); - SetHistoryOutputValue("AVG_TEMPERATURE", heat_solver->GetTotal_AvgTemperature()); - SetHistoryOutputValue("RMS_TEMPERATURE", log10(heat_solver->GetRes_RMS(0))); - SetHistoryOutputValue("MAX_TEMPERATURE", log10(heat_solver->GetRes_Max(0))); - if (multiZone) - SetHistoryOutputValue("BGS_TEMPERATURE", log10(heat_solver->GetRes_BGS(0))); + LoadHistoryDataImpl(config, geometry, solver, this); SetHistoryOutputValue("LINSOL_ITER", heat_solver->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(heat_solver->GetResLinSolver())); - SetHistoryOutputValue("CFL_NUMBER", config->GetCFL(MESH_0)); ComputeSimpleCustomOutputs(config); @@ -94,20 +101,23 @@ void CHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver SetCustomAndComboObjectives(HEAT_SOL, config, solver); } +void CHeatOutput::SetHistoryOutputFieldsImpl(CConfig *config, COutput* output) { -void CHeatOutput::SetHistoryOutputFields(CConfig *config){ + output->AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL); + output->AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL); + output->AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "Block-Gauss-Seidel residual of the temperature", HistoryFieldType::RESIDUAL); - AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); - AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); + output->AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total 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"); +} - AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root mean square residual of the temperature", HistoryFieldType::RESIDUAL); - AddHistoryOutput("MAX_TEMPERATURE", "max[T]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the temperature", HistoryFieldType::RESIDUAL); - AddHistoryOutput("BGS_TEMPERATURE", "bgs[T]", ScreenOutputFormat::FIXED, "BGS_RES", "Block-Gauss-Seidel residual of the temperature", HistoryFieldType::RESIDUAL); +void CHeatOutput::SetHistoryOutputFields(CConfig *config) { - AddHistoryOutput("TOTAL_HEATFLUX", "HF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Total heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("MAXIMUM_HEATFLUX", "MaxHF", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Maximum heatflux on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("AVG_TEMPERATURE", "AvgTemp", ScreenOutputFormat::SCIENTIFIC, "HEAT", "Average temperature on all surfaces defined in MARKER_MONITORING", HistoryFieldType::COEFFICIENT); - AddHistoryOutput("CFL_NUMBER", "CFL number", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current value of the CFL number"); + SetHistoryOutputFieldsImpl(config, this); + + AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); + AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 9d19aa7c5a05..ec11a36889d1 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -30,6 +30,7 @@ #include "../../include/numerics/elasticity/CFEAElasticity.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/solvers/CHeatSolver.hpp" #include using namespace GeometryToolbox; @@ -563,6 +564,11 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, const bool disc_adj_fem = (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_FEM); const bool topology_mode = config->GetTopology_Optimization(); + /*--- Set the pointer to the heat solver so we can access temperatures. ---*/ + if (config->GetWeakly_Coupled_Heat()) { + heat_nodes = solver_container[HEAT_SOL]->GetNodes(); + } + /* * For topology optimization we apply a filter on the design density field to avoid * numerical issues (checkerboards), ensure mesh independence, and impose a length scale. @@ -646,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. ---*/ @@ -688,6 +695,9 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, element->SetRef_Coord(iNode, iDim, val_Coord); element->SetCurr_Coord(iNode, iDim, val_Sol); } + if (heat_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the stiffness. ---*/ @@ -740,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. ---*/ @@ -795,6 +806,9 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri de_elem->SetRef_Coord(iNode, iDim, val_Coord); } } + if (heat_nodes) { + fea_elem->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the stiffness. ---*/ @@ -878,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(); @@ -918,6 +933,9 @@ 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_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the mass, @@ -965,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(); @@ -999,6 +1018,9 @@ 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_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the mass, @@ -1044,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. ---*/ @@ -1093,6 +1116,9 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric element->SetCurr_Coord(iNode, iDim, val_Sol); element->SetRef_Coord(iNode, iDim, val_Coord); } + if (heat_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the stiffness ---*/ @@ -1137,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]; @@ -1209,6 +1236,9 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, element->SetCurr_Coord(iNode, iDim, val_Sol); element->SetRef_Coord(iNode, iDim, val_Coord); } + if (heat_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- In topology mode determine the penalty to apply to the stiffness ---*/ @@ -1410,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 @@ -1449,6 +1481,9 @@ 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_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } /*--- Penalize the dead load, do it by default to avoid unecessary "ifs", since it @@ -2990,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; @@ -3021,6 +3057,9 @@ 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_nodes) { + element->SetTemperature(iNode, heat_nodes->GetSolution(indexNode[iNode], 0) * t_ref); + } } // Avoid double-counting elements: diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index dbbb7191f172..58e53f03421e 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -35,7 +35,7 @@ template class CScalarSolver; CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CScalarSolver(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++) { diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index 743775ad139b..7c055241624b 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -173,6 +173,9 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon break; case MAIN_SOLVER::FEM_ELASTICITY: solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); + if (config->GetWeakly_Coupled_Heat()) { + solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); + } break; case MAIN_SOLVER::DISC_ADJ_FEM: solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel); diff --git a/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg new file mode 100644 index 000000000000..7c5b49dd6d8f --- /dev/null +++ b/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 + diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref index f633bb537397..ac3cc62550ad 100644 --- a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref +++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/of_grad_findiff.csv.ref @@ -1,2 +1,2 @@ -"VARIABLE" , "AVG_DENSITY[0]", "AVG_ENTHALPY[0]", "AVG_NORMALVEL[0]", "DRAG[0]" , "EFFICIENCY[0]" , "FORCE_X[0]" , "FORCE_Y[0]" , "FORCE_Z[0]" , "LIFT[0]" , "MOMENT_X[0]" , "MOMENT_Y[0]" , "MOMENT_Z[0]" , "SIDEFORCE[0]" , "SURFACE_MACH[0]", "SURFACE_MASSFLOW[0]", "SURFACE_MOM_DISTORTION[0]", "SURFACE_PRESSURE_DROP[0]", "SURFACE_SECONDARY[0]", "SURFACE_SECOND_OVER_UNIFORM[0]", "SURFACE_STATIC_PRESSURE[0]", "SURFACE_STATIC_TEMPERATURE[0]", "SURFACE_TOTAL_PRESSURE[0]", "SURFACE_TOTAL_TEMPERATURE[0]", "SURFACE_UNIFORMITY[0]", "AVG_TEMPERATURE[1]", "MAXIMUM_HEATFLUX[1]", "TOTAL_HEATFLUX[1]", "FINDIFF_STEP" -0 , 0.0 , -1000000.0242143869, -9.992000000081167e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -1.4199999995301305, 0.0 , -36.519999999096164 , 0.0 , 39.07999999996914 , 73.22000000198337 , 6360.000003269306 , -279.99999474559445 , 109.99999915384251 , -290.0000026784255 , -34.39999998189336 , -110.00000199601345, 0.0 , 260.00000161729986, 1e-08 +"VARIABLE" , "AVG_DENSITY[0]", "AVG_ENTHALPY[0]", "AVG_NORMALVEL[0]", "DRAG[0]" , "EFFICIENCY[0]" , "FORCE_X[0]" , "FORCE_Y[0]" , "FORCE_Z[0]" , "LIFT[0]" , "MOMENT_X[0]" , "MOMENT_Y[0]" , "MOMENT_Z[0]" , "SIDEFORCE[0]" , "SURFACE_MACH[0]", "SURFACE_MASSFLOW[0]", "SURFACE_MOM_DISTORTION[0]", "SURFACE_PRESSURE_DROP[0]", "SURFACE_SECONDARY[0]", "SURFACE_SECOND_OVER_UNIFORM[0]", "SURFACE_STATIC_PRESSURE[0]", "SURFACE_STATIC_TEMPERATURE[0]", "SURFACE_TOTAL_PRESSURE[0]", "SURFACE_TOTAL_TEMPERATURE[0]", "SURFACE_UNIFORMITY[0]", "AVG_TEMPERATURE[1]", "TOTAL_HEATFLUX[1]", "FINDIFF_STEP" +0 , 0.0 , -1000000.0242143869, -9.992000000081167e-08, 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -1.4199999995301305, 0.0 , -36.519999999096164 , 0.0 , 39.07999999996914 , 73.22000000198337 , 6360.000003269306 , -279.99999474559445 , 109.99999915384251 , -290.0000026784255 , -34.39999998189336 , -110.00000199601345, 260.00000161729986, 1e-08 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 1ef6d9ecf7bb..46eb842f808b 100644 --- a/TestCases/parallel_regression.py +++ b/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.140220, -5.842734, -5.972391, -8.091358, 262, -8.246755, 81, -8.298569, 135620, 144.65] + 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" @@ -1314,8 +1323,8 @@ def main(): solid_periodic_pins.cfg_dir = "solid_heat_conduction/periodic_pins" solid_periodic_pins.cfg_file = "configSolid.cfg" solid_periodic_pins.test_iter = 750 - solid_periodic_pins.test_vals = [-15.878977, -14.569206, 300.900000, 425.320000, 0.000000, 5.000000, -1.672737] - solid_periodic_pins.test_vals_aarch64 = [-15.879016, -14.569206, 300.900000, 425.320000, 0.000000, 5.000000, -1.672666] + solid_periodic_pins.test_vals = [-15.878977, -14.569206, 300.900000, 425.320000, 5.000000, -1.672737] + solid_periodic_pins.test_vals_aarch64 = [-15.879016, -14.569206, 300.900000, 425.320000, 5.000000, -1.672666] test_list.append(solid_periodic_pins) # ###############################