Accepted Manuscript

A peridynamic plastic model based on von Mises criteria with

isotropic, kinematic and mixed hardenings under cyclic loading

Hossein Pashazad , Mahsa Kharazi

PII: S0020-7403(18)34321-2
DOI: https://doi.org/10.1016/j.ijmecsci.2019.03.033
Reference: MS 4840

To appear in: International Journal of Mechanical Sciences

Received date: 29 December 2018

Revised date: 17 March 2019
Accepted date: 24 March 2019

Please cite this article as: Hossein Pashazad , Mahsa Kharazi , A peridynamic plastic model based
on von Mises criteria with isotropic, kinematic and mixed hardenings under cyclic loading, International
Journal of Mechanical Sciences (2019), doi: https://doi.org/10.1016/j.ijmecsci.2019.03.033

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service
to our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please
note that during the production process errors may be discovered which could affect the content, and
all legal disclaimers that apply to the journal pertain.
 ACCEPTED MANUSCRIPT

Highlights
 A Peridynamic plasticity model with various hardenings is proposed for cyclic load.
 Equivalent plastic stretch is used as internal variable for different hardenings.
 The modified kinematic hardening is based on a scalar-valued internal variable.
 Numerical method is proposed to solve the elasto-plastic Peridynamic equations.
 The results of Peridynamic model are in good agreement with corresponding FEM ones.

T
IP
CR
US
AN
M
ED
PT
CE
AC

1
 ACCEPTED MANUSCRIPT

A peridynamic plastic model based on von Mises criteria with

isotropic, kinematic and mixed hardenings under cyclic loading

Hossein Pashazad, Mahsa Kharazi

Mechanical Engineering Department, Sahand University of Technology, Sahand Newtown, 51335-1996 Tabriz,
Iran

T
Abstract

IP
In this study, an elastoplastic model based on ordinary state-based peridynamic theory is presented.
The von Mises yield criteria is used to describe plastic yielding and the equivalent plastic stretch is

CR
utilized as internal variable for the general form of isotropic hardening, the modified form of
kinematic hardening and the mixed of isotropic and kinematic hardening. Like classical approach to
plasticity, a plastic flow rule is proposed based on the yield function. The proposed plastic model is
described in the thermodynamic framework and it is shown that the proposed plastic flow and

US
hardening of the peridynamic model are satisfied the requirements of the second law of
thermodynamics. The proposed plastic model is rate-independent and the quasi-static analysis is
considered. A numerical approach is proposed for the elastoplastic model and Newton Raphson
AN
algorithm is used to solve the nonlinear peridynamic equations. The numerical validation of proposed
peridynamic models are done by comparing the results of 2-D peridynamic models under cyclic
loading with the results of finite element methods based on the classical local continuum mechanic
assumptions. The numerical results show that the proposed PD elastoplastic model can predict plastic
M

yielding and linear/nonlinear hardening of materials beyond initial yield stress accurately. Moreover,
it is observed that the presented method it is capable of modeling the kinematic hardening and
Bauschinger effect in cyclic loading too.
ED

Keywords: Peridynamics, Plasticity, Isotropic hardening, Kinematic hardening, Numerical solution

method, Cyclic loading.
PT
CE
AC


Corresponding author.
Email addresses: [email protected] )M. Kharazi)

2
 ACCEPTED MANUSCRIPT

Nomenclature

, , material parameters
cross-sectional area of a beam in a 1-D peridynamic model
material parameter
hardening thermodynamical force
body force vector of a material point
̂ material parameter
̅

T
equivalent back-stress
a positive factor

IP
volume correction factor

CR
strain energy density correction factor
dilatation correction factor
distance between material points
̂ material parameter
extension state
US
AN
̂ unit vector
linear elastic modulus
internal force vector
M

function of incremental equivalent plastic stretch

function of equivalent plastic stretch
ED

the thickness of a plate in a 2-D peridynamic model

family of a material point
PT

invariant of the stress deviator tensor

isotropic hardening modulus
CE

, Material parameters
tangent stiffness matrix
kinematic hardening modulus
AC

boundary line in a continuum mechanic model

̂ normal unit vector of loading
number of material points located in the family region
ratio of the horizon to the distance between material points
uniform distributed loading
norm of residual vector
⃗ residual vector
fictitious boundary layer

3
 ACCEPTED MANUSCRIPT

stretch between two material points

̅ equivalent plastic stretch
sign of equivalent back-stress
force density vector
scalar part of force density vector
force density state
force vector state
displacement of a material point

T
̈ acceleration of a material point

IP
volume of a material point
strain energy density

CR
dilatational part of the strain energy density
distortional part of the strain energy density

US
coordinate of a material point in undeformed configuration
coordinate of a material point in deformed configuration
deformation vector state
AN
Greek
M

internal variable of hardening

material parameters
ED

horizon
width of fictitious boundary layer
strain
PT

dilatation term of a material point

bulk modulus
CE

small positive perturbation

shear modulus
Poisson ratio
AC

mass density of a material point

effect factor of hardening
stress
̅ equivalent stress of a material point
yield stress
yield function
free energy density

4
 ACCEPTED MANUSCRIPT

subscripts

initial value of a parameter

j a material point
k a material point
value of a parameter in previous load step
trial value of a parameter
, , Cartesian coordinate axes

T
IP
superscripts

́ [prime] parameters of neighbor material points

CR
CM a parameter based on the continuum mechanics
elastic part of a parameter
iteration number
iteration step US
plastic part of a parameter
AN
PD a parameter based on the peridynamics
load step
M

1. Introduction
ED

Elastoplastic model plays a significant role in the structural analysis of ductile material. In these
types of materials, when the stress level of material in loading process reaches the yield stress,
permanent (plastic) deformation occurs in the material. The yield stress may be constant, known as
PT

perfect plasticity, or may change accompanying the evolution of plastic deformation which is called
hardening. In plasticity with hardening, the material resists loading when the stress exerts beyond the
initial yield stress, however this resistance is less than that occurs in elastic domain. Hardening is
CE

defined by a function that depends on the history of plastic deformation and could be a linear or a
nonlinear function. In some types of material models, the plastic deformation doesn’t depend on the
rate of loading and this behavior is called rate-independent plasticity.
AC

The plasticity is one of the most interesting field in the engineering researches and the elastoplastic
theory and its applications were studied from many years ago. The definition of perfect plasticity and
the hardening were presented in [1] and the kinematic hardening models and its mechanical behavior
in plastic region were discussed in [2]. In a general study on the plasticity [3], the different yield
criteria and plasticity flow rule as well as different hardening models were presented. Simo and
Hughes [4] proposed the finite element solution of plasticity and viscoplasticity. In their work, the
numerical solutions for different rate-independent plastic models and viscoplastic models were
presented. An extensive mathematical theory of the plasticity was presented in [5]. In this study, finite
element method and numerical algorithm of different elastoplastic models were proposed.
furthermore, the thermodynamic aspect of plasticity was presented. It should be noted that all of
mentioned plasticity studies and their proposed models are based on the continuum mechanics.
Continuum mechanics assumes that material is distributed continuously (without gap or discontinuity)

5
 ACCEPTED MANUSCRIPT

and remains continuous even after deformation. In mathematical elastoplastic models based on this
assumption, partial differential equations are used and these models are just defined in continuous
media.
When a ductile material is subjected to a high value loading, especially cyclic loading, damage
(deterioration) develops together with plastic deformation [6]. As damage develops, microcavities and
microcracks initiate, propagate and coalesce together to create a crack. When a discontinuity appears
in the material, such as crack or void, the spatial derivatives of constitutive model cannot be defined.
Therefore, continuum elastoplastic models have inherent limitation to solve the problem in which
crack growth or fracture occurs with the evolution of plastic deformation.
Peridynamics (PD) is an innovative theory in solid mechanics introduced by Silling [7]. In this -
theory, instead of partial differential equations used in continuum mechanics, integro-differential
equations are used in developing the governing equations, in which the displacements appear in the

T
integral form and the differential operator defines respect to time. This feature makes PD’s

IP
constitutive relations remain valid encounter with discontinuities. Another important feature of PD is
the nonlocal nature of this framework. To apply the nonlocal effects in modeling the medium, PD
represents a novel length scale notion, called horizon , which plays a significant role in the PD

CR
models. A material particle, which is called material point, interacts with all of material points located
in distance less than from it. The interaction was defined in different forms in peridynamic theory.
In one form, interaction between two material points only depends on the deformation of these

US
material points alone and it is referred as bond-based theory [7]. Another form which is referred as
state-based theory is proposed by Silling et al. [8] and in this form, a material point is affected by the
deformation of all material points within its horizon. State-based PD is divided into ordinary and non-
ordinary types. In ordinary state-based peridynamics, the force vector between two material points is
AN
collinear with their relative position vector but in non-ordinary state-based PD, the force vector might
have different direction from the relative position vector.
Concerning the elastic model, the extensive studies in PD have been introduced in [9] and [10]. In
these studies, a response function was introduced to determine the force vector interacted between
M

material points and the brittle damage models are presented to eliminate this interaction in special
conditions. A coupled peridynamic and thermomechanic model based on thermodynamic
consideration was proposed in [11] to consider the thermal effects in fracture analysis. Some PD
ED

studies has been considered damage initiation and crack propagation in brittle and ductile materials. A
modified bond-based peridynamic model was proposed in [12] to simulate the linear elastic behavior
and predict crack nucleation and propagation in brittle materials. In this study, a continuous function
of distance between material points is used as bond stiffness between them to increase the accuracy of
PT

results. In [13], the bond-based peridynamic theory was utilized to simulate the bifurcation and
coalescence of flaws in the rock-like materials under tensile loads. In this work, the propagation of
micro-flaws and macro-flaws in brittle materials were studied. It should be noted that one of the
CE

limitations of bond-based models is constant poisson’s ratio for elastic materials. For overcoming this
limitation, an extended bond-based model was presented in [14]. In this study, the interactions among
material points depend on not only length of bonds between material points but also the rotation
angles between pair bonds. The proposed model also is employed a bond failure criterion to
AC

investigate crack propagation or coalescence in brittle materials such as rock specimens. The
mentioned extended bond-based PD model referred as bond-pair-based peridynamic model, were used
in various fracture studies. Using a 3-D bond-pair-based model, the fracture of rocks with intermittent
fissures under compressive-shear loading were studied in [15]. In [16], the critical stretch and the
critical shear energy density of a bond were used in a bond-pair-based PD model to simulate different
types of fracture modes in brittle materials. In another field study, a coupled thermo-mechanical bond-
based model was proposed in [17] to simulate crack propagation in brittle materials under thermal
shock. In addition to elastic deformation, heat transfer was considered in the proposed model. In [18],
a thermo-mechanical fracture model was proposed in the peridynamic framework to simulate
thermally-driven crack growth in the thin glass plates.
Non-ordinary state-based model is an alternative approach in the peridynamic framework to study
fracture problems of material. In this model, the force vector between material points is obtained from

