The essential work of fracture in peridynamics

In this work, the essential work of fracture (EWF) method is introduced for a peridynamic (PD) material model to characterize fracture toughness of ductile materials. First, an analytical derivation for the path-independence of the PD J-integral is provided. Thereafter, the classical J-integral and PD J-integral are computed on a number of analytical crack problems, for subsequent investigation on how it performs under large scale yielding of thin sheets. To represent a highly nonlinear elastic behavior, a new adaptive bond stiffness calibration and a modified bond-damage model with gradual softening are proposed. The model is employed for two different materials: a lower-ductility bainitic-martensitic steel and a higher-ductility bainitic steel. Up to the start of the softening phase, the PD model recovers the experimentally obtained stress–strain response of both materials. Due to the high failure sensitivity on the presence of defects for the lower-ductility material, the PD model could not recover the experimentally obtained EWF values. For the higher-ductility bainitic material, the PD model was able to match very well the experimentally obtained EWF values. Moreover, the J-integral value obtained from the PD model, at the absolute maximum specimen load, matched the corresponding EWF value.


Introduction
The peridynamic (PD) theory is a nonlocal formulation of solid mechanics, originally introduced for handling crack initiation, extension and final failure of a body, without the need of supplementary methods (Silling 2000;Silling and Askari 2005). PD is based upon integro-differential equations, thereby avoiding spatial derivatives, which are not defined at discontinuities, such as crack surfaces. For brittle and small scale yielding fracture, in which methods from linear and nonlinear elastic fracture mechanics apply, fracture characterization can be achieved with a single parameter, the J -integral value, for example. In certain cases of large scale yielding, ductile failure may be similarly characterized by the essential work of fracture (EWF) method. We next briefly review these concepts from the classical theory and, when available, their perdynamic-based versions.

