Modeling of the Split-Hopkinson-Pressure-Bar experiment with the explicit material point method

Abstract

The material point method (MPM) represents an alternative discretization method for numerical simulations. It aims to combine the benefits of a Lagrangian representation of bodies and an Eulerian numerical solution approach. Therefore, especially at high material deformations the method is not prone to mesh distortions such as the finite element method (FEM). For this reason, the MPM is used to a great extent for modeling granular materials as in geo-mechanics. However, high deformations occur in many industrial processes on metallic materials. The Split-Hopkinson-Pressure-Bar (SHPB) experiment is used to characterize material properties at high deformation rates. Although widely used, this experiment is not yet standardized and shows a variety of sensitivities, e.g. to friction. Inter alia for this reason, simulations are conducted with the experiment to allow for a better evaluation of the measured data. The purpose of this work from an engineering point of view is to analyze the performance of the MPM on an SHPB experiment. In order to validate the experimental results for the material characterization under dynamic loading conditions we introduce frictional contact. We use arbitrary tri-linear brick domains in a 3D CPDI1 scheme, instead of originally used parallelepipeds. This allows for a more flexible geometry approximation using standard meshes. The results of the method are analyzed with respect to discretization sensitivity and discussed in the context of the experimental results for a 42CrMo4 steel. We were able to show that the method is capable to reproduce the SHPB experiment. Additionally the method shows convergency in the results with finer discretizations. Thus, the MPM has underlined its importance as an alternative simulation technique for problems with high deformation.

1 Introduction

The material point method (MPM) represents a mesh-based particle method. Its origins were made with the particle in cell method which proposed the idea of combining the advantages of a Eulerian computational frame with a Lagrangian representation of physical bodies, see Sulsky et al. [1]. In contrast to Lagrangian mesh-based methods, the MPM does not suffer from problems associated with distorted meshes which usually occur at very large deformations. The actual balance equations are solved on a computational background grid in an Eulerian form. However, conservation of mass is inherently satisfied as mass is carried with the particles. The kinematic and constitutive properties are attached to the material points constituting the Lagrangian representation of the body. In an MPM simulation, the material points are moving in space, independent of the computational grid. Therefore, the most common choice is to reset the computational grid at each time/load step. That means the position of grid nodes is the same at the beginning of all time steps and mesh distortion is completely avoided. Furthermore, a regular grid pattern leads to a performance advantage, respectively.

In general, an MPM algorithm has three main steps. In the first, we map the relevant quantities from the material points onto the computational grid, dependent on their location. Essentially this assembles the required vectors and matrices on the grid. In the second step the balance equations on the grid are solved, defined by the grid vectors and matrices, assembled before. The actual solution procedure, as well as the required vectors and matrices, rely on the choice of the solution scheme. In the third step of the algorithm, the solution of the grid is mapped down to the particles, which move in space, and update their kinematical and constitutive states based on these grid information.

This generic scheme is subject to several different approaches. In a broad scope, there exist implicit solution schemes such as for example in [2,3,4,5], but also explicit solution schemes are widely applied, see, e.g. [6, 7] or [8]. Different update strategies for the material points can be found, such as the Update-Stress-Last (USL), or the contrary Update-Stress-First (USF) scheme, to name just a few. A more detailed overview can be found in, e.g. [7, 9]. Time-dependent, dynamic problems require a temporal discretization, where the MPM uses standard techniques, such as implicit/explicit Euler forms, or central difference schemes such as the Leapfrog or Velocity-Verlet time integration. Impact problems (similar to the SHPB experiment) represent a well-suited use case for MPM because of the high deformations that occur, as shown, e.g. in [6, 10].

A research focus has been the mapping of material point properties onto the computational grid and back to material points. Stability problems associated with poor mapping operators were referred to as the cell-crossing noise and also unphysical discontinuities under tension in solid mechanics which set the motivation for these investigations. To date, the greatest coverage is set by the generalized interpolation material point (GIMP) method, introduced by Bardenhagen and Kober [8] and varieties of convected particle domain interpolation (CPDI) methods as proposed in Sadeghirad et al. [11]. A more detailed overview of this subject is given, e.g. in Nguyen et al. [12].

Right from the beginning, special interests were given to contact modeling using MPM. As all material points on the same computational grid share the velocity field in standard MPM, the non-penetration condition is fulfilled a priori, as early reported by Sulsky et al. [1]. As no dedicated contact search is needed, the resulting implementations are very efficient. This approach however satisfies a perfect stick condition on body interfaces. This give reason to further development as applications often require constitutive behaviour of contact which is not possible with the standard MPM. Several approaches have been studied to overcome this issue, following the idea of only mapping separating force components at a contact interface, see, e.g. [7, 13, 14] to name just a few. Most challenging in this approaches is the computation of the contact normal vector, as the bodies are not represented by an explicit geometry but rather with their properties projected on the computational grid. This ansatz has been followed intensively for modeling contact between MPM bodies, with recent advances in Nairn et al. [15]. Specialised contact models for geotechnical applications using a penalty function in order to reduce oscillation were discussed in Ma et al. [16]. Another approach uses grid boundary conditions to establish contact against rigid surfaces, ranging from first applications as in Sulsky and Schreyer [17] to more complex cases, e.g. inclined surfaces with physical contact models as in Li et al. [18].