6
 ACCEPTED MANUSCRIPT

stress tensor. An extended non-ordinary state-based PD model were employed in brittle fracture
studies in which the maximum tensile stress criterion and the Mohr-Coulomb criterion were utilized
to determine the tensile and the shear failure of the bonds among material points, respectively. Using
the mentioned model, the crack propagation process were simulated in brittle materials under
compressive loads in [19] and under dynamic loads in [20] as well as in pre-cracked Brazilian disks
under compressive line loads in [21]. In another research [22], a non-ordinary state-based PD model
was developed to simulate fracture in quasi-brittle materials. In this study, the modified Johnson
Holmquist damage model was implemented into the peridynamic state-based formulation. Tupek et
al. [23] incorporated classical continuum damage model in peridynamic theory by modifying the
peridynamic influence function. In the proposed model, the specific case of the Johnson–Cook
damage model was used for the analysis of ductile material behavior. In this study, the behavior of
ductile materials was considered and a damage model presented based on PD theory. In [24], the

T
phase field constitutive equations was incorporated in a non-ordinary state-based model. The
proposed approach was employed to simulate damage evolution in solids and delamination in

IP
composites.
Despite the advantages of peridynamic theory, the computational cost of PD simulation is more than
the continuum mechanic model. Therefore, the coupled models of PD and continuum mechanic theory

CR
were developed to utilize the advantages of these two theories. According to this point a coupled
model of continuum damage mechanics and PD was proposed in [25] to simulate a complete material
failure that involves damage nucleation, crack formation and its propagation. In their study, the

US
continuum damage mechanics simulates material degradation and damage nucleation before a
macrocrack formation and after that, the peridynamic modeling is employed to simulate the crack
propagation. In another approach, a coupling technique between PD grids and finite element method
meshes was proposed in [26] to solve 2-D static problems. In this technique, the critical region of the
AN
domain where a crack growth could be happened, are analysed by PD and in other regions the finite
element method are applied. A novel technique was presented in [27] to couple local and non-local
meshless methods. This technique couples the Meshless Exponential Basis Functions method to the
bond-based PD model to solve static problems. In [28], a different approach was employed to
M

decrease the computational cost of PD models. In this study, a coupling of finite point method and PD
meshless method was proposed to simulate crack propagation in brittle material. In the proposed
technique, the region near a discontinuity, is analyzed by bond-based PD model and finite point
method is employed in other regions. Shojaei et al. [29] proposed a novel numerical technique to
ED

reduce computational cost of PD models. In their proposed method, the domain of problem is
analyzed by different grid size in which the grids of critical region are finer than other regions. In their
study, a coupling method between the parts of domain with different grid size were presented to
PT

simulate dynamic crack propagation accurately.

Besides the coupling of PD and classical continuum mechanics, the peridynamic model can be
coupled with the multiscale models such as Molecular dynamics (MD) too. A coupled model of PD
and MD was proposed in [30] in which the material is divided to atomistic region, macroscale region
CE

and transition region that transfers information between two first regions.
Another area, that researchers applied PD theory is plasticity. According to the fact that plastic
deformation depends on loading history, the relation between force density vector and deformation,
AC

which is used in the elastic PD models, should be replaced by relation between the increment form of
force density and deformation. For the first tries, Silling et al. [8] decomposed deformation state to
volumetric and deviatoric parts and defined yield surface to model the plasticity. In this study a plastic
flow rule was proposed for plastic deformation evolution. In another study, a non-ordinary state-based
PD model was presented to simulate the viscoplasticity by Foster et al. [31]. The proposed model is
able to reproduce the experimental results of Taylor impact tests over a wide range of impact
velocities. Mitchell [32] presented a perfect plasticity model according to ordinary state-based PD
theory and described a nonlocal yield criterion which depends on the yield stress and PD horizon. An
ordinary state-based PD model was proposed in [33] to simulate pressure-dependent plasticity. In
their work, the nonlocal PD variables were related to the local parameters of Drucker-Prager
plasticity. In addition, a bond failure criterion was defined to simulate the elimination of the
interactions which would finally predict the fracture of material. This method was used to model the

7
 ACCEPTED MANUSCRIPT

plasticity flow and the fracture of quasi-brittle heterogeneous materials, such as concrete. Sun and
Sundararaghavan [34] presented a model for the simulation of crystal plasticity for first time in PD. In
this work, a quasi-static implementation of PD was used and PD forces between material points were
computed from stress tensors obtained from crystal plasticity. A formulation for the fracture analysis
under thermomechanical loadings using non-ordinary state-based PD, was presented in [35]. In their
work, Johnson-cook model was employed for plastic hardening and thermal softening in damage
analysis. However, it should be noted that their model was restricted to adiabatic problems. Madenci
and Oterkus [36] presented constitutive relations for plastic deformation according to von Mises yield
criteria with linear isotropic hardening. The suggested model was based on PD J-integral to predict
crack growth. In addition, in their work a new method was proposed to impose the nonlocal boundary
condition. In a different work, a non-ordinary state-based PD model was proposed in [37] to analyze
strain localization in geomaterials. In this study, a non-local elastoplastic constitutive model for

T
geomaterials was implemented in peridynamic state-based formulation. In one of the latest works in
this field, a non-ordinary state-based PD plasticity model was presented in [38] to simulate thermo-

IP
viscoplasticity. In this study the micro-inertial term in plastic flow rule was considered which could
affect the plastic deformation and heat generation.
According to the literature review, there are different studies based on PD theory which each of

CR
them applied this theory in diverse problems and applications. The plasticity is one of the most
important fields in which the PD theory is applied. The PD elastoplastic models can be coupled with
damage models to simulate degradation of mechanical properties of ductile materials in which

US
damage develops with plastic deformation. The damage developing leads to initiation and propagation
of a crack and the coupled PD models could simulate all of the failure steps. For developing the
mentioned PD model, the PD elastoplastic models and damage models must be coupled in order to
simulate degradation and fracture in ductile materials. The PD elastoplastic model could prepare a
AN
framework to develop a PD ductile damage model to study damage growth together with the plastic
deformation. Therefore, to study the damage of ductile materials based on PD model, developing a
complete elastoplastic PD model considering different hardenings, accordance with the real
M

elastoplastic behavior of ductile materials, is the essential stage in damage analysis. However, in
plasticity models based on the PD theory, there are few studies which are focused on different
hardening models. As, the damage in metals is often coupled with plasticity, therefore to model the
ED

damage in PD models, it is required to interpret the plasticity in PD formulation accurately which is

consistent with the real plasticity behavior of materials. As regards, the isotropic and kinematic
hardening happens in real elastoplastic behavior of ductile materials, it is necessary to develop a
complete PD model which includes both of isotropic and kinematic and also the mixed of these two
PT

hardenings in formulation. According to the authors’ knowledge, only the isotropic hardening is
considered in the PD plasticity models in previous works. Therefore the main purpose and originality
of this study is developing a PD plasticity model with isotropic, kinematic and mixed hardenings.
CE

In non-ordinary state-based PD plasticity models, the yield criteria and hardening internal
variables are based on the stress or strain tensors. In this study, similar to the approach proposed in
[36] for isotropic hardeing, an ordinary state-based peridynamic is presented in which the yield
criterion and hardening internal variable are based on PD parameters. It should be noted that in
AC

modeling the kinematic hardening, a tensor type internal variable is defined to model the movement
of the yield locus in principal stress space in continuum based plasticity models. The originality of the
presented ordinary state-based peridynamic, is introducing a suitable internal variable based on PD
parameters to formulate the kinematic and mixed hardenings. Furthermore, the proposed PD plasticity
model is considered from the thermodynamic aspect as well.
In this study, ordinary state-based PD constitutive relations are proposed to simulate plastic
deformation with linear/nonlinear isotropic, kinematic and mixed hardening. Similar to classical
approach to plasticity, the von Mises hypothesis is used for yield criteria which means that in the
loading path, when the equivalent stress of each material point reaches the critical value, the plastic
deformation occurs. To this purpose, a scalar valued, equivalent plastic stretch, is used as an internal
variable for both isotropic and kinematic hardening. In this work, a plastic flow rule is considered to

8
 ACCEPTED MANUSCRIPT

determine the rate of plastic deformation as well as a consistency condition is adopted to evaluate the
loading/unloading criterion. In addition, according to the thermodynamic framework the plastic
deformation and hardening dissipation obtained from the proposed model don’t violate the dissipation
inequality.
In this paper, the rate-independent plasticity model is studied and the quasi-static loading is chosen
for simulation and a numerical algorithm is proposed to solve the obtained equations. It should be
noted that in the studied cases, the loading is applied as a cyclic displacement control load and the
obtained results of PD model in this work are compared with those obtained from the Abaqus
software based on continuum mechanic assumption.

T
2. Model description

IP
In this section, a review of peridynamic theory, its features and PD based equations of motion are
presented. Afterward, in the second part of this section, an approach of elastic formulation to obtain
the force density existed between material points is presented.

CR
2.1. Peridynamic theory

US
Peridynamic theory is a nonlocal continuum theory in solid mechanics described in [7]. It is
thought that this theory is a continuum version of molecular dynamics [39]. In this theory, it is
assumed that each material point interacts with other material points located in its neighborhood. In an
undeformed body, material point is identified by its coordinate and has an infinity small
AN
volume, , and mass density, . The material point may experience a displacement, , and
after that its location in deformed configuration is identified by, . All material points within a
specific distance, which is called horizon , of are called the family of point and defined by
M

. As shown in Fig. 1, material points, and , interact with their own families of material points,
and , respectively.
ED
PT
CE
AC

Fig. 1. The uniform distribution of peridynamic material points and the family region of the material points and

9
 ACCEPTED MANUSCRIPT

Two material points and experience displacement and during the deformation,
respectively, and relative position vector between them in the deformed configuration is shown in the
form of . All relative position vectors associated with a material point, are stored in a
deformation vector state, , which is defined in [8]. Deformation vector state plays an analogous role
in peridynamic theory as deformation gradient does in continuum mechanics. Moreover, in
peridynamic theory, stretch between two material points and can be defined as [10]

| | | |
(1)
| |

In the state-based peridynamics, a material point is influenced collectively by the deformation of

T
all material points within its family region and it is interacted with all of them [10]. These interactions