The classical and the PD J -integral
The J -integral is an expression for calculating the strain energy release rate of a cracked body, or the amount of energy available at the tip of a crack to form new crack surfaces as the crack extends (Rice 1968). It is a firmly established parameter that applies both to linear and nonlinear elastic materials. The J -integral is used as a toughness measure and failure criterion, for ranking of materials and for dimensioning, but it can also be used for testing the accuracy of material models, such as PD models. Furthermore, since the J -integral can be expressed and discretized in various ways, there is also a need for comparison of different approaches. In the PD literature, both Rice's classical J -integral (J R ) and the nonlocal PD formulation of Rice's J -integral (J PD ) have been studied for testing the expressions and for testing PD numerical implementations, in comparison with analytical values, experimental data and the finite element method (FEM) of the classical model.
The J PD has been derived by Silling and Lehoucq (2010) for state-based PD, based on an energy balance approach. Later, Hu et al. (2012) presented a bondbased J PD , using an infinitesimal virtual crack extension method. The bond-based J PD expression is a special case of the more general state-based J PD expression.
The displacements are the principal unknowns in PD, from which other quantities subsequently can be obtained. A general expression of the contour J R as a function of displacement derivatives was derived by Bruck (1989) and Hart and Bruck (2021). A rectangular integration treatment of the J R has been given by Stenström and Eriksson (2019) for use in PD.
Both the J R and J PD expressions have for separate studies shown results in line with analytical expressions, experimental results and FEM calculations of the classical model. Hu et al. (2012) evaluated the J PD on single and double edge notched PD models, and with refined discretization, achieved results approaching J R of FEM solutions of the classical model. In the formulation of Hu et al. (2012), the stress tensor is considered through computation of a set of force interactions over inner and outer regions associated with the contour integral. The width of these regions depends on the degree of nonlocality of the J PD . Panchadhara and Gordon (2016) used the method of Hu et al. (2012), but approximated the stress tensor of J PD from force interactions at the contour. Details were left out, but the number of interactions to account for is stated by the authors to be fewer. Nevertheless, it need to be noticed that in Hu et al.'s expression, the force interactions is counted only from one region to the other and not from both. Breitenfeld et al. (2014) computed the J R by using a material point method (particle-based) and the displacement gradient of state-based PD. The displacement boundary condition along a circular boundary of a mode I PD specimen corresponded to the classical analytical K-field solution. Using a rectangular J R , the K-field solution and the corresponding PD model were compared, while varying the lattice orientation. Stenström and Eriksson (2019) also used the J R expression by writing it as a function of displacement derivatives and approximating the displacement derivatives by using the central difference scheme on the material model's displacement field, as used by Hu et al. (2012) for the J PD . The specimen was a central crack with mode I stress boundary condition along a rectangular boundary, extracted from an exact analytical stress-strain-displacement solution. The J R was then applied on the classical exact analytical displacement field and on the corresponding PD model and compared, while studying the convergence.
The above studies were for mode I loading. Panchadhara and Gordon (2016) studied the J PD for mode I edge notch and mode II Kalthoff-Winkler PD models, with results in agreement to published experimental results. Jung and Seok (2017) studied fatigue crack modeling for mixed-mode loading and decomposed the J PD into mode I and II, with comparison to experimental results. Furthermore on mixed mode, Imachi et al. (2018) applied the J R to inclined cracks of rectangular and L-shaped PD specimens, with close resemblance to experimental results. Imachi et al. derived the displacement derivatives based on the moving least-squares approximation.
Considering the area/domain formulation of J , the area J R has been studied within PD by Imachi et al. (2018), Fallah et al. (2020 and Stenström and Eriksson (2021). Imachi et al. used the area formulation in the above previously mentioned mixed mode study. Fallah et al. derived the bond-based PD stress tensor and used it for calculating the area J R and compared it to the J PD of Hu et al. (2012). Lastly, Stenström and Eriksson carried out a comparison of the contour J R and the area J R . Nevertheless, the studies does not compare the J PD and J R integrals.
The various J expressions of the above mentioned studies show in general errors less than one percent when applied on PD displacement fields and compared with the analytical expressions of a classical model or FEM approximations of the J R . With finer PD discretization, the error can be much less than one percent, and vice versa for coarser discretization. When the J expressions are applied to classical exact analytical displacement solutions, the error is less than a per mille from J 0 = K 2 I /E (Breitenfeld et al. 2014;Stenström and Eriksson 2019).
In the present study, we will use the following notation of J : -Classical J computed on displacements of a classical continuum mechanics material model is termed J R . -PD J computed on displacements of a PD material model is termed J PD . -J R and J PD values computed on a different model than their natural one is termed J * R and J * PD , respectively.
-Critical J are denoted by adding a 'c', e.g. J * Rc . -Lastly, the classical expression K 2 I /E = J 0 .
The reason why we would like to compute J on a different model than their natural one is that the contour J R is less affected by PD boundary effects of narrow ligaments of double edge notched tension (DENT) specimens typically used in the EWF method. As J PD is computed along a strip, it comes closer to the boundaries of a coarsely discretized ligament. However, an advantages of J PD is that it can average out local disturbances of the displacement field, similar to the classical area J R .

The EWF method to evaluate fracture toughness
The EWF method considers the fracture energy with two parameters; the energy required to form new crack surfaces and the energy dissipated during plastic deformation. The former one can be approximated to the critical J c as long as a set of requirements are met. Since J only considers crack initiation resistance, the EWF method may be preferable over J for studying overall fracture resistance, if there is a significant energy contribution from plastic yielding (Frómeta et al. 2020). Broberg (1968) laid the foundation for the determination of fracture toughness of thin sheet material. The region in front of a crack tip can be separated into an inner and outer zone. The inner zone is the fracture process zone (FPZ). For thin ductile materials, it is possible to identify the FPZ by necking. In plane stress, the necking region is dependent on material thickness. The energy dissipated from the FPZ can be seen as a thickness dependent pseudo material constant, while the outer much larger region is characterized by plastic deformation. Cotterell and Reddel (1977) defined the conditions for the determination of EWF for thin sheets with regard to ligament length and thickness. An experimental method for EWF to be applied to thin sheet metal was proposed by Marchal and Delannay (1996). The EWF method have been used successfully on a wide range of materials including polymers (Bárány et al. 2010), aluminium (Vendra et al. 2017), and steel (Efthymiadis et al. 2017). Casellas et al. (2017) and Frómeta et al. (2020) showed the possibility of fracture toughness determination using the EWF method for several advanced high-strength steel (AHSS) grades like, dual and complex phase, press hardened, TRIP (transformation induced plasticity) and high manganese steels.

Peridynamic material models and their calibration
In ductile failure problems in which the EWF can be useful, the material response is nonlinear. To represent that behavior, a PD model requires calibration of the parameters that determine the force-strain relationship of a PD bond. Calibration for linear-elastic material response is simple and it is based on enforcing a match between the strain energy density of a classical and of a PD material. For a nonlinear material response, this calibration can be accomplished in several ways.
PD constitutive models can be categorized into bond-based, ordinary state-based, and non-ordinary state-based, with the bond-based type being a special case of the ordinary state-based type. In bond-based PD, the pairwise force between two material points depends only on the two material points and no other points (Silling 2000). In state-based PD, deformation at each bond depends on the collective deformation in its vicinity ). In ordinary statebased PD, the pairwise force acts in the direction of the corresponding bond, a restriction that is relaxed for non-ordinary state-based PD. Classical constitutive models can be "directly embedded" in PD formulations via the so-called correspondence models (initiated in Silling et al. 2007), which tend to be of the non-ordinary kind but suffer from zero energy modes instabilities in regions of high strain gradients, such as around crack tips. Various stabilization methods have been proposed to resolve this issue, among them being the bondassociated PD models, which provide accurate and stable correspondence-based PD solutions (Breitzman and Dayal 2018;Chen 2018;Gu et al. 2019;Behzadinasab and Foster 2020;Behzadinasab et al. 2021).
The original constitutive law is the prototype microelastic brittle material (Silling 2000), accompanying a linear elastic bond law up to a sudden failure, to mimic brittle-type failure. Subsequently, in order to better capture quasi-brittle damage and failure, models in which bond force-strain features a softening phase before full bond failure was developed (Macek and Silling 2007). These models utilize, for example, a piecewise linear representation for bond force in terms of bond relative elongation, with "bilinear" (linear up to max load, then linear for the softening phase to zero force) and "trilinear" being the most popular. Multi-parameter models, like the piecewise linear ones, are needed in order to be able to match both the critical fracture energy and strength of a material when using an arbitrary PD horizon size. With a single-parameter bond-damage model, like the prototype microelastic brittle model, one can match, for example, the material's fracture toughness when using an arbitrary horizon size, but it can also match the material's strength only for a specific value of the horizon size. That value, in many instances, ends up being too small to be of practical use, since it leads to very large-scale computational models.
Bilinear bond force-strain models have been studied previously for quasi-brittle materials, considering crack nucleation ) and crack propagation (Zaccariotto et al. 2015;Niazi et al. 2021). In these studies, curve fitting was used to match the bilinear model to an experimental stress-strain response, at the point where the stress-strain departs from linearity and crack nucleation/propagation takes place. Bilinear fitting for quasi-brittle materials have also been accomplished by use of material properties (Yang et al. 2018) and by analytical formulas (Niazi et al. 2021). Bilin-ear models have also been employed to model elastoplastic behavior of metals, considering high speed loading (Macek and Silling 2007;Lee et al. 2016) and crack mouth opening displacement (Yolum et al. 2016). Material softening of quasi-brittle materials has as well been studied, by using a trilinear law (Yang et al. 2018;Niazi et al. 2021).
In state-based PD, the original constitutive model is the linear peridynamic solid , later accompanied by perfect plasticity (Mitchell 2011b), visco-elasticity (Mitchell 2011a) and correspondence models, reviewed in Mitchell et al. (2015).
Work on bond force-strain curves to match global elasto-plastic behavior with hardening is sparse. One study can be found, on state-based PD, where an automated calibration method, together with manually adjusted piecewise linear hardening, is employed to model experimental tensile test of 304 L stainless steel (Littlewood et al. 2012;Littlewood 2015a). However, since the purpose of these studies was to introduce the reader to the PD theory, the bond calibration method was not provided.

