1 Introduction

Nonlocality or a characteristic length scale in modelling of multiphase materials in the continuum theory approach can be associated with the physical interactions of material points within finite distance. The effects of these nonlocal constituents are usually not incorporated in conventional constitutive models for multiphase materials in classical solid mechanics. In the local theory, it is assumed that material points only interact with their immediate neighbours. This limitation of the local continuum theory for modelling multiphase materials is not considered in a nonlocal mechanical model. To overcome the above-mentioned limitations, state-based peridynamics can be applied. The direct precursors of the peridynamics are the non-local descriptions of crystal mechanics, see Kunin [1] and Rogula [2]. The state-based peridynamics is a continuum theory that employs a non-local model to describe material properties. The essence of peridynamics is that spatial integration, rather than spatial differentiation, is used to compute the force state on a material point (Silling et al. [3]). In the state-based peridynamics, the momentum balance equation is represented by an integral formulation instead of the partial differential equations. A peridynamic material model corresponds to a classical material model when the strain energy density of both the classical and the peridynamic material are equal under affine deformations. The constitutive laws can be formulated based on a non-local version of classical non-linear strain measures directly computed from an approximate non-local deformation gradient (Silling et al. [3]). The form of such laws provides a description of deformation near a local point based on a weighted average of the deformation of all the neighbouring bonds. The resulting force vector-state for constitutive correspondence is in the form of the force vector-state and leads to a peridynamic formulation that conserves angular momentum provided that the classical constitutive model used to ensure conservation of angular momentum and is non polar. It happens if the resulting Cauchy stress is symmetric (Silling et al. [3]).

In literature, coupled formulations involving peridynamics have been developed to analyse thermal, chemical, diffusion, and electrical problems as well as for porous flow [4,5,6,7,8,9,10,11,12,13]. As the peridynamics formulation can inherently account for discontinuities in the displacement field, peridynamics has been predominantly used to solve fracture and damage problems [14, 15]. The extent of peridynamics is not restricted to fracture mechanics. The concise categorization of the peridynamics applications and related contributions in the literature, in particular, are presented in the review paper by Javili et al. [16]. There are described applications of gradient plasticity methods, fractional methods, micropolar plasticity methods, combinations of molecular dynamics and finite element method with peridynamics. An alternative coupling of peridynamics and finite element method is discussed in [17, 18]. Particular developments of the fractional methods are given in Sumelka [19].

In paper by Song and Khalili [20], the elastoplastic constitutive model for geomaterials by adapting a representative model for unsaturated soils via the constitutive correspondence principle is formulated. The total stress is adjusted by a constant value for its possible extension to coupled poromechanical problems. With the effective force vector state, a version of an elastoplastic model for geomaterials formulated in the framework of critical state soil mechanics through the constitutive correspondence principle in the state-based peridynamics. The numerical results show that the proposed constitutive model can be used to model the inception and propagation of shear banding with no need for regularization. Furthermore, the numerical results demonstrate that the variable can be used to characterize the thickness of shear banding zone.

Foster et al. [21] incorporated a viscoplasticity model within the state-based peridynamics framework. Sun and Huang [22] used the peridynamics theory to study impact damage in composite materials. The plastic deformation due to high-velocity impact using ordinary state-based peridynamic theory is presented in the paper by Kazemi [23]. In the paper by Wang et al. [24] a non-ordinary state-based peridynamic formulation for thermo-visco-plastic deformation and impact fracture is applied.

A novel continuum-kinematics-inspired approach for peridynamics and revisited peridynamics thermodynamic foundations are presented by Javili et al. [16]. The related approach is bridging the gap between classical continuum thermodynamics and peridynamics. They propose the peridynamics formulation whose underlying concepts are reminiscent of classical continuum mechanics. In particular, they propose to decompose the interaction potentials into three parts corresponding to one-neighbour interactions, two-neighbour interactions and three-neighbour interactions within the horizon. One-neighbour interactions are identical to bond-based interactions in the peridynamics formulation of Piola (dell’Isola et al. [25]) and Silling [26]. Javili et al. [27] derive the balance of linear and angular momentum corresponding to our interaction potentials and identify the fundamental properties of these potentials such that angular momentum balance is a priori fulfilled and propose thermodynamically consistent constitutive laws.

The classical bond-based peridynamic model is limited due to a fixed Poisson’s ratio. To overcome this limitation, an improved bond-based peridynamic model is proposed in the paper by Zheng et al. [28]. Zheng et al. analyse the deformation and crack propagation of microelastic brittle materials emphasising varying Poisson’s ratios. In the proposed model, the bond is subjected to axial and transverse pairwise forces, and the particle rotation angle is added to eliminate the additional bending moment caused by transverse forces, which is a key factor to satisfy the balance of angular momentum exactly. As a result, the bond has not only axial displacement but also has transverse displacement and particle rotation. This extension in the displacement mode overcomes the limitation of the fixed Poisson’s ratio in the classical peridynamic model. Lee et al. [29] develop a numerical scheme to model contact and impact between a peridynamic domain and a non-peridynamic domain such as conventional finite elements and rigid bodies. In this paper, the contact is modelled as a node-to-surface type, and a penalty method is used for transient analyses by the explicit time integration. Therefore, the computational burden of the proposed numerical scheme might be reduced. A quantitative study by Lee et al. [29] of ballistic perforation using the proposed numerical scheme and the simulation results show good agreement with the experimental measurements and the analytical calculation of residual velocities.

The subject of the study presented in the paper is the deformation of open-cell copper foams as an example of cellular materials under dynamic compression. The skeleton structure of the cellular material is considered as virtual material with geometry and topology generated by tomographic images of polyurethane foams, produced by Strȩk [30]. The geometry of the virtual material skeleton was obtained with use of the reconstruction methodology presented in Nowak et al. [31] and applied for metallic foams by Pęcherski et al. [32]. The skeleton material is assumed as the oxygen-free high conductivity OFHC copper and can be rescaled according to requirements. The OFHC copper powder can be applied in additive manufacturing to produce open-cell multifunctional structures, e.g., crush-resistant heat exchangers, heat capacitors, etc. In the presented peridynamic computations the foam skeleton material is described using an elastic-plastic model with isotropic hardening. The dynamic process of compression and crushing during impact is simulated. In Pęcherski et al. [32], the elastic-viscoplastic model with Huber–Mises–Hencky quasi-static yield criterion was used. The experimental data of the strain-rate dependency of yield strength reveal the plateau within the range of strain-rates between \(10^{-4} s^{-1}\) and \(10^{3} s^{-1}\) of rather small influence of strain-rate with the rapid growth of yield limit for higher strain-rates, cf. Fig. 7 in Pęcherski et al. 2017. Such a behaviour is characteristic for \({fcc}\) metals, and in particular for copper. In the paper by Postek et al. [33], the elastic-plastic rate-independent model is assumed to simplify the peridynamic simulations of crushing processes and reduce the computational costs.

The paper is organized as follows. Section 2 presents a problem statement with the presentation of the constitutive law adjusted to peridynamics format. Section 3 presents the example of Taylor impact test of solid circular bar with the assumed viscoplasticity model and material properties, and Sect. 4 describes copper open-cell foam sample analysis. Section 5 provides conclusions.

2 Problem statement

An outline of the peridynamics method is shown below. It is presented state-based peridynamics [3, 34, 35]. The application is narrowed to the Huber–Mises–Hencky plasticity [36]. It is applied to the rate-dependent model [21].

The dynamic equation of equilibrium reads:

$$\begin{aligned} \rho {\ddot{\mathbf{u}}}\left[ {\mathbf{{x}},t} \right] = \nabla \cdot {{\varvec{\sigma }}}^\mathrm{{I}} \left[ {\mathbf{{x}},t} \right] + \mathbf{{b}}\left[ {\mathbf{{x}},t} \right] \end{aligned}$$
(1)

where \(\rho\) is the material density, \(\mathbf {u}\) is the displacement vector, \({\varvec{\sigma }}^\mathrm{{I}}\) is the first Piola–Kirchhoff stress tensor, and \(\mathbf {b}\) is the body force vector. Equation  (1) is determined at the position \(\mathbf {x}\) of the body \(\Omega\) at time t, Fig. 1a.

Fig. 1
figure 1

Peridynamic body; a continuous representation; b discretized form

Full size image

The state-based peridynamic equation of motion reads [3]:

$$\begin{aligned} \rho {\ddot{\mathbf{u}}}[\mathbf{{x}},t] = \int \limits _H {\left\{ {\mathbf{{T}}[\mathbf{{x}},t]\left\langle {\mathbf{{x}}' - \mathbf{{x}}} \right\rangle } \right. } - \mathbf{{T}}[\mathbf{{x}}',t]\left. {\left\langle {\mathbf{{x}} - \mathbf{{x}}'} \right\rangle } \right\} dV_{x'} + \mathbf{{b}}[\mathbf{{x}},t] \end{aligned}$$
(2)

In the equation above, Eq. (2), \(\mathbf {T}\) is the force-vector state, \(\mathbf {x}^{'}\) is the position of a point within the body \(\Omega\) with respect to \(\mathbf {x}\) within the subdomain H of radius h. The brackets 〈 〉, denote the vector on which the state acts.

The discrete form of Eq. (2) reads:

$$\begin{aligned}&\rho {\ddot{\mathbf{u}}}\left[ {\mathbf{{x}}_i ,t} \right] \nonumber \\&\quad = \sum \limits _{j = 1}^k {\left\{ {\mathbf{{T}}\left[ {\mathbf{{x}}_i ,t} \right] } \right. } \left\langle {\mathbf{{x}}_j - \mathbf{{x}}_i } \right\rangle \nonumber \\&\qquad - \mathbf{{T}}\left[ {\mathbf{{x}}_j ,t} \right] \left. {\left\langle {\mathbf{{x}}_i - \mathbf{{x}}_j } \right\rangle } \right\} V_j \nonumber \\&\qquad +\mathbf{{b}}\left[ {\mathbf{{x}}_i ,t} \right] \quad i = 1,2...n \end{aligned}$$
(3)

The integral is replaced with a sum of k points within the subdomain H, and n is the number of points in the body \(\Omega\). The \(V_j\) is a volume associated with point j.

Fig. 2
figure 2

The deformation state; a continuous form; b discrete form

Full size image

Following Fig. 2, the relative position vector between the points in the undeformed configuration of the body \(\Omega\) is denoted by \(\mathbf {{\varvec{\xi }}}\) . It is given in continuous and discrete forms, Fig. 2a, b, respectively:

$$\begin{aligned} {{\varvec{\xi }}} = \mathbf{{x}}' - \mathbf{{x}} , \quad {{\varvec{\xi }}} = \mathbf{{x}}_j - \mathbf{{x}}_i \end{aligned}$$
(4)

The displacement state in its continuous and discrete form is as follows:

$$\begin{aligned}&{{\varvec{\eta }}} = \mathbf{{u}}\left[ {\mathbf{{x}}',t} \right] - \mathbf{{u}}\left[ {\mathbf{{x}},t} \right] , \nonumber \\&{{\varvec{\eta }}} = \mathbf{{u}}\left[ {\mathbf{{x}}_j ,t} \right] - \mathbf{{u}}\left[ {\mathbf{{x}}_i ,t} \right] \end{aligned}$$
(5)

The deformation state reads:

$$\begin{aligned}&\mathbf{{Y}}\left[ {\mathbf{{x}},t} \right] \left\langle {{\varvec{\xi }}} \right\rangle \nonumber \\&\quad = \mathbf{{y}}(\mathbf{{x}}_j ,t) - \mathbf{{y}}(\mathbf{{x}}_i ,t) = {{\varvec{\xi }}} + {{\varvec{\eta }}} \end{aligned}$$
(6)

Further on, we will deal with the finite deformation and Huber–Mises–Hencky type of plasticity. The gradient definition is of the following form:

$$\begin{aligned} \mathbf{{F}} = \mathbf{{I}} + \mathbf{{u}}\nabla _\mathbf{{x}} \end{aligned}$$
(7)

The discrete form of the gradient definition is defined by introducing an influence function \(\omega\) [37].

$$\begin{aligned} \mathbf{{F}}\left[ {\mathbf{{x}}_i ,t} \right] \approx \left[ {\sum \limits _{j = 1}^k {\omega \left( {\left| {{\varvec{\xi }}} \right| } \right) \left( {\mathbf{{Y}}\left\langle {{\varvec{\xi }}} \right\rangle \otimes {{\varvec{\xi }}}} \right) V_j } } \right] \mathbf{{K}}^{-1} \end{aligned}$$
(8)

The \(\mathbf {K}\) is the shape tensor:

$$\begin{aligned}&\mathbf{{K}}\left[ {\mathbf{{x}}_i ,t} \right] \nonumber \\&\quad \approx \sum \limits _{j = 1}^k {\omega \left( {\left| {{\varvec{\xi }}} \right| } \right) } \left( {{{\varvec{\xi }}} \otimes {{\varvec{\xi }}}} \right) V_j \end{aligned}$$
(9)

Since the yield condition depends on deviatoric stress, we need both the stress deviator and the deviatoric strain rate evaluated at the unrotated configuration. The deviatoric strain rate is obtained subtracting the volumetric part of the unrotated rate of deformation tensor from the total rate of deformation \(\mathbf {d}\) as follows:

$$\begin{aligned} {\dot{\mathbf{e}}} = \mathbf{{d}} - \frac{1}{3}tr(\mathbf{{d}})\mathbf{{I}} \end{aligned}$$
(10)

The stress deviator depends on unrotated Cauchy stress tensor \({\varvec{\tau }}\):

$$\begin{aligned} \mathbf{{S}} = {{\varvec{\tau }}} - \frac{1}{3}tr\left( {{\varvec{\tau }}} \right) \mathbf{{I}} \end{aligned}$$
(11)

Our goal is to show the relationship between the peridynamic force state \(\mathbf {T}\) and the first Piola-Kirchhoff stress tensor. The left and right gradient decomposition into rotation and stretch tensor are of the form [38]

$$\begin{aligned} \mathbf{{F}} = \mathbf{{RV}} = \mathbf{{VR}} \end{aligned}$$
(12)

The unrotated rate of deformation is of the form:

$$\begin{aligned} \mathbf{{d}} = \mathbf{{R}}^\mathrm{{T}} \mathbf{{DR}} \end{aligned}$$
(13)

In Eq. (13), \(\mathbf {D}\) is obtained from the symmetric decomposition of the spatial velocity gradient \(\mathbf {L}\) as follows:

$$\begin{aligned}&\mathbf{{L}} = {\dot{\mathbf{FF}}}^{ - 1} , \nonumber \\&\mathbf{{D}} = \frac{1}{2}\left( {\mathbf{{L}} + \mathbf{{L}}^\mathrm{{T}} } \right) \end{aligned}$$
(14)

The Cauchy stress tensor \({\varvec{\tau }}\) is rotated with rotation tensor \(\mathbf {R}\):

$$\begin{aligned} {{\varvec{\tau }}}_{rot} = \mathbf {R} \,{\varvec{\tau }}\, \mathbf{{R}}^\mathrm{{T}} \end{aligned}$$
(15)

The first Piola–Kirchhoff stress is related to the rotated Cauchy stress tensor by the following relationship:

$$\begin{aligned} {{\varvec{\sigma }}}^\mathrm{{I}} = \det (\mathbf{{F}}) \, {{\varvec{\tau }}}_{rot} \, \mathbf{{F}}^{ - \mathrm{{T}}} \end{aligned}$$
(16)

The force state is expressed in terms of the first Piola–Kirchhoff stress tensor as follows, Silling et al [3]:

$$\begin{aligned}&\mathbf{{T}}\left\langle {{\varvec{\xi }}} \right\rangle \nonumber \\&\quad = \omega \left( {\left| {{\varvec{\xi }}} \right| } \right) \, {{\varvec{\sigma }}}^\mathrm{{I}} \, \mathbf{{K}}^{ - 1} {{\varvec{\xi }}} \end{aligned}$$
(17)

The deviatoric strain rate is decomposed into its elastic and plastic parts:

$$\begin{aligned} {\dot{\mathbf{e}}} = {\dot{\mathbf{e}}}^e + {\dot{\mathbf{e}}}^p \end{aligned}$$
(18)

The elastic and plastic parts of the deviatoric strain rate read:

$$\begin{aligned} {\dot{\mathbf{e}}}^e = \frac{{{\dot{\mathbf{S}}}}}{{2\mu }} , \, {\dot{\mathbf{e}}}^p = \dot{\lambda }\mathbf{{Q}} \end{aligned}$$
(19)

where \(\mu\) is the shear modulus. \(\mathbf {Q}\) is the unit tensor normal to the yield surface that comes from the normality condition and \(\dot{\lambda }\) is the plastic multiplier. It is assumed that the plastic potential function is equal to the yield function. The assumption defines the associativity of the plasticity condition.

The viscoplastic constitutive law is given in the form that is formulated in [39, 40]. The flow stress exhibit a strain and strain rate sensitivity and reads:

$$\begin{aligned}&\overline{\sigma }({\bar{\epsilon }}, \dot{\bar{\epsilon }}) \nonumber \\&\quad = \sigma _y \left( 1+\frac{\bar{\epsilon }}{{\bar{\epsilon }}_0}\right) ^{1/n}\left( 1+\frac{\dot{\bar{\epsilon }}}{{\dot{\bar{\epsilon }}}_0}\right) ^{1/m} \end{aligned}$$
(20)

where \(\overline{\sigma } \doteq \sqrt{3/2 ~\mathbf {S} : \mathbf {S}}\) is the equivalent stress with \(\mathbf {S}\) being the deviatoric part of the Cauchy stress tensor, \(\dot{\bar{\epsilon }} \doteq \sqrt{2/3 ~{\dot{\mathbf{e}}}^p : {\dot{\mathbf{e}}}^p }\) is the equivalent plastic strain rate, \(\dot{\bar{\epsilon }}_0\) is the reference equivalent plastic strain rate, \({\bar{\epsilon }} \doteq \int \limits _t {\left\{ ~{\dot{\mathbf{e}}}^p : {\dot{\mathbf{e}}}^p \right\} dt}\) is the equivalent plastic strain and \({\bar{\epsilon }}_0\) is the reference equivalent plastic strain. The \(\sigma _y\) is the initial yield limit. The exponents m and n are the strain rate coefficient and strain hardening coefficient, respectively. In the analysis of the foam, the strain hardening term is neglected.

According to the energy-based equivalence hypothesis presented in Nowak et al. [41], one can conclude that in two kinds of cellular structures: the periodic one, e.g. of fcc type, and the random one corresponding to foams, the same elastic energy density cumulated in both structures leads to the following assessment of the failure strength of real foam under compression

$$\begin{aligned}&\left( \sigma ^C_f\right) _{real} \nonumber \\&\quad = \left( \sigma ^C_f\right) _{fcc}\sqrt{\left( \frac{E_{real}}{E_{fcc}}\right) } \end{aligned}$$
(21)

where the failure strength, \((\sigma ^C_f)\), corresponds in this case to the elastic limit, while \(E_{real}\) and \(E_{fcc}\) denote the Young moduli for real and periodic skeleton, respectively. Therefore, accordingly, the subsequent discussion of numerical examples appears valid for a broader class of cellular materials.

3 Example of Taylor impact test

3.1 Model and material properties

To demonstrate the rate hardening effects, we consider an example of Taylor impact test of solid round bar. The scheme of the structure is given in Fig. 3. The copper cylindrical sample hits a rigid obstacle with an initial velocity.

The length of the bar is 0.0324 m, and the diameter is 0.00032 m. The bar is discretized with 21,000 points. The material parameters of the bar are as follows: Young’s modulus is 110 GPa, Poisson’s ratio is 0.296, \(\dot{\bar{\epsilon }}_0\) = 1.0E−03 1/s, \({\bar{\epsilon }}_0\) = 0.00218, yield limit is 50 MPa, and material density is 8960 kg/m\(^3\). In this case, the rigid obstacle is modelled, assuming fixed zero displacements along the axis of the bar (z-displacement) at the impacting surface of the bar. The horizon is assumed 0.0065 m. The structure in Fig. 1a is presented as a set of calculation points surrounded by spheres of radius calculated assuming the volume of the sphere is equal to the hexahedral volumes that are shown in Fig. 1b, c. The results are presented on the points-spheres representation.

Fig. 3
figure 3

Scheme of the impact of the bar (dimensions are in m); a bar with boundary and initial conditions; b hexahedral grid, distance of the observed point from the impacting edge; c impacting edge view, position of the observed point O

Full size image

We consider four cases of material properties. Case A stands for the elastic ideally plastic material. The next three cases concern the increasing strain rate hardening. For B, C and D, the strain rate hardening is 0.001, 0.01 and 0.1. The strain hardening coefficient is 0.1 for all of the cases. The strain rate hardening values are chosen to demonstrate the effect of high, medium and low value of the coefficient. The values of their range are given in industrial applications; see Peirce et al. [42].

The simulations are performed with an Open Access program Peridigm [43,44,45,46]. The program is massively parallelized with MPI. It uses the set of libraries Trilinos [47] that provides tools for the equation solving and parallelization. The geometry is prepared using Patran system [48]. Then, the neutral file is transformed into Exodus/Genesis ASCII format, [49, 50]. The ASCII file is converted into binary one. The binary file is decomposed into calculation domains before the Peridigm execution. The problem definition is given in an xml file. The results have to be merged into one file using a python script. Then, the results are converted into ASCII file. To visualize the results the GiD program is used [51]. The analysis is done with the explicit integration method. The stable time step is 5.28E−07 s. However, it is employed a much lower fixed time step of 6.609E−08 s to observe the process at the same time instants for different testing parameters. The analysis takes 120 s on a Linux cluster with Intel/Xeon processors E5-2680 v3. It has been applied 2 processors of 24 cores each.

3.2 Numerical results

We compare the shapes of the bar at the end of the process. The displacement fields for all of the cases are shown in Figs. 4 and 5.

Fig. 4
figure 4

Shape of the bar, impact velocity 50 m/s; a Case A; b Case B; c Case C; d Case D

Full size image
Fig. 5
figure 5

Shape of the bar, impact velocity 80 m/s; a Case A; b Case B; c Case C; d Case D

Full size image

The series of pictures in Figs. 4 and 5 starts from the elastic-plastic case with no strain rate hardening and ends up with the highest value of the strain rate hardening for Case D. We note that the shape of the sample strongly depends on the strain rate hardening coefficient. In the case of its highest value, there is no visible mushroom-like foot in contrast to the elastic-plastic case and the lowest value of the coefficient. The series of shapes are qualitatively similar for both impact velocities, namely 50 m/s and 80 m/s. The values of the maximum displacements are collected in Tables 1 and 2. The difference between the case without rate hardening and with a low value of rate hardening for lower impact velocity is 59.6\(\%\) and 50.8\(\%\) for higher impact velocity. The equivalent plastic strain distribution at the end of the process is shown in Figs. 6 and 7. The distributions that are shown in both Figs. 6 and 7 are qualitatively similar. It means they are similar for both impact velocities. However, the picture is different for particular cases of the strain rate hardening coefficient. In the case of elastic-plastic example without strain hardening, the base of the sample shows significant effect of mushrooming, Figs. 6a and 7a. The equivalent plastic strain can be considered as smeared on the surface of of the bar’s deformed base but the equivalent plastic strain is the highest close to the centre of the base. The cases of high and moderate values of the coefficient, namely Figs. 6b, c and 7b, c are qualitatively similar to Case A. They are more similar to Case A when the sample undergoes higher impact velocity when observing the shape of the developed foot of the sample and the equivalent plastic strain distribution. The Case D is qualitatively different from Cases A, B and C. Namely, the foot is practically not developed, and the equivalent plastic strain is high in the middle of the base of the sample and close to the outer surface of the sample, Figs. 6d and 7d.

Fig. 6
figure 6

Equivalent plastic strain distribution, impact velocity 50 m/s; a Case A; b Case B; c Case C; d Case D

Full size image
Fig. 7
figure 7

Equivalent plastic strain distribution, impact velocity 50 m/s; a Case A; b Case B; c Case C; d Case D

Full size image

Further on, we analyse the variation in time of the equivalent plastic strain at the observed point O. The point is located in the middle of the cross-section of the bar at 3.17E−03 m from its base. We note that for both impact velocities, the curves are qualitatively similar, Fig.8. When concerning the elastic-plastic without strain hardening (Case A), the variation is practically linear throughout the process. In the cases with strain hardening, the variations tend to stabilize at certain levels indicated in Tables 1 and 2.

Table 1 Values at the end of the process, impact velocity 50 m/s
Full size table
Table 2 Values at the end of the process, impact velocity, impact velocity 80 m/s
Full size table
Fig. 8
figure 8

Equivalent plastic strain variation in time at the observed point; a impact velocity 50 m/s; b impact velocity 80 m/s

Full size image
Fig. 9
figure 9

Maximum equivalent plastic strain in the sample at the end of the process, the dependence on strain rate hardening coefficient; a impact velocity 50 m/s; b impact velocity 80 m/s

Full size image

The maximum equivalent plastic strain dependence on the strain hardening coefficient is given in Fig. 9. We find that the dependence is not linear. The difference between the high and low values of the equivalent plastic strain for the two impact velocities is much higher for the low values of the coefficient than for the high one. The percentage difference reads 154\(\%\) and 127\(\%\), respectively.

The maximum Huber–Mises–Hencky stress at the end of the process strongly depends on the rate hardening coefficient. The values are collected in Tables 1 and 2. The highest stress is reached when the strain rate hardening coefficient is the highest, namely, in the Case of D, for both impact velocities. The differences between Cases B and D for the impact velocity 50 m/s an 80 m/s are 44\(\%\) and 39\(\%\), respectively. To summarize, the rate hardening coefficient is vital for the behaviour of the sample. Its value affects the shape of sample strongly, maximum equivalent plastic strain and Huber–Mises–Hencky stress.

4 Copper foam sample analysis

4.1 Model and material properties

We analyse the impact of a copper foam sample into an ideally elastic block. The scheme of the analysed system is given in Fig. 10. The foam is represented by tetrahedral volumes while the foam is represented by hexahedrals. It is shown in Fig. 10b, c. The copper foam sample hits the elastic obstacle with an initial velocity. Similarly to the considered example of Taylor impact test of solid bar the numerical model consists of set of points surrounded by spheres. It is shown in Fig. 11. The results are presented on points-spheres representation.

Fig. 10
figure 10

Impacting foam sample, basic volume discretized geometry; a scheme; b detail A; c detail B

Full size image
Fig. 11
figure 11

Impacting foam sample, points-spheres representation; a entire structure; b detail A; c detail B

Full size image

The dimensions of the foam sample are 0.0025\(\times\)0.0025\(\times\)0.0025 m. The plate is of the dimensions 0.0035\(\times\)0.0035\(\times\)0.001 m. The sample hits the plate with initial velocities of 20 m/s and 40 m/s, respectively. The plate is discretized with 1,600,000 points, and the foam skeleton is discretized with 2,651,427 points. The material properties of the foam skeleton are the same as those of the cylindrical bar considered in the Taylor test. The plate is ideally elastic of the material properties as follows: Young’s modulus is 210 GPa, Poisson’s ratio is 0.296 and material density is 8960 kg/m\(^3\). The process is followed up to 0.095E−08 s. The explicit time stepping is used. The stable time increment is 2.8E−09 s. However, the fixed time incrementation is used. The time step is 2.0E−09 s. The plate is fixed at its bottom surface. The general contact conditions are assumed. It means that the searching for contact is between all outer surfaces of the bodies. It includes the self-contact as well. The penalty formulation is employed. The contact force is computed based on a short-range formulation [43, 45, 52, 53]. The contact force reads:

$$\begin{aligned}&{{f_{cont}}}\nonumber \\&\quad = {\left\{ {\begin{array}{*{20}{c}} {{K_{cont}}\left( {{r_{cont}} - d} \right) {V_1}{V_2}}\\ 0 \end{array} \begin{array}{*{20}{c}} {\text{ for } ~\,{d} \le {r_{cont}}}\\ \text{ for } ~\,{d > r_{cont}} \end{array}} \right. } \end{aligned}$$
(22)

In Eq. (22), \(K_{cont} = 18K_p/(\pi h^5)\) where \(K_p\) is the penalty number (spring stiffness) that replaces the bulk modulus in the spring constant [52], \(r_{cont}\) is the contact radius, and d is the distance between the calculation points. The penalty number is assumed 1.0E+12, and the contact radius is 55.0E−06 m. \(V_1\) and \(V_2\) are the volumes of the corresponding points in contact and d is the distance between them. The horizon is assumed 50.0\({\text E}-\)06 m.

In non-local methods, the boundary conditions appear non-local as well. The problem for gradient methods has been discussed in [54, 55]. As in the gradient methods that are described above, the effect of boundary is propagated into the considered domain [3]. The problem is solved by the PALS model (Position Aware Linear Solid) [56] employing the surface correction. However, in the report [57], it has been shown that the effect of the free surface diminishes when employing fine discretization. The discretization of each branch of the foam, Fig. 11, is close to or finer than the case N = 7 (slide 4 in [57]). The effect of the boundary diminishes.

As in the Taylor impact test example, we assume four cases of the material properties. Case A stands for the elastic, ideally plastic material. As in the previous consideration, we take into account three further cases, namely, cases B, C and D of the strain rate hardening 0.001, 0.01 and 0.1, respectively. The strain hardening effect is neglected. As previously, the strain rate hardening is chosen to demonstrate the effect of the low, medium and high value of the coefficient.

In this case, the geometry of the foam is taken from CT (Computed Tomography) scans that are converted into finite element mesh. The CT scans are transformed into a basic tetrahedral mesh with ScanIP+FE [58]. Then, the basic mesh is processed using GiD [51] and GMSH [59] programs. Employing the GiD program, an outer triangular surface mesh is taken, so-called ’skin’. Further, using the GMSH program, a new refined tetrahedral mesh is obtained. The new mesh is converted into the Exodus/Genesis ASCII format as in the Taylor test example. The rest steps are the same as in the latter case. The GiD program is used for the post-processing of the results.

The analysis is carried out on the same cluster as the Taylor test example. It is used 40 processors with 24 cores each. A run requires about 5100 s. The value varies slightly due to the level of engagement of the cluster with other tasks.

4.2 Numerical results

We compare the equivalent plastic strain and HMH stress distributions. The equivalent plastic strain distributions for all of the considered cases are shown in Figs. 14 and 15. The highest equivalent plastic strain appears when the material is assumed elastic-plastic, Figs. 12a and 13a. When the strain rate coefficient increases what is represented by Cases B, C and D, the maximum equivalent plastic strain decreases. The relevant values are given in Tables 3 and 4. From the qualitative analysis of the subsequent results in Figs. 12 and 13, we note that even though the maximum equivalent plastic strain decreases, the extent of the equivalent plastic strain in the structure increases. It happens in the case of both impact velocities.

Fig. 12
figure 12

Equivalent plastic strain (logarithmic scale), impact velocity 20 m/s; a Case A; b Case B; c Case C; d Case D

Full size image
Fig. 13
figure 13

Equivalent plastic strain (logarithmic scale), impact velocity 40 m/s; a Case A; b Case B; c Case C; d Case D

Full size image

Further on, we follow the HMH stress distribution. The maximum values for both impact velocities are collected in Tables 3 and 4. The maximum HMH stress depends on the strain rate coefficient, significantly, reaching 10 times the yield stress. When observing Figs. 14 and 15, we note that the elastic-plastic case (Case A) is very similar to Case B that is the lowest strain rate coefficient, Figs. 14a, b (lower impact velocity) and Figs. 15a, b (higher impact velocity). The picture is different in Cases C and D for both impact velocities, Figs 14c, d and 15c, d. We observe that the extent of the high HMH stress is similar, but in Case C, the region of higher stress (red colour and colours closer to upper part of the scale) that value is closer to the maximum stress is much more extended than in Case D. In Case D, the pictures are dominated by the colours from midscale of the palette. It means that the highest stress is concentrated.

Fig. 14
figure 14

Huber–Mises–Hencky stress distribution, impact velocity 20 m/s; a Case A; b Case B; c Case C; d Case D

Full size image
Fig. 15
figure 15

Huber–Mises–Hencky stress distribution, impact velocity 40 m/s; a Case A; b Case B; c Case C; d Case D

Full size image
Table 3 Values at the end of the process, impact velocity 20 m/s
Full size table
Table 4 Values at the end of the process, impact velocity 40 m/s
Full size table

The points where the highest equivalent plastic strain appears are located at the bottom of the foam sample, Fig. 16. The positions are different in the case of different impact velocities. The locations depend on the complex shape of the structure. In the case of the impact velocity 20 m/s, we identify two such points (K, L), and in the case of the higher impact velocity, one point is identified (M).

Fig. 16
figure 16

Positions of the observed points at the bottom of the sample (dimensions are multiplied by 1.0E–04 m); a impact velocity 20 m/s; b impact velocity 40 m/s

Full size image
Fig. 17
figure 17

Equivalent plastic strain dependence on time at the observed points (in Fig. 16), impact velocity 20 m/s; a Point K; b Point L

Full size image

In the case of the lower impact velocity, the plots of the equivalent plastic strain variation in time are given in Fig. 17. We note that the strain rate hardening coefficient’s influence starts to be significant for its higher values, namely, in Case D. In Fig. 17, Case B is omitted since it is close to the elastic-plastic case. The curves for Cases A and B are practically overlapped. In the case of higher strain hardening coefficient, the increase of the equivalent plastic strain stabilizes that is in contrast to Cases A, B and C, where the equivalent plastic strain grows monotonically. When considering the higher impact velocity, the observations are the same as in the case of the lower impact velocity. However, the values of the maximum equivalent plastic strain are higher, Fig. 18a.

Fig. 18
figure 18

Equivalent plastic strain dependence on time at the observed point M (in Fig. 16), impact velocity 40 m/s (a); Maximum equivalent plastic strain in the sample at the end of the process, the dependence on the strain rate hardening coefficient, impact velocity 20 m/s and 40 m/s (b)

Full size image

The dependence of maximum values of the equivalent plastic strain at the end of the process on the strain hardening coefficient is depicted in Fig. 18b. Approximately, we note a linear decrease of the maximum equivalent plastic strain in the structure for both impact velocities and the strain rate hardening coefficient’s interval.

5 Conclusions

We investigated the influence of strain rate hardening coefficient on the behaviour of a copper foam that is under impact conditions. The analysis precedes an introductory example of the Taylor test impact bar. The observations could be collected in the following points:

  • the strain rate hardening effect is vital for the behaviour of the structure,

  • an increase of the strain hardening coefficient causes a decrease of maximum equivalent plastic strain in the sample, however the extent of plastic zone is larger,

  • the maximum HMH stress increases in the sample when the strain rate hardening coefficient increases,

  • when the strain rate hardening effect is stronger the HMH stress is more concentrated than in the opposite case.

The presented study contains preliminary results of the virtual cellular material analysis with the assumed skeleton of OFHC copper. The investigations were conducted as numerical simulations of quasi-static and dynamic loading processes. The peridynamics computations provide a convenient and versatile method of repeating numerical experiments. It is in accord with the point of view presented by Zohdi and Wriggers [60]. They stress that making predictions of a new materials behaviour by numerical simulations reduces laboratory expense and accelerates the trial and error laboratory development of new materials.

The elaborated virtual cellular material analysis provides particular studies on demand preparing the experimental investigations program for functional cellular copper materials. Such a goal finds industrial confirmation for copper foam has excellent function of conducting heat and electricity. It can be used as the electrode substrate of lithium-ion battery or fuel cell and electric double layer capacitor electrode material. On the other hand, copper foam can also be used as heat sink materials, heat absorption materials, chemical catalyst carriers and electromagnetic shielding materials, as well as damping or deadening materials.