IP
can be described in the form of a force density vector, , which shows the interaction force that
the material point exerts on the material point . All of the force density vectors associated with a
material point, are stored in force vector state, , similar to the deformation vector state. A difference

CR
between force vector state and deformation vector state is that the deformation vector state is
independent, while the force vector state depends on the deformation vector state [10].
Considering the force vector state, the integro-differential equation of motion of state-based
peridynamics is described in [8] as

̈ ∫ [ ]〈 ́ 〉 [ ́ ]〈
US ́〉 ́ (2)
AN
in which , ̈ and represent the family region, the acceleration and the body force vector of a
material point, respectively, and superscript prim shows the parameters of neighbor material points.
M

As presented in [10], the discretized form of equation of motion of state-based peridynamics can be
written as
ED

̈ ∑[ ( )
(3)
( )]
PT

in which and ̈ represent body force vector and acceleration at material point , respectively,
CE

and is the number of material points located in , the family region of material point , and
is the volume of these material points. According to Eq. (2), the spatial derivatives of displacement
doesn’t appear in integrand, therefore the PD equations of motion are valid even where the
discontinuities appear in the medium.
AC

2.2. Strain energy density and force density vector

For an ordinary state-based PD, interaction between two material points can be obtained from a
scalar-valued strain energy density developed in [10] for an isotropic elastic material. Defined strain
energy density, depends on both material properties and the deformation vector state. Peridynamic
strain energy density at a material point by ignoring the thermal effects, was proposed in [10] as

̂ ∑( ) | | (4)

10
 ACCEPTED MANUSCRIPT

where and ̂ in Eq. (4), are material parameters and , which is called dilatation term, can be
identified [10] in the form of:

̂ ∑ (5)

in which

(6)
| | | |

T
In Eq. (5), ̂ is a material parameter. The material parameters ̂ , and ̂ are determined in the form of

IP
engineering parameters in [10] as below:

̂ ̂

CR
, and (for 1-D) (7a)

, ̂ and ̂ (for 2-D) (7b)

, ̂ and ̂ US (for 3-D) (7c)

AN
in which , , and represent bulk modulus, shear modulus, linear elastic modulus and Poisson
ratio, respectively. The parameters and used in Eqs. (7a) and (7b), represent the cross-sectional
area of a beam in a 1-D model and the thickness of a plate in a 2-D model, respectively. The
parameters ̂ and ̂ depend on horizon size and they are determined for a material point which is
M

located far from the boundaries of the domain and its horizon is completely embedded in the domain.
For other material points whose horizon are not completely embedded in the domain and their family
ED

members are not complete, a correction factor is required. The correction factors are achieved by
comparing strain energy density and dilatation obtained from PD theory to their corresponding values
obtained from classical continuum theory in the same conditions. For a loading condition such as an
applied displacement , the strain energy density correction factor for material point can be
PT

obtained [10] as
CE

(8)

in which, represents the strain energy density computed in the continuum mechanics and
AC

is the strain energy density obtained from Eq. (4). The dilatation correction factor
for material point due to the applied displacement , can be calculated [10] as

(9)

in which, represents the dilatation based on the continuum mechanics and is the
dilatation computed from Eq. (5). These correction factors are employed to calculate the strain energy
density and the dilatation of material points whose family members are not complete in the domain.
The general calculation of these correction factors due to the applied displacement loading in different

11
 ACCEPTED MANUSCRIPT

directions and the use of them in the computing of interaction between material points are explained
in [10].
The strain energy density is divided to dilatational and distortional parts and therefore, by
partitioning the strain energy density of Eq. (4), the dilatational and distortional parts can be shown as
below:

(10a)

and

̂ ∑( ) | | (10b)

T
IP
where and represent dilatational and distortional parts of the strain energy density. In
addition, and are material parameters obtained from dividing parameter . Considering Eqs.

CR
(10a) and (10b), the material parameters and determine the contributions of dilatation ( )
and distortional ( ) strain energy densities which are derived from dilatation ( ).
In the continuum mechanics, dilatational ( ) and distortional ( ) parts of the strain energy

US
density, by ignoring the thermal effects, were presented [40] as:

(11a)
AN
and
M

*( ) ( ) ( )+ (11b)

In 3-D peridynamic models, the contribution of dilatational part ( ) of strain energy density is
ED

same as its contribution in the continuum mechanics. Therefore, in 3-D formulation the parameters
and are obtained as
PT

, (for 3-D) (12)

CE

The general form of dilatation ( ) in the continuum mechanics is the sum of normal strains in 3
directions. With assumption of plane stress, out of plane normal strain, , depends on the in-plane
normal strains, and . Therefore, the and dilatation ( ) are obtained as below [40]:
AC

(13)

and

(14)

Considering Eq. (14), the dilatational part of strain energy density for the plane stress formulation
in the continuum mechanic framework is calculated in the form of

12
 ACCEPTED MANUSCRIPT

(15)

In 2-D peridynamic models, two displacement degrees of freedom are defined and the
displacement in out of plane direction is ignored. Therefore, the dilatation obtained from 2-D PD
model, , is equaled the sum of in-plane normal strain obtained from continuum mechanics as

(16)

Considering Eq. (16) and the equality of dilatational part of strain energy obtained from PD and

T
continuum mechanics in the same conditions, the parameters and for 2-D PD model are
obtained in the form of

IP
, (for 2-D) (17)

CR
It should be noted that as and for 1-D PD model, therefore and in 1-D
peridynamic model are equal to zero.
Force density vector
US
which is described in previous section, can be derived from strain
energy density. The relation between the force density vector and the strain energy density is
expressed as below [10]:
AN
* + (18)
(| |) | |
M

in which, the terms in bracket show the direction of force density vector. Therefore, by substituting
Eq. (4) into Eq. (18), the scaler part of force density vector, , can be obtained as
ED

̂ ̂ (19)
| |
PT

3. Theory of plasticity in peridynamics

CE

According to plasticity definition in the classical continuum mechanics, if the stress level of a
material, which is considered under loading/unloading condition, exceeds the yield stress, permanent
(or plastic) deformation is sustained in material when completely unloaded. Unlike elastic
AC

deformation, plastic deformation depends on the history of loading. Therefore, unlike the elastic
deformation, the incremental form of the force density and the stretch relations are required in plastic
deformation. In this study, it is assumed that the plastic deformation does not depend on the rate of
loading and the applied formulation can predict rate-independent elastoplastic behavior of materials.
Same as the classical formulation of plasticity in continuum mechanics, a description for
equivalent stress is required in plasticity analysis by peridynamics. This is because in plastic analysis,
the elastic region is restricted by the yield surface. Furthermore, the growth and the movement of the
yield surface should be determined in every load step, therefore a description for equivalent plastic
stretch is required as a hardening internal variable.

3.1. Incremental form of PD parameters

13
 ACCEPTED MANUSCRIPT

As mentioned above, the incremental forms of parameters are required for PD plasticity
formulation. It is noted that each applied load step leads to incremental stretch between
material points and , and also incremental force density . Similar to the classical plasticity
theory, the incremental stretch can be decomposed into elastic and plastic parts as

(20)

in which and are the elastic and plastic parts of incremental stretch, respectively.
Similar to the type plasticity models in classical continuum mechanics [5], it is assumed that plastic
deformation doesn’t depend on volumetric deformation in presented peridynamic theory and therefore
dilatation hasn’t any contribution in plastic stretch expansion, .

T
Considering the incremental force density and also by decomposing the incremental stretch, the

IP
incremental form of strain energy can be expressed as below:

(21)

CR
where and , are the incremental form of elastic strain energy density and plastic strain
energy density, respectively and their incremental forms can be written as [10]

∑( US ) (22a)
AN
and
M

∑( ) (22b)
ED

3.2. Yield criterion and equivalent stress

PT

According to the classical approach to plasticity, the elastic domain is restricted by the yield
surface. The yield function, , is defined as

(̅ ) |̅ |
CE

(23)

in which ̅ is the equivalent stress of material point and is the yield stress. In Eq. (23), plastic
deformation occurs when yield function becomes to zero, where isn’t zero, and in other
AC

conditions, elastic deformation occurs in the material. Therefore, the elastic and plastic
loading/unloading conditions can be written as following:

(̅ ) and elastic loading (24a)

(̅ ) and elastic unloading (24b)

(̅ ) and plastic loading (24c)

According to von Mises model, plastic deformation begins when the distortional part of strain
energy density reaches a critical value [5]. This critical value equals the distortional strain energy of

14
 ACCEPTED MANUSCRIPT

material when the stress reaches to initial yield stress , in uniaxial tension test. Therefore, the
initial yield stress can be written as [3]

√ √ (25)

in which is the second invariant of the stress deviator tensor. In order to utilize this equation in
general loading cases, it is required to define the equivalent stress. The equivalent stress can be
defined as

̅ √ ( ) (26)

T
IP
By substituting Eq. (10b) into Eq. (26), the equivalent stress can be written as following:

CR
̅ √ ( ̂ ∑( ) | | ) (27)

3.3. Plastic flow rule

US
Similar to the definition of incremental strain energy density, the incremental form of the yield
function can be written in the terms of incremental force density as
AN

∑( ) (28)
M

where in loading path is non-positive ( ) and the case of in general shows

the neutral loading. However, in the perfect plasticity (without hardening), shows loading
ED

beyond elastic region. According to classical approach to plasticity [1,36], the plastic part of the
incremental strain energy density due to the development of plastic deformation is non-negative,
. In the case of plasticity models with hardening, plastic strain energy density is positive,
PT

and in perfect plasticity, it is zero, . In the material with the assumption of

perfect plasticity, the incremental form of and can be rewritten as
CE

∑( ) (29)
AC

and

∑( ) (30)

By considering the Eqs. (29) and (30) and vanishing the incremental force density, , the following
equation is obtained as

(31)

15
 ACCEPTED MANUSCRIPT

in which is a positive factor. By substituting Eqs. (23) and (26), in Eq. (31), the incremental form
of plastic stretch can be obtained. As derived in [36], incremental plastic stretch can be expressed as

̂
( | | ( ) ) (32)
̂ ̂

and by substituting force density, , which is expressed in Eq. (19), the incremental plastic
stretch can be rewritten as

| | (33)

T
Eq. (33) plays the role of plastic flow rule in the peridynamic theory. Therefore, when the material
point experiences the plastic yielding, by considering the total plastic stretch could be

IP
obtained as

CR
(34)

in which, the and represent the plastic stretch at the end of the last load step and the
previous load step, respectively.
US
AN
3.4. Equivalent plastic stretch and hardening law