The objective and paper organization
The goal of this paper is to establish the use of the EWF method with a PD model that employs nonlinear elasticity for bond deformation, up to maximum bond-load, followed by bond-softening, to represent the deformation theory of plasticity and ductile failure. In this way, the J PD can be used to characterize fracture for brittle or elasto-plastic fracture problems in which small scale yielding assumptions holds, and the EWF for elastoplastic fracture under large scale yielding that can be reasonably approximated by the deformation theory of plasticity.
The paper is organized as follows: we first introduce the bond-based PD, J PD , J R and give a proof for the PD J -integral path independence. In the section thereafter, we describe the EWF and then give the algorithm (procedure) for the PD linear material model, the numerical computation of the J -integrals, and their verification for a number of linear example problems with available analytical solutions for the classical model: an infinite center cracked tension (CCT) specimen, a finite CCT specimen, and a DENT specimen. We then extend the model to include adaptive bond calibration for nonlinear elastic behavior of AHSS grade speci- Fig. 1 Uniformly discretized peridynamic body Ω, showing point x's interaction with one of its neighbors x , within the region H x of radius δ. ξ = x − x is the relative position of the points, η = u − u is the relative displacement and f is the pairwise force mens (therefore implementing the deformation theory of plasticity). An AHSS grade is chosen because of its Poisson's ratio close to 1 /3 fits well to a bond-based constitutive equation. Using the so-calibrated nonlinear model with brittle failure, the EWF is computed for several ligament lengths of DENT specimens. In order to better represent plastic dissipation, gradual softening is added. The last two sections present an in-depth discussion and conclusions.
The major contributions are, therefore, the introduction of a new adaptive bond micromodulus calibration method, and the extension of PD modelling to include the evaluation of the EWF and its validation against experiments.

Bond-based peridynamic theory
The peridynamic equation of motion of the material point at position x at time t is given as (Silling 2000;Silling and Askari 2005) where Ω is the domain of the body, u is the displacement vector field, ρ is the mass density and b is a prescribed body force density. f is the pairwise force function (a vector) per unit volume squared, denoting the force the material point at x exerts on the material point at x. This interaction between pairs of material points is called bond. The integral is defined over a region H x of radius δ, called the horizon, Fig. 1, which can be seen as a sphere, disk or range, for 3D-, 2D-and 1D-models, respectively. A suitable horizon size is chosen and the material body discretized in accordance with problem geometry, loading and desired accuracy of the results (Silling 2016;Xu et al. 2018). Convergence studies may be performed to justify the selection of horizon and discretization. The horizon factor (Bobaru et al. 2009), m = δ/Δx, where Δx is the uniform grid spacing, should in a plane square lattice arrangement have a ratio of at least 3 (Silling and Askari 2005; Madenci and Oterkus 2014) and in many cases 4 or higher (Ha and Bobaru 2010;Henke and Shanbhag 2014;Dipasquale et al. 2016), to provide grid independent crack growth patterns. Studying the effect of increasing m is a socalled m-convergence study (Bobaru et al. 2009). A material is called microelastic if the pairwise force between material points is derivable from a micropotential ω (Silling 2000): where ξ = x − x is the relative position of two material points of the reference configuration, and η = u x , t −u (x, t) = u −u is the corresponding relative displacement of the deformed configuration. A linear microelastic material results in the micropotential where (ξ +η)/||ξ +η|| = e, is a unit vector along a line through the two points of a bond in the deformed configuration. As we assume that a material point x does not interact with material points outside its horizon, f = 0 for ||ξ || > δ. The particular kernel of the integrand in Eq.

Micromodulus
The elastic stiffness of a bond is determined by the micromodulus function c(||ξ ||), which is found by calibrating the peridynamic strain energy density against the classical strain energy density, for a homogeneous body under a) isotropic (dilatational) deformation and b) pure shear (distortional) deformation. The micropotential ω is the energy of a single bond with dimension 'energy per unit volume squared'. The strain energy density of a single point is therefore The factor 1 /2 appears as the points for a bond shares the bond energy between them equally. For a 2D body where t is the thickness of the body, and c 1 comes from assuming a constant micromodulus c(||ξ ||) = c 1 . Isotropic deformation and plane stress give the classical strain energy Setting W = W 0 yields the corresponding micromodulus: The isotropic and pure shear deformations must result in the same c 1 , which for 2D plane stress restricts Poisson's ratio to 1 /3, and to 1 /4 for 2D plane strain and 3D models (Silling 2000;Gerstle et al. 2005). This is because the forces within a bond depends only on the two material points of a bond (and no other points). This restriction has been overcome for state-based PD ). However, an advantage of bondbased PD is that it is computationally less expensive.
Other types of micromoduli (triangular, conical, higher degree polynominal) are derived in a similar fashion as the constant one and are available for 1D (Bobaru et al. 2009), 2D (Ha and Bobaru 2010) and 3D (Silling and Askari 2005). Notice that these micromoduli are calculated assuming a continuous body, i.e. an infinite number of points inside the horizon. For a finite number of points in this region, the micromodulus is, depending upon the type, more or less dependent upon the number of points there; see Eriksson and Stenström (2020a, b) for 1D and (Stenström and Eriksson (2021), App. A) for 2D. This dependence must usually be taken into account.
Since we are using relatively low grid density in the present study, approximate values of c(m = 3) and c(m = 5) have been derived (see appendix) and are c m=3 = 1.1507c ∞ and c m=5 = 1.0697c ∞ . An alternative approach for correcting c, with a similar effect is to apply PD surface effect correction procedures to the bonds of the entire material body.

Peridynamic surface effect
Each material point x interacts with its neighboring points, the family of x, within a radius δ. Therefore, material points at a distance smaller than δ from a free surface will experience a truncated horizon, resulting in a smaller domain of integration. This PD surface effect (Ha and Bobaru 2011) leads to a softer material and larger strains near boundaries. Various correction methods exist, and have conveniently been reviewed and compared by Le and Bobaru (2018).
In this study, 'the energy method' method by Madenci and Oterkus (2014) will be used for PD surface effect correction (Le and Bobaru 2018, p. 507).