Particle methods are a great alternative to Lagrangian methods like the FEM where the meshes are prone to numerical failures. Using particle methods it is possible to model the behaviour of laser melting, as shown, e.g. in Zohdi [19] or Fürstenau et al. [20] and the references therein. To date, the material point method is used in several applications. In solid mechanics, the method is widely used for modeling granular materials such as soil in geo-mechanics, see, e.g. Wang et al. [5], but also for snow as in Li et al. [18]. Also for metallic materials, several works exist. Sulsky and Schreyer [17] investigated the method in axisymmetric problems such as a Taylor-Bar impact experiment and a steady-state variant of a compression test. In this case, the authors rely on the inherit MPM contact or grid boundary conditions, which do not follow an actual contact law. Lian et al. [6] also used the MPM on metallic materials for impact problems.

The Split-Hopkinson-Pressure-Bar (SHPB) is an experimental technique to determine material parameters and to analyze the deformation behavior at high strain rates. Also the effect of varying temperatures can be observed. Usually the gathered data are used to determine stress–strain curves which are used to fit the material model of choice. Multiple variants are known where the specimen is subject to rapid tension or compression as shown in Al-Mousawi et al. [21]. Also, different material types can be tested such as brittle and ductile materials, e.g. concrete as in Kveton [22]. In a nutshell, the experimental setup consists of a specimen which is clamped between two bars of material which behave elastically under the expected loads. For the test, one of the elastic bars is subject to an impact shock, e.g. driven by a striker bar, accelerated in a gas gun. The resulting stress wave will move through this incident bar until a fraction is reflected at the interface to the specimen. Another fraction of the incident wave is transferred into the specimen, deforming it elasto-plastically. On the other side, the wave is further transmitted into the transmission bar. These waves are measured in both bars, using strain gauges allowing to apply simple mechanics, to recover the stress–strain curve of the specimen. This is mainly possible through the simplifying assumption of a one-dimensional problem, neglecting effects such as friction, wave dispersion, and three-dimensional constitutive phenomena such as material anisotropy. Although this test is often applied in material testing, see Chen and Song [23], it is not yet standardized for a variety of materials. This leads to an undesired variation when material parameters are determined from the test signal, as recently investigated by Kariem et al. [24]. Reason of these variations might stem from the filtering techniques which are applied to the raw signal or from the neglected phenomena as, e.g. friction, see Hartley et al. [25]. Numerical simulations are usually conducted to detect these influences, to minimize the experimental expenses. Particularly in the context of inverse methods, simulations are used to fit a specific material law, as shown, e.g. in Afdahl et al. [26] or Manes et al. [27].

In this paper, an MPM implementation for solid mechanics in 3D is presented. An explicit solution scheme based on a Leapfrog time integration scheme is used to solve the balance of momentum on a regular computational grid composed of tri-linear cuboid cells. The material points carry a convected particle domain of first-order (CPDI). A contact formulation for MPM discretized bodies and moving, rigid surfaces is used to establish frictional contact based on Coulomb’s law. A finite J$_2$ von Mises plasticity law is used on the material point level as a constitutive equation.

With this computational framework, a numerical simulation of an SHPB experiment is presented. The numerical model is analyzed for the sensitivity of selected results with respect to friction and discretization density. Furthermore, the simulation results are compared to experimental data, showing the applicability of the MPM for the SHPB. To the author’s knowledge this is the first report on using the MPM to simulate metallic specimens in a SHPB. Also the computational framework, using a 3D CPDI in combination with the rigid body penalty contact is unique to this point. The aim of this work is to analyze the MPM for its applicability to engineering problems. The investigations made, implies that the MPM is a reliable tool for supporting material analysis under difficult conditions, such as observed in an SHPB experiment.

Throughout this paper, the subscript $(\bullet )_{{\texttt {MP}}}$ is used to indicate a quantity at a material point, itself indicated by ${\texttt {MP}}$. The superscript $(\bullet )^{I}$ denotes nodes I of the computational grid. Tensorial quantities are typeset in bold-face, while scalars are depicted by italic characters. Furthermore we denote a single contraction tensor multiplication by a single dot, i.e. $\varvec{c} = \varvec{A} \cdot \varvec{b}$ and a double contraction with the two dot operator $c = \varvec{A}:\varvec{B}$. Lower case spatial differential operators indicate differentiation with respect to the current placement. The current placement by means of a Lagrangian description of the solid is defined at different time steps t. Time discrete positions are denote by $(\bullet )^{t}$ for the current time state, $(\bullet )^{t+1}$ or $(\bullet )^{t+1/2}$ for next or intermediate states, where $t+1$ characterizes time $t + \Delta t$. The delta is hence used to denote any increment within a time step, also in multiplicative updates such as $\varvec{F}^{t+1} = \Delta \varvec{F} \cdot \varvec{F}^{t}$.

2 MPM for solid mechanics

