diff --git a/README.md b/README.md index f8aecea3..968aafe6 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,7 @@ example can be run with: ./CabanaPD/build/install/bin/KalthoffWinkler CabanaPD/examples/inputs/kalthoff_winkler.json ``` -The third example is crack branching in soda-lime glass [3]. The example can be +The third example is crack branching in a pre-notched soda-lime glass plate due to traction loading [3]. The example can be run with: ``` @@ -138,13 +138,18 @@ run with: ### Thermomechanics Examples which demonstrate temperature-dependent mechanics and fracture are with `examples/thermomechanics`. -The first example demonstrates a thermoelastic problem without fracture with a homogeneous plate -under linear thermal loading [4]. The example can be run with: +The first example is thermoelastic deformation in a homogeneous plate due to linear thermal loading [4]. The example can be run with: ``` ./CabanaPD/build/install/bin/ThermalDeformation CabanaPD/examples/thermomechanics/thermal_deformation.json ``` +The second example is crack initiation and propagation in an alumina ceramic plate due to a thermal shock caused by water quenching [5]. The example can be run with: + +``` +./CabanaPD/build/install/bin/ThermalCrack CabanaPD/examples/thermomechanics/thermal_crack.json +``` + ## References [1] P. Seleson and D.J. Littlewood, Numerical tools for improved convergence @@ -161,6 +166,8 @@ Kunze, and L.W. Meyer, eds., Vol 1, DGM Informationsgesellschaft Verlag (1988) [4] D. He, D. Huang, and D. Jiang, Modeling and studies of fracture in functionally graded materials under thermal shock loading using peridynamics, Theoretical and Applied Fracture Mechanics 111 (2021): 102852. +[5] C.P. Jiang, X.F. Wu, J. Li, F. Song, Y.F. Shao, X.H. Xu, and P. Yan, A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock, Acta Materialia 60 (2012): 4540–4550. + ## Contributing We encourage you to contribute to CabanaPD! Please check the diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt index 84788a90..73a932b3 100644 --- a/examples/thermomechanics/CMakeLists.txt +++ b/examples/thermomechanics/CMakeLists.txt @@ -1,3 +1,7 @@ add_executable(ThermalDeformation thermal_deformation.cpp) target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD) -install(TARGETS ThermalDeformation DESTINATION ${CMAKE_INSTALL_BINDIR}) + +add_executable(ThermalCrack thermal_crack.cpp) +target_link_libraries(ThermalCrack LINK_PUBLIC CabanaPD) + +install(TARGETS ThermalDeformation ThermalCrack DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/inputs/thermal_crack.json b/examples/thermomechanics/inputs/thermal_crack.json new file mode 100644 index 00000000..4427f6c6 --- /dev/null +++ b/examples/thermomechanics/inputs/thermal_crack.json @@ -0,0 +1,16 @@ +{ + "num_cells" : {"value": [501, 101, 11]}, + "system_size" : {"value": [0.05, 0.01, 0.001], "unit": "m"}, + "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}"}, + "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"}, + "output_frequency" : {"value": 200}, + "output_reference" : {"value": true} +} diff --git a/examples/thermomechanics/thermal_crack.cpp b/examples/thermomechanics/thermal_crack.cpp new file mode 100644 index 00000000..28db3ef9 --- /dev/null +++ b/examples/thermomechanics/thermal_crack.cpp @@ -0,0 +1,169 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 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 +#include + +#include "mpi.h" + +#include + +#include + +// Simulate crack initiation and propagation in a ceramic plate under thermal +// shock caused by water quenching. +void thermalCrackExample( 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 G0 = inputs["fracture_energy"]; + double delta = inputs["horizon"]; + double alpha = inputs["thermal_coefficient"]; + + // Problem parameters + double temp0 = inputs["reference_temperature"]; + double temp_w = inputs["background_temperature"]; + double t_ramp = inputs["surface_temperature_ramp_time"]; + + // ==================================================== + // Discretization + // ==================================================== + // FIXME: set halo width based on delta + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array 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::TemperatureDependent; + + // ==================================================== + // Particle generation + // ==================================================== + // Does not set displacements, velocities, etc. + auto particles = std::make_shared< + CabanaPD::Particles>( + exec_space(), low_corner, high_corner, num_cells, halo_width ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles->sliceDensity(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) { rho( pid ) = rho0; }; + particles->updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Force model + // ==================================================== + auto force_model = + CabanaPD::createForceModel( model_type{}, CabanaPD::Fracture{}, + *particles, delta, K, G0, alpha, temp0 ); + + // ==================================================== + // Create solver + // ==================================================== + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model ); + + // -------------------------------------------- + // Thermal shock + // -------------------------------------------- + auto x = particles->sliceReferencePosition(); + auto temp = particles->sliceTemperature(); + + // Plate limits + double X0 = low_corner[0]; + double Xn = high_corner[0]; + double Y0 = low_corner[1]; + double Yn = high_corner[1]; + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + // Define a time-dependent surface temperature: + // An inverted triangular pulse over a 2*t_ramp period starting at temp0 + // and linearly decreasing to temp_w within t_ramp, then linearly + // increasing back to temp0, and finally staying constant at temp0 + double temp_infinity; + if ( t <= t_ramp ) + { + // Increasing pulse + temp_infinity = temp0 - ( temp0 - temp_w ) * ( t / t_ramp ); + } + else if ( t < 2 * t_ramp ) + { + // Decreasing pulse + temp_infinity = + temp_w + ( temp0 - temp_w ) * ( t - t_ramp ) / t_ramp; + } + else + { + // Constant value + temp_infinity = temp0; + } + + // Rescale x and y particle position values + double xi = ( 2.0 * x( pid, 0 ) - ( X0 + Xn ) ) / ( Xn - X0 ); + double eta = ( 2.0 * x( pid, 1 ) - ( Y0 + Yn ) ) / ( Yn - Y0 ); + + // Define profile powers in x- and y-directions + double sx = 1.0 / 50.0; + double sy = 1.0 / 10.0; + + // Define profiles in x- and y-direcions + double fx = 1.0 - Kokkos::pow( Kokkos::abs( xi ), 1.0 / sx ); + double fy = 1.0 - Kokkos::pow( Kokkos::abs( eta ), 1.0 / sy ); + + // Compute particle temperature + temp( pid ) = temp_infinity + ( temp0 - temp_infinity ) * fx * fy; + }; + auto body_term = CabanaPD::createBodyTerm( temp_func, false ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init( body_term ); + cabana_pd->run( body_term ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + thermalCrackExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 4dea6829..80a6309c 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -81,9 +81,8 @@ void thermalDeformationExample( const std::string filename ) // ==================================================== // Force model // ==================================================== - auto force_model = - CabanaPD::createForceModel( - *particles, delta, K, alpha, temp0 ); + auto force_model = CabanaPD::createForceModel( + model_type{}, CabanaPD::Elastic{}, *particles, delta, K, alpha, temp0 ); // ==================================================== // Create solver diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 2f614e8e..7799435e 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -62,6 +62,7 @@ #include +#include #include namespace CabanaPD @@ -128,6 +129,36 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, template class Force; +template +class Force +{ + protected: + bool _half_neigh; + + Timer _timer; + Timer _energy_timer; + + public: + Force( const bool half_neigh ) + : _half_neigh( half_neigh ) + { + } + + template + void computeWeightedVolume( ParticleType&, const NeighListType&, + const ParallelType ) const + { + } + template + void computeDilatation( ParticleType&, const NeighListType&, + const ParallelType ) const + { + } + + auto time() { return _timer.time(); }; + auto timeEnergy() { return _energy_timer.time(); }; +}; + /****************************************************************************** Force free functions. ******************************************************************************/ diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 2e60ad8d..fb3dc2d9 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -211,21 +211,19 @@ class SolverElastic void init( BoundaryType boundary_condition, const bool initial_output = true ) { + // Add non-force boundary condition. if ( !boundary_condition.forceUpdate() ) - { - // Add boundary condition. boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. - if constexpr ( std::is_same::value ) - comm->gatherTemperature(); - } + // Communicate temperature. + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); // Force init without particle output. init( false ); - // Add boundary condition. + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); @@ -248,14 +246,11 @@ class SolverElastic // Add non-force boundary condition. if ( !boundary_condition.forceUpdate() ) - { boundary_condition.apply( exec_space(), *particles, step * dt ); - if constexpr ( std::is_same< - typename force_model_type::thermal_type, - TemperatureDependent>::value ) - comm->gatherTemperature(); - } + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); // Update ghost particles. comm->gatherDisplacement(); @@ -442,7 +437,7 @@ class SolverElastic }; template + class ForceModel> class SolverFracture : public SolverElastic { @@ -459,12 +454,12 @@ class SolverFracture using force_model_type = ForceModel; using force_type = typename base_type::force_type; using neigh_iter_tag = Cabana::SerialOpTag; - using prenotch_type = PrenotchType; using input_type = typename base_type::input_type; + template SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, prenotch_type prenotch ) + force_model_type force_model, PrenotchType prenotch ) : base_type( _inputs, _particles, force_model ) { init_mu(); @@ -516,12 +511,19 @@ class SolverFracture void init( BoundaryType boundary_condition, const bool initial_output = true ) { + // Add non-force boundary condition. if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); + // Communicate temperature. + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + // Force init without particle output. init( false ); + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); @@ -544,7 +546,11 @@ class SolverFracture // Add non-force boundary condition. if ( !boundary_condition.forceUpdate() ) - boundary_condition.apply( exec_space{}, *particles, step * dt ); + boundary_condition.apply( exec_space(), *particles, step * dt ); + + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); // Update ghost particles. comm->gatherDisplacement(); @@ -567,27 +573,32 @@ class SolverFracture // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); - // Print output. - if ( step % output_frequency == 0 ) - { - auto W = computeEnergy( *force, *particles, *neighbors, mu, - neigh_iter_tag() ); - - particles->output( step / output_frequency, step * dt, - output_reference ); - _step_timer.stop(); - this->step_output( step, W ); - } - else - { - _step_timer.stop(); - } + output( step ); } // Final output and timings. this->final_output(); } + void output( const int step ) + { + // Print output. + if ( step % output_frequency == 0 ) + { + auto W = computeEnergy( *force, *particles, *neighbors, mu, + neigh_iter_tag() ); + + particles->output( step / output_frequency, step * dt, + output_reference ); + _step_timer.stop(); + this->step_output( step, W ); + } + else + { + _step_timer.stop(); + } + } + using base_type::dt; using base_type::num_steps; using base_type::output_frequency; @@ -621,14 +632,25 @@ auto createSolverElastic( InputsType inputs, inputs, particles, model ); } +template +auto createSolverFracture( InputsType inputs, + std::shared_ptr particles, + ForceModel model ) +{ + return std::make_shared< + SolverFracture>( + inputs, particles, model ); +} + template auto createSolverFracture( InputsType inputs, std::shared_ptr particles, ForceModel model, PrenotchType prenotch ) { - return std::make_shared>( + return std::make_shared< + SolverFracture>( inputs, particles, model, prenotch ); } diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 3380de37..7f920b24 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -90,6 +90,13 @@ struct ForceModel s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } + + KOKKOS_INLINE_FUNCTION + bool criticalStretch( const int, const int, const double r, + const double xi ) const + { + return r * r >= bond_break_coeff * xi * xi; + } }; template <> @@ -168,16 +175,82 @@ struct ForceModel } }; -template -auto createForceModel( ParticleType particles, const double delta, +template +auto createForceModel( PMB, Elastic, ParticleType particles, const double delta, const double K, const double alpha, const double temp0 ) { auto temp = particles.sliceTemperature(); using temp_type = decltype( temp ); - return ForceModel( + return ForceModel( delta, K, temp, alpha, temp0 ); } + +template +struct ForceModel + : public ForceModel, + BaseTemperatureModel +{ + using base_type = ForceModel; + using base_temperature_type = BaseTemperatureModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using thermal_type = TemperatureDependent; + + using base_type::c; + using base_type::delta; + using base_type::K; + + // Does not use the base bond_break_coeff. + using base_type::G0; + using base_type::s0; + + // Thermal parameters + using base_temperature_type::alpha; + using base_temperature_type::temp0; + using base_temperature_type::temperature; + + // Explicitly use the temperature-dependent stretch. + using base_temperature_type::thermalStretch; + + // ForceModel(){}; + ForceModel( const double _delta, const double _K, const double _G0, + const TemperatureType _temp, const double _alpha, + const double _temp0 = 0.0 ) + : base_type( _delta, _K, _G0 ) + , base_temperature_type( _temp, _alpha, _temp0 ) + { + set_param( _delta, _K, _G0, _alpha, _temp0 ); + } + + void set_param( const double _delta, const double _K, const double _G0, + const double _alpha, const double _temp0 ) + { + base_type::set_param( _delta, _K, _G0 ); + base_temperature_type::set_param( _alpha, _temp0 ); + } + + KOKKOS_INLINE_FUNCTION + bool criticalStretch( const int i, const int j, const double r, + const double xi ) const + { + double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0; + double bond_break_coeff = + ( 1.0 + s0 + alpha * temp_avg ) * ( 1.0 + s0 + alpha * temp_avg ); + return r * r >= bond_break_coeff * xi * xi; + } +}; + +template +auto createForceModel( PMB, Fracture, ParticleType particles, + const double delta, const double K, const double G0, + const double alpha, const double temp0 ) +{ + auto temp = particles.sliceTemperature(); + using temp_type = decltype( temp ); + return ForceModel( + delta, K, G0, temp, alpha, temp0 ); +} + } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 54991207..cd5acce7 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -63,6 +63,7 @@ #include #include +#include #include #include #include @@ -71,36 +72,27 @@ namespace CabanaPD { template class Force> + : public Force { public: using exec_space = ExecutionSpace; using model_type = ForceModel; + using base_type = Force; protected: - bool _half_neigh; + using base_type::_half_neigh; model_type _model; - Timer _timer; - Timer _energy_timer; + using base_type::_energy_timer; + using base_type::_timer; public: Force( const bool half_neigh, const model_type model ) - : _half_neigh( half_neigh ) + : base_type( half_neigh ) , _model( model ) { } - template - void computeWeightedVolume( ParticleType&, const NeighListType&, - const ParallelType ) const - { - } - template - void computeDilatation( ParticleType&, const NeighListType&, - const ParallelType ) const - { - } - template void computeForceFull( ForceType& f, const PosType& x, const PosType& u, @@ -181,22 +173,19 @@ class Force> _energy_timer.stop(); return strain_energy; } - - auto time() { return _timer.time(); }; - auto timeEnergy() { return _energy_timer.time(); }; }; template class Force> - : public Force> + : public Force { public: using exec_space = ExecutionSpace; - using model_type = ForceModel; + using model_type = ForceModel; protected: - using base_type = - Force>; + using base_model_type = typename model_type::base_type; + using base_type = Force; using base_type::_half_neigh; model_type _model; @@ -205,7 +194,7 @@ class Force> public: Force( const bool half_neigh, const model_type model ) - : base_type( half_neigh, model ) + : base_type( half_neigh ) , _model( model ) { } @@ -220,7 +209,6 @@ class Force> _timer.start(); auto model = _model; - auto break_coeff = _model.bond_break_coeff; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); @@ -247,7 +235,7 @@ class Force> model.thermalStretch( s, i, j ); // Break if beyond critical stretch unless in no-fail zone. - if ( r * r >= break_coeff * xi * xi && !nofail( i ) && + if ( model.criticalStretch( i, j, r, xi ) && !nofail( i ) && !nofail( j ) ) { mu( i, n ) = 0; @@ -329,14 +317,14 @@ class Force> template class Force> - : public Force> + : public Force { public: using exec_space = ExecutionSpace; - using model_type = ForceModel; + using model_type = ForceModel; protected: - using base_type = Force>; + using base_type = Force; using base_type::_half_neigh; model_type _model; @@ -345,7 +333,7 @@ class Force> public: Force( const bool half_neigh, const model_type model ) - : base_type( half_neigh, model ) + : base_type( half_neigh ) , _model( model ) { }