Constraint conditions and stresses
The PD equation of motion is a nonlinear integrodifferential equation in time and space, without spatial derivatives, and therefore does not involve any classical boundary constraints. Instead volume constraints can be applied through a nonzero volume rather than on a surface (Silling and Askari 2005;Du et al. 2012), commonly introduced within a layer of thickness Δx or δ (Hu et al. 2012).
The treatment of stresses at boundaries also differ from that of classical continuum mechanics, for that boundary stresses are introduced to PD as body forces b within a layer Δx (Madenci and Oterkus 2014) or δ (Silling and Askari 2005). In practice, to prescribe a stress σ to a material point x, one divides σ with the grid spacing Δx, to obtain an equivalent force per unit volume.

Partial volume scheme
Numerical integration of a material point x over its region H x including the entire volume of each material point x within the horizon δ (Fig. 20), results in an integration area somewhat larger than the area of the discshaped (2D) integration region. A correction factor was introduced, commonly varying linearly between 1 and 1 /2 (Parks et al. 2010), for elements intersected by the horizon: where r = Δx/2, and a uniform cubic lattice is assumed. Equation (8), which is called the LAMMPS method, will be used in this study. Other correction methods exist and can be found in Seleson (2014) and Seleson and Littlewood (2018).

Evaluation of the peridynamic equation of motion
For meshfree numerical implementation (Silling and Askari 2005), the integral for the PD equation of motion, Eq. (1), is replaced by a summation (Niazi et al. 2021): where Fam(i) is the family of nodes j of the horizon region of node i. For dynamic or quasi-static simulation, the common implementation in practice is by time integration through a finite number of small time or load steps and looping over the equilibrium equations of the material points. Detailed algorithms are found in Littlewood (2015b) and Madenci and Oterkus (2014), and a supportive flowchart can be found in Javili et al. (2019). For static modeling, matrix solutions are found in Bobaru et al. (2009 and Eriksson and Stenström (2020a, b). For further reading, reviews are offered by Javili et al. (2019) and Diehl et al. (2019). For peridynamic modeling of Mode I loading, which is used in the present study, see Diehl et al. (2016) for investigation on convergence and crack initiation. For the state-based peridynamic theory, the reader is referred to Silling et al. (2007).

The bond-based nonlocal peridynamic J-integral ( J PD )
J is a contour integral evaluated counterclockwise along an arbitrary path ∂R, enclosing the tip of a straight through crack along the x-axis. The bond-based peridynamic J , derived by Hu et al. (2012), is given as: where W is the strain energy density, Eq. (6), n 1 is the x-component of the outward unit normal vector of ∂R, ds is an element of arc length along ∂R, R 1 is the region inside the contour, and R 2 is the region outside the contour (see Fig. 2). J PD was derived using the infinitesimal virtual crack extension method. The identity of Eqs. (12) and (18) of Hu et al. (2012) can likewise be derived from the left and right sides of the energy balance Eq. (31) of Silling and Lehoucq (2010).

Path independence
So far the J PD has been shown through numerical calculations to be path independent within an error of a Fig. 2 Illustration of the J PD contour, here made rectangular, enclosing a crack of length a. Ω is the domain of the body and includes the regions R 0, 1, 2 and the contour ∂R. R 0 + R 1 = R few percents (Hu et al. 2012;Panchadhara and Gordon 2016). However, like the original, path independence can be derived analytically.
Let R be a subregion enclosing a crack tip of a plane body Ω, in static equilibrium. Path independence is proved by showing that J PD vanishes for any regular region, i.e. a region not containing a crack.
We start from the first integral in Eq. (10), apply the divergence theorem and substitute for the energy density W , Eq. (6): Here, the integral Ω f · ∂u ∂ x 1 dV x vanishes because of static equilibrium. Thus, where the integration boundary of the inner integral is different from that of the corresponding part of the J PD expression. To adjust the integration boundary and the integrand accordingly, Eq. (12) is expanded to where Ω \ R is Ω but not R. The identity Eq. (31), (see appendix), is invoked in the first integral, Eq. (13) then becomes 1 2 Thus, the left side of Eq. (11) equals Eq. (14), which implies that J PD vanishes for any region R not containing a crack tip. This statement proves path independence of J PD . J R and discretization For J R formulated as a function of displacement derivatives, and the discretization of both the two J formulations, see appendix.

The EWF method
The EWF method follows the approach developed by Cotterell and Reddel (1977) and Marchal and Delannay (1996). Performing a test series on geometrically similar specimens but of different ligament lengths, allows the energy spent in the fracture process to be split into two terms. The energy dissipated from the FPZ is called essential work of fracture (W e ), and the energy spent at the region outside the FPZ is called non-essential plastic work (W p ). A DENT specimen is commonly used for thin sheets as its symmetry prevents buckling and ensures mode I crack propagation during the whole test. The total work of fracture where L is the DENT ligament length, i.e. the distance between the two edge notches, and t is the specimen thickness. W e is the specific essential work of fracture per unit area. Following this definition W e is the fracture toughness, equivalent to the J -integral value. W p is the specific non-essential plastic work per unit volume. β is a dimensionless shape factor, depending on the shape of the plastic zone. Normalizing Eq. (15) by the cross-sectional area, allows determination of the specific total work of fracture W f : By tensile testing, with load F as a function of the displacement s, W f is obtained. The ligament length L is determined using, e.g., light optical microscopy.
The value W f is calculated by integration of the loaddisplacement curve: where s f is displacement at fracture.
By plotting W f against the ligament length L, a straight line is obtained and its intercept at the limit L → 0 is taken as the W e . A schematic representation of the evaluation of the EWF method is shown in Fig. 3. However, some restrictions must be met for the validity of the EWF method. The entire ligament must be in a state of plane stress and completely yielded before crack initiation. Cotterell and Reddel (1977) found that if the shortest ligament length is three to five times the sheet thickness and the longest ligament length not greater than one third of the width, W , of the specimen, or two times the radius of the plastic zone, r p , in plane stress: valid results are obtained, and W e ≈ J c . This relationship has been studied for example by Rink et al. (2014).

Algorithm for the linear PD material model and verification of J PD
The algorithm for the linear PD model and evaluation of J PD and J * R is as follows: 1. Specify input parameters: -Material body width and height -Material properties ν and E -Grid spacing Δx and horizon factor m -Remote loading σ 0 -Integration contour number, e.g. contour 1 is the contour through the four points next to the crack tip 2. Discretize the material body into spatial coordinates x and y.  36) and (37) The J expressions are going around the crack tip, from one crack face to the other. Therefore, to avoid peridynamic boundary effects where J crosses the x-axis, i.e. a symmetry line, the peridynamic models include the material body of both sides of the crack.
In the set-up of the overall PD model, we apply a constant micromodulus c as per Eq. (7) (see appendix), PD surface effect correction by using the energy method (Le and Bobaru 2018, p. 507), constraints and stresses within a layer of thickness Δx (Sect. 2.3), volume correction as per Parks et al. (2010), and adaptive dynamic relaxation of Madenci and Oterkus (2014) for obtaining quasi-static solutions.
A constant c gives a linearly increasing bond force with increasing strain, and thus a linear elastic material model. Therefore, for modeling non-linear materials and EWF, c will be calibrated accordingly.
The numerical implementation, including EWF, is done using Matlab software.