One of the greatest appeals in the MPM is the discretization of the equation of motion. The strong form of the balance of momentum is given by $- \text {div}\, \varvec{\sigma } - \rho \, \left( \varvec{b} - \ddot{\varvec{x}} \right) = \varvec{0}$ $\ \forall \, \varvec{x} \in {\mathcal B}$. The boundary of the domain of the solid body $\mathcal B$ is decomposed into $\partial \mathcal B_{u}$ and $\partial \mathcal B_{\sigma }$ with outward-pointing normal vector $\varvec{n}$. On $\partial \mathcal B_{\sigma }$ we describe the surface traction $\bar{\varvec{t}}=\varvec{\sigma }\cdot \varvec{n}$. The symbol $\varvec{\sigma }$ denotes the symmetric Cauchy stress tensor, $\rho$ the current density, $\rho \, \varvec{b}$ a body force, $\varvec{x}$ the position of a material point in the current configuration and $\ddot{\varvec{x}}$ the acceleration. In order to solve the balance of momentum we transform it into a Galerkin weak form with an appropriate test function $\delta \varvec{u}$, i.e.

$$\begin{aligned} G= & {} \int _{\mathcal B} \, \varvec{\sigma } : \text {grad}\, \delta \varvec{u} \, \text {d}v - \int _{\mathcal B} \, \rho \, \varvec{b} \cdot \delta \varvec{u} \, \text {d}v\nonumber \\&+ \int _{\mathcal B} \, \rho \, \ddot{\varvec{x}} \cdot \delta \varvec{u} \, \text {d}v\nonumber \\&- \int _{\partial \mathcal B_{\sigma }} \, \bar{\varvec{t}} \cdot \delta \varvec{u} \, \text {d}a = 0. \end{aligned}$$

(1)

In the following we neglect the boundary traction. As the domain of the solid is in fact discretized not by the computational grid, but by material points, a Riemann integration scheme is introduced. The discrete counterpart of the weak form of balance of momentum appears as

$$\begin{aligned} G \approx \sum _{{\texttt {MP}}= 1}^{\text {No}{\texttt {MP}}} \, \left[ \varvec{\sigma } : \text {grad}\, \delta \varvec{u} \,v - \varvec{b} \cdot \delta \varvec{u} \, \rho \, v + \ddot{\varvec{x}} \cdot \delta \varvec{u} \, \rho \, v \right] _{{\texttt {MP}}} = 0, \end{aligned}$$

(2)

where the total number of material points is denoted by $\text {No}{\texttt {MP}}$ and $v_{{\texttt {MP}}}$ the associated volume to a material point. In the MPM the equation is solved for a discrete-time step on a computational background grid (CBG).

In this contribution, we rely on an explicit solution scheme for equation (2), hence considering the acceleration at the nodes of the computational grid as the primary unknown. Following the Galerkin principle we use the approximations

$$\begin{aligned} \delta \varvec{u}(\varvec{x})= & {} \sum _{I=1}^{\text {No}I} \, N^{I}(\varvec{x}) \,\delta \varvec{u}^{I}, \nonumber \\ \text {grad}\, \delta \varvec{u}\left( \varvec{x}\right)= & {} \sum _{I=1}^{\text {No}I} \, \varvec{B}^{I}(\varvec{x}) \, \delta \varvec{u}^{I},\nonumber \\ \ddot{\varvec{x}}(\varvec{x})= & {} \sum _{I=1}^{\text {No}I} \, N^{I}(\varvec{x}) \, \ddot{\varvec{x}}^{I}, \end{aligned}$$

(3)

where the matrix $\varvec{B}^{I}$ contains the spatial derivatives of the nodal shape functions $N^{I}(\varvec{x})$. Using these, the global system is assembled and solved explicitly for nodal accelerations ${{\varvec{a}}}^{I}$ and velocities ${{\varvec{v}}}^{I}$, using a Velocity-Verlet time discretization. With nodal solutions on the grid for the current time step, these are mapped to the material points, using

$$\begin{aligned} \varvec{a}_{{\texttt {MP}}}(\varvec{x})= & {} \sum _{I=1}^{\text {No}I} \, N^{I}(\varvec{x}) \, \varvec{a}^{I}, \quad \nonumber \\ \varvec{v}_{{\texttt {MP}}}(\varvec{x})= & {} \sum _{I=1}^{\text {No}I} \, N^{I}(\varvec{x}) \, \varvec{v}^{I}, \nonumber \\ \text {grad}\, \varvec{v}_{{\texttt {MP}}}(\varvec{x})= & {} \sum _{I=1}^{\text {No}I} \, \varvec{B}^{I}(\varvec{x}) \, \varvec{v}^{I}, \end{aligned}$$

(4)

and then update the kinematical state of the material points

$$\begin{aligned} \varvec{x}^{t+1}_{{\texttt {MP}}} = \varvec{x}^{t}_{{\texttt {MP}}} + ( \varvec{v}^{t+1/2} )_{{\texttt {MP}}}\, \Delta t, \quad \varvec{v}^{t+1}_{{\texttt {MP}}} = \varvec{v}^{t}_{{\texttt {MP}}} + ( \varvec{a}^{t} )_{{\texttt {MP}}}\, \Delta t. \end{aligned}$$

(5)

With the spatial velocity gradient, we update the particles deformation state

$$\begin{aligned} \varvec{F}^{t+1}_{{\texttt {MP}}} = \Delta \varvec{F}_{{\texttt {MP}}} \cdot \varvec{F}^{t}_{{\texttt {MP}}}, \quad \text {and} \quad \Delta \varvec{F}_{{\texttt {MP}}} = \varvec{I} + ( \text {grad}\, \varvec{v}^{t+1/2})_{{\texttt {MP}}} \, \Delta t, \end{aligned}$$

