Skip to content

Commit

Permalink
Merge pull request #108 from streeve/heat_transfer
Browse files Browse the repository at this point in the history
Add heat transfer
  • Loading branch information
streeve authored Nov 6, 2024
2 parents bd587d2 + 0be8643 commit ec71208
Show file tree
Hide file tree
Showing 15 changed files with 609 additions and 37 deletions.
1 change: 1 addition & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion examples/thermomechanics/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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})
16 changes: 16 additions & 0 deletions examples/thermomechanics/inputs/heat_transfer.json
Original file line number Diff line number Diff line change
@@ -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}
}
3 changes: 2 additions & 1 deletion examples/thermomechanics/inputs/thermal_crack.json
Original file line number Diff line number Diff line change
Expand Up @@ -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}
}
25 changes: 14 additions & 11 deletions examples/thermomechanics/inputs/thermal_deformation.json
Original file line number Diff line number Diff line change
@@ -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}
}
Original file line number Diff line number Diff line change
@@ -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}
}
2 changes: 1 addition & 1 deletion examples/thermomechanics/thermal_crack.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"];
Expand Down
2 changes: 1 addition & 1 deletion examples/thermomechanics/thermal_deformation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"];
Expand Down
158 changes: 158 additions & 0 deletions examples/thermomechanics/thermal_deformation_heat_transfer.cpp
Original file line number Diff line number Diff line change
@@ -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();
}
40 changes: 40 additions & 0 deletions src/CabanaPD_ForceModels.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#ifndef FORCE_MODELS_H
#define FORCE_MODELS_H

#include <CabanaPD_Constants.hpp>
#include <CabanaPD_Types.hpp>

namespace CabanaPD
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit ec71208

Please sign in to comment.