Verification of J PD
The J PD and J R integrals are given in Eqs. (36) and (37). We start out the modeling with the infinite CCT specimen.

Infinite CCT specimen
For the classical exact analytical stress-strain-displacement solution of the infinite geometry problem, see appendix. From the remote loading σ 0 of the exact solution, we can calculate the boundary stresses using Eq. (38a), of the specimen shown in Fig. 4a. The boundary stresses are then introduced as input to the PD model. The J methods can now be evaluated for the exact solution using Eq. (38c) and on the modeled PD solution' displacement field.

Problem set-up
To facilitate comparison with previous results of others, e.g. Hu et al. (2012), Eriksson (2019, 2021), we will use a Young's modulus of 72 GPa and a Poisson's ratio of 1 /3 (plane stress) for modeling. The dimensions of the specimen is set to 10 cm by 10 cm (W × 2h), the half crack length a = 5 cm and σ 0 = 1 MPa.
Our reference value for the problem is the classical analytical J 0 = K 2

Peridynamic model discretization
We will keep the horizon factor m constant at three, i.e. δ = 3Δx, and decrease the horizon δ for δconvergence. The error of the PD model is expected to decrease as the number of material points of the model is increased and the PD surface effect diminishes. Note that for the examples shown in Sects. 6-8 below, we use a grid factor m of 4 and 5 for the DENT and EWF examples to achieve a better accuracy for the energy flow. In certain problems, especially for cracks with curving or branching paths, using a larger m-value may be warranted (see Dipasquale et al. 2016). In the examples shown in this work, the crack paths are simple: straight. For such problems, as long as the grid is aligned with the propagation direction, using a large m-value is generally not needed.

Evaluation of J
The resulting J values are shown in Fig. 5. The relative difference is calculated as (J − J 0 )/J 0 . Using the finer δ = 0.6 mm, the difference is less than one percent for both J methods, on both the analytical solution and the PD solution. The difference of the coarser δ is due to PD surface effect and treatment of corner points, and a reason for the difference at the finer discretization is numerical error.

Finite CCT specimen
The finite CCT specimen is loaded at the upper and lower edges by a constant σ y as shown in Fig. 4b. For estimating the stress intensity factor, Isida's (1971) classical solution is commonly regarded as the most accurate one (ASM 1996): where the dimensionless geometry function f 1 can be found in standard books, e.g. (Anderson 2005), i.e. different from the components of PD f. The problem set-up and discretization is the same as for the previous infinite CCT specimen. Thus, a W and h W both becomes 0.5, giving f 1 = 1.967 (Isida 1971), given σ y = σ 0 = 1 MPa. The relative difference is then calculated based on the reference value J 0 = K 2 I /E = 8.44 N/m.
The results are presented in Fig. 6 and show differences of about one percent or less at the finer δ. We notice that the difference at the coarser δ is less than of the infinite model results (Fig. 5).

Half DENT specimen
In a similar fashion as for the finite CCT specimen, we test the J methods on a half DENT specimen, Fig. 4c. Our classical reference solution is that of Wu and Carlsson (1991): The problem set-up differs from the previous case only in height; h = 0.2 m, leading to h/W = 2 and f 2 (2) = 1.169 (Wu and Carlsson 1991). With σ y = σ 0 = 1 MPa, the relative difference is based on the reference value J 0 = K 2 I /E = 2.981 N/m. Since the verification on the CCT models shown good agreement with their classical analytical J values, we use only one discretization here; δ = 3 mm, leading to 100×400 materials points (W × h). This δ gives the differences e J R = 2.2% and e J PD = 2.7%. The δ is the same as for the previous 100 2 models, but we have less PD surface effect in the y-direction due to 400 material points. Thus, the differences are similar to the differences at the 100 2 and 250 2 models in Figs. 5 and 6.

Adaptive calibration of PD model to a given nonlinear stress-strain hardening curve
For evaluating EWF, the PD model needs to be calibrated for non-linear materials. The EWF method is going to be applied to a PD model calibrated to an experimentally obtained stress-strain response. Therefore, an adaptive calibration is carried out with data of a standard (defect free) tensile test specimen. The calibration is carried out on the bond force level so as to obtain a global stress-strain behavior similar to the experimentally measured stress-strain response. Once calibrated, a DENT specimen and EWF is modeled. The calibration procedure includes five constants; c cal , c Δ , c t ,ê c and c E . c cal is a vector of calibration constants of unit interval. Each constant of c cal corresponds to a local bond strain segment (subinterval) of vector c Δ . Since multiplying each bond with a calibration factor is computationally intensive, calibration of local bonds and checking of the resulting global error e c = σ PD /σ 0 , are carried out at load steps evenly divisible with c t . σ PD is the PD global stress and σ 0 is the experimental stress.ê c is the error limit when a new calibration constant is calculated and added to c cal . The fifth and optional constant, c E , is a global strain value between the true elastic limit and the proportionality limit; this is the strain where the calibration starts. (1) as a summation -At each load step, add v 0 to the outermost material points of the specimen's short edges (strips with thickness of Δx or mΔx) -At load steps that are evenly divisible with c t : Calculate the bond force strain of all bonds and multiply each bond force with its corresponding c cal factor. Also, calculate σ PD through the gauge coordinates (Fig. 7). -At load steps that are evenly divisible with c t : If the global strain is larger than c E and the model error e c >ê c , calculate the next calibration constant c i+1 cal = c i cal /e c that corresponds to the strain subinterval c E + (i + 1)c Δ (c 0 cal = 1, i = 0, 1, 2, ...). The obtained c cal vector can now be used in the general four steps algorithm (Sect. 5) by calibrating all bonds at given bond strains.
σ PD is the stress through a point interface by summing the x-components along m points (Fig. 7): where the first summation are the m points along a line along the stress direction, and Fam + (i) are the family members of point i that are located in front of the m points dividing line. Multiplication with the area of the In the above algorithm, a calibration constant is added when e c >ê c . Alternatively, the calibration constant can be added at certain (evenly distributed) load steps, bond strains or global strains. However, we have found that using aê c limit gives a PD stress-strain curve with less oscillations and better necking compared to the other criteria.

