Plasticity

examples/mechanics/CMakeLists.txt

@@ -13,4 +13,7 @@ target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD)
 add_executable(RandomCracks random_cracks.cpp)
 target_link_libraries(RandomCracks LINK_PUBLIC CabanaPD)
 
-install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DESTINATION ${CMAKE_INSTALL_BINDIR})
+add_executable(CompactTensionTest compact_tension_test.cpp)
+target_link_libraries(CompactTensionTest LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks CompactTensionTest DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/mechanics/compact_tension_test.cpp

@@ -0,0 +1,203 @@
+/****************************************************************************
+ * 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 crack propagation in a compact tension test.
+void crackBranchingExample( 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 parameters
+    // ====================================================
+    double rho0 = inputs["density"];
+    double E = inputs["elastic_modulus"];
+    double nu = 0.25; // Use bond-based model
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double G0 = inputs["fracture_energy"];
+    double delta = inputs["horizon"];
+    delta += 1e-10;
+
+    // ====================================================
+    //                  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"];
+
+    // ====================================================
+    //                    Pre-notch
+    // ====================================================
+    double height = inputs["system_size"][0];
+    double thickness = inputs["system_size"][2];
+    double L_prenotch = height / 2.0;
+    double y_prenotch1 = 0.0;
+    Kokkos::Array<double, 3> p01 = { low_corner[0], y_prenotch1,
+                                     low_corner[2] };
+    Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
+    Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
+    Kokkos::Array<Kokkos::Array<double, 3>, 1> notch_positions = { p01 };
+    CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions );
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type =
+        CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::ElasticPerfectlyPlastic>;
+    model_type force_model( delta, K, G0 );
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Note that individual inputs can be passed instead (see other examples).
+    auto particles = CabanaPD::createParticles<memory_space, model_type>(
+        exec_space{}, inputs );
+
+    // ====================================================
+    //    Custom particle generation and initialization
+    // ====================================================
+
+    // Rectangular prism containing the full specimen: original geometry
+    CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
+        low_corner[0], high_corner[0], low_corner[1], high_corner[1],
+        low_corner[2], high_corner[2] );
+    std::vector<CabanaPD::RegionBoundary<CabanaPD::RectangularPrism>> planes = {
+        plane };
+
+    // Geometric parameters of specimen
+    double L = inputs["system_size"][1];
+    double W = L / 1.25;
+    double a = 0.45 * W;
+
+    // Grid spacing in y-direction
+    double dy = particles->dx[1];
+
+    // Remove particles from original geometry
+    auto init_op = KOKKOS_LAMBDA( const int, const double x[3] )
+    {
+        // Thin rectangle
+        if ( x[0] < low_corner[1] + 0.25 * W + a &&
+             Kokkos::abs( x[1] ) < 0.5 * dy )
+        {
+            return false;
+        }
+        // Thick rectangle
+        else if ( x[0] < low_corner[1] + 0.25 * W &&
+                  Kokkos::abs( x[1] ) < 25e-4 )
+        {
+            return false;
+        }
+        else
+        {
+            return true;
+        }
+    };
+
+    particles->createParticles( exec_space(), init_op );
+
+    auto rho = particles->sliceDensity();
+    auto x = particles->sliceReferencePosition();
+    auto v = particles->sliceVelocity();
+    auto f = particles->sliceForce();
+    auto nofail = particles->sliceNoFail();
+
+    // Pin radius
+    double R = 4e-3;
+    // Pin center coordinates (top)
+    double x_pin = low_corner[0] + 0.25 * W;
+    double y_pin = 0.37 * W;
+    // Pin velocity magnitude
+    double v0 = inputs["pin_velocity"];
+
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+        // Density
+        rho( pid ) = rho0;
+
+        auto xpos = x( pid, 0 );
+        auto ypos = x( pid, 1 );
+        auto distsq =
+            ( xpos - x_pin ) * ( xpos - x_pin ) +
+            ( Kokkos::abs( ypos ) - y_pin ) * ( Kokkos::abs( ypos ) - y_pin );
+        auto sign = Kokkos::abs( ypos ) / ypos;
+
+        // pins' y-velocity
+        if ( distsq < R * R )
+            v( pid, 1 ) = sign * v0;
+
+        // No-fail zone
+        if ( distsq < ( 2 * R ) * ( 2 * R ) )
+            nofail( pid ) = 1;
+    };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    auto cabana_pd =
+        CabanaPD::createSolver<memory_space>( inputs, particles, force_model );
+
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+
+    // Create BC last to ensure ghost particles are included.
+    f = particles->sliceForce();
+    x = particles->sliceReferencePosition();
+    // Create a symmetric force BC in the y-direction.
+    auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
+    {
+        auto xpos = x( pid, 0 );
+        auto ypos = x( pid, 1 );
+        if ( ( xpos - x_pin ) * ( xpos - x_pin ) +
+                 ( Kokkos::abs( ypos ) - y_pin ) *
+                     ( Kokkos::abs( ypos ) - y_pin ) <
+             R * R )
+            f( pid, 1 ) = 0;
+    };
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes,
+                                       true );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    cabana_pd->init( prenotch );
+    cabana_pd->run( bc );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    crackBranchingExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

examples/mechanics/inputs/compact_tension_test.json

@@ -0,0 +1,15 @@
+{
+    "num_cells"              : {"value": [187, 181, 15]},
+    "system_size"            : {"value": [0.12625, 0.1212, 0.01], "unit": "m"},
+    "density"                : {"value": 2440,    "unit": "kg/m^3"},
+    "elastic_modulus"        : {"value": 72.6e+9, "unit": "Pa"},
+    "fracture_energy"        : {"value": 16160,   "unit": "J/m^2"},
+    "horizon"                : {"value": 0.00203, "unit": "m"},   
+    "yield_stretch"          : {"value": 304e6,   "unit": "Pa"},
+    "pin_velocity"           : {"value": 7.15e6,  "unit": "m/s"},
+    "final_time"             : {"value": 43e-6,   "unit": "s"},
+    "timestep"               : {"value": 4.5e-8,  "unit": "s"},
+    "timestep_safety_factor" : {"value": 0.85},
+    "output_frequency"       : {"value": 5},
+    "output_reference"       : {"value": true}
+}

src/CabanaPD_Types.hpp

@@ -36,6 +36,9 @@ struct is_fracture<Fracture> : public std::true_type
 struct Elastic
 {
 };
+struct ElasticPerfectlyPlastic
+{
+};
 
 // Contact and DEM (contact without PD) tags.
 struct NoContact

src/force/CabanaPD_ForceModels_PMB.hpp

@@ -56,13 +56,60 @@ struct ForceModel<PMB, Elastic, NoFracture, TemperatureIndependent>
     }
 };
 
+template <>
+struct ForceModel<PMB, ElasticPerfectlyPlastic, NoFracture,
+                  TemperatureIndependent>
+    : public ForceModel<PMB, Elastic, NoFracture, TemperatureIndependent>
+{
+    using base_type = ForceModel<PMB, Elastic, NoFracture>;
+    using base_model = PMB;
+    using fracture_type = NoFracture;
+    using plasticity_type = ElasticPerfectlyPlastic;
+    using thermal_type = TemperatureIndependent;
+
+    using base_type::c;
+    using base_type::delta;
+    using base_type::K;
+    // FIXME: hardcoded
+    const double s_Y = 0.0014;
+
+    ForceModel( const double delta, const double K )
+        : base_type( delta, K )
+    {
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    auto forceCoeff( const double s, const double vol ) const
+    {
+        if ( s < s_Y )
+            return c * s * vol;
+        else
+            return c * s_Y * vol;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    auto energy( const double s, const double xi, const double vol ) const
+    {
+        double stretch_term = 0.0;
+        if ( s < s_Y )
+            stretch_term = s * s;
+        else
+            stretch_term = s_Y * ( 2 * s - s_Y );
+
+        // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
+        // the integrand (pairwise potential).
+        return 0.25 * c * stretch_term * xi * vol;
+    }
+};
+
 template <>
 struct ForceModel<PMB, Elastic, Fracture, TemperatureIndependent>
     : public ForceModel<PMB, Elastic, NoFracture, TemperatureIndependent>
 {
     using base_type = ForceModel<PMB, Elastic, NoFracture>;
     using base_model = typename base_type::base_model;
     using fracture_type = Fracture;
+    using mechanics_type = Elastic;
     using thermal_type = base_type::thermal_type;
 
     using base_type::c;
@@ -88,13 +135,65 @@ struct ForceModel<PMB, Elastic, Fracture, TemperatureIndependent>
     }
 };
 
+template <>
+struct ForceModel<PMB, ElasticPerfectlyPlastic, Fracture,
+                  TemperatureIndependent>
+    : public ForceModel<PMB, Elastic, Fracture, TemperatureIndependent>
+{
+    using base_type = ForceModel<PMB, Elastic>;
+    using base_model = typename base_type::base_model;
+    using fracture_type = Fracture;
+    using mechanics_type = ElasticPerfectlyPlastic;
+    using thermal_type = base_type::thermal_type;
+
+    // Purposely not using the (static) base class bond_break_coeff
+    using base_type::c;
+    using base_type::delta;
+    using base_type::G0;
+    using base_type::K;
+    using base_type::s0;
+
+    // FIXME: hardcoded
+    const double s_Y = 0.0014;
+
+    ForceModel( const double delta, const double K, const double G0 )
+        : base_type( delta, K, G0 )
+    {
+    }
+
+    // FIXME: avoiding multiple inheritance..
+    KOKKOS_INLINE_FUNCTION
+    auto forceCoeff( const double s, const double vol ) const
+    {
+        if ( s < s_Y )
+            return c * s * vol;
+        else
+            return c * s_Y * vol;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    auto energy( const double s, const double xi, const double vol ) const
+    {
+        double stretch_term = 0.0;
+        if ( s < s_Y )
+            stretch_term = s * s;
+        else
+            stretch_term = s_Y * ( 2 * s - s_Y );
+
+        // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
+        // the integrand (pairwise potential).
+        return 0.25 * c * stretch_term * xi * vol;
+    }
+};
+
 template <>
 struct ForceModel<LinearPMB, Elastic, NoFracture, TemperatureIndependent>
     : public ForceModel<PMB, Elastic, NoFracture, TemperatureIndependent>
 {
     using base_type = ForceModel<PMB, Elastic, NoFracture>;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
+    using mechanics_type = Elastic;
     using thermal_type = base_type::thermal_type;
 
     using base_type::base_type;
@@ -111,6 +210,7 @@ struct ForceModel<LinearPMB, Elastic, Fracture, TemperatureIndependent>
     using base_type = ForceModel<PMB>;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
+    using mechanics_type = Elastic;
     using thermal_type = base_type::thermal_type;
 
     using base_type::base_type;
@@ -187,6 +287,7 @@ struct ForceModel<PMB, Elastic, Fracture, TemperatureDependent, TemperatureType>
     using base_temperature_type = BaseTemperatureModel<TemperatureType>;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
+    using mechanics_type = Elastic;
     using thermal_type = TemperatureDependent;
 
     using base_type::c;