(6)

which is then argument to the constitutive routines that results in an update of material history and stresses. Performing the stress update as the final action within a time step, the algorithm represents a Update Stress Last (USL) scheme. This treatment shows some amount of numerical damping as shown early by Bardenhagen [28] and is thus not fully energy conserving. This choice is intended however, as dissipation may increase numerical stability, also suggested in Nairn [29]. Some details of the algorithm are given in Appendix, see Table 3.

2.1 MPM contact against rigid surfaces

In the introduction we discuss some aspects of contact mechanics in MPM. In this contribution, we need to model contact between rigid surfaces of a tool i.e. the pressure bars, and a deformable body, modeled by the MPM. Furthermore, the contact conditions need to be modeled with an actual contact law, such as Coulomb. As shown by Hartley et al. [25], friction has an impact to the results in an SHPB.

To account for this requirement an implementation of the model proposed in Ding and Schroeder [30] was chosen. Their approach of a penalty contact is compatible with the utilized explicit scheme and allows for a straight forward implementation. This scheme allows to model rigid bodies as geometric objects with exact contact normals. Also the rigid body does not need to share nodes with the CBG. Rather this contact is resolved directly between the material points and the geometric shape of the rigid body. The resulting penalty force is mapped to the CBG to contribute to the solution. This does however lead to some amount of penetration.

Only single surfaces are considered in this contribution as contact partners. For a detailed insight into the algorithm, we refer to the work of Ding and Schroeder [30], while we give a brief overview of the strategy in to following, respectively.

Based on a hierarchical contact search, we cover the case of a material point actually in contact with a rigid surface. Contact is defined as a penetration of the material points position $\varvec{x}_{MP}$ into the region covered by a plane, with outward-facing normal $\varvec{n}$ and tangent vectors $\varvec{t}$. At first detection of this contact, an attachment point ${{\varvec{x}}}_{A}$ is initialized for this contact using the closes point projection ${{\varvec{x}}}_{C}$. The contact penalty force $\varvec{f}^{t}_{pen} = \kappa \, \left( \varvec{x}_{A}^{t} - \varvec{x}_{{\texttt {MP}}}^{t}\right)$ will be computed from the distance of the particle to the attachment point throughout future time steps, until the contact is resolved. The parameter $\kappa$ represents the penalty parameter (Fig. 1).

Fig. 1

Contact over different time steps. a A particle at position $\varvec{x}_{{\texttt {MP}}}$ penetrates the surface of a body with surface normal $\varvec{n}$ at time t and b particle exceeds the Coulomb penalty force

In case that the trial penalty force violates Coulombs law $\varvec{f}_{pen} \cdot {{\varvec{t}}}\le \mu \, \varvec{f}_{pen} \cdot {{\varvec{n}}}$, the attachment point $\varvec{x}_{A}^{t}$ is relaxed towards the closest point until the friction cone is satisfied. Hereby the friction cone is defined using the friction coefficient parameter $\mu$.

2.2 Convected particle domain interpolation

Since its advent, several modification were developed to the mapping operations in MPM. An overview is given in the introduction. For the proposed scheme, the convected particle domain interpolation (CPDI) is chosen. The idea is that the particle domain is represented by a geometric domain defined by a set of vertices, that represent a brick element. This domain follows the material points motion and deforms according to its deformation state. In CPDI1 the deformation is assumed to be constant over the material point domain. The interested reader is referred, e.g. to Wang et al. [2] where an overview along with example simulation can be found.

The main consequences of this procedure and therefore the motivation of use in this contribution are:

• The possibility to use finite element meshes for the initial placement of vertices for the material point domains allows for a more accurate tracking (and visualization) of the deformation of bodies.

• This approach improves the numerical stability of the MPM especially in tensile, as reported in [11] while also reducing mesh dependency.

• Implementation is straight forward via a modification of the grid shape function, taking into account the extend and deformation of material point domain.

In this contribution a three dimensional, regular computational grid, composed of 8-node hexahedral cells is used. Also in compliance with the CPDI reported in Sadeghirad et al. [11] the same cell type is used as convected particle domain, but not restricted to any regularity. Also, we do not restrict the domains to be parallelepipeds in 3D but treat them as actual tri-linear brick elements. However, in contrast to CPDI2 as proposed in Sadeghirad et al. [31], the domain vertices are updated by the constant deformation gradient, not by the global velocity field. The actual improvement in this ansatz is by using arbitrary meshes that cover the cylindrical geometries in this research better than any mesh of parallelepipeds could possibly do. Furthermore by using a CPDI1 update, we allow for some material discontinuity in contrast to CPDI2.

To compute the actual shape functions used for the mapping in equation (3) and (4) we introduce two sets of standard Lagrangian shape functions $\bar{N}^{I}$ and $Q^{\alpha }_{{\texttt {MP}}}$ which are well known and can be found in many standard finite element text books , see for example [32, 33]. Hereby $\bar{N}^{I}$ is constructed with respect to each node of the computational grid and $Q^{\alpha }_{{\texttt {MP}}}$ are constructed with respect to each vertex $\alpha$ of a particle domain, illustrated in Fig. 2.

The shape functions that are actually used for constructing the global vector/matrix are computed from these via

