A Stress-State based Peridynamics model for elastio-plastic material modeling

A Stress-State based PD (SSPD) model using a well-known yield criteria is proposed in this paper and tested on the modeling of two dimensional bars under different loading levels as a first step for further development. SSPD is based upon peridynamics (PD) which utilize temporal spatial integro-differential equation of motion and formulates continuum problems in terms of integral equations, which are capable of modeling discontinuities such as cracks. The proposed bond strength not only depends on the bond stretch, but on the current state of all bonds connected to a particle as well. Thus, a stress-based peridynamics model is obtained. The tensile simulation compared to conventional FEM shows promising performance with an error of 5%. Compression simulations, however, need more investigation to include the effect of contact forces.


I INTRODUCTION
Peridynamics (PD) is a recently-developed novel continuum mechanics that strives to unify the distinction between "discrete" and "continuous" media currently present in the classical continuum mechanics theory [1] Since its development in the landmark paper [1], PD was used to model problems where classical continuum mechanics, which rely on partial differential equations, fail to do so. The majority of research in classical continuum mechanics tries to amend the theory by "adding-on" terms or equations, which usually results in a modified theory with limited application for a certain class of problems [2].
One of the most prominent of such additions is eXtended Finite Element Method (XFEM) presented in Moes et al. [3] and Fries and Belytschko [4] as an extension to conventional Finite Element Method (FEM) to capture discontinuities. PD, on the other hand, inherently provides a suitable framework that predicts material failure and discontinuities. A comparison between PD and XFEM can be found in the work of Agwai et al. [5].
Nevertheless, PD has been successfully implemented in various problems investigating a continuum with discontinuities.
The first bond-based PD model was presented in the landmark paper by Silling [1] and further explained in [6]. Yet, bond-based PD model suffer from the limitations imposed on incorporating the Poisson's ratio. The introduction of a state-based PD model by Silling et al. [7] and Silling and Lehoucq [8] presented a method of dealing with such limitation [9,10].
PD was also successfully applied in both material and geometric non-linear analysis. The behavior of elastic, plastic and damaged models based on PD was investigated by Gerstle [11] Most recently, an ordinary state-based peridynamic model was proposed for the analysis of models with geometrical non-linearities [12]. The bond stretch was defined using logarithmic functions suitable for large deformations.
The flexibility and simplicity of PD enabled it to be couploed within a finite element analysis. The work of Kilic and Madenci [13] used PD to model the regions of discontinuities within the displacement field and coupled those with FEM in an attempt to take advantage of both methods. Similar research can be found in [14,15,16]. PD was also implemented in a FEA framework in [17,18,19,20].
Interested readers may consult extensive literature surveys provided by Madenci and Oterkus [21] and Javili [22]. PD is, however, still not a flawless theory. As will be shown in section (II), a continuum is modeled via discrete particles connected via bonds that can be modeled as springs. The derivation of the spring constant relying on the bulk modulus seems to be heuristic [21].
Engineers are used to using the Young's Modulus, also known as the Elastic Modulus, in describing tensile tests and its resulting stress-strain relationships. The research presented in this contribution derives bond forces within a PD state based on a well-established yield criteria, which enables a deeper understanding of how bonds "perform" and simplifies the implementation of advanced analysis via PD such as fatigue.
The research presented in this paper continues with a brief introduction of the general PD approach in section (II), followed by the derivation of the governing equations for PD using different yield criteria in section (III). The numerical procedure of analysis is presented in section (IV) followed by the case studies in section (V). Finally, conclusions and recommendations are drawn in section (VI).

II BOND-BASED PERIDYNAMICS
Peridynamics relies on a temporal spatial integrodifferential equation of motion [1]. Thus, in a PD model, the body being modeled is discritized into a set of particles with differential volumes. The response of a body , shown in Figure 1, to external forces is assumed to depend on the displacement of particles relative to their initial positions in the reference configuration [1].
The internal forces affecting a particle in PD are derived from a number of particles in its vicinity [22]. Thus, the particle cannot interact with other particles beyond a horizon i.e. a particle in the reference configuration can only interact with another particle ′ that lies within a neighborhood of defined as The basis of PD is the integral equation , which relates to the force per unit volume of a particle at time resulting from the interaction of all other particles ∈ [1]. This integral equation is defined as As indicated in eq. (2), a bond force exists between two particles and that is defined by the relative position and the relative displacement , , t where is the displacement field at time .
The PD equation of motion [1] is defined as where , represents the body forces on particle at time .
One of the assumptions that can be used when modeling the bond force in PD is -in its simplest form -a linear spring.

III STRESS-STATE BASED PD (SSPD)
The PD theory formulates continuum problems in terms of integral equations unlike its classical counterparts, which relies heavily on partial differential equations, which -as previously mentioned -face difficulties whenever discontinuities are present in the continuum. One of the most common discontinuities is crack propagation and growth.
Cracks are considered in PD as bond breakage [2], which requires the definition of a limiting bond stretch that, when exceeded, results in a bond failure.
However, the definition of PD bonds as linear -or even non-linear -springs represents a somehow heuristic approach. Thus, the proposed model in this research formulates the bond forces using well-known stress-strain relationships and models the material behavior using established yield criteria.
The proposed bond strength depends not only on the bond stretch, called strain in this research, but on the current state of all bonds connected to the particle as well. Thus, a Stress-State based PD (SSPD) is realized.
The bond strength between two particles in an SSPD, can be calculated as where is the modulus of elasticity (young's modulus), is the strain at yield and is the effective area allocated to the bond. The bond described in eq. (4) is subjected to the strain | | | | | | not exceeding the rupture strain which is a characteristic value in a uni-axial tension test. If the strain between two particles exceeds , the bond is broken. Broken bonds in the presented work cannot be healed.
The implication of considering strain and young's modulus in eq. (4) is that the behavior of materials subjected to tensile stresses are to be incorporated in the SSPD. Such behavior is usually characterized by a linear PD2-figure0.pdf part, followed by a non-linear part. A typical stress-strain behavior for steel is shown in Figure 2. To simplify such behavior, strain hardening is neglected and a bi-linear stress-strain curve shown in Figure 3 is adopted.
The stresses are obtained from the bond force by dividing it by the "tributary area". In the research presented in this contribution, the area is assumed to be where is the lattice spacing and is the length of a side of a regular polygon inscribed in a circle of radius and is calculated as where is the number of bonds for the particle in question.
A particle is usually bound to multiple particles in its vicinity, as shown in Figure 4. The bond forces described in eq. (4) can reach the forces obtained when stresses within the material reach yielding stresses .
Having multiple bonds with stresses that could reach leads to the conclusion that a particle could be subjected to stresses that produce an equivalent stress that is well beyond , even if eq. (4) is limited to .
The equivalent stress resulting from a two dimensional stress state can be calculated to be √ After calculating for a particle, two cases can emerge: Case 1: does not cause the material to yield. Case 2: causes the material to yield.
The question whether causes the material to yield or not is based on the selected yielding criteria. For example, using the yielding criteria according to Tresca [23] would result in a yield-surface for two dimensional stresses as shown in Figure 5.

If
is outside the envelope defined by the yield criterion, a return-mapping should be performed to correct back to the yielding surface, as no stresses higher than the yield stress are acceptable, since strain hardening is not considered in the research presented here. Such return-mapping is shown in Figure 6.

IV NUMERICAL IMPLEMENTATION
The simulation of a material with SSPD starts by generating a reference configuration using a lattice. In the research presented here, an equally spaced lattice is adopted with a spacing defined to be , resulting in particles with a volume of .    After defining the reference configuration, the material is defined by selecting , , and the yield criterion to be considered.
As for integrating the equation of motion in eq. (3), the Verlet Integration scheme is used.
The pseudo-code of SSPD is shown in algorithm (1), whereas the return mapping is shown in algorithm (2).

V CASE STUDIES
Three study cases are presented to investigate the perfor-mance of SSPD for 2D problems and gain first insights. Of interest are the force displacement behavior and failure patterns. The former is important to check that the suggested stress-state peridynamics is able to capture the behavior of the test specimen and compare it with the well-known stressstrain behavior obtained in practice, whereas the later is important to check whether the failure patterns do meet the expectations based on sound engineering judgment.
For all the cases stated below, the horizon was set to be √ , i.e. for an equidistant lattice, a particle is bound to its neighbors from 8 directions. The radius required to calculate the "tributary area" of the bonds was taken to be . . The analysis is performed under a gradually increasing external forces. This is realized using the function [ ] where is the targeted stress level.
A typical stress history is shown in Figure 7.

2D Bar subjected to Tension
To start with, a simple bar, shown in Figure 8, is subjected to tensile forces and analyzed.
PD2-figure0.pdf  The bar dimensions were 60 mm x 21 mm x 3mm and was assumed to be 3 mm. The discretized bar is shown in Figure 9 in which the center of each element is shown as nodes and the bonds are shown as links.

Case 1 -
To test the accuracy of SSPD, a targetted stress level well below the yield stress is applied in eq. (8). Thus, = 250 MPa is selected.
The displacements resulting from the SSPD analysis of the aforementioned bar are shown in Figure 10 and Figure  11 for the x-direction and y-direction respectively. For comparison purposes, the resulting displacements using a conventional finite element analysis software (Autodesk Robot) is shown in Figure 12 and Figure 13, respectively.
A qualitative comparison of the displacement fields shown in Figure 10 through 13 shows that both methods provide displacement fields that are similar.
Furthermore, a quantitative comparison shows that the maximum displacement is 0.70mm and 0.73mm for SSPD and FEM respectively. This indicates an error of 3mm or 4.3%.
Similarly, the maximum displacement is 0.042mm and 0.040mm for SSPD and FEM, indicating an error of 0.002mm or 5%.

Case 2 -
For the material characteristics defined in this research -= 200 GPa, = 0.002 -the yield stress is expected to be 400 MPa. Thus, a stress level of 500 MPa is selected.
Multiple strain gauges are defined along the length of the bar at the centerline of the bar, namely at = [16.5, 31.5, 46.5] mm.
The resulting stress strain diagram measured at the strain gauges is shown in Figure 14. This stress-strain diagram was obtained just before failure, which is evident from the strain of the sensor at x = 46.5mm which approaches a value of 0.025, which was the value set for .
After failure, the stress strain diagram becomes chaotic, as the bar starts vibrating violently. The linear portion of the stress-strain diagram has a slope of 200 GPa, which corresponds to the provided young's modulus. The inelastic part follows a profile corresponding to a force-driven tensile test with no strain hardening.
The bar just after exceeding is shown in Figure 15. It can be seen that failure occurs near model non-discontinuities, i.e. near the boundary conditions and force application particles.

2D Bar subjected to Compression
The same bar is used in a compression test. Here, a targeted stress level of 500 MPa is assumed. Furthermore, identical SSPD parameters is used in this test.
An identical stress strain and failure diagram as shown in Figure 14 and 15 is obtained. This, however, is not compatible with usual stress strain diagrams obtained from the usual compression tests, which indicates an area for future investigation.
The SSPD fails to capture compression correctly mainly due to the absence of contact forces particles used in the PD analysis. In the current proposed SSPD, particles are allowed to get "crushed", i.e. to have a distance of zero between its neighbors. Even more, the distance can actually become negative, which indicates that a particle has been crushed and even pierced through its neighbor. Such behavior is non-physical and requires further investigation. A recommendation that can be given here is to include contact forces within the framework of SSPD [22]. This is going to be the subject of future investigations.

2D Bar with notches subjected to Tension
The bar described previously is notched from both sides and used in a tensile test. The bar is shown in Figure 16. The targeted stress level of 500 MPa is applied.
The resulting stress strain diagram is shown in Figure 17, obtained just before failure. The linear portion of the stressstrain diagram has a slope of 200 GPa, which corresponds to the provided young's modulus. The inelastic part follows a profile corresponding to a force-driven tensile test with no strain hardening.
The bar just after exceeding is shown in Figure 18. It can be seen that failure occurs near model non-discontinuities, i.e. near the crack as well as boundary conditions and force application particles.

CONCLUSION
A newly developed stress-based peridynamics model (SSPD) was presented in the research work of this contribution and applied on three numerical cases:  The first numerical case applies stresses well below the yielding stresses and shows that the proposed SSPD conforms well with FEM results with an error of around $5\%$. Furthermore, it shows the application of stresses well beyond the yielding stresses, resulting in stress-strain diagrams that conform with the expected stress strain diagram when a bilinear stress-strain relation ship is assumed.  The second numerical case shows a bar under compression. The results are identical to those obtained from tension stresses.  The SSPD fails to capture compression correctly mainly due to the absence of contact forces particles used in the PD analysis. A recommendation that can be given here is to include contact forces within the framework of SSPD. This is going to be the subject of future investigations.  The third numerical case show a bar with notches to encourage failure near those notches. The failure does indeed occur near those notches, but still occur near the force application points too.
For future research, it is planned to investigate the inaccuracies with regard to compression simulation and derive ways of including contact forces to prevent the ``crushing'' of the mesh.
Also, implementing a stress-strain diagram that features strain hardening, as well as cyclic load analysis and fatigue.
Furthermore, tests with 3D bar models are planned to be implemented next and compared to the results obtained form 2D bars.