Problem set-up and discretization
For model calibration and evaluation, we will use experimental data of a bainite-martensite sheet metal (AHSS) provided by Golling et al. (2019). The engineering stress-strain from Fig. 4 of Golling et al. is reproduced in Fig. 8a. The stress-strain data is obtained using a standard specimen of 12.5 mm in width and 2.0 mm in thickness (ISO 2019; ASTM 2016), and gives an E 0 of 200 GPa. Thus, the PD model constants E and ν are set to 200 GPa and 1 /3, respectively, and the problem geometry used for the calibration is set to 12.5 mm × 62.5 mm.
The geometry is discretized into 40×200 material points, with m = 5, thus δ = 5 · 0.3125 mm. The opposite ends of the specimen is loaded with a constant velocity v 0 = 1.6 · 10 −7 m/s, and a virtual strain gauge is set 6.4 mm (20.5Δx) from the lengthwise mid. The loading rate is chosen small enough to avoid specimen oscillations.
The PD recovered E modulus before calibration is 0.8 % below E 0 , measured at the strain gauge, and 0.1 % Fig. 8 a Calibrated PD global stress-strain and experimental engineering stress-strain (Golling et al. 2019). b Corresponding local PD bond force. c c cal versus c Δ . d Displacement at a point for steady check below when measured as a mean across the specimen width. The approach used to find the stress, and thus E, is presented in Fig. 7 and Eq. (21).

Adaptive calibration
The adaptive calibration is performed during a single specimen loading. The calibration procedure begins at c E = 0.003 global strain and the first calibration constant is of unity. The global stress is checked against the experimentally obtained stress-strain response for each c t = 5 steps. As the specimen is loaded, a new calibration constant is supplied each time that the global stress error reaches the limitê c = 2%. Calibration of all bonds are carried out for each c t = 5 steps.
The calibrated bond force-strain model and global PD stress-strain response are shown in Fig. 8. The measurement is carried out at the virtual strain gauge point and the bond force-strain is of one bond in the global load direction, Fig. 8a, b. c cal and c Δ amounts to 107 calibration constants/parameters. The total number of load steps, at the strain where the standard sheet metal specimen fractured, is 10 540.
Since material bonds reach their respective yield strain at different times/loads, the global stress-strain response is "smoothed". Not least for a bi-linear law, this leads to an apparent hardening and/or necking effect (Silling 2016, p. 44). The noise in Fig. 8b comes from oscillations due the quasi-static model and numerical approximation.

Recovered stress-strain curve
Since the calibration factors are adaptively supplied one at a time during the calibration, based on monitoring at one gauge point, the recovered stress-strain curve will differ when all factors (c cal and c Δ ) are supplied initially. The recovered stress-strain with a loading rate of v 0 = 0.2 · 10 −7 m/s is shown in Fig. 9. 'No dmg' is without bond break and 'With dmg' is with bond break at the last calibration constant of 0.0465 bond strain. The green cross in Fig. 9a is were it fractured. s 1 in Fig. 9b indicates the failure bond law used here and in the next section.
To notice is that the calibration values are more or less independent of the number of material points and grid density m, as seen in Fig. 10.

EWF on full DENT specimen
We found from the literature review that bond force calibration studies are presently going on and that EWF modeling for PD is new. Therefore, this preliminary EWF study is limited regarding influencing parameters and applications. Moreover, the EWF modeling includes estimation of the J -integral and to avoid generation of duplicate results, we will apply the classical J R formulation, Eq. (37), only, in particular because of the small number of points of narrow ligaments.  The full DENT specimen of Golling et al. (2019) is 45 mm wide and 2.0 mm thick. The corresponding tensile test was monitored with an extensometer with a gauge length of 50 mm. In the PD model, we use a geometry of 45 mm × 225 mm (2W × 2h; see Fig. 4d). A virtual gauge for logging the displacement is set 25 mm from the two notches, i.e. 50 mm gauge length. For logging the force, another gauge is set 10 mm from the two notches. This is to avoid a delay/bump of the force-elongation response that appears if the force is logged at 25 mm. The grid density m is set to 5. Discretization and ligament lengths (L) are given in Table 1. The loading rate v 0 is 0.4 · 10 −7 m/s for δ = 2.8, 1.9 and 1.4, and 0.6 · 10 −7 m/s for δ = 1.1.
We will use two different crack tip definitions that we call sharp and fuzzy tips; see explanation in Fig. 12.

Evaluation of EWF and J
The adaptive calibration of Fig. 8c is used for the full DENT specimen without further adjustments. PD bonds break at the bond strain s 1 = 0.0465 and drops suddenly to zero, Fig. 9b, and the PD model is of a nonlinear-elastic type. The simulation result for δ = 1.1 mm is shown in Fig. 11a. Force is the global cross sectional force. The dashed lines indicate the amount of damage, or fractured ligament, between the two notches. 100% damage means fully fractured ligament. For calculation of W e , the upper integration limit is set to 90% damage; see the dotted vertical lines.
Integration of the force-elongation curves using Eq. (17) gives W e . See Fig. 11b. Experimental data (crosses) of Golling et al. (2019) gives W e = 78±6 kN/m while PD gives W e = 182 kN/m and lack plastic work.  Fig. 11c. Here, the vertical lines of the maximum global force have been set manually and is therefore somewhat subjective.
For studying δ-convergence, we complement with three more discretizations. δ-convergence is performed keeping the grid density m constant while decreasing the horizon radius δ, i.e. increasing the total number of material points of the model (Bobaru et al. 2009). See Fig. 15. By increasing the number of material points in the model, the PD surface effect is reduced. Only J c is traced as W e lack plastic work.
Lastly, damage maps are presented in Figs. 13 and 14.

