Merge pull request #102 from streeve/thermal_cracking

Thermal cracking

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

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})

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}
+}

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 <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// 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<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::TemperatureDependent;
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Does not set displacements, velocities, etc.
+    auto particles = std::make_shared<
+        CabanaPD::Particles<memory_space, model_type, thermal_type>>(
+        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<memory_space>(
+        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();
+}

examples/thermomechanics/thermal_deformation.cpp

@@ -81,9 +81,8 @@ void thermalDeformationExample( const std::string filename )
     // ====================================================
     //                    Force model
     // ====================================================
-    auto force_model =
-        CabanaPD::createForceModel<model_type, CabanaPD::Elastic, thermal_type>(
-            *particles, delta, K, alpha, temp0 );
+    auto force_model = CabanaPD::createForceModel(
+        model_type{}, CabanaPD::Elastic{}, *particles, delta, K, alpha, temp0 );
 
     // ====================================================
     //                   Create solver

src/CabanaPD_Force.hpp

@@ -62,6 +62,7 @@
 
 #include <cmath>
 
+#include <CabanaPD_ForceModels.hpp>
 #include <CabanaPD_Particles.hpp>
 
 namespace CabanaPD
@@ -128,6 +129,36 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i,
 template <class ExecutionSpace, class ForceType>
 class Force;
 
+template <class ExecutionSpace>
+class Force<ExecutionSpace, BaseForceModel>
+{
+  protected:
+    bool _half_neigh;
+
+    Timer _timer;
+    Timer _energy_timer;
+
+  public:
+    Force( const bool half_neigh )
+        : _half_neigh( half_neigh )
+    {
+    }
+
+    template <class ParticleType, class NeighListType, class ParallelType>
+    void computeWeightedVolume( ParticleType&, const NeighListType&,
+                                const ParallelType ) const
+    {
+    }
+    template <class ParticleType, class NeighListType, class ParallelType>
+    void computeDilatation( ParticleType&, const NeighListType&,
+                            const ParallelType ) const
+    {
+    }
+
+    auto time() { return _timer.time(); };
+    auto timeEnergy() { return _energy_timer.time(); };
+};
+
 /******************************************************************************
   Force free functions.
 ******************************************************************************/