As mentioned previously, hardening is a phenomenon that shows the dependence of yield surface
on plastic flow rule. Therefore, internal variables are required to represent the growth and movement
M

of the yield surface. In the proposed hardening model, equivalent plastic stretch is chosen for the
hardening internal variable and yield surface is considered as a function of equivalent plastic stretch.
Madenci and Oterkus [36] defined the incremental form of equivalent plastic stretch as
ED

̅ √ ( ) (35)
PT

in which, ( ) is the distortional part of strain energy density obtained from the
incremental plastic stretch and is a material parameter. Considering Eqs. (10b) and (35),
CE

the incremental form of equivalent plastic stretch can be calculated as

(̂ ∑ (
AC

̅ ) | | ) (36)

In the case of the simple uniaxial tension, without dilatation term, the distortional strain energy
density due to incremental plastic stretch can be determined and accordingly, parameter could be
extracted. Parameter have been determined in [36] for different body dimensions as

(for 1-D) (37a)

√ ̂

16
 ACCEPTED MANUSCRIPT

√
(for 2-D) (37b)
√ ̂

√
(for 3-D) (37c)
√ ̂

Like plastic stretch, total equivalent plastic stretch can be obtained as

̅ ̅ ̅ (38)

T
in which, the ̅ and ̅ represent the equivalent plastic stretch at the last load step and the

IP
previous load step, respectively.
In Eq. (23), in perfect plasticity (without hardening) is constant during all loadings/unloadings.
However, for isotropic hardening, new yield surface is a function of equivalent plastic stretch, ̅

CR
and the yield function can be rewritten as

|̅ | ( ̅ ) (39a)

where US
AN
( ̅ ) ( ̅ ) (39b)

in which is a function of equivalent plastic stretch. For nonlinear isotropic hardening, can be
defined a nonlinear function of equivalent plastic stretch. However, in the case of linear isotropic
M

hardening, is defined as ̅ , in which represents isotropic hardening modulus which is

constant during loading/unloading (independent from plastic stretch).
ED

The experimental studies in plasticity field show that many materials after loading (hardening) in
one direction, show decrease in resistance to plastic yielding in opposite direction. This phenomenon
is called as Bauschinger effect and in classical continuum mechanics, it is modeled by description of
kinematic hardening [2]. When kinematic hardening takes place in loading path, the behavior of
PT

material in unloading/reloading is not isotropic anymore. In kinematic hardening modeling with

classical approach to plasticity, a symmetric stress tensor, which is known as back-stress, is defined.
In contrast to isotropic hardening which requires the scalar valued parameter (equivalent plastic
CE

stretch), the kinematic hardening requires the tensor-valued parameter (back-stress) to identify the
changes in the yielding surface location in principle stress space. The most challenging point in
modeling the kinematic hardening coincide with peridynamics is that, unlike the classical continuum
mechanics, the back-stress concept in the tensorial form is not suitable to utilize in PD. Therefore, an
AC

equivalent form of back-stress which could be applicable in PD theory should be defined. The
modified yield function can be written as

|̅ ̅ | (40)

in which ̅ is equivalent back-stress at the material point . Note that the sign of von Mises stress
̅ , used in the yield function of kinematic hardening model, Eq. (40), must be determined. The von
Mises stress is positive under tension loading and it is negative in compressive loading. In Eq. (40),
the equivalent back-stress, ̅ , can be defined as the function of equivalent plastic stretch and for
general kinematic hardening, ̅ can be obtained as following:

17
 ACCEPTED MANUSCRIPT

̅ ( ̅ ) ( ̅ ) (41)

in which is a function of incremental equivalent plastic stretch. Similar to isotropic hardening, for
nonlinear kinematic hardening, is a nonlinear function of incremental equivalent plastic stretch and
for linear kinematic hardening, is defined as ̅ , in which is kinematic hardening
modulus. In addition, in Eq. (41), represents sign of equivalent back-stress
obtained from loading condition of the problem and it is negative under compressive loading and
positive in tension loading. The total equivalent back-stress at the last load-step, ̅ , can be obtained
as

T
̅ ̅ ̅ (42)

IP
in which ̅ represents equivalent back-stress at the previous load-step.
As mentioned for monotonic loading under uniaxial stress in [5], the stress-strain curve of material

CR
in kinematic hardening model with initial yield stress and initial hardening ̅ , could
coincide with the stress-strain curve of isotropic hardening model with the same initial yield stress
and initial equivalent plastic stretch ̅ . Therefore, the yield functions of these two models is

US
equal with mentioned condition in monotonic tension loading. The equality of these two yield
functions can be expressed as below:
AN
|̅ ̅ | |̅ | ( ̅ ) (43)

For monotonic uniaxial tension, Eq. (43) can be rewritten as

M

(̅ ̅ ) (̅ ) ( ̅ ) (44)
ED

By Simplifying the Eq. (44), the relation between equivalent back-stress and isotropic hardening term
can be found as

̅ ̅ (45)
PT

Considering Eqs. (41) and (45) and ̅ (for linear kinematic hardening), it is revealed that the
kinematic hardening modulus equals to isotropic hardening modulus, .
CE

In addition, by considering these two hardening models, the mixed hardening formulation can be
defined in the general form as

̅
AC

|̅ ( ̅ )| ( ̅ ) (46)

The yield function of linear mixed hardening model can be written as below [4,41]:

|̅ ̅ | ̅ (47)

in which is the effect factor of each hardening and its value can be varied between 0 and 1. When
is zero, material hardened just with isotropic hardening and in this case, the initial yield surface
expands uniformly [5]. When is 1, material deforms with kinematic hardening in plastic region and
the yield surface translates in the stress space [5]. A real-life material under uniaxial cyclic loading
could be hardened with two mechanisms [4,41], therefore the effect factor of this kind of materials

18
 ACCEPTED MANUSCRIPT

should be chosen between 0 and 1 and this factor for a real-life materials can be obtained from
uniaxial tests.

3.5. Thermodynamic consideration

In this section, it is shown that the proposed PD plasticity model satisfies the requirements of the
second law of thermodynamics. For this purpose, it should be demonstrated that the plastic dissipation
of the model is positive and the inequality of Classius-Duhem is not violated. As a whole, the
thermodynamic restrictions on constitutive PD models were studied in [9] and particularly for
plasticity models (without hardening), were studied in [32] and [33]. The 2nd law of thermodynamics
in absence of the thermal effects, which was proposed in [9], can be written in the form of

T
̇ ̇ (48)

IP
in which first term, ̇ , is defined absorbed power density in PD theory [9] and is free energy
density. By defining the extension state as | | | |, the inequality, Eq. (48), can be rewritten in

CR
the form of

̇ ̇ (49)

extension state , plastic extension state

US
in which ̇ is the rate of extension state. It can be assumed that the free energy density is a function of
and the internal variable of hardening , and it can be
AN
divided to two parts as below:

( ) ( ) (50)
M

in which and are elastic and hardening contributions of free energy density, respectively. By
ignoring thermal part, the elastic portion of free energy density can be equal to the strain energy
density [33] therefore,
ED

( ) (51)

By substituting Eqs. (50) and (51) into Eq. (49), the inequality of Classius-Duhem is rewritten in the
PT

form of

( ) ̇ ̇ ̇ (52)
CE

in which the first term, in parenthesis, implies the Eq. (18) which gives the force density as
AC

(53)
| |

In Eq. (52), term ⁄ , is the thermodynamical force associated with plastic extension which can
be obtained as following:

(54)

By considering Eq. (54) and the plastic flow rule, Eq. (31), the second term in Eq. (52) can be
rewritten as following:

19
 ACCEPTED MANUSCRIPT

| |
̇ (55)

in which is a positive coefficient. In Eq. (55), if is considered as a non-negative convex function

and zero-valued at the origin, ⁄ would be always positive [5]. It should be noted that the
proposed yield function possesses the mentioned features and therefore, the plasticity portion of
dissipation is satisfied the inequality.
The term ⁄ in Eq. (52) can be defined as the hardening thermodynamical force associated
with internal variable [5] therefore the third term of Eq. (52) could be rewritten as

T
̇ ̇̅ (56)

IP
in which and ̇ ̅ are the hardening thermodynamical force and equivalent plastic stretch rate,

CR
respectively. According to classical approach to plasticity, if hardening law obtained from a flow
potential which is a non-negative convex function and zero-valued at the origin, the hardening
dissipation is always positive and the inequality of Classius-Duhem is not violated.

4. Boundary conditions
US
AN
In PD models, boundary conditions, such as applied tractions and displacements, may be exerted
on the domain by describing a non-zero volume layer along original boundary of the region. In a
numerical experience, the width of the layer is suggested to be equal with horizon in [42] to ensure
that the boundary conditions are correctly exerted on the actual region in the nonlocal PD model. As
M

shown in Fig. 2, a fictitious layer, , with width , is extended along the boundary of real region, ,
in 2-D model.
ED
PT
CE
AC

Fig. 2. Extending a fictitious layer with width , along the boundary of real region , in 2-D model

20
 ACCEPTED MANUSCRIPT

The traction boundary conditions cannot be imposed directly in PD equations of motion, Eq. (2),
therefore, describing of the traction condition is implemented by exerting of a specific body force
density to all material points of the fictitious layer. The uniform distributed loading in each
timestep, , which is exerted on a boundary in a continuum mechanic model, can be imposed by
defining the body force density in the fictitious boundary layer, , as [10]

̂ (57)

in which ̂ is the normal unit vector of loading and is the width of boundary layer equaled with
horizon. The subscript in Eq. (57) represents the coordinates of material points located in the

T
fictitious boundary layer, .
Similar to the traction boundary conditions, the displacement boundary conditions can be

IP
described in a non-zero volume fictitious layer along the real region boundary with specific
displacement. The displacement condition , which is exerted on a boundary in a continuum

CR
mechanic model, is imposed on material points located in the fictitious boundary layer with an
extrapolation of . This is because the coordinate of the material points within the layer are different
from the coordinate of boundary where displacement condition is exerted in continuum mechanic
model. Therefore, the displacement boundary conditions of the fictitious layer at each timestep can be
written as

( )
US (58)
AN
in which is the coordinate of the material points located in the fictitious boundary layer, , and
is the coordinate of real boundary in the continuum mechanic model. It should be noted that, the
coordinate chosen in both of real boundary and fictitious boundary layer is coincided with the
M

direction of displacement boundary condition . As shown in Fig. 3, in a 2-D PD model with

horizon , in which d is the length of square material points, the fictitious layer, , is extended
in top boundary of the basic region and is the top boundary line where displacement boundary
ED

