Add heat transfer

README.md

@@ -115,6 +115,7 @@ CabanaPD currently includes the following:
    - Elastic (no failure)
    - Brittle fracture
  - Thermomechanics (bond-based only)
+   - Optional heat transfer (elastic only)
  - Time integration
    - Velocity Verlet
  - Pre-crack creation

examples/thermomechanics/CMakeLists.txt

@@ -4,4 +4,7 @@ target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD)
 add_executable(ThermalCrack thermal_crack.cpp)
 target_link_libraries(ThermalCrack LINK_PUBLIC CabanaPD)
 
-install(TARGETS ThermalDeformation ThermalCrack DESTINATION ${CMAKE_INSTALL_BINDIR})
+add_executable(ThermalDeformationHeatTransfer thermal_deformation_heat_transfer.cpp)
+target_link_libraries(ThermalDeformationHeatTransfer LINK_PUBLIC CabanaPD)
+
+install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/thermomechanics/inputs/heat_transfer.json

@@ -0,0 +1,16 @@
+{
+    "num_cells"               : {"value": [49, 49, 49]},
+    "system_size"             : {"value": [0.1, 0.1, 0.1], "unit": "m"},
+    "density"                 : {"value": 8915,  "unit": "kg/m^3"},
+    "elastic_modulus"         : {"value": 115e+9, "unit": "Pa"},
+    "thermal_expansion_coeff" : {"value": 17e-6, "unit": "oC^{-1}"},
+    "thermal_conductivity"    : {"value": 387, "unit": "W/(m.K)"},
+    "specific_heat_capacity"  : {"value": 385, "unit": "J/(kg.K)"},
+    "reference_temperature"   : {"value": 100.0, "unit": "oC"},
+    "horizon"                 : {"value": 0.00615, "unit": "m"},  
+    "final_time"              : {"value": 8, "unit": "s"},
+    "timestep"                : {"value": 0.02, "unit": "s", "note": "This timestep is too large for mechanics to be stable. Use thermal_deformation_heat_transfer.json instead if coupled thermomechanics is desired."},
+    "thermal_subcycle_steps"  : {"value": 1},
+    "output_frequency"        : {"value": 10},
+    "output_reference"        : {"value": true}
+}

examples/thermomechanics/inputs/thermal_crack.json

@@ -4,13 +4,14 @@
     "density"                       : {"value": 3980, "unit": "kg/m^3"},
     "elastic_modulus"               : {"value": 370e+9, "unit": "Pa"},
     "fracture_energy"               : {"value": 24.32, "unit": "J/m^2"},
-    "thermal_coefficient"           : {"value": 7.5E-6, "unit": "oC^{-1}"},
+    "thermal_expansion_coeff"       : {"value": 7.5E-6, "unit": "oC^{-1}"},
     "reference_temperature"         : {"value": 300.0, "unit": "oC"},
     "background_temperature"        : {"value": 20.0, "unit": "oC"},
     "surface_temperature_ramp_time" : {"value": 0.001, "unit": "s"},
     "horizon"                       : {"value": 3.0e-4, "unit": "m"},  
     "final_time"                    : {"value": 7.5e-4, "unit": "s"},
     "timestep"                      : {"value": 7.5E-9, "unit": "s"},
+    "thermal_subcycle_steps"        : {"value": 100},
     "output_frequency"              : {"value": 200},
     "output_reference"              : {"value": true}
 }

examples/thermomechanics/inputs/thermal_deformation.json

@@ -1,13 +1,16 @@
 {
-    "num_cells"             : {"value": [100, 30, 3]},
-    "system_size"           : {"value": [1.0, 0.3, 0.03], "unit": "m"},
-    "density"               : {"value": 3980,  "unit": "kg/m^3"},
-    "elastic_modulus"       : {"value": 370e+9, "unit": "Pa"},
-    "thermal_coefficient"   : {"value": 7.5E-6, "unit": "oC^{-1}"},
-    "reference_temperature" : {"value": 20.0, "unit": "oC"},
-    "horizon"               : {"value": 0.03, "unit": "m"},  
-    "final_time"            : {"value": 0.01, "unit": "s"},
-    "timestep"              : {"value": 7.5E-7,  "unit": "s"},
-    "output_frequency"      : {"value": 100},
-    "output_reference"      : {"value": true}
+    "num_cells"               : {"value": [101, 31, 3]},
+    "system_size"             : {"value": [1.0, 0.3, 0.03], "unit": "m"},
+    "density"                 : {"value": 3980,  "unit": "kg/m^3"},
+    "elastic_modulus"         : {"value": 370e+9, "unit": "Pa"},
+    "thermal_expansion_coeff" : {"value": 7.5E-6, "unit": "oC^{-1}"},
+    "thermal_conductivity"    : {"value": 31, "unit": "W/(m.K)"},
+    "specific_heat_capacity"  : {"value": 880, "unit": "J/(kg.K)"},
+    "reference_temperature"   : {"value": 20.0, "unit": "oC"},
+    "horizon"                 : {"value": 0.03, "unit": "m"},  
+    "final_time"              : {"value": 0.01, "unit": "s"},
+    "timestep"                : {"value": 7.5E-7,  "unit": "s"},
+    "thermal_subcycle_steps": {"value": 100},
+    "output_frequency"        : {"value": 100},
+    "output_reference"        : {"value": true}
 }

examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json

@@ -0,0 +1,16 @@
+{
+    "num_cells"               : {"value": [49, 49, 49]},
+    "system_size"             : {"value": [0.1, 0.1, 0.1], "unit": "m"},
+    "density"                 : {"value": 8915,  "unit": "kg/m^3"},
+    "elastic_modulus"         : {"value": 115e+9, "unit": "Pa"},
+    "thermal_expansion_coeff" : {"value": 17e-6, "unit": "oC^{-1}"},
+    "thermal_conductivity"    : {"value": 387, "unit": "W/(m.K)"},
+    "specific_heat_capacity"  : {"value": 385, "unit": "J/(kg.K)"},
+    "reference_temperature"   : {"value": 100.0, "unit": "oC"},
+    "horizon"                 : {"value": 0.00615, "unit": "m"},  
+    "final_time"              : {"value": 8, "unit": "s"},
+    "timestep"                : {"value": 4.2e-07, "unit": "s"},
+    "thermal_subcycle_steps"  : {"value": 5e+4},
+    "output_frequency"        : {"value": 10},
+    "output_reference"        : {"value": true}
+}

examples/thermomechanics/thermal_crack.cpp

@@ -43,7 +43,7 @@ void thermalCrackExample( const std::string filename )
     double K = E / ( 3 * ( 1 - 2 * nu ) );
     double G0 = inputs["fracture_energy"];
     double delta = inputs["horizon"];
-    double alpha = inputs["thermal_coefficient"];
+    double alpha = inputs["thermal_expansion_coeff"];
 
     // Problem parameters
     double temp0 = inputs["reference_temperature"];

examples/thermomechanics/thermal_deformation.cpp

@@ -41,7 +41,7 @@ void thermalDeformationExample( const std::string filename )
     double nu = 0.25;
     double K = E / ( 3 * ( 1 - 2 * nu ) );
     double delta = inputs["horizon"];
-    double alpha = inputs["thermal_coefficient"];
+    double alpha = inputs["thermal_expansion_coeff"];
 
     // Problem parameters
     double temp0 = inputs["reference_temperature"];

examples/thermomechanics/thermal_deformation_heat_transfer.cpp