Softening
PD materials with a linear bond force-strain relationship contains a single parameter (critical bond strain) to simulate both the strength and toughness. For a cer-tain horizon size, one is then able to match both, but the likelihood is that the particular horizon is too small to make it practical for computations. In order to allow matching of both strength and toughness for any horizon size, one needs to introduce a two-parameter bond force-strain model. For crack nucleation, propagation and strongly nonlinear materials, bi-and trilinear constitutive laws have been developed (see Sect. 1.3). PD trilinear laws for brittle and quasi-brittle materials have been studied by Yang et al. (2018) and Niazi et al. (2021). The adaptively calibrated constitutive law of the present study, resulted in a piecewise linear approximation of the experimental stress-strain response. In the preset model, for a more pronounced softening part, we also assume a linear behavior for the bond-response, for simplicity. See Fig. 16. Softening starts at s 1 and continues linearly to s 2 .
We keep s 1 = 0.0465 and manually calibrate s 2 = 2 for capturing the plastic work. We decrease the horizon δ to 0.28 mm, and for computational efficiency set m = 4. With the added softening, the PD model is still able to recover the experiment tensile specimen (Fig. 17). For the PD DENT, v 0 = 0.2 · 10 −7 m/s, bond calibration every c t = 5, sample size is 45 mm × 63 mm for computational efficiency, and the global force gauge is shifted from 10 mm to 5 mm from the mid for better force response. Lastly, we set the ligament lengths closer to the experimental, to 6.61, 8.72, 11.25 and 12.80 mm (this was not possible before with δ = 2.8 mm). The resulting DENT model is shown in Fig. 18. Ligament damage, Fig. 18a, is the fraction of bonds that have passed s 1 . The PD DENT model shows improved plastic work but appears too strong, i.e. the bainitic-martensitic material appears too brittle. We therefore try a more ductile material of Golling et al. (2019) called 'Bainite'; see calibrated and recovered standard tensile test in Fig. 17.
The values s 1 and s 2 are calibrated manually to 0.03 and 1.2, respectively. The loading speed is lowered to v 0 = 0.1 · 10 −7 m/s, and for computational efficiency, the horizon δ is increased to 1.4 mm, with m = 5. The ligament lengths are set as close as possible to the experimental ones.
The resulting PD DENT model now better match the maximum global force and the EWF-line (Fig. 19). PD obtained W e = 369 kN/m, experimental W e = 349 kN/m ≈ J c , and PD obtained J * R = 365 kN/m. If the force-elongation curves are integrated to the crack initiation, an alternative toughness measure is obtained, called fracture toughness at crack initiation (W * e ), and is the average of the resulting EWF data points, with the slope assumed zero ). In Fig. 19b, W * e = (128+154+170+182+195) /5 · 10 3 = 166 kN/m.
The bainite standard tensile test response has a narrow bend from the elastic region to the hardening phase (Fig. 17). This caused slightly more noise to the calibration constants compared to the constants of the bainite-martensite material (Fig. 10). It did not affect the recovered tensile test, but caused more oscillations in the adaptive dynamic relaxation solver. To remove the noise, a curve was fitted to the calibration constants. Curves tested includes polynomial, two terms exponential, Weibull, three terms Gaussian and two terms power function. Several fittings gave an accurate recovered tensile test, but a tenth order polynomial gave the least amount of noise for the DENT, and where therefore chosen. By curve fitting, the number of calibration constants can be chosen. 30, 100 and 180 constants were tested. More constants produced smoother DENT  Golling et al. (2019). For the s 1 values used, softening occurred at the cross symbols test and therefore the polynomial function was supplied directly to the PD code instead of calibration constants. However, 30 constants did work as well.

Discussion
The martensitic-bainitic material turned out to have too low ductility for the PD material model. Experimental data showed high ultimate strength (1730 MPa) and sudden fracture of DENT specimens, at the point where the DENT experimental curves ends in Fig. 18a. This is due to microscale effects discussed later.
The bainitic PD material model was able to match both the standard tensile test specimen strength and the DENT specimens strength. By studying the softening part, Fig. 19a, the PD DENT curves meet at a single point (beyond 1 mm elongation) while the experimental curves ends at different elongations, e.g. 0.7, 0.8 and 0.9 mm. To obtain similar EWF, we integrated the PD DENT curves to a remaining strength of 6.1 kN. The reason to stop at this value is the difference between the actual last stages of the ductile failure and the deformation theory (nonlinear elasticity) approximation of it in our model. Note that if one chooses a lower force level (on the softening side) to integrate to, the resulting EWF line (in Fig. 18b) only slightly shifts downward. Changing the S 1 or S 2 values affects both the maximum force and ductility of the DENT.
Since we are using the deformation theory (nonlinear elasticity) plasticity instead of an incremental plasticity model, the part where unloading happens (monotonic loading is no longer satisfied) will not be very accurate. DENT softening is a challenge experienced for FEM as well (Sandin et al. 2021), but due to reasons inherent to FEM.
The disturbance of the PD DENT force elongation responses seen in Fig. 19a) comes from bonds entering the linear softening (discontinuous function). The linear softening can be replaced by a curved softening (bell-shaped), to better match the curvature of the experimental softening part. Here we use linear for computational simplicity.
The PD DENT responses makes slightly sharper bend from the elastic-to-hardening region compared to the experimental results (Fig. 19a). This can be due to microscale effects.

