[ThermalDEM] Heat generation features

applications/ThermalDEMApplication/CMakeLists.txt

@@ -69,8 +69,11 @@ set( KRATOS_THERMAL_DEM_APPLICATION_CORE_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/radiation/radiation_continuum_zhou.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_model.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_lu.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_morris.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_zhou.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_morris_area.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_morris_area_time.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_rangel_area.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/real_contact/real_contact_rangel_area_time.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_elements/thermal_spheric_particle.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_elements/thermal_spheric_continuum_particle.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_elements/sintering_spheric_continuum_particle.cpp

applications/ThermalDEMApplication/README.md

@@ -2,14 +2,14 @@
 
 This application is an extension of the [DEM Application](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/DEMApplication) to include **thermal effects** such as:
 
-- Heat transfer between particle-particle, particle-rigid wall, and particle-surrounding fluid.
+- Heat transfer between particle-particle, particle-wall, and particle-surrounding fluid.
 - Heat transfer mechanisms by conduction, convection, and radiation.
 - Heat generation by energy dissipation and internal sources.
 - Temperature dependent material properties.
 
-Theoretical information on thermal DEM analysis can be found [here](./ThermalDEMTheory.pdf).
+Theoretical information on thermal DEM analysis can be found [here](https://doi.org/10.5281/zenodo.13846126).
 
-A [Matlab version](https://gitlab.com/rafaelrangel/demlab) of this application is also available.
+A [Matlab version](https://github.com/rlrangel/DEMLab) of this application is also available.
 
 ## Table of Contents
 - [Authorship](#authorship)
@@ -124,17 +124,20 @@ Add **thermal settings** with desired options:
 
 Add **post options** with desired options:
 
-	"PostTemperature"                : true or false,
-	"PostHeatFlux"                   : true or false,
-	"PostGraphParticleTempMin"       : true or false,
-	"PostGraphParticleTempMax"       : true or false,
-	"PostGraphParticleTempAvg"       : true or false,
-	"PostGraphParticleTempDev"       : true or false,
-	"PostGraphModelTempAvg"          : true or false,
-	"PostGraphHeatFluxContributions" : true or false,
-	"PostGraphHeatGenContributions"  : true or false,
-	"PostGraphEnergyContributions"   : true or false,
-	"PostMapHeatGeneration"          : true or false
+	"PostTemperature"                      : true or false,
+	"PostHeatFlux"                         : true or false,
+	"PostGraphParticleTempMin"             : true or false,
+	"PostGraphParticleTempMax"             : true or false,
+	"PostGraphParticleTempAvg"             : true or false,
+	"PostGraphParticleTempAvgVol"          : true or false,
+	"PostGraphParticleTempDev"             : true or false,
+	"PostGraphMechanicalEnergy"            : true or false,
+	"PostGraphDissipatedEnergy"            : true or false,
+	"PostGraphThermalEnergy"               : true or false,
+	"PostGraphHeatFluxContributions"       : true or false,
+	"PostGraphHeatGenerationValues"        : true or false,
+	"PostGraphHeatGenerationContributions" : true or false,
+	"PostHeatMapGeneration"                : true or false
 
 ### Materials (json file)
 
@@ -209,35 +212,35 @@ Add **SubModelPartData** to sub model parts with desired options:
 - *"thermal_integration_scheme"*:\
   Selected scheme for time integration of thermal problem.\
   Default: "forward_euler"
-  
+
 - *"numerical_integration_method"*:\
   Selected metohd for solving integration expressions numerically.\
   Default: "adaptive_simpson"
-  
+
 - *"thermal_solve_frequency"*:\
   Number of steps in which thermal problem is solved.\
   Default: 1
-  
+
 - *"voronoi_tesselation_frequency"*:\
   Number of steps in which Voronoi diagram is built, in case it is required.\
   Default: 1000
-  
+
 - *"porosity_update_frequency"*:\
   Number of steps in which porosity is computed, in case it is required.\
   Default: 1000
 
 - *"automatic_solve_frequency"*:\
   Boolean for automatically setting the thermal solve frequency based on the maximum allowed time step (it overrides the value set for thermal_solve_frequency).\
   Default: false
-  
+
 - *"compute_motion"*:\
   Boolean for solving mechanical problem.\
   Default: true
 
 - *"compute_direct_conduction"*:\
   Boolean for computing heat transfer between elements by direct conduction.\
   Default: true
-  
+
 - *"compute_indirect_conduction"*:\
   Boolean for computing heat transfer between elements by indirect conduction.\
   Default: false
@@ -277,7 +280,7 @@ Add **SubModelPartData** to sub model parts with desired options:
 - *"heat_generation_model"*:\
   List of selected models for simulating heat generation by energy dissipation.\
   Default: ["sliding_friction"]
-  
+
 - *"adjusted_contact_model"*:\
   Selected model for adjusting contact geometry.\
   Default: "zhou"
@@ -301,7 +304,7 @@ Add **SubModelPartData** to sub model parts with desired options:
 - *"conduction_radius"*:\
   Radius of cylindrical conductive region (ratio of particles radii) required for conduction model "batchelor_obrien_complete" and "batchelor_obrien_modified".\
   Default: 1.0
-  
+
 - *"fluid_layer_thickness"*:\
   Thickness of particle fluid layer (ratio of particles radii) required for conduction model "surrounding_layer".\
   Default: 0.4
@@ -337,7 +340,7 @@ Add **SubModelPartData** to sub model parts with desired options:
 - *"heat_map_subdivisions"*:\
   Number of subdivisions in X,Y,Z directions to defined the grid of the heat map.\
   Default: [10,10,10]
-  
+
 - *"global_fluid_properties"*:\
   Prescribed values for the properties of the interstitial fluid, assumed as constant throughout all the analysis (fluid behavior does not change as it is not simulated).
 
@@ -377,47 +380,59 @@ Add **SubModelPartData** to sub model parts with desired options:
 - *"PostTemperature"*:\
   Boolean for showing elements temperature in post processing.\
   Default: false
-  
+
 - *"PostHeatFlux"*:\
   Boolean for showing elements heat flux in post processing.\
   Default: false
-  
+
 - *"PostGraphParticleTempMin"*:\
   Boolean for writing a graph with the minimum particle temperature.\
   Default: false
-  
+
 - *"PostGraphParticleTempMax"*:\
   Boolean for writing a graph with the maximum particle temperature.\
   Default: false
-  
+
 - *"PostGraphParticleTempAvg"*:\
   Boolean for writing a graph with the average temperature of all particles.\
   Default: false
-
+
+- *"PostGraphParticleTempAvgVol"*:\
+  Boolean for writing a graph with the volume-weighted average temperature of all particles.\
+  Default: false
+
 - *"PostGraphParticleTempDev"*:\
   Boolean for writing a graph with the standard deviation of the temperature of all particles.\
   Default: false
-  
-- *"PostGraphModelTempAvg"*:\
-  Boolean for writing a graph with the weighted average of the temperature of all particles (taking the particles volume into account).\
+
+- *"PostGraphMechanicalEnergy"*:\
+  Boolean for writing a graph with the mechanical energy components of all partilces.\
   Default: false
-
+
+- *"PostGraphDissipatedEnergy"*:\
+  Boolean for writing a graph with the accumulated energy dissipation of all partilces.\
+  Default: false
+
+- *"PostGraphThermalEnergy"*:\
+  Boolean for writing a graph with the accumulated thermal energy (from heat generation by energy dissipation) of all partilces.\
+  Default: false
+
 - *"PostGraphHeatFluxContributions"*:\
   Boolean for writing a graph with the contribution of each heat transfer mechanism to the total heat transfer.\
   Default: false
 
-- *"PostGraphHeatGenContributions"*:\
-  Boolean for writing a graph with the contribution of each heat generation mechanism to the total heat generation.\
+- *"PostGraphHeatGenerationValues"*:\
+  Boolean for writing a graph with the total values of  each heat generation mechanism.\
   Default: false
-  
-- *"PostGraphEnergyContributions"*:\
-  Boolean for writing a graph with the energy composition (conservative and accumulated dissipative components) of all partilces.\
+
+- *"PostGraphHeatGenerationContributions"*:\
+  Boolean for writing a graph with the contribution of each heat generation mechanism to the total heat generation.\
   Default: false
 
-- *"PostMapHeatGeneration"*:\
-  Boolean for assemblying and writing the density map of heat generation.\
+- *"PostHeatMapGeneration"*:\
+  Boolean for assemblying and writing the heat map of heat generation.\
   Default: false
-  
+
 **Material properties**
 - *"materials.Variables.THERMAL_CONDUCTIVITY"*:\
   Thermal conductivity of material (always required).

applications/ThermalDEMApplication/custom_constitutive/generation/generation_dissipation.cpp

@@ -38,6 +38,7 @@ namespace Kratos {
     const double coeff      = conversion * partition;
 
     // Initialize contribution from different sources of heat generation
+    double thermal_energy;
     double heat_gen;
     double heat_gen_damping_pp = 0.0;
     double heat_gen_damping_pw = 0.0;
@@ -48,49 +49,58 @@ namespace Kratos {
 
     // Damping thermal power
     if (r_process_info[GENERATION_DAMPING_OPTION]) {
-      heat_gen = coeff * contact_params.viscodamping_energy / time;
+      thermal_energy = coeff * contact_params.viscodamping_energy;
+      heat_gen = thermal_energy / time;
 
       if (particle->mNeighborType & PARTICLE_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_damp_particle += thermal_energy;
+        particle->mGenerationHeatFlux_damp_particle      += heat_gen;
         heat_gen_damping_pp = heat_gen;
-        particle->mGenerationHeatFlux_damp_particle += heat_gen;
       }
       else if (particle->mNeighborType & WALL_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_damp_wall += thermal_energy;
+        particle->mGenerationHeatFlux_damp_wall      += heat_gen;
         heat_gen_damping_pw = heat_gen;
-        particle->mGenerationHeatFlux_damp_wall += heat_gen;
       }  
     }
 
     // Sliding friction thermal power
     if (r_process_info[GENERATION_SLIDING_OPTION]) {
-      heat_gen = coeff * contact_params.frictional_energy / time;
+      thermal_energy = coeff * contact_params.frictional_energy;
+      heat_gen = thermal_energy / time;
 
       if (particle->mNeighborType & PARTICLE_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_slid_particle += thermal_energy;
+        particle->mGenerationHeatFlux_slid_particle      += heat_gen;
         heat_gen_sliding_pp = heat_gen;
-        particle->mGenerationHeatFlux_slid_particle += heat_gen;
       }
       else if (particle->mNeighborType & WALL_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_slid_wall += thermal_energy;
+        particle->mGenerationHeatFlux_slid_wall      += heat_gen;
         heat_gen_sliding_pw = heat_gen;
-        particle->mGenerationHeatFlux_slid_wall += heat_gen;
       }
     }
 
     // Rolling friction thermal power
     if (r_process_info[GENERATION_ROLLING_OPTION] && particle->Is(DEMFlags::HAS_ROTATION) && particle->Is(DEMFlags::HAS_ROLLING_FRICTION)) {
-      heat_gen = coeff * contact_params.rollresist_energy / time;
+      thermal_energy = coeff * contact_params.rollresist_energy;
+      heat_gen = thermal_energy / time;
 
       if (particle->mNeighborType & PARTICLE_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_roll_particle += thermal_energy;
+        particle->mGenerationHeatFlux_roll_particle      += heat_gen;
         heat_gen_rolling_pp = heat_gen;
-        particle->mGenerationHeatFlux_roll_particle += heat_gen;
       }
       else if (particle->mNeighborType & WALL_NEIGHBOR) {
+        particle->mGenerationThermalEnergy_roll_wall += thermal_energy;
+        particle->mGenerationHeatFlux_roll_wall      += heat_gen;
         heat_gen_rolling_pw = heat_gen;
-        particle->mGenerationHeatFlux_roll_wall += heat_gen;
       }
     }
 
     // Fill heat map
     if (r_process_info[HEAT_MAP_GENERATION_OPTION])
-      FillDensityMap(r_process_info, particle, time, heat_gen_damping_pp, heat_gen_damping_pw, heat_gen_sliding_pp, heat_gen_sliding_pw, heat_gen_rolling_pp, heat_gen_rolling_pw);
+      FillHeatMap(r_process_info, particle, time, heat_gen_damping_pp, heat_gen_damping_pw, heat_gen_sliding_pp, heat_gen_sliding_pw, heat_gen_rolling_pp, heat_gen_rolling_pw);
 
     return heat_gen_damping_pp + heat_gen_damping_pw + heat_gen_sliding_pp + heat_gen_sliding_pw + heat_gen_rolling_pp + heat_gen_rolling_pw;
 
@@ -109,7 +119,7 @@ namespace Kratos {
   }
 
   //------------------------------------------------------------------------------------------------------------
-  void GenerationDissipation::FillDensityMap(const ProcessInfo& r_process_info,
+  void GenerationDissipation::FillHeatMap(const ProcessInfo& r_process_info,
                                              ThermalSphericParticle* particle,
                                              const double time,
                                              const double heat_gen_damping_pp,

applications/ThermalDEMApplication/custom_constitutive/generation/generation_dissipation.h

@@ -32,15 +32,15 @@ namespace Kratos
       // Public methods
       double ComputeHeatGeneration (const ProcessInfo& r_process_info, ThermalSphericParticle* particle) override;
       double ComputePartitionCoeff (ThermalSphericParticle* particle);
-      void   FillDensityMap        (const ProcessInfo& r_process_info,
-                                    ThermalSphericParticle* particle,
-                                    const double time,
-                                    const double heat_gen_damping_pp,
-                                    const double heat_gen_damping_pw,
-                                    const double heat_gen_sliding_pp,
-                                    const double heat_gen_sliding_pw,
-                                    const double heat_gen_rolling_pp,
-                                    const double heat_gen_rolling_pw);
+      void   FillHeatMap (const ProcessInfo& r_process_info,
+                          ThermalSphericParticle* particle,
+                          const double time,
+                          const double heat_gen_damping_pp,
+                          const double heat_gen_damping_pw,
+                          const double heat_gen_sliding_pp,
+                          const double heat_gen_sliding_pw,
+                          const double heat_gen_rolling_pp,
+                          const double heat_gen_rolling_pw);
 
       // Clone
       HeatGenerationMechanism* CloneRaw() const override {

applications/ThermalDEMApplication/custom_constitutive/real_contact/real_contact_morris.cpp → applications/ThermalDEMApplication/custom_constitutive/real_contact/real_contact_morris_area.cpp

@@ -11,15 +11,15 @@
 // External includes
 
 // Project includes
-#include "real_contact_morris.h"
+#include "real_contact_morris_area.h"
 
 namespace Kratos {
   //-----------------------------------------------------------------------------------------------------------------------
-  RealContactMorris::RealContactMorris() {}
-  RealContactMorris::~RealContactMorris() {}
+  RealContactMorrisArea::RealContactMorrisArea() {}
+  RealContactMorrisArea::~RealContactMorrisArea() {}
 
   //------------------------------------------------------------------------------------------------------------
-  void RealContactMorris::AdjustContact(const ProcessInfo& r_process_info, ThermalSphericParticle* particle) {
+  void RealContactMorrisArea::AdjustContact(const ProcessInfo& r_process_info, ThermalSphericParticle* particle) {
     KRATOS_TRY
 
     // Parameters
@@ -35,7 +35,7 @@ namespace Kratos {
     const double correction_area = pow(hertz_force * eff_radius / eff_young_real, 1.0 / 3.0);
 
     // Time correction
-    double correction_time = 1.0; // TODO: Compute time correction from collision time
+    double correction_time = 1.0;
 
     // Adjusted value of contact radius
     particle->mContactRadiusAdjusted = correction_area * correction_time;