@@ -0,0 +1,158 @@
+/****************************************************************************
+ * Copyright (c) 2022 by Oak Ridge National Laboratory                      *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of CabanaPD. CabanaPD is distributed under a           *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// Simulate heat transfer in a pseudo-1d cube.
+void thermalDeformationHeatTransferExample( const std::string filename )
+{
+    // ====================================================
+    //             Use default Kokkos spaces
+    // ====================================================
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    // ====================================================
+    //                   Read inputs
+    // ====================================================
+    CabanaPD::Inputs inputs( filename );
+
+    // ====================================================
+    //            Material and problem parameters
+    // ====================================================
+    // Material parameters
+    double rho0 = inputs["density"];
+    double E = inputs["elastic_modulus"];
+    double nu = 0.25;
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double delta = inputs["horizon"];
+    double alpha = inputs["thermal_expansion_coeff"];
+    double kappa = inputs["thermal_conductivity"];
+    double cp = inputs["specific_heat_capacity"];
+
+    // Problem parameters
+    double temp0 = inputs["reference_temperature"];
+
+    // ====================================================
+    //                  Discretization
+    // ====================================================
+    // FIXME: set halo width based on delta
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+    std::array<int, 3> num_cells = inputs["num_cells"];
+    int m = std::floor( delta /
+                        ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) );
+    int halo_width = m + 1; // Just to be safe.
+
+    // ====================================================
+    //                Force model type
+    // ====================================================
+    using model_type = CabanaPD::PMB;
+    using thermal_type = CabanaPD::DynamicTemperature;
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Does not set displacements, velocities, etc.
+    auto particles =
+        std::make_shared<CabanaPD::Particles<memory_space, model_type,
+                                             typename thermal_type::base_type>>(
+            exec_space(), low_corner, high_corner, num_cells, halo_width );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    auto rho = particles->sliceDensity();
+    auto temp = particles->sliceTemperature();
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+        // Density
+        rho( pid ) = rho0;
+        // Temperature
+        temp( pid ) = temp0;
+    };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    auto force_model = CabanaPD::createForceModel(
+        model_type{}, CabanaPD::Elastic{}, *particles, delta, K, kappa, cp,
+        alpha, temp0 );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
+        inputs, particles, force_model );
+
+    // ====================================================
+    //                   Boundary condition
+    // ====================================================
+    // Temperature profile imposed on top and bottom surfaces
+    double dy = particles->dx[1];
+    using plane_type = CabanaPD::RegionBoundary<CabanaPD::RectangularPrism>;
+
+    // Top surface
+    plane_type plane1( low_corner[0], high_corner[0], high_corner[1] - dy,
+                       high_corner[1] + dy, low_corner[2], high_corner[2] );
+
+    // Bottom surface
+    plane_type plane2( low_corner[0], high_corner[0], low_corner[1] - dy,
+                       low_corner[1] + dy, low_corner[2], high_corner[2] );
+
+    // This is purposely delayed until after solver init so that ghosted
+    // particles are correctly taken into account for lambda capture here.
+    temp = particles->sliceTemperature();
+    auto temp_bc = KOKKOS_LAMBDA( const int pid, const double )
+    {
+        temp( pid ) = 0.0;
+    };
+
+    std::vector<plane_type> planes = { plane1, plane2 };
+    auto bc = CabanaPD::createBoundaryCondition(
+        temp_bc, exec_space{}, *particles, planes, false, 1.0 );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    cabana_pd->init( bc );
+    cabana_pd->run( bc );
+
+    // ====================================================
+    //                      Outputs
+    // ====================================================
+    // Output temperature along the y-axis
+    int profile_dim = 1;
+    auto value = KOKKOS_LAMBDA( const int pid ) { return temp( pid ); };
+    std::string file_name = "temperature_yaxis_profile.txt";
+    createOutputProfile( MPI_COMM_WORLD, num_cells[1], profile_dim, file_name,
+                         *particles, value );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    thermalDeformationHeatTransferExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

src/CabanaPD_ForceModels.hpp

@@ -12,6 +12,7 @@
 #ifndef FORCE_MODELS_H
 #define FORCE_MODELS_H
 
+#include <CabanaPD_Constants.hpp>
 #include <CabanaPD_Types.hpp>
 
 namespace CabanaPD
@@ -71,6 +72,45 @@ struct BaseTemperatureModel
     }
 };
 
+// This class stores temperature parameters needed for heat transfer, but not
+// the temperature itself (stored instead in the static temperature class
+// above).
+struct BaseDynamicTemperatureModel
+{
+    double delta;
+
+    double thermal_coeff;
+    double kappa;
+    double cp;
+    bool constant_microconductivity;
+
+    BaseDynamicTemperatureModel( const double _delta, const double _kappa,
+                                 const double _cp,
+                                 const bool _constant_microconductivity = true )
+    {
+        set_param( _delta, _kappa, _cp, _constant_microconductivity );
+    }
+
+    void set_param( const double _delta, const double _kappa, const double _cp,
+                    const bool _constant_microconductivity )
+    {
+        delta = _delta;
+        kappa = _kappa;
+        cp = _cp;
+        const double d3 = _delta * _delta * _delta;
+        thermal_coeff = 9.0 / 2.0 * _kappa / pi / d3;
+        constant_microconductivity = _constant_microconductivity;
+    }
+
+    KOKKOS_INLINE_FUNCTION double microconductivity_function( double r ) const
+    {
+        if ( constant_microconductivity )
+            return thermal_coeff;
+        else
+            return 4.0 * thermal_coeff * ( 1.0 - r / delta );
+    }
+};
+
 template <typename ModelType, typename DamageType,
           typename ThermalType = TemperatureIndependent, typename... DataTypes>
 struct ForceModel;