$$\begin{aligned} \begin{array}{ll} N^{I}(\varvec{x}_{{\texttt {MP}}}) &{} = \dfrac{1}{v_{{\texttt {MP}}}} \, \displaystyle \sum _{\alpha } \, \left( \, \displaystyle \int _{{\mathcal B}_{{\texttt {MP}}}} \, Q^{\alpha }_{{\texttt {MP}}}(\varvec{x}) \, \text {d}\varvec{x} \right) \, \bar{N}^{I}(\varvec{x}_{{\texttt {MP}}}^{\alpha }),\\ \text {grad}\, N^{I}(\varvec{x}_{{\texttt {MP}}}) &{} = \dfrac{1}{v_{{\texttt {MP}}}} \, \displaystyle \sum _{\alpha } \, \left( \, \displaystyle \int _{{\mathcal B}_{{\texttt {MP}}}} \, \text {grad}\, Q^{\alpha }_{{\texttt {MP}}}(\varvec{x}) \, \text {d}\varvec{x} \right) \, \bar{N}^{I}(\varvec{x}_{{\texttt {MP}}}^{\alpha }). \end{array} \end{aligned}$$

(7)

The integral in equation (7) of the local ansatz function $Q^{\alpha }_{{\texttt {MP}}}$ over the material points domain are computed numerically and are material point wise operations. The positions of the vertices $\varvec{X}^{\alpha }$ are updated with the material points deformation gradient as shown in Fig. 2. The gradient operator on the particle domain is computed by means of the chain rule as in standard finite element procedure.

Fig. 2

CPDI1: 8-noded hexahedral particle domain deforms with the particles deformation gradient

3 SHPB experiment and numerical analysis

3.1 Experimental setup

Fig. 3

Experimental setup of the Split-Hopkinson-Pressure-Bar a Exemplary specimen consisting of 42CrMo4 before and after the experiment. b Schematic of experimental setup including all components. c Depiction of the placement of the specimen between incident and transmitter bar with aligned pyrometer

Fig. 4

Measurement data from the strain gauges after initial impulse

In order to investigate deformations at high strain-rates, specimen consisting of the low-alloyed 42CrMo4 steel were used. The samples were austenitized at an austenitizing temperature of $T={850}{^{\circ }}$, quenched in oil and subsequently tempered for two hours at a temperature of $T={660}{^{\circ }}$. Both processes were conducted in an inert gas atmosphere. The specimen, originally cylindric with a diameter $d_0={4}\,{\mathrm{mm}}$ and height $l_0={4}\,{\mathrm{mm}}$, were deformed, see Fig. 3a using a SHPB. The experimental setup, see Fig. 3b consisted of a striker bar, two steel bars (incident and transmitter bar) and a specimen, which is placed between the two bars, see Fig. 3c as shown in Zabel et al. [34]. As described by Follansbee and Frantz [35], the striker bar as well as the incident and transmitter bar need specific length to diameter ratios in order to assume one-dimensional wave propagation. For this investigation, the incident and transmitter bars had a Young’s modulus $E={210}{\hbox { GPa}}$, density $\rho ={7850}{\hbox { kg}\,\hbox {m}^{-3}}$, a length of $l={1.5}\,{\mathrm{m}}$ and diameter of $d={14}\,{\mathrm{mm}}$, assuring sufficient length to diameter ratios. The positioning of specimen was applied manually using a laser beam referencing point. A gas gun was utilized with a pressure of $p={2.5}\,\mathrm{bar}$ to accelerate the striker bar. The specimen was additionally heated to the temperature of $T={400}{^{\circ }}$ using an induction unit, which was verified using a digital infrared pyrometer. To obtain boundary conditions for the conduction of simulations, the strain signals were measured using strain gauges, which were placed centrally on the incident and transmitter bars, see Fig. 3b.

3.2 SHPB experimental analysis

The strain signals were measured using strain gauges placed on the incident and transmitter bars of the setup. In order to calculate appropriate boundary conditions, the signals were analyzed. The experimental setup leads to a considerable amount of plastic deformations resulting in corresponding impulses, see Fig. 4, where the first two impulses of the incident signal correspond to the incident $\varepsilon _I$ and the reflected impulse $\varepsilon _R$, whereas the first impulse of the transmitter signal corresponds to the transmitted impulse $\varepsilon _T$.

As seen in Fig. 5a, suitable sections of the signals were selected for analysis. In this contribution, the signals were processed using a low-pass filter based on Savitzky and Golay [36], which applies a polynomial regression of sixth order for time-domain signals. This ensures the removal of distortions such as noise or dispersion effects while maintaining a close approximation of the measured signal, see Fig. 5b.

Following the methodology of Al-Mousawi et al. [21], stress and strain values were calculated as

$$\begin{aligned} \varepsilon \cong \frac{2\, c_0}{l_0} \int _0^t \varepsilon _R \texttt { and } \sigma \cong E \cdot \left( \frac{A}{A_0} \right) \cdot \varepsilon _T, \end{aligned}$$

(8)

with $c_0 = \sqrt{E/\rho }$ , $A_0$ is the initial specimen cross-sectional area and A correspond to the cross-sectional areas of the pressure bars. The associated force signals were approximated for further evaluation with $F=A_0\cdot \sigma$ compare Fig. 5c.

Fig. 5

Processing of experimental data. a Selected signal sections of measured data. b Filtered signals to reduce distortions. c Computed force signals of the experiment. d Stress-strain-diagram for the analyzed specimen including the calibrated material model

