Add HertzianContact with normal spring and damping forces;

includes new test with physically realistic values for 316L

src/CabanaPD.hpp

@@ -33,7 +33,9 @@
 #include <force/CabanaPD_ForceModels_LPS.hpp>
 #include <force/CabanaPD_ForceModels_PMB.hpp>
 #include <force/CabanaPD_Force_Contact.hpp>
+#include <force/CabanaPD_Force_HertzianContact.hpp>
 #include <force/CabanaPD_Force_LPS.hpp>
 #include <force/CabanaPD_Force_PMB.hpp>
+#include <force/CabanaPD_HertzianContact.hpp>
 
 #endif

src/force/CabanaPD_ContactModels.hpp

@@ -20,6 +20,24 @@
 
 namespace CabanaPD
 {
+/******************************************************************************
+  Contact model helper functions
+******************************************************************************/
+template <class VelType>
+KOKKOS_INLINE_FUNCTION void
+getRelativeNormalVelocity( const VelType& vel, const int i, const int j,
+                           const double rx, const double ry, const double rz,
+                           const double r, double& vx, double& vy, double& vz,
+                           double& vn )
+{
+    vx = vel( i, 0 ) - vel( j, 0 );
+    vy = vel( i, 1 ) - vel( j, 1 );
+    vz = vel( i, 2 ) - vel( j, 2 );
+
+    vn = vx * rx + vy * ry + vz * rz;
+    vn /= r;
+};
+
 /******************************************************************************
   Contact model
 ******************************************************************************/

src/force/CabanaPD_Force_HertzianContact.hpp

@@ -0,0 +1,163 @@
+/****************************************************************************
+ * 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                                    *
+ ****************************************************************************/
+
+#ifndef CONTACT_HERTZIAN_H
+#define CONTACT_HERTZIAN_H
+
+#include <cmath>
+
+#include <CabanaPD_Force.hpp>
+#include <CabanaPD_Input.hpp>
+#include <CabanaPD_Output.hpp>
+#include <force/CabanaPD_HertzianContact.hpp>
+
+namespace CabanaPD
+{
+/******************************************************************************
+  Normal repulsion forces
+******************************************************************************/
+template <class MemorySpace>
+class Force<MemorySpace, HertzianModel>
+    : public Force<MemorySpace, BaseForceModel>
+{
+  public:
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+
+    template <class ParticleType>
+    Force( const bool half_neigh, const ParticleType& particles,
+           const HertzianModel model )
+        : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(),
+                     particles.n_local, particles.ghost_mesh_lo,
+                     particles.ghost_mesh_hi )
+        , _model( model )
+    {
+        for ( int d = 0; d < particles.dim; d++ )
+        {
+            mesh_min[d] = particles.ghost_mesh_lo[d];
+            mesh_max[d] = particles.ghost_mesh_hi[d];
+        }
+    }
+
+    template <class ForceType, class PosType, class ParticleType,
+              class ParallelType>
+    void computeForceFull( ForceType& fc, const PosType& x, const PosType& u,
+                           const ParticleType& particles, const int n_local,
+                           ParallelType& neigh_op_tag )
+    {
+        auto Rc = _model.Rc;
+        auto radius = _model.radius;
+        auto Es = _model.Es;
+        auto Rs = _model.Rs;
+        auto beta = _model.beta;
+
+        const double coeff_h_n = 4.0 / 3.0 * Es * std::sqrt( Rs );
+        const double coeff_h_d = -2.0 * sqrt( 5.0 / 6.0 ) * beta;
+
+        const auto vol = particles.sliceVolume();
+        const auto rho = particles.sliceDensity();
+        const auto y = particles.sliceCurrentPosition();
+        const auto vel = particles.sliceVelocity();
+
+        _neigh_timer.start();
+        _neigh_list.build( y, 0, n_local, Rc, 1.0, mesh_min, mesh_max );
+        _neigh_timer.stop();
+
+        auto contact_full = KOKKOS_LAMBDA( const int i, const int j )
+        {
+            using Kokkos::abs;
+            using Kokkos::min;
+            using Kokkos::pow;
+            using Kokkos::sqrt;
+
+            double fcx_i = 0.0;
+            double fcy_i = 0.0;
+            double fcz_i = 0.0;
+
+            double xi, r, s;
+            double rx, ry, rz;
+            getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz );
+
+            // Contact "overlap"
+            const double delta_n = ( r - 2.0 * radius );
+
+            // Hertz normal force coefficient
+            double coeff = 0.0;
+            double Sn = 0.0;
+            if ( delta_n < 0.0 )
+            {
+                coeff =
+                    min( 0.0, -coeff_h_n * pow( abs( delta_n ), 3.0 / 2.0 ) );
+                Sn = 2.0 * Es * sqrt( Rs * abs( delta_n ) );
+            }
+
+            coeff /= vol( i );
+
+            fcx_i = coeff * rx / r;
+            fcy_i = coeff * ry / r;
+            fcz_i = coeff * rz / r;
+
+            fc( i, 0 ) += fcx_i;
+            fc( i, 1 ) += fcy_i;
+            fc( i, 2 ) += fcz_i;
+
+            // Hertz normal force damping component
+            double vx, vy, vz, vn;
+            getRelativeNormalVelocity( vel, i, j, rx, ry, rz, r, vx, vy, vz,
+                                       vn );
+
+            double ms = ( rho( i ) * vol( i ) ) / 2.0;
+            double fnd = coeff_h_d * sqrt( Sn * ms ) * vn / vol( i );
+
+            fcx_i = fnd * rx / r;
+            fcy_i = fnd * ry / r;
+            fcz_i = fnd * rz / r;
+
+            fc( i, 0 ) += fcx_i;
+            fc( i, 1 ) += fcy_i;
+            fc( i, 2 ) += fcz_i;
+        };
+
+        _timer.start();
+
+        // FIXME: using default space for now.
+        using exec_space = typename MemorySpace::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, n_local );
+        Cabana::neighbor_parallel_for(
+            policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(),
+            neigh_op_tag, "CabanaPD::Contact::compute_full" );
+
+        _timer.stop();
+    }
+
+    // FIXME: implement energy
+    template <class PosType, class WType, class ParticleType,
+              class ParallelType>
+    double computeEnergyFull( WType&, const PosType&, const PosType&,
+                              ParticleType&, const int, ParallelType& )
+    {
+        return 0.0;
+    }
+
+  protected:
+    HertzianModel _model;
+    using base_type::_half_neigh;
+    using base_type::_neigh_list;
+    using base_type::_timer;
+    Timer _neigh_timer;
+
+    double mesh_max[3];
+    double mesh_min[3];
+};
+
+} // namespace CabanaPD
+
+#endif