condition is applied on it in continuum mechanic model. In this PD model, if exerted displacement

conditions in the fictitious boundary layer are in y direction, the components of and in Eq.
(58) will be the component of y direction.
PT
CE
AC

21
 ACCEPTED MANUSCRIPT

T
IP
CR
US
Fig. 3. A 2-D fictitious boundary layer with width and uniform distribution of material points
AN
5. Numerical solution method

As mentioned in PD theory, the equations of motion, Eq. (2), are integro-differential type and in
the discretized form, the spatial integral are replaced with finite summations. For solving these types
M

of equations, the domain is discretized to the material points and for each of these material points, a
specific volume is dedicated. In quasi-static analysis, mass inertia effects are negligible and the
acceleration term in Eq. (3) is set to zero. Since PD equations of motion for quasi-static analysis can
ED

be written as:

∑[ ] (59)
PT

in which is the volume of each material point , located in the family region of material point . It
CE

should be noted that the volume of these material points, , may be not embedded completely in the
horizon of material point . Therefore, a correction factor is required to adjust and remove the extra
parts of volumes placed out of the horizon. Different volume correction factors are presented by
AC

Seleson [43] which were compared together.

For a uniform distribution of material points in 2-D plate, each material point is a cubic with the
volume of , in which represents distance between material points and is thickness of
material points which is equal to thickness of the plate. As shown in Fig. 4, the entire volume of
material point located in the range of | | , are not embedded completely in
the horizon of material point . According the proposed method in [43], the volume of these material
points could be corrected with a linear approximation as

( | |) (60)

22
 ACCEPTED MANUSCRIPT

T
IP
CR
Fig. 4. Material point US
and its family members located in different distances
AN
It should be noted that the volume correction factor for the family members of material point
which is located in the range | | , is set to 1. Considering the volume correction
factor, Eq. (59) can be rewritten as:
M

∑[ ] (61)
ED

A quasi-static analysis consists of several load steps and in each load step, the displacement or
loading boundary conditions should be updated. Updating the boundary conditions makes the static
equilibrium to be violated and therefore, Eq. (61) for the material points of the domain are not
PT

satisfied. A vector, which is called residual vector ⃗ , is defined and the size of vector is set to the total
number of degrees of freedom of all material points. In unequilibrated condition, the residual of Eq.
(61) for each material point is stored in the corresponding element of residual vector. Note that in
CE

residual vector, the elements related to displacement boundary conditions are set to zero. In an ideal
equilibrated system, all the elements of this vector should be equal to zero, but in numerical approach,
residual vector must be minimized to approach to zero vector which makes the equilibrium condition
AC

to be established. For this purpose, the norm of vector ⃗ , must be less than a specific value to enforce
the static equilibrium to be stablished. This has significant effect on convergence rate of nonlinear
solution and in this study, this tolerance value is considered .
In linear elastic peridynamic model, however the internal forces are linearly related to stretch,
nonlinear relation exists between the displacements of material points and the internal forces.
Therefore, a nonlinear solution procedure is required to solve the peridynamic equations of motion
[44]. In this study, the Newton Raphson method is utilized to solve these nonlinear equations. In order
to implement Newton Raphson algorithm, it is required that the tangent stiffness matrix be defined.
The tangent stiffness matrix is the stiffness of a domain due to small displacement exerted on the
current configuration of the domain. The tangent stiffness matrix can be defined in indicial notation as

23
 ACCEPTED MANUSCRIPT

(62)

in which, represents the component i of the internal force vector related to the displacement
vector , and is the component j of the displacement vector. First-order derivation, which is
imposed on Eq. (62), is numerically obtained from a finite differences technique (Taylor series
expansion). Different techniques were expressed for the calculation of first-order derivation to obtain
tangent stiffness matrix in [45]. By applying the central difference formulation on Eq. (62), it can be
written as

T
⃗ ̂ ⃗ ̂
(63)

IP
in which, represents a small positive perturbation and ̂ is jth degree of freedom (component) of unit

CR
vector ̂ . The accuracy of this numerical equality depends on the perturbation size and it should be
chosen very smaller than material points spacing. In this study, the perturbation size is chosen
times as much as material points spacing. In the following sections, the solution procedures for elastic
analysis and also for elastoplastic analysis are represented.

5.1. Numerical solution procedure for elastic analysisUS

AN
The solution procedure for the quasi-static elastic analysis utilizing Newton Raphson algorithm
could be explained as:

1) Set load step to zero, . Initialize the boundary conditions and the displacements of material
M

points and set both of them to zero, and .

2) Update the load step , and set iteration step to zero, . Update boundary conditions
ED

and initialize the trial value of displacement in the iteration , (or set an initial guess to
).

3) Update the iteration step and calculate the residual vector ⃗ and the norm of residual vector, .
PT

4) If go to step 7 and if go to step 5.

CE

5) Obtain the tangent stiffness matrix regarding to the current configuration, . Determine
by solving equation, ⃗.
AC

6) Set and return to step 3.

7) Set the displacement of material points with trial displacement, .

8) If exerting boundary conditions are not complete, go to step 2.

9) End.

The flowchart of the solution procedure for elastic analysis is demonstrated in Fig. 5. It should be
noted that in any box of the flowchart (see Fig. 5), the corresponding step numbers of the solution
procedure are written in parentheses.

24
 ACCEPTED MANUSCRIPT

T
IP
CR
US
AN
M
ED
PT
CE

Fig. 5. Flowchart of solution procedure for quasi-static elastic analysis

AC

5.2. Numerical solution for elastoplastic analysis

As the elastoplastic analysis is a nonlinear analysis, solving the equations of plasticity is required
an extra iterative procedure. The implemented algorithm for elastoplastic analysis could be
represented as following:

1) Set load step to zero, . Initialize boundary conditions and the displacement of material points
and set both of them to zero, and .

2) Set plastic stretch, equivalent plastic stretch and equivalent back-stress to zero,
, ̅ ̅ .

25
 ACCEPTED MANUSCRIPT

3) Update the load step , and set iteration step to zero, . Update incrementally
boundary conditions and set an initial guess to trial displacement, .

4) Set plastic stretch, equivalent plastic stretch and equivalent back-stress equal to their value in
previous load step as

(64a)

and

̅ ̅ (64b)

T
and

IP
̅ ̅ (64c)

CR
5) Update the iteration step. Obtain incremental stretch between material points, ,
according to incremental displacement,

6) Assume that plastic deformation is not occurred,

as following: US , and set the trial elastic stretch
AN
( ) (65)

7) Find the trial dilatation and the trial equivalent stress as following:
M

̂ ∑ (66a)
ED

and
PT

̅ √ (̂ ∑ ( ) | |

(66b)
CE

( ) )
AC

8) Using the equivalent stress, obtain the value of trial yield function, .

8a) If , the assumption taken in step 6, is correct and

and go to step 16.
8b) If , the assumption taken in step 6, is not correct ( )
and go to 9.

9) Obtain the elastic stretch as

26
 ACCEPTED MANUSCRIPT

(67)

10) Calculate the dilatation and the equivalent stress as

̂ ∑( ) (68a)

and

T
̅ √ (̂ ∑ ( ) | |

IP
(68b)
( ) )

CR
11) Update ̅
hardening model as
̅ ̅
US
and yield function for isotropic, kinematic and mixed
AN
| ̅ | ( ( ̅ ̅ )) (69a)

and
M

| ̅ ̅ ( ̅ ̅ )| (69b)

and
ED

| ̅ ̅ ( ̅ ̅ )|
( ̅ ̅ ) (69c)
PT

12) and ̅ , are functions of as following:

CE

| | (70a)
AC

and

̅ (̂ ∑ ( ) (| |) ) (70b)

13) By using the Newton Raphson method, obtain as

27
 ACCEPTED MANUSCRIPT

( ) (superscript m is iteration number used in an

iterative algorithm)
|
(71)

14) By utilizing obtained in previous step (step 13), calculate the incremental plastic stretch, the
incremental equivalent plastic stretch and the incremental equivalent back-stress (in kinematic and
mixed hardening model) as

| | (72a)

T
and

IP
̅ (̂ ∑ ( ) | | ) (72b)

CR
and

̅ ( ̅ US ) (72c)
AN
15) Obtain the plastic and the elastic stretch of each material point considering incremental plastic
stretch calculated in previous step (step 14) as

(73a)
M

and
ED

(73b)

16) Calculate residual vector ⃗ and norm of residual vector, . If go to step 19 and else go
to step 17.
PT

17) Obtain the tangent stiffness matrix considering as current configuration. Determine
from solving equation, ⃗.
CE

18) Set and return to step 5.

AC

19) Set and update the plastic stretch, the equivalent plastic stretch and the
equivalent back-stress as following:

(74a)

and

̅ ̅ ̅ (74b)

and

̅ ̅ ̅ (74c)

28
 ACCEPTED MANUSCRIPT

20) If exerting boundary conditions are not complete, go to step 3.

21) End.

The flowchart of the algorithm for elasto-plastic analysis is shown in Fig. 6. It should be noted that in
any box of the flowchart (see Fig. 6), the corresponding step numbers of the algorithm are written in
parentheses.

T
IP
CR
US
AN
M
ED
PT
CE
AC

29
 ACCEPTED MANUSCRIPT

T
IP
CR
US
AN
M
ED
PT
CE
AC

Fig. 6. Flowchart of algorithm for elastoplastic analysis

30
 ACCEPTED MANUSCRIPT

6. Results

In this section, different numerical case studies are considered to validate the proposed PD
plasticity models and also to show the capability of this theory in plasticity analysis. For this purpose,
at first, square plates with linear hardening models under uniaxial and biaxial cyclic loading are
considered. Afterward, square plates with nonlinear isotopic and kinematic hardening are analyzed
under cyclic loading and finally, a plate with linear isotropic hardening model containing a hole,
subjected to monotonic tension is considered. For all plates, the results of PD theory are compared
with finite element method results based on local continuum mechanics obtained from Abaqus
software. For comparison purpose, von Mises stress and equivalent stretch are employed in both PD
and FEM models.

T
In all of studied plates, Young's modulus and Poisson's ratio are set as and ,
respectively and the thickness of models is considered as .

IP
To insure the convergence of the presented numerical method in this study, before performing the
2-D elastoplastic numerical studies, the numerical convergence study of the proposed PD model in

CR
elastic domain has been investigated. For this purpose, as presented in [46], -convergence and -
convergence of the PD model are studied in a 2-D benchmark example. In the -convergence, the
horizon is kept fixed and the aim of the analysis is to find the optimal ratio ⁄ , in which
represents distance between material points in the uniform distribution. In -convergence study, as a

US
finial point, the numerical model converges to the precise non-local peridynamic solution for the
given . In -convergence, the best horizon for the PD model is chosen by keeping the ratio
fixed. In this case, the numerical model lastly converges to the local classical solution.
AN
In the benchmark example for the mentioned convergence analyses, a square plate with length of
under a uniaxial displacement loading is considered. The displacement in y direction
, is exerted on the top boundary of the plate and the bottom boundary of the plate is fixed
in y direction as well as the midpoints of the top and bottom boundaries are fixed in x direction. It
M

should be noted that, in the numerical convergence case studies, the von Mises stress is chosen as the
parameter for comparison purpose. According to the analytical solution, the distribution of von Mises
stress in the domain will be uniformly ̅ , due to the displacement loading .
ED

For the -convergence case study, the horizon is kept fixed and the domain is
uniformly discretized by the ratio . The variation of von Mises stress of PD material
points (in circular shapes) in the three discretized domains due to the displacement loading
PT

are demonstrated in Fig. 7(a)-(c). These results show that, by increasing the ratio , the
distribution of von Mises stress becomes smother and von Mises stress of most material points
approach to ̅ . As demonstrated in these figures, the variation of von Mises stress of
CE

material points located near to boundaries are more than other material points. The reason of this
variation is that the horizons of these material points are not completely embedded in the domain.
Although, the strain energy density correction factor and the dilatation correction factor are employed
in these kinds of material points, the differences of results don’t become zero. It should be noted that
AC

by increasing the ratio , the computational cost of the numerical model increases especially for
multistep elastoplastic models and it is crucial to choose a suitable ratio . As shown in Fig. 7(b) and
(c), there aren’t noticeable differences between the results of the discretized domains by the
. Accordingly, the ratio is chosen for the numerical studies of this paper.

31
 ACCEPTED MANUSCRIPT

a b

T
IP
CR
c

US
AN
M
ED

Fig. 7. The distribution of von Mises stress (MPa) due to displacement loading in three discretized domains with
PT

the horizon and the ratio (a) (b) (c)

CE

In the -convergence case study, by keeping the ratio fixed, the domain of benchmark example
is uniformly discretized with horizons . The
variation of von Mises stress of PD material points (in circular shapes) in the four discretized domains
AC

due to the displacement loading are shown in Fig. 8(a)-(d). Considering the results, by
decreasing the horizon size, the distribution of von Mises stress of the domain becomes smother and
approaches to the exact value, ̅ . As demonstrated in Fig. 8(a)-(d), the results of
discretized domain with have considerable differences with others. Therefore, when
the ratio is chosen, the maximum acceptable horizon size for the presented PD model is
in which the errors of von Mises stresses in the domain (in comparison to the exact
solution) don’t exceed .

32
 ACCEPTED MANUSCRIPT

a b

T
IP
CR
c d

US
AN
M
ED
PT

Fig. 8. The distribution of von Mises stress (MPa) due to displacement loading in four discretized domains with
the ratio and the horizon size (a) (b) (c) (d)
CE

6.1. Plates with different hardening under uniaxial cyclic loading

AC

Three square plates with length of are considered with linear kinematic, isotropic and
mixed hardening models. All of models have constant hardening modulus, and
initial yield stress . In addition, the effect factor of mixed hardening model is set as
. In each model, the domain is uniformly discretized by material points (by choosing
the ratio and the horizon ) and a uniaxial displacement in y direction is
cyclically exerted on the top boundary of the domain, as shown in Fig. 9. The displacement in y
direction at the bottom boundary of the domain is set to zero and the midpoints of the top and bottom
boundaries are fixed in x direction.

33
 ACCEPTED MANUSCRIPT

2.5
first c ycle
second cycle
2

1.5

Displacement (mm)
0.5

0
t0 t1 t2 t3 t4 t5 t6
-0.5

-1

-1.5

T
-2

IP
-2.5
T ime step
Fig. 9. A cyclic displacement loading path in y direction

CR
Considering the midpoint of the plate as the origin of Cartesian coordinate, the results of the
material point located in coordinate and are demonstrated in Fig. 10(a)-(d).
The results of the mentioned point in three plates with linear kinematic hardening (see Fig. 10(a)),

US
isotropic hardening (see Fig. 10(b)) and mixed hardening (see Fig. 10(c)), obtained from PD models
are compared with the results of the plates obtained from FEM considering same features and
boundary conditions.
AN
As shown in Fig. 10(a)-(c), the von Mises stress vs. equivalent stretch for three different hardening
models obtained from the proposed PD plasticity formulation for uniaxial cyclic loading are in good
agreement with the FEM results.
The differences of results of PD models with three different hardenings under first cycle of loading
M

are demonstrated in Fig. 10(d). Furthermore, the results of a PD plate with mixed hardening model
and the effect factors and are shown in Fig. 10(d). It should be noted the only
difference between two PD models with mixed hardening is the effect factor. As shown in Fig. 10(d),
ED

all of graphs are coincided each other up to the end of tensile loading, according to their same
hardening modulus, same initial yield stress as well as their initial hardening which is set to zero,
̅ ̅ . However, in unloading and reloading path the models behave diversely, as
PT

expected, and the curves are not coincided.

CE
AC

34
 ACCEPTED MANUSCRIPT

a b
300 400

300
200

200
Von Mises stress (MPa)

Von Mises stress (MPa)

100
100

0 0
-0.003 -0.002 -0.001 0 0.001 0.002 0.003 -0.003 -0.002 -0.001 0 0.001 0.002 0.003

-100

T
-100

-200

IP
-200 PD- kinem atic
-300 PD-isotropic

FEM-kinem atic

CR
FEM-isotropic
-300 -400
Equivalent stretch (mm/mm) Equivalent stretch (mm/mm)

c d
300

200
US 300
AN
200
Von Mises stress (MPa)

100
Von Mises stress (MPa)

100

0
M

-0.003 -0.002 -0.001 0 0.001 0.002 0.003

0
-0.002 -0.001 0 0.001 0.002
-100

-100 PD-kinema tic

ED

-200
PD-isotropic

-300 PD-mixed -200

P D-mixed ξ = 0.7

FEM-m ixed P D-mixed ξ = 0.3

PT

-400 -300
Equivalent stretch (mm/mm) Equivalent stretch (mm/mm)

Fig. 10. The von Mises stress vs. equivalent stretch of the material (integration) point located in coordinate ,
CE

, obtained from PD and FEM models with (a) linear kinematic hardening (b) linear isotropic hardening (c) linear
mixed hardening ( ); (d) obtained from PD models with linear isotropic, kinematic and two different mixed
hardenings ( )
AC

As mentioned in section (3.4), in the mixed models, the expansion and movement of initial yield
surface are happened simultaneously in the stress space and the results include the combination of two
isotropic and kinematic hardenings. The effect factor in the linear mixed hardening model
determines the proportion of isotropic and kinematic hardening. In the mixed model with effect factor
, the proportion of isotropic hardening is more than kinematic hardening and as expected, the
curve of mixed model is closer to the isotropic hardening. However, for the mixed model with effect
factor , the effect of translation of yield surface is more than isotropic expansion of the yield
locus and the behavior of mixed model is closer to the kinematic hardening.
In Fig. 11 the displacement components of material points in x and y directions are compared for
PD and FEM models. For this purpose, the plate is considered to be under the applied displacement

35
 ACCEPTED MANUSCRIPT

loading (t1 in first cycle). It should be noted that the PD material points are shown in
circular shape in Fig. 11 which are different from FEM elements.
As shown in Fig. 11(a) and (b), the variation of displacements in y direction for FEM and PD
domains are the same. Similarly, as shown in Fig. 11(c) and (d), the variation of displacements in x
direction demonstrated in FEM and PD domains are very close to each other too. In addition, similar
to FEM model, in PD model the volume of domain under plastic deformation is not changed. This is
because of assuming the incompressibility condition under plastic stretch in PD plasticity model.

a b

T
IP
CR
US
AN

c d
M
ED
PT
CE
AC

Fig. 11. The variation of displacement (mm) due to at the top of the plate: (a) FEM displacement in y
direction (b) PD displacement in y direction (c) FEM displacement in x direction (d) PD displacement in x direction

6.2. A plate with linear kinematic hardening under biaxial cyclic loading

To validate the PD plasticity model with kinematic hardening, a square plate with length of
is cyclically loaded under two different biaxial displacement loading paths. The plate has
the constant hardening modulus and initial yield stress . The domain of
plate is uniformly discretized by material points with same volumes (by choosing the ratio

36
 ACCEPTED MANUSCRIPT

and the horizon ). The displacement loading in y-direction, , and in x-direction,

, are exerted on top and right boundaries of plate, respectively. Fig. 12(a) and (b), show two
different paths for displacement loadings and . In both of loading paths, the left and bottom
boundaries of the domain are fixed in x and y direction, respectively.
2 b
2
1.5

1.5
1

1
Displacement (mm)

0.5

Displacement (mm)
0.5

T
0
t0 t1 t2 t3 0

IP
-0.5 t0 t1 t2 t3
-0.5
-1

CR
-1
-1.5 ux
u_x -1.5 ux
-2
u_x

US
T ime step
-2
a T ime step

Fig. 12. Variation of displacement boundary conditions in x direction and y direction in (a) first loading path (b)
AN
second loading path

As shown in Fig. 13(a) and (b), in both loading paths, the curve of von Mises stress vs. equivalent
stretch in the midpoint of plate obtained from PD and FEM models are very close and almost
coincided to each other. These results show that the proposed PD plasticity model with kinematic
M

hardening, properly and accurately estimates the movement of yield surface and the hardening of
medium in loading/unloading processes beyond initial yield surface.
ED

a b
300 300
PT

200 200
Von Mises stress (MPa)
CE
Von Mises stress (MPa)

100 100

0 0
-0.002 -0.001 0 0.001 0.002 -0.003 -0.002 -0.001 0 0.001 0.002 0.003
AC

-100 -100

-200 -200 PD
PD

FEM FEM
-300 -300
Equivalent stretch (mm/mm) Equivalent stretch (mm/mm)

Fig. 13. The von Mises stress vs. equivalent stretch of midpoint of the plate obtained from PD and FEM models under: (a)
first loading path (b) second loading path

37
 ACCEPTED MANUSCRIPT

6.3. Plates with different nonlinear hardenings under uniaxial cyclic loading

In this section, the case studies are square plates with nonlinear hardening and length of
which are subjected to uniaxial cyclic loading path. For this purpose, an exponential hardening model,
known as Voce hardening, is chosen for the nonlinear isotopic hardening model. The general form of
Voce hardening yield surface is proposed by Simo and Hughes [4] as

( ̅ ) ̅ * ( ̅ )+ (75)

in which , and are material parameters.

For kinematic hardening, the same nonlinear hardening model is chosen and the incremental
equivalent back-stress, ̅ , in any time step is defined as following:

T
IP
̅ [ ( ̅ )] (

CR
(76)
( ̅ ))

, , , and
The domains of the models are uniformly discretized by
. US
In both of isotropic and kinematic models, the initial yield stress and the material parameters are set as

material points (by choosing the

AN
ratio and the horizon ) and a cyclic displacement in y direction (see Fig. 14), is
applied on top boundary of the domain. The bottom boundary of the domain is fixed in y direction as
well as the midpoints of top and bottom boundaries are fixed in x direction.
M

2.5

2
ED

1.5
Displacement (mm)

0.5
PT

0
t0 t1 t2 t3
-0.5

-1
CE

-1.5

-2

-2.5
AC

T ime step

Fig. 14. The displacement load path in y direction applied on the top boundary of domain

Considering the midpoint of the plate as the origin of Cartesian coordinate, the results of material
point located in coordinate , , obtained from PD and FEM models are shown
in Fig. 15(a)-(c). As shown in Fig. 15(a), the curves of von Mises stress vs. equivalent stretch for one
cycle obtained from PD and FEM models with nonlinear isotropic hardening are coincided. Fig. 15(b)
shows the curves of von Mises stress vs. equivalent stretch obtained from PD and FEM models for the
mentioned kinematic hardening model. As demonstrated in this figure, the PD and FEM results are
coincided for two cycles of loadings.

38
 ACCEPTED MANUSCRIPT

The result of PD models with nonlinear isotropic and kinematic hardening are compared in Fig.
15(c). As shown in this figure, the curves of von Mises stress vs. equivalent stretch for isotropic and
kinematic hardening are coincided up to end of first loading path (note that their initial hardening is
set to zero, ̅ ̅ ), however in unloading and reloading paths, the models behave in
different ways. Therefore, the proposed PD plasticity models with the nonlinear isotropic or kinematic
hardening could predict the expected results which agreed with the classical theory of plasticity in
cyclic manner.

a b
300 300

T
200
200
Von Mises stress (MPa)

IP
100

Von Mises stress (MPa)

100

CR
0
-0.003 -0.002 -0.001 0 0.001 0.002 0.003
0
-0.003 -0.002 -0.001 0 0.001 0.002 0.003
-100

-100
-200

-300 PD

FEM
US -200 PD

FEM
AN
-400 -300
Equivalent stretch (mm/mm) Equivalent stretch (mm/mm)

c
M

300

200
ED
Von Mises stress (MPa)

100

0
-0.003 -0.002 -0.001 0 0.001 0.002 0.003
PT

-100

-200
CE

isotr opic
hardening
-300
kinematic
hardening
-400
AC

Equivalent stretch (mm/mm)

Fig. 15. The von Mises stress vs. equivalent stretch of point and : obtained from PD and FEM
models with (a) nonlinear isotropic hardening (b) nonlinear kinematic hardening; (c) obtained from PD models with
nonlinear isotropic and kinematic hardening models

6.4. A plate with a hole under uniaxial loading

A square plate with length of containing a circular hole with radius of at the
center of the plate, is subjected to a uniaxial displacement loading symmetrically exerted on the top
and bottom boundaries of the plate. The material model of the plate is considered as linear isotropic

39
 ACCEPTED MANUSCRIPT

hardening with constant hardening modulus and initial yield stress . The
applied displacement in y direction, , is incrementally increased from zero to and the
midpoints of top and bottom boundaries are fixed in x direction. Using Abaqus software, the FEM
domain is discretized to 9889 non-uniform elements with different element size. Same as the FEM
mesh generation, the non-uniform material point distribution is used in PD model discretization. It
should be noted that in the discretized PD domain, the average space between material points is about
and by choosing the ratio , the horizon is fixed for all of material points.
̂ ̂
As mentioned in section (2.2), the parameters and depend on horizon and for a non-uniform
distribution of material points, the family members located in the horizon differs for each material
point. Therefore, in a non-uniform distribution, the values of parameters ̂ and ̂ presented in section
(2.2), are not suitable for all of material points and all the PD parameters depended on parameters ̂

T
and ̂ need to be corrected. For this purpose, the strain energy density correction factor and the

IP
dilatation correction factor are employed for all of material points to correct the PD parameters such
as the strain energy density, von Mises stress, the force density vector based on the strain energy
density and the dilatation.

CR
The variation of von Mises stress of PD material points and FEM elements are demonstrated in
Fig. 16(a) and (b), due to the displacement loading . The PD material points are shown
in circular shapes (see Fig. 16(b)), which their radiuses are varied according to their volume size and

US
the elements in FEM models are shown by quadrilateral elements. It should be noted that the von
Mises stress in these figures is not reached the initial yield stress of the material.
AN
a b
M
ED
PT
CE
AC

Fig. 16. The distribution of von Mises stress (MPa) due to displacement loading in (a) FEM elements (b) PD
material points

By increasing the displacement load up to , the stress in some material points exceeds
the initial yield stress and these material points are experienced plastic deformation. For this loading,
the distribution of von Mises stress for two FEM and PD models are demonstrated in Fig. 17(a) and
(b). It is obvious from these figures that, the PD plasticity model can predict the same distribution for
von Mises stress of the domain as FEM model.

40
 ACCEPTED MANUSCRIPT

a b

T
IP
CR
Fig. 17. The distribution of von Mises stress (MPa) due to displacement loading in (a) FEM elements (b) PD
material points

US
Similar to the von Mises stress, under displacement loading , the variation of
equivalent plastic stretch of PD material points and FEM elements are shown in Fig. 18(a) and (b).
According to these figures, it could be concluded that the proposed PD plasticity model is capable of
AN
distinguishing the material points which experience the plastic deformation during the loading. In
addition, the total equivalent plastic stretch of the material points in PD model is in good agreement
with those obtained from the FEM.
M

a b
ED
PT
CE
AC

Fig. 18. The total equivalent plastic stretch (mm/mm) due to displacement loading in (a) FEM elements (b) PD
material points

Fig. 19(a) and (b) demonstrate the variation of displacement in y direction for FEM domain and
PD material points. The variation of displacement in x direction for both FEM and PD models are
shown in Fig. 20(a) and (b). By comparing the FEM results and PD ones in Figs. 19 and 20, it is
obtained that the proposed PD model is capable of predicting the displacement fields analogous and

41
 ACCEPTED MANUSCRIPT

near the displacement filed in FEM model which is formulated based on classical local continuum
mechanic assumptions.

a b

T
IP
CR
PD material points US
Fig. 19. The variation of displacement (mm) in y direction due to displacement loading (a) FEM elements (b)
AN
To compare the results quantitatively in FEM and PD models, two material points near the hole,
where the stress gradient is noticeably great, are chosen (see Fig. 21). Material point 1 is located at
and and material point 2 is located at and
.
M

The von Mises stress for applied displacement loadings and and equivalent
plastic stretch due to applied displacement loading obtained from PD and FEM models
for material points 1 and 2 are presented in Table 1. According to the results given in this table, it is
ED

obvious that the non-uniform PD model can predict the plasticity deformation in domain accurately,
in comparison to the FEM results.

a b
PT
CE
AC

Fig. 20. The variation of displacement (mm) in x direction due to displacement loading (a) FEM elements (b)
PD material points

42
 ACCEPTED MANUSCRIPT

a b

T
IP
CR
Table 1
US
Fig. 21. Two chosen material points (a) in the domain (b) in the magnified view
AN
The results of two mentioned points obtained from PD and FEM models

Material point number

Results Model
Point 1 Point 2
M

von Mises stress (MPa) PD 483.627 354.439

for applied
(elastic domain) FEM 483.667 349.025
ED

von Mises stress (MPa) PD 652.736 579.606

for applied
(beyond elastic domain) FEM 659.996 580.546
PT

Equivalent plastic stretch PD 0.003060 0.001593

(mm/mm) for applied
CE

FEM 0.003204 0.001614

7. Conclusions
AC

In this paper, an elastoplastic model based on ordinary state-based peridynamic theory is proposed
for cyclic loading analyses. The von Mises yield criteria is used to describe plastic yielding and the
equivalent plastic stretch is used as an internal variable for different hardening models.
The general form of isotropic hardening, the modified form of kinematic hardening and the mixed
of kinematic and isotopic hardenings are considered as the hardening of material in PD constitutive
relations. The proposed plastic model is also described in the thermodynamic framework and it is
shown that the proposed plastic flow and hardening of the PD model are not violated the second law
of thermodynamics.
The equations of elastoplastic PD models with different kinds of hardening are linearized and the
obtained equations are solved numerically by the proposed numerical solution approach.

43
 ACCEPTED MANUSCRIPT

To validate the presented PD plasticity model, different benchmark problems are considered. The
2-D plate models with linear/nonlinear isotopic and kinematic hardening as well as linear mixed
hardening under uniaxial cyclic loading are analyzed by the presented PD plasticity model and FEM
too. In addition, the 2-D plate model with linear kinematic hardening under biaxial cyclic loading is
investigated by using PD and FEM methods. Finally, a plate with hole and linear isotropic hardening
are considered as benchmark problem and the results of PD models with non-uniform distribution of
material points are compared with FEM ones.
By comparison the obtained results, it could be concluded that the proposed PD elastoplastic
models can predict plastic yielding and hardening of materials beyond initial yield stress accurately
which are very close to the FEM results.
Although, defining the kinematic hardening, which depends on back stress tensor in classical

T
theory of plasticity, is complicated in PD formulated, the presented model can handle this challenge
too. The obtained results show that the presented PD formulation is capable of modeling the

IP
kinematic hardening in cyclic loading very well and the movement of the yield surface in the stress
space could be captured accurately by the presented PD model.

CR
According to the results it could be concluded that, the developed PD constitutive relations is
capable of modeling the general form of elastoplastic model used in continuum mechanics and
classical theory of plasticity. Therefore, by employing this model in conjunction with the damage
assumptions, it would be possible to predict the damage growth and crack nucleation in ductile
materials by PD models in future studies.
US
AN
References
[1] Drucker DC. A more fundamental approach to plastic stress-strain relations. Proc. First U.S.
Natl. Congr. Appl. Mech. ASME, 1951, p. 487–91.
[2] Prager W. The Theory of Plasticity: A Survey of Recent Achievements. Proc Inst Mech Eng
M