Microscale effects
Microscale properties that can have an influence on the result are: Different materials The experimental DENT specimens have fatigue cracks and is therefore in a different damage phase compared to the PD DENT that does not have any fatigued bonds. Also, there are different failure mechanisms taking place; the standard tensile specimen used for calibration is polished (Golling et al. 2019) and has a fracture zone with a radius approaching infinity while the DENT specimen has a crack tip approaching zero in radius.
Microstructure The PD model lacks a microstructure. In the experimental DENT specimens, following Crack tip shielding The crack tip of the experimental DENT specimens is under compression when the loading starts due to asperities mismatch, thus, leading to higher strain at fracture compared to the PD model.
Discretization The horizon radius δ is much larger than an actual crack tip radius; thus, the not enough refined horizon leads to a blunt PD crack tip.
Surface correction The PD model includes a surface correction factor, which increase the surface stiffness but leave some stiffness fluctuation.
Strain hardening The experimental DENT speci-men may experience strain hardening, and therefore a strengthening effect.

J -integral
Before adding softening, the J -integral of the PD DENT model converges to the experimental J , when J is taken at maximum global force (Fig. 15).
With added softening to the bainite-martensite material, J takes of about where the softening starts, before the maximum force. J shows an increase, or elbow, around 100-150 kN/m, Fig. 18c. This elbow becomes more pronounced and closer to the experimental J when δ is shrunken. At δ = 0.14 mm, the elbow is around 100 kN/m, but was not studied further due to computational time. An elbow may also be distinguished in Fig. 11 at 120-140 kN/m, which explains the slower convergence of the model with softening. The elbow can also be seen by plotting J versus bond damage.
With the bainite material and softening, at the absolute maximum force, the mean J * Rc of 365 kN is close to the value obtained by Golling et al. (2019) W e = 349 kN/m ≈ J c . An elbow cannot be seen.

Experimental data
The experimental DENT results are somewhat uncertain. The narrow temperature-time process window of bainite formation under quenching makes it difficult to obtain specimens with small variation in phase volume fractions (Golling et al. 2016). Another source of uncertainty is that the crack length of the experimental DENT was measured using a stereomicroscope. Recommendation is, e.g. by ASTM, to heat tint the specimen surface and measure the crack length after a fracture test, as the actual crack tip is often not visible in a stereomicroscope and the crack tip front is convex.

Conclusions
A peridynamic (PD) material model calibrated for highly nonlinear elastic behavior has been adopted to investigate, computationally, the use of the EWF method in characterizing fracture toughness of ductile materials.
The J PD -integral was shown analytically to be path independent, in analogy with its classical counterpart.
First, we compared the classical and peridynamic forms of the J integral value on a number of example problems. On the classical analytical displacement solution of the infinite center cracked tension specimen, the relative differences between the J values were less than one percent compared to the classical analytical J 0 value (Fig. 5).
Thereafter, to match nonlinear constitutive behavior, the PD model was calibrated using a novel auto-mated procedure. Under small scale yielding assumptions, the classical J -integral, applied on the PD double edge notched tension (DENT) model, converged to the experimentally obtained J value (Fig. 15).
For capturing the energy dissipated during plastic deformation, a linear bond softening model was adopted (Fig. 16). The model was employed for two very different materials; a lower-ductility bainiticmartensitic steel and a higher-ductility bainitic steel. Up to the start of the softening phase, the model recovers excellently the material response, of both materials (Fig. 17).
The lower-ductility material has an ultimate strength of 1730 MPa, while the higher-ductility material has a strength of 970 MPa. A similar large difference in strength is not seen in the experimental DENT response due to high sensitivity to defects of the lower-ductility material. This sensitivity to defects is seen in the results, as we could not recover the specific essential work of fracture (W e ) value of the lower-ductility martensiticbainitic steel with the nonlinear elastic model (Fig. 18).
While we only used a nonlinear elastic model (equivalent to the deformation theory of plasticity), we were able to match very well the experimentally obtained W e value for the higher-ductility bainitic steel (Fig. 19). Also, the J -integral value obtained from the PD model, at the absolute maximum specimen load, matched the corresponding W e value.
While being aware of the model limitations, future work will investigate the use of incremental plasticity models on this topic. to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

The classical J-integral as a function of displacement derivatives
For a rectangular integration path, Rice's J -contour integral on strain components formulation can be written as (Stenström and Eriksson 2019): On the left side is the general path independent formulation, where W for a linear elastic material is 1 2 σ i j i j , and σ i j n j is (the components of) the traction vector. On the right side, the factor 2 accounts for the rectangular contour to be symmetric across the straight crack. I, II and III are the right, top and left sides, respectively, of the contour above the symmetry line (Fig. 2). By substitution of stress-displacements derivatives relationships of W for a linear elastic material, we obtain: w i corresponds to the second term of the general J expression. With similar substitution of stressdisplacements derivatives, and by taking the contour as a rectangle around the crack tip, the second term of the general J expression develops into three parts; right, top and left hand side contours: Thus, w I , w II and w III are calculated on the right, top and left sides of the contour, respectively.

Discretization of the J-integral expressions
The material domain is discretized uniformly to allow for a mid-point integration scheme (one-point Gaussian quadrature) (Silling and Askari 2005). The strain components of W Eqs. (10) and (32) can be approximated by applying the central difference scheme on the displacement field. For a material point i, the strains are given as follows (Hu et al. 2012): The J -integrals are approximated by using the trapezoidal rule. The J PD of Eq. (10) is then discretized as follows: n ∂R , n R1 and n R2 are the number of material points within the contour path, inner region and outer region, respectively. The subscripts 1 and 2 refer to the xand y-components, t is the material thickness and A j/i refers to the area. The classical J of Eq. (32) is discretized as: where n i is the number of material points along the three regions I-III.

Exact analytical solution of the CCT specimen
Classical exact analytical solutions of stresses and displacements for the infinite center cracked tension (CCT) specimen have been derived by Unger et al. (1983) under plane strain conditions. We have in an analogue manner derived the stresses and displacements under plane stress conditions (Stenström and Eriksson 2019): where X , Y and A-C are dimensionless variables, α and β are constants for plane stress and plane strain, respectively. The dimensionless variables are given as α and β are given by: (3 − ν)/(1 + ν) for plane stress 3 − 4ν for plane strain (40c) x and y in Eqs. (39a) and (39b) are the spatial coordinates. With Eq. (38) the stresses and displacements can be calculated at any point of the infinite plate. The benefit of this exact solution is that the stresses or displacements can be used as input to peridynamic modeling, or more precisely, as boundary conditions for finite models of the infinite geometry problem.