src/force/CabanaPD_HertzianContact.hpp

@@ -0,0 +1,70 @@
+/****************************************************************************
+ * 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                                    *
+ ****************************************************************************/
+
+#ifndef CONTACTMODEL_HERTZIAN_H
+#define CONTACTMODEL_HERTZIAN_H
+
+#include <cmath>
+
+#include <CabanaPD_Force.hpp>
+#include <CabanaPD_Input.hpp>
+#include <CabanaPD_Output.hpp>
+
+namespace CabanaPD
+{
+struct HertzianModel : public ContactModel
+{
+    // FIXME: This is for use as the primary force model.
+    using base_model = PMB;
+    using fracture_type = Elastic;
+    using thermal_type = TemperatureIndependent;
+
+    using ContactModel::Rc; // Contact horizon (should be > 2*radius)
+
+    double nu;     // Poisson's ratio
+    double radius; // Actual radius
+    double Rs;     // Equivalent radius
+    double Es;     // Equivalent Young's modulus
+    double e;      // Coefficient of restitution
+    double beta;   // Damping coefficient
+
+    HertzianModel( const double _Rc, const double _radius, const double _nu,
+                   const double _E, const double _e )
+        : ContactModel( 1.0, _Rc )
+    {
+        set_param( _radius, _nu, _E, _e );
+    }
+
+    void set_param( const double _radius, const double _nu, const double _E,
+                    const double _e )
+    {
+        using Kokkos::log;
+        using Kokkos::pow;
+        using Kokkos::sqrt;
+
+        nu = _nu;
+        radius = _radius;
+        Rs = 0.5 * radius;
+        Es = _E / ( 2.0 * pow( 1.0 - nu, 2.0 ) );
+        e = _e;
+        double ln_e = log( e );
+        beta = -ln_e / sqrt( pow( ln_e, 2.0 ) + pow( pi, 2.0 ) );
+    }
+};
+
+template <>
+struct is_contact<HertzianModel> : public std::true_type
+{
+};
+
+} // namespace CabanaPD
+
+#endif

unit_test/CMakeLists.txt

@@ -1,3 +1,8 @@
+configure_file(
+  ${CMAKE_CURRENT_SOURCE_DIR}/inputs/hertzian_contact.json
+  ${CMAKE_CURRENT_BINARY_DIR}/hertzian_contact.json
+  COPYONLY
+)
 ##--------------------------------------------------------------------------##
 ## On-node tests
 ##--------------------------------------------------------------------------##
@@ -48,6 +53,6 @@ macro(CabanaPD_add_tests)
   endforeach()
 endmacro()
 
-CabanaPD_add_tests(NAMES Force Integrator)
+CabanaPD_add_tests(NAMES Force Integrator Hertz)
 
 CabanaPD_add_tests(MPI NAMES Comm)

unit_test/inputs/hertzian_contact.json

@@ -0,0 +1,16 @@
+{
+    "num_cells"              : {"value": [1, 1, 1]},
+    "system_size"            : {"value": [1.0e-3, 1.0e-3, 1.0e-3], "unit": "m"},
+    "density"                : {"value": 7.95e3,  "unit": "kg/m^3"},
+    "volume"                 : {"value": 5.236e-13,  "unit": "m^3"},
+    "elastic_modulus"        : {"value": 195.6e9, "unit": "Pa"},
+    "poisson_ratio"          : {"value": 0.25, "unit": ""},
+    "restitution"            : {"value": 0.5},  
+    "horizon"                : {"value": 1e-4 },  
+    "radius"                 : {"value": 5e-5 },
+    "final_time"             : {"value": 1e-5, "unit": "s"},
+    "timestep"               : {"value": 1e-8,  "unit": "s"},
+    "timestep_safety_factor" : {"value": 0.8},
+    "output_frequency"       : {"value": 10000},
+    "output_reference"       : {"value": false}
+}