1955;169:41–57. doi:10.1243/PIME_PROC_1955_169_015_02.
[3] Chen WF, Han DJ. Plasticity for Structural Engineers. New York, NY: Springer New York;
1988. doi:10.1007/978-1-4612-3864-5.
ED

[4] Simo JC, Hughes TJ. Computational Inelasticity. vol. 7. New York: Springer-Verlag; 1998.
doi:10.1007/b98904.
[5] de Souza Neto EA, Peri D, Owen DRJ. Computational Methods for Plasticity. Chichester, UK:
John Wiley & Sons, Ltd; 2008. doi:10.1002/9780470694626.
PT

[6] Lemaitre J. A Course on Damage Mechanics. Berlin, Heidelberg: Springer Berlin Heidelberg;
1992. doi:10.1007/978-3-662-02761-5.
[7] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J
CE

Mech Phys Solids 2000;48:175–209. doi:10.1016/S0022-5096(99)00029-0.

[8] Silling SA, Epton M, Weckner O, Xu J, Askari E. Peridynamic States and Constitutive
Modeling. J Elast 2007;88:151–84. doi:10.1007/s10659-007-9125-1.
[9] Silling SA, Lehoucq RB. Peridynamic Theory of Solid Mechanics. Adv. Appl. Mech., vol. 44,
AC

Elsevier Inc.; 2010, p. 73–168. doi:10.1016/S0065-2156(10)44002-8.

[10] Madenci E, Oterkus E. Peridynamic Theory and Its Applications. vol. c. New York, NY:
Springer New York; 2014. doi:10.1007/978-1-4614-8465-3.
[11] Oterkus S, Madenci E, Agwai A. Fully coupled peridynamic thermomechanics. J Mech Phys
Solids 2014;64:1–23. doi:10.1016/j.jmps.2013.10.011.
[12] Huang D, Lu G, Qiao P. An improved peridynamic approach for quasi-static elastic
deformation and brittle fracture analysis. Int J Mech Sci 2015;94–95:111–22.
doi:10.1016/j.ijmecsci.2015.02.018.
[13] Zhou X-P, Gu X-B, Wang Y-T. Numerical simulations of propagation, bifurcation and
coalescence of cracks in rocks. Int J Rock Mech Min Sci 2015;80:241–54.
doi:10.1016/j.ijrmms.2015.09.006.
[14] Wang Y, Zhou X, Shou Y. The modeling of crack propagation and coalescence in rocks under

44
 ACCEPTED MANUSCRIPT

uniaxial compression using the novel conjugated bond-based peridynamics. Int J Mech Sci
2017;128–129:614–43. doi:10.1016/j.ijmecsci.2017.05.019.
[15] Wang Y-T, Zhou X-P, Kou M-M. Three-dimensional numerical study on the failure
characteristics of intermittent fissures under compressive-shear loads. Acta Geotech 2018;7.
doi:10.1007/s11440-018-0709-7.
[16] Wang Y, Zhou X, Wang Y, Shou Y. A 3-D conjugated bond-pair-based peridynamic
formulation for initiation and propagation of cracks in brittle solids. Int J Solids Struct
2018;134:89–115. doi:10.1016/j.ijsolstr.2017.10.022.
[17] Wang Y, Zhou X, Kou M. An improved coupled thermo-mechanic bond-based peridynamic
model for cracking behaviors in brittle solids subjected to thermal shocks. Eur J Mech -
A/Solids 2019;73:282–305. doi:10.1016/j.euromechsol.2018.09.007.
[18] Xu Z, Zhang G, Chen Z, Bobaru F. Elastic vortices and thermally-driven cracks in brittle

T
materials with peridynamics. Int J Fract 2018;209:203–22. doi:10.1007/s10704-017-0256-5.
[19] Wang Y, Zhou X, Xu X. Numerical simulation of propagation and coalescence of flaws in

IP
rock materials under compressive loads using the extended non-ordinary state-based
peridynamics. Eng Fract Mech 2016;163:248–73. doi:10.1016/j.engfracmech.2016.06.013.
[20] Zhou X, Wang Y, Qian Q. Numerical simulation of crack curving and branching in brittle

CR
materials under dynamic loads using the extended non-ordinary state-based peridynamics. Eur
J Mech - A/Solids 2016;60:277–99. doi:10.1016/j.euromechsol.2016.08.009.
[21] Zhou X-P, Wang Y-T. Numerical simulation of crack propagation and coalescence in pre-

[22] US
cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics. Int J Rock
Mech Min Sci 2016;89:235–49. doi:10.1016/j.ijrmms.2016.09.010.
Lai X, Liu L, Li S, Zeleke M, Liu Q, Wang Z. A non-ordinary state-based peridynamics
modeling of fractures in quasi-brittle materials. Int J Impact Eng 2018;111:130–46.
AN
doi:10.1016/j.ijimpeng.2017.08.008.
[23] Tupek MR, Rimoli JJ, Radovitzky R. An approach for incorporating classical continuum
damage models in state-based peridynamics. Comput Methods Appl Mech Eng 2013;263:20–
6. doi:10.1016/j.cma.2013.04.012.
M

[24] Roy P, Pathrikar A, Deepu SP, Roy D. Peridynamics damage model through phase field
theory. Int J Mech Sci 2017;128–129:181–93. doi:10.1016/j.ijmecsci.2017.04.016.
[25] Han F, Lubineau G, Azdoud Y. Adaptive coupling between damage mechanics and
peridynamics: A route for objective simulation of material degradation up to complete failure.
ED

J Mech Phys Solids 2016;94:453–72. doi:10.1016/j.jmps.2016.05.017.

[26] Galvanetto U, Mudric T, Shojaei A, Zaccariotto M. An effective way to couple FEM meshes
and Peridynamics grids for the solution of static equilibrium problems. Mech Res Commun
PT

2016;76:41–7. doi:10.1016/j.mechrescom.2016.06.006.
[27] Shojaei A, Zaccariotto M, Galvanetto U. Coupling of 2D discretized Peridynamics with a
meshless method based on classical elasticity using switching of nodal behaviour. Eng Comput
2017;34:1334–66. doi:10.1108/EC-03-2016-0078.
CE

[28] Shojaei A, Mudric T, Zaccariotto M, Galvanetto U. A coupled meshless finite

point/Peridynamic method for 2D dynamic fracture analysis. Int J Mech Sci 2016;119:419–31.
doi:10.1016/j.ijmecsci.2016.11.003.
AC

[29] Shojaei A, Mossaiby F, Zaccariotto M, Galvanetto U. An adaptive multi-grid peridynamic

method for dynamic fracture analysis. Int J Mech Sci 2018;144:600–17.
doi:10.1016/j.ijmecsci.2018.06.020.
[30] Tong Q, Li S. Multiscale coupling of molecular dynamics and peridynamics. J Mech Phys
Solids 2016;95:169–87. doi:10.1016/j.jmps.2016.05.032.
[31] Foster JT, Silling SA, Chen WW. Viscoplasticity using peridynamics. Int J Numer Methods
Eng 2009;81:n/a-n/a. doi:10.1002/nme.2725.
[32] Mitchell JA. A Nonlocal , Ordinary , State-Based Plasticity Model for Peridynamics.
Technical Report SAND2011-3166, Sandia National Laboratories, Albuquerque, New Mexico,
2011.
[33] Lammi CJ, Vogler TJ. A Nonlocal Peridynamic Plasticity Model for the Dynamic Flow and
Fracture of Concrete. SAND2014-1825, Sandia National Laboratories, Albuquerque, New
Mexico, 2014.

45
 ACCEPTED MANUSCRIPT

[34] Sun S, Sundararaghavan V. A peridynamic implementation of crystal plasticity. Int J Solids

Struct 2014;51:3350–60. doi:10.1016/j.ijsolstr.2014.05.027.
[35] Amani J, Oterkus E, Areias P, Zi G, Nguyen-Thoi T, Rabczuk T. A non-ordinary state-based
peridynamics formulation for thermoplastic fracture. Int J Impact Eng 2016;87:83–94.
doi:10.1016/j.ijimpeng.2015.06.019.
[36] Madenci E, Oterkus S. Ordinary state-based peridynamics for plastic deformation according to
von Mises yield criteria with isotropic hardening. J Mech Phys Solids 2016;86:192–219.
doi:10.1016/j.jmps.2015.09.016.
[37] Song X, Khalili N. A peridynamics model for strain localization analysis of geomaterials. Int J
Numer Anal Methods Geomech 2019;43:77–96. doi:10.1002/nag.2854.
[38] Rahaman MM, Roy P, Roy D, Reddy JN. A peridynamic model for plasticity: Micro-inertia
based flow rule, entropy equivalence and localization residuals. Comput Methods Appl Mech

T
Eng 2017;327:369–91. doi:10.1016/j.cma.2017.07.034.
[39] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics.

IP
Comput Struct 2005;83:1526–35. doi:10.1016/j.compstruc.2004.11.026.
[40] Saad MH. Elasticity. 2nd ed. Burlington, MA: Elsevier; 2009. doi:10.1016/B978-0-12-
374446-3.X0001-6.

CR
[41] Key SW, Stone CM, Krieg RD. Dynamic Relaxation Applied to the Quasi-Static, Large
Deformation, Inelastic Response of Axisymmetric Solids. Nonlinear Finite Elem. Anal. Struct.
Mech., Berlin, Heidelberg: Springer Berlin Heidelberg; 1981, p. 585–620. doi:10.1007/978-3-

[42]

[43]
642-81589-8_30.

US
Macek RW, Silling SA. Peridynamics via finite element analysis. Finite Elem Anal Des
2007;43:1169–78. doi:10.1016/j.finel.2007.08.012.
Seleson P. Improved one-point quadrature algorithms for two-dimensional peridynamic
AN
models based on analytical calculations. Comput Methods Appl Mech Eng 2014;282:184–217.
doi:10.1016/j.cma.2014.06.016.
[44] Littlewood DJ. Roadmap for Peridynamic Software Implementation. SAND2015-9013, Sandia
National Laboratories, Albuquerque, New Mexico, 2015.
M

[45] Brothers MD, Foster JT, Millwater HR. A comparison of different methods for calculating
tangent-stiffness matrices in a massively parallel computational peridynamics code. Comput
Methods Appl Mech Eng 2014;279:247–67. doi:10.1016/j.cma.2014.06.034.
[46] Bobaru F, Yang M, Alves LF, Silling SA, Askari E, Xu J. Convergence, adaptive refinement,
ED

and scaling in 1D peridynamics. Int J Numer Methods Eng 2009;77:852–77.

doi:10.1002/nme.2439.
PT
CE
AC

46