3.3 Numerical model of the SHPB

The numerical model for the simulation of the SHPB experiment is based on an MPM discretization of the specimen only. The surfaces of incident and transmitter bar that are in contact with the specimen are modelled as rigid surfaces. The geometry is illustrated, alongside with the used dimensions in Fig. 6, respectively. Notice that the computational background grid (CBG) is technically not restricted in space.

Fig. 6

Geometry of the idealized numerical model of the SHPB

The model is free of any grid boundary conditions, contact mechanics between the specimen and both plates are active during the whole simulation. The loading conditions for the specimen are introduced by a prescribed movement of the incident surface $s_{i}$. The transmission surface $s_{t}$ does not move at all. The contact forces on each surface are assumed to be the model equivalent of the bar forces of the experiment, as derived in Fig. 5c. A comparison of these forces derived from the experiment and captured from the simulation is chosen as a macroscopic benchmark property. According to Fig. 5c a time frame of ${150}\upmu \hbox {s}$ is chosen. Throughout this contribution the numerical simulations are conducted with a constant time step size of $\Delta t= 5\cdot 10^{-9}\hbox {s}$.

Table 1 Fitted material parameters for 42CrMo4 (h+At) at $400^{\circ }\hbox {C}$ for a strain rate of $\dot{\varepsilon } \approx {2661.9}\,{\hbox {s}^{-1}}$

Fig. 7

SHPB initial CPDI configurations and CPDI centers, representing the material points

The incident surface $s_{i}$ will move at a constant velocity, while we aim to overlap the simulation and experimental time span, used for the material fitting. For the designated time segment, the induced strain rate can be approximated to be constant. From this we calculate the relative velocity of the surfaces $s_{i}$ and $s_{t}$ which is required for an analogue macroscopic deformation. Considering the measured rate of $\dot{\varepsilon }_{IRT} = 2661.9$ $\hbox {s}^{-1}$ and an initial height of the specimen $h={4}\,{\mathrm{mm}}$, we compute the relative velocity to $v_{y}= {10.6476}\,\hbox {m}\,\hbox {s}^{-1}$. In the simulation this velocity is applied constantly to the upper surface $s_{i}$ as depicted in Fig. 6 while surface $s_{t}$ remains fixed.

For the material model, a finite strain J$_2$ plasticity law based on the multiplicative decomposition $\varvec{F} = \varvec{F}^{e} \cdot \varvec{F}^{p}$ is chosen, using a von Mises yield criterion

$$\begin{aligned} \phi= & {} \sqrt{\dfrac{3}{2}\, \text {dev}\, \varvec{\tau }:\text {dev}\, \varvec{\tau }}- y(\alpha )\nonumber \\&\quad \begin{array}{l} \phi \ge 0 \quad \rightarrow \quad \text {elastic behavior}\\ \phi < 0 \quad \rightarrow \quad \text {plastic behavior.} \end{array} \end{aligned}$$

(9)

with the yield stress

$$\begin{aligned} y (\alpha ) = y_{\infty } + \left( y_{0}-y_{\infty }\right) e^{-\eta \, \alpha } + h \, \alpha , \end{aligned}$$

(10)

as a function of the accumulated plastic arc length $\alpha$. The material parameters are given in Table 1.

The Kirchhoff stresses $\varvec{\tau }$ as used for the yield criterion and the first Piola-Kirchhoff stress tensor are computed from a free energy function $\Psi$, by means of the elastic part of the deformation gradient:

$$\begin{aligned} \varvec{\tau } = \partial _{\varvec{F}} \Psi \left( \varvec{F}^{e}\right) \cdot \varvec{F}^{T}, \quad \varvec{P} = \partial _{\varvec{F}} \Psi \left( \varvec{F}^{e}\right) . \end{aligned}$$

(11)

Here we chose a Neo-Hooke type free energy function

$$\begin{aligned} \Psi \left( \varvec{F}^{e}\right) = \dfrac{\lambda }{4}\, \left( J^{2} -1 - \text {ln}\,J^{2} \right) + \dfrac{\mu }{2}\, \left( \text {tr}\left( \varvec{F}^{e\,T}\cdot \varvec{F}^{e}\right) -3 - \text {ln}\,J^{2} \right) , \end{aligned}$$

(12)

where J denotes the determinant of the deformation gradient. In the case of plastic loading the update of the internal variables is computed by means of the associative flow rule. After some algebraic manipulations we obtain the update scheme

$$\begin{aligned} \varvec{F}^{e}-\text {exp}\left( -\Delta \alpha \, \sqrt{\dfrac{3}{2}} \, \varvec{n}\right) \cdot \varvec{F}\cdot \varvec{F}^{p-1}_{n} = \varvec{0} \, , \end{aligned}$$

(13)

where $\varvec{n}$ denotes the associative flow direction $\varvec{n} = \partial _{\varvec{\tau }} \phi$. The plastic evolution equation is time-integrated locally (at the end of each time step for each material point) with an implicit scheme.

The resulting algorithms to implement the plastic constitutive equations are the most computational expensive part of the proposed MPM code. In this contribution we used the automatic code generation tool AceGen to build the constitutive routine. AceGen offers automatic differentiation as well as code optimization, see [37, 38].

The material parameters for 42CrMo4 at the temperature of 400$^{\circ }$C and for the observed strain rate regime $\dot{\varepsilon } \approx 2661.9$, where calibrated according to the calculated stress-strain curve of the experiment within the area after initial yield strength (see Fig. 5d). The calibration was done using a bounded BFGS optimization scheme based on Zhu et al. [39].

There are other material models available, such as the Johnson-Cook model, accounting for temperature-driven softening, strain rate hardening, damage, and other phenomena. Such a more elaborate model will be subject to future investigations.

The specimen is discretized by material points in a CPDI scheme. This scheme allows to utilize finite element meshes as shown in Fig. 7 with parameters given in Table 2. Each finite element from these meshes is used as a particle domain in the proposed scheme. Notice that during the simulation, each domain will deform individually, according to the accumulated deformation gradient at the material point.

Upfront we also introduce an array of CBG that we used for the SHPB simulation. All grids are composed of cuboid cells with constant grid spacing for each spatial direction. The actual spacings are given in Table 2.

Table 2 CPDI parameters for specimen discretization (left) and computational background grids (right) used for the numerical MPM simulation and directional grid spacings

From a computational engineering point of view, it is important to investigate the sensitivity of results with respect to the discretization for any computational model. In order to validate this aspect of the MPM-SHPB problem, various combinations of particle and mesh density where simulated.

In the following we discuss results obtained by using MPM discretization B5 (compare Fig. 7) on a computational grid with a constant grid spacing of $4\times 10^{-4}$m (G4, compare Table 2). This discretization represents a good balance between computational effort and approximation quality as found in the convergency study discussed later. In Fig. 8 we show the force amplitudes of the experimental dataset from the previous chapter (see Fig. 5d) alongside the contact force amplitudes of the numerical model. Hereby $f_{i}$ is recorded at the upper surface $s_{i}$ and $f_{t}$ at the lower surface $s_{t}$, respectively. Note that in the numerical models surface $s_{i}$ is moving and hence introduces the compression wave. For a better comparison, the experimental forces are given with an offset in time to meet the moment of impact on the numerical simulation. The friction coefficient was set to $\mu = 0.1$. A good agreement can be observed for the forces gathered from the experiment and the contact forces of the numerical model. However, some offset in force intensity is visible. Motivated by the post-processing of the raw experimental data we introduce another quantity for comparison. Here we introduce $f_{a}$, as the volume average over the current specimen volume V, of the vertical stress component multiplied with the initial specimen cross sectional area

$$\begin{aligned} f_{a} = A_{0}\, \dfrac{1}{V} \, \int _{\mathcal B} \, \sigma _{yy} \, \text {d}v = A_{0}\, \dfrac{1}{V} \, \sum _{{\texttt {MP}}= 1}^{\text {No}{\texttt {MP}}} \, {\sigma _{yy}}_{{\texttt {MP}}} \, v_{{\texttt {MP}}}. \end{aligned}$$

(14)

It can be seen in Fig. 8 that this average force is in better agreement to the experimental data that the contact forces $f_i$, $f_t$ which show an offset. The reason for the better agreement might be that the averaging introduces smoothing in time as it computes an average over the whole volume of the specimen. Furthermore, the actual response stress of the specimen is reduced to its center, which actually corresponds to the assumptions of the SHPB experiment as a 1-dimensional process. Hence $f_{a}$ includes very similar assumptions as where used for computing the experimental forces $f_I$ and $f_T$.

Fig. 8

SHPB numerical results for B5 on G4 with $\mu = 0.1$. Contact force amplitudes of the rigid surfaces and averaged force amplitude in comparison with the experimental data

Figure 9 shows the deformed configuration (after ${150}\,\mu \hbox {s}$) of the CPDI discretized Specimen (B5) in 3D and in 2D an exemplary comparison to a photograph of specimen from the set of conducted SHPB experiments, which shows the typical deformation. It can be concluded that the numerical deformation mode is physically reasonable. Furthermore, the comparison with a specimen from the experiment shown in Fig. 9b confirms this assessment.

Fig. 9

SHPB numerical results. a Deformed configuration of the specimen in 3D after ${150}\,\upmu \hbox {s}$, clamped between the two rigid surfaces. b 2D visualisation of the same configuration in direct comparison with a deformed specimen from the experimental set. The deformation mode is reasonable and aligns with the observed pattern from experiments

As a last observation we consider the stress distribution, plotted over the deformed configuration (after ${150}\upmu \hbox {s}$) in Fig. 10. The plot shows the von Mises stresses which reflects the expected shear band pattern.

Fig. 10

Plot of von Mises stresses over the deformed configuration after ${150}\upmu \hbox {s}$. Results show stress patterns expected for this experiment

Fig. 11

Sensitivity analysis of SHPB simulation with respect to the friction coefficient. The final system force $f_i$ of $\mu = 0.05$ deviates less than $5\%$ from the force at the highest friction coefficient $\mu = 0.35$

Fig. 12

Convergence study. a Simulation forces using different CPDI meshes for the MPM discretization with time offset to align the moment of contact on constant CBG G4. b Simulation forces on varying CBG for constant MPM discretization B4. The studies show a stronger influence of the CBG that for the MPM discretization but an overall result convergency with a finer discretization in both patterns

3.4 Convergence study

So far we have presented and discussed the results for CPDI discretization scheme B5 on the CBG G4. In the following, we discuss the numerical models sensitivities. The friction behavior between the specimen and the bars is a very interesting aspect of the SHPB experiment. With actual frictional contact, the 1-dimensional approximation of the SHPB is less accurate, which gets worse with higher friction coefficients. However, previous investigations have shown that it is challenging to determine a friction coefficient from the physical experiment, see Hartley et al. [25]. In Fig. 11 the force amplitudes $f_{i}$ and $f_{a}$ of numerical analysis are given with respect to different friction coefficients.

Fig. 13

Failure mode of SHPB simulation using CPDI discretization B6 from different angles (a, b). The modes indicate the failure to be driven by a strong localization effect during the simulation

It can be seen that a higher friction coefficient does lead to higher stresses in the specimen. Nevertheless, the effect is very subtle as the maximum final value of the system force $f_i$ at a friction coefficient of $\mu = 0.05$ and $\mu = 0.35$ deviate from each other in less than $5\%$. This conclusion is of course restricted to the test conducted here, given the material and process parameter like, e.g. the deformation rate. Regarding the discretization it is shown that the method behaves predictable, rendering higher stresses on more sticky contact.

The numerical method must be validated with respect to discretization sensitivity as well. In Fig. 12 the system forces obtained in a convergence study are given. The study spans around the CPDI discretizations depicted in Table 2 and illustrated in Fig. 7 and also with respect to varying computational background meshes given in Table 2. For the analysis of results with respect to the number of material points, the CBG is kept constant at G4. Additionally, the forces in Fig. 12 are shifted in time to compensate for discretization dependent differences with respect to the contact. In Fig. 1 we have shown the contact law to be handled for the material points and not for the CPDI vertices which results in slightly different contact initiations in time, which is compensated in the plot. This offset was only applied to this plot for a better comparison.

In the convergence study for the CBG a constant CPDI variant (B4) was chosen, respectively. The sensitivity of the results with respect to the applied CBG is best visible in the contact force magnitude of the upper surface $f_{i}$. However one can see that the force trends do converge with finer grids. The authors chose grid G4 as the preferred grid for the results discussed above and also for the study on varying MPM discretization for representing the desired tradeoff between computational expense and stability, but also computational accuracy as shown by this analysis. When finding stable parameters for this analysis, the well-known trend could be confirmed that finer meshes do require smaller time steps. However, for all numerical results, a constant time step could be realized with computational background grids used in this contribution.

In Fig. 12a. the deviation of force trends are relatively low, indicating a generally weak sensitivity of the result with respect to the CPDI discretization. The same conclusion could be drawn from analyzing the deformation modes obtained in these simulations. Yet, a finer discretization of the specimen allows for more precise tracking of these deformations as well as the resulting stress distributions.

It is worth putting some attention to the failure mode of the finest MPM discretization B6 which catches a numerical failure in the simulation. The literature state on MPM does indicate a certain relation between material points and grid density to result in stable results. This observation is generally supported by the results of the present analysis. However, with the finest discretization of B6 one would intentionally expect a stable simulation on CBG G4. This gave reason to analyze the deformation mode where the numerical failure occurs, which is illustrated in Fig. 13.

The deformation mode shows the degeneration of the numerically modeled specimen in a concentrated region. The effect of temporal discretization was eliminated by conducting analysis with smaller time step size, for which we do not show the results here, though the same behavior could be observed. The author’s hypothesis for the failure is a strong localization, driven by an unfortunate combination of MPM discretization and CBG. This behavior of the MPM require further investigations in the future.

4 Conclusion and outlook

In this contribution, we have introduced a 3-dimensional numerical simulation setup using the material point method (MPM). The computational framework combines the convected particle domain interpolation (CPDI) discretization of the MPM with a state of the art frictional contact formulation for particles against rigid surfaces. A 42CrMo4 heat treatable steel was used in a Split-Hopkinson-Pressure-Bar (SHPB) experiment. Experimental data were used to simulate the experiment with the MPM framework. Several numerical simulations were conducted to analyze the sensitivities of macroscopic forces, which were also compared with the experimental data.

The numerical model has proven to be capable to approximate physical experiments to a satisfactory degree. The applied numerical framework based on the MPM shows interesting potential as a tool for numerical material analysis. For the SPHB model, the sensitivities are in a reasonable range and we were able to show that viable results can be achieved for inexpensive discretizations even in a 3D simulation. The proposed variant of CPDI1 based on meshes allows for an accurate geometry approximation, also at high deformations. From a computational engineering point of view, this CPDI variant can be concluded to represent an efficient way for the simulation of the SHPB in 3D.

Further research is required to specify stability criteria in MPM regarding the relation of material point density and the computational grid. This is especially of interest as the used CPDI variant is supposed to enhance this stability.

Implementation details

Table 3 MPM algorithm with leapfrog time stepping and a USL scheme

Table 3 gives a brief overview of the MPM algorithm, specifically used in this contribution. Besides the material point mass matrix $\varvec{m}_{{\texttt {MP}}}$ and force vector $\varvec{f}_{{\texttt {MP}}}$, the local momentum vector $\varvec{p}_{{\texttt {MP}}}$ is required in the Velocity-Verlet scheme. Further, global vectors/matrices are build by classical assembly operators . As tribute to performance, we use a lumped mass matrix.