The Use of Peridynamic Virtual Fibres to Simulate Yielding and Brittle Fracture

The forces in the ‘arms’ joining the particles in a peridynamic analysis depend upon the state of stress in the equivalent continuum and the orientation, length and density of the arms. Short and long arms carry less force than medium length arms as controlled by the weighting kernel. We introduce an intermediate step of imagining a mat of long fibres in which the fibre forces only depend upon the stress, the fibre orientation and the length of fibres per unit volume without the added complexity of the arm lengths. The effect of the arm lengths can then be considered as a separate exercise, which does not involve the continuum properties. The arm length is proportional to size of the particles and the separation of length from the state of stress allows for modelling of variable particle density in the discretisation of a problem domain, which enables computationally efficient accurate analysis. We then introduce the concept of arm elongation to fracture in order to model surface energy in fracture mechanics. This means that shorter arms have a larger strain to fracture than longer arms. Numerical implementation demonstrates that this produces a fracture stress that is inversely proportional to the square root of the crack length as predicted by the Griffith theory [1, 2].


Backgound
In this section, we give a brief historical overview of the development of peridynamics, as far as possible using the notation of previous authors. However, in Section 2, which follows this background, we assume that the reader is not already familiar with this theoretical domain and more comprehensive definitions and derivations are presented. A summary of the symbols, variable names and nomenclature that is used in the paper is compiled in Appendix 1, Table 2.
Peridynamics is a non-local continuum theory that was developed for the simulation of fracture phenomenon for brittle materials and was introduced by Silling in 2000 [3]. As opposed to other techniques, peridynamics does not require cracks to be predefined but rather to appear in a spontaneous fashion from an initial intact domain. Unlike classical continuum mechanics which rely on the evaluation of partial derivatives, peridynamics works by replacing the differentiation with integration which remains valid also in the presence of discontinuities. The peridynamics continuum is approximated with a finite set of material points with intermediate bonds. Each material point is connected by bonds to the neighbouring material points within a distance which defines a horizon H which is circular in 2D and spherical in 3D. The force in a bond is related to the bond deformation and the choice of constitutive model. The peridynamics theory was developed on the basis of regular material point distribution which may lead to computationally costly simulations. Additions have been made to enable variable material point distribution for homogeneous materials in [4,5] and [6] with the motivation to reduce computational cost (see Fig. 1). However, these additions complicates the theory and an argument is presented here for a potentially simpler approach where the bond force is based on the theory of Smooth Particle Hydrodynamics (SPH).
The equation of motion of a material point in classical continuum mechanics according to Newton's second law of motion is, in which the gradient of the stress tensor ∇ is replaced in the peridynamic formulation with an integration of force density, which does not rely on assumptions of spatial continuity, and thus enables simulation of phenomena such as fracture. Here is the material density, ü( , t) is the acceleration of a material point at position at time t, is a bond force density acting between two particles at positions and ′ , dV x ′ is the volume associated with particle ′ and is the external loading.
(1) ( )ü( , t) = ∇ + ( , t), The calculation of the force density varies with the different peridynamic formulations, which are typically classified in terms of bond based, ordinary state-based and non-ordinary state-based theory as described in [7]. Silling introduced the bond based theory in [3], and following the notation in [7], where and represents initial and displaced position of a material point, the bond force is formulated as where ′ is the displaced position of a particle with initial position ′ . The constitutive model is introduced with the micro modulus c, which is calculated from, where E is the Young's modulus, is the radius of the horizon H x . The bond stretch s in Eq. (3) is calculated from, Ordinary state-based peridynamics (OSB-PD) was introduced by Silling et. al in 2007 [8] and can be understood as a generalisation of the bond based theory where the limitation of fixed Poisson's ratio is overcome. The OSB-PD builds on the concept of a state which can be described as a mapping. For the application in solid mechanics the deformation state and the force state are the most central concepts for which definitions and derivation are given in [8]. The equivalent expression for the force density in Eq. (1) as outlined for bondbased theory in Eq. (3) but instead for the OSB-PD theory is according to [7] given as, where ⟨⟩ is the force state acting on the vector x � − x within the angle brackets. For a linear elastic isotropic material the force state is according to [7] given by, where a, b and d are peridynamic parameters and ( , t) is the dilatation term.
Smooth particle hydrodynamics is another meshless methods that has be used for similar problems. It was introduced in 1977 by Monaghan et al [9] and the first use of SPH for the simulation of solids was carried out by Libersky and Petschek in 1991 [10]. Just like with peridynamics, the SPH continuum is discretised with a set of particles and each particle is influenced by all neighbours within reach of a kernel function. The similarities between peridynamics and SPH is explored by Ganzenmuller et al in [11]. Following the notation in [12] the equation for the motion of a particle in a continuum with the SPH formulation can be written as, where is the material density, v i is the i th component of the velocity, ij is the stress tensor, x j is the j th cartesian component of the position vector and g i denotes the i th component of the body force per unite mass. This expression is developed in [12] Eq.(3.1) into, where W ab is the kernel for a bond between particles a and b, m b is the mass of particle b, Π ab and ij are terms related to viscosity. The summation acts over all particles b that are within reach of the kernel function, i.e the SPH-horizon. If we simplify Eq. (9) by ignoring body load and the viscosity term (which is mainly relevant for fluid problems) and we multiply both sides with the mass of particle a and focus the attention on the bond between the a and b, the i th component of the force ab acting on particle a can be interpreted as The same bond generates a force of opposite direction and magnitude on particle b, such that Equation (10) seems to lend it self for variable size of particles a and b since the masses, stress, and densities are separate entities. However, the kernel function W ab − , h which in [13] is defined as a function of , and h is the same for both particles and need to be separated into W ab and W ba , allowing for different values of the particle size parameter h.

Introduction to the Fibre Model
In meshless methods such as aforementioned smoothed particle hydrodynamics (SPH) [9,13,14] and peridynamics [3], a fluid or solid continuum is approximated by a system of particles in which each particle is joined by arms to its near neighbours. The forces in the arms are used to model the state of stress in the continuum and it is usual to assume only axial tensile and compressive forces in the arms, although in non-ordinary state based peridynamics [8] non-axial forces in the arms are allowed, effectively producing bending moments in the arms. We will limit our discussion to the case when there are only axial forces.
The left-hand image in Fig. 2 shows a collection of particles and arms. The particles are arranged randomly, but with a limit on the minimum spacing. The arms are shown grayscale with the darkness proportional to the weighting ab in Eq. (41).
The axial force in a particular arm depends upon: 1. the local state of stress in the continuum that we are modelling, 2. the orientation of the arm relative to the state of stress, 3. the size of the particles at its ends and 4. the length of the arm. Usually, longer arms and shorter arms will carry less force than those somewhere in the middle.
The reasons for criterion 4 are that we want the forces to die away smoothly as we consider further away particles and large forces in short arms cause problems due to the change in force as an arm rotates due to transverse displacements. Or, to put it another way, the geometric stiffness of an arm is the tension divided by the length, which is large for short members, unless the tension is small. A compression produces a negative geometric stiffness. In the following, we shall use the word 'tension' to mean tension or compression, if its value is negative.
The assumption is that we have a sufficiently large number of particles and arms to model a continuum with sufficient accuracy. The problem is compounded by the fact that the system of arms and particles is highly statically indeterminate and many different distributions of arm tension will be equivalent to the same stress. We therefore have to make some further assumptions to make the problem determinate. However, if all the arms are in tension, or all the arms are in compression then the sum of the product of the arm forces times their lengths is independent of the arm arrangement [15,16]. This can be easily shown using virtual work with a virtual uniform dilatation.
In SPH and peridynamics [11] it is usual to derive expressions for the arm tensions by focusing attention on the divergence of the stress tensor, that is how the stress changes with position and hence causes acceleration of the continuum. This approach produces some mathematical complexity which is beyond most engineers who are happier with the ideas behind the finite element as first described by Turner et al. [17].
We will therefore use simpler ideas based on resolution of forces, which we can do if we separate consideration of criteria 1 and 2 listed above from consideration of criteria 3 and 4. The approximation usually introduced in finding the derivative of functions in SPH and peridynamics is replaced by the approximation in replacing summations by integrals. However, the resulting expressions are identical, at least for simple states of stress corresponding to elasticity, fluid pressure or viscosity. Regardless of how arm tensions are derived in SPH and peridynamics, it is always assumed that the change of stress is small over distances comparable with the particle spacing. We will therefore consider the situation of a continuum under stress that varies with direction, but not position. But rather than consider a system of particles and arms, we will begin by considering a mat of long random continuous fibres in 2D as shown in the central image in Fig. 2, or the three-dimensional equivalent in which the fibres are randomly arranged in 3D.
Even though the fibres are arranged randomly, we assume that on average the number within a cone of a certain solid angle is the same in all directions.
We assume that there are a large number of fibres, sufficient to model a constant state of stress. At this stage we make no assumption as to what causes the stress. We could use the fibres to represent an elastic or plastic solid or a fluid with pressure and viscosity. This theory is therefore equally applicable to SPH and peridynamics.
We will be discussing stress in 2 and 3 dimensions, and it is important to remember that stress has the units force per unit length in 2D and force per unit area in 3D.
Our model is similar to a fibre reinforced composite, but without the matrix. We do not need a matrix because we are assuming that the stress does not vary with position, only with orientation and that the fibres do not stop and restart. See Hashin [18] for a discussion of composites. When we introduce particles later, the stress will then be able to vary with position, exactly like any other pin jointed framework.
The force in a fibre in a particular direction will depend upon the state of stress and the orientation of the fibre relative to the stress, that is criteria 1 and 2 above. It will also depend upon the number of fibres, as expressed by the total length of fibres per unit area in 2D or volume in 3D. Let us imagine that we double the number of fibres. The total length of fibres will also be doubled, but the force in each fibre will be halved. Thus the sum of the length of each fibre times the force it is carrying per unit area in 2D or volume in 3D will be unchanged. Force times length per unit volume or area has the same units as stress in 3D and 2D respectively.
The arrows in Fig. 1 refer to the fact that a knowledge of the peridynamic arm forces implies a unique set of fibre tensions which implies a unique stress in the continuum. But the reverse does not apply. The peridynamic arms forces depend upon the kernel which gives arm forces dependent upon the arm length. In the case of the fibre mat, different variations of fibre tension with orientation can give rise to the same state of stress in the equivalent continuum.

Force Flux Density
Let us introduce a quantity which we shall call the force flux density, S, which is a function of the state of stress and the fibre orientation relative to that state of stress. In 2D it is given by where T is the tension in each fibre in a particular direction and L is the total length of the fibres per unit area. In 3D it is where L is now the total length of the fibres per unit volume. L has unit one over length in 2D and one over length squared in 3D. Therefore, in both cases, S has the dimensions of stress. In 2D, the fibres contained within the angle from to + d would be considered to be the same as the fibres contained in the angle from + to + + d . Therefore all the fibres would be contained within an angle of radians. Similarly, all the fibres in 3D would be contained within a hemisphere, that is a solid angle of 2 steradians. Thus, in general the tension in a fibre is given by Note that S and T will in general vary with direction, whereas L is not associated with a direction. We will see that the reason for choosing the constant in this way is so that when S is isotropic, then its value is equal to the mean stress.
If we want to calculate S for a given state of stress, the problem is statically indeterminate, as we have already noted, and we would have to make certain assumptions. The reverse procedure requires no assumptions since a given S as a function of direction corresponds to a unique state of stress, which we can find by resolving an area into a particular direction and then resolving the corresponding force. In 2 dimensions, this double resolving produces a product of cosines and/ or sines in exactly the same way as in producing the Mohr's circle for stress [19] to give Where , and are unit vectors in the directions of the coordinate axes and indicates the dyadic product, also known as the outer product or tensor product of and , and sometime written as ⊗ . In 3 dimensions, the resolving is slightly more complicated, where is the angle from the z direction in spherical polar coordinates. The reason for the sin between the d and the d is that d sin d is an element of solid angle, that is surface area on the unit sphere. Thus, for example, taking the and components from both sides, The force crossing a cut in 3D specified by the vector whose magnitude is equal to the area of the cut and whose direction is normal to the cut is given by We first resolve in the direction specified by and and then resolve the force so produced into the coordinate directions. The integral is required to sum the forces from the fibres in different directions. The force flux density S associates a scalar with each direction in space, but it is clearly more than a vector or tensor field since the values of the scalar associated with each direction are completely independent. When we come to consider plasticity we will find that some directions will be yielding while other directions are still elastic.

Isotropic Stress
If S is independent of orientation in Eq. (15) we obtain in 2D because the average value of cos 2 and sin 2 is 1 2 , whereas the average value of cos sin is zero. If S is independent of orientation in Eq. (16), in 3D. Thus, in both 2 and 3D, we find that the mean stress is equal to S, when S is isotropic. Even when S is not isotropic, the mean value of S is equal to the mean stress, which follows from the use of virtual work with an isotropic virtual expansion, as noted above.

Uniaxial Stress
When we have a uniaxial stress the force flux density must vary with direction. In 2D let us consider a uniaxial stress in the x direction, corresponding to = 0 . The simplest force flux density distribution to choose is This satisfies S( ) = S( + ) , which is a requirement since each fibre points in two opposite directions. In addition, the mean force flux density is corresponding to half the isotropic value and the stress in the x direction is as required. Eq. (21) tells us that S = 3 2 when = 0 and S = − 2 when = 2 . This means that if we treat the fibres as simply elastic bars the Poisson's ratio is equal to 1 3 in twodimensional plane stress.
We can repeat this process in 3D, now with a uniaxial stress in the z direction, because then the mean force flux density, corresponding to one third the isotropic value and the stress in the z direction is as required. Equation (24) tells us that S = 2 when = 0 and S = − 2 when = 2 . This means that if we treat the fibres as simply elastic bars the Poisson's ratio is equal to 1 4 in 3D, as we would have expected [3].

General Stress State
We know that we can generate any stress state by superimposing the two or three orthogonal uniaxial principal stresses. Thus if is a unit vector we can generalise Eq. (21) to give

Calculation of Equivalent Peridynamic Inter-Particle Arm Tensions
We can modify Eq. (35) to find the tensions in arms joining particles for SPH or peridynamics, thus introducing criteria 3 and 4 that were introduced in Section 2. We can write the tension in the arm joining particles a and b as where ab is a non-dimensional weighting function which will depend upon the length of the arm and the size of the particles a and b. Again, is the deviatoric stress and ̄ is the mean stress.
is the unit vector from a to b. Particle a is currently located at a and is the actual length of the arm. The weighted length of the arm is, and L is now the total weighted length of the arms per unit area or unit volume, If we multiply all the weightings by the same constant, it makes no difference to T ab since T ab depends upon the ratio ab L .

Choice of Arm Length Weighting Function
There is no 'correct' arm weighting function since any weighting function should in principal be able to model a constant stress. However some functions will work better than others when it comes to numerical work. In choosing a weighting function, we want to be able to estimate the the value of L from Eq. (40) by replacing the summation by an integral and we want to be able to produce an equation for T ab which is consistent with the results from conventional SPH and peridynamics. Let us assume we have a large number of randomly but evenly distributed particles of variable mass. The mass of particle a is m a . The average density is with the units of mass per unit area in 2D and mass per unit volume in 3D.
We will will write the weighting of the arm ab as where C is a constant with the units of length and N = 2 or N = 3 is the number of dimensions. As noted above the value of C does not influence T ab , but we have included it to make Eq. (41) non-dimensional. We have chosen this rather complicated expression for ab because we shall see that [13], except we assume that the 'size' of the particle b as represented by h b may vary from particle to particle. Equation (38) of Ganzenmüller et al. [11], in which their X ij is our r ab , relates W() to the function ⟨ ⟩ used in peridynamics [8].
ing upon how we define h b . Thus r ab W r ab h a , h a is negative and R ab is positive, as we would expect. Figure 3 shows a typical kernel, and note that the length r ab is always taken as positive.
Clearly particles of larger mass should have a larger dimension, so we could write If we do this r ab h b will almost certainly not be equal to 1.0 where the kernel becomes equal to 0.
W is a scalar function and in writing W , h b we are assuming that W is radially symmetric and depends only on the ratio and h b . It is possible to have kernels which are not radially symmetric, but that would require additional information, for example a symmetric second order tensor to define the dimensions and orientations of the axes of an ellipsoid. Following Monaghan [13] will write may not be equal to W ba . To remember which is which is the weighting of particle b at particle a. In 3D, the kernel has the property where is the distance from particle b to a typical point . Therefore in 3D. The upper limit of integration is taken as infinity, but in practice kernels will always be of finite size and be equal to zero beyond a certain horizon.
where the non-dimensional function The reason for 1 , rather than the 1 h N b that one might have expected is that We can see that Eq. (55) is dimensionally consistent with the units of length on both sides.

Total Weighted Length of Arms Per Joining Particles
Let us associate only the first part of R ab with particle a and the second part with b. If we add the contributions for all the particles to which a is connected we obtain the sum of the weighted lengths associated with a, This is the total weighted length associated with particle a and to find the total weighted length per unit area or per unit volume as appropriate we need to multiply by m a to give Thus, finally, which is independent of whatever value we have chosen for the constant length C.

Constant Stress and Density
The formula for the tension in the arm ab is given by Eq. (36) and upon substituting in Eq. (63), Thus the total force applied to particle a divided by its mass is Following Monaghan [13] we can write the gradient so that In the case of an hydrostatic pressure P then S = −P so that

Variable Stress and Density
In Eq. (67), we assumed constant stress and density. But in general the stress and density will vary, but slowly compared to the particle spacing. Let us write a for the density associated with particle a and S a ab for the force flux density associated with particle a in the direction of the unit vector ab . Remember that reversing the direction of does not change S( ) . So let us now write The corresponding arm tension would be In the case of isotropic pressure P a at a and P b at b, Finally, if all the particles have the same size, so that h a = h b , and therefore ignoring any influence of how variable density modifies the size of particles, then W ab = W ba producing This is identical to Equation (3.3) of Monaghan [13].

Linear Elasticity
Let us assume that the fibres are linear elastic and we know the strain tensor , which is assumed to be small from the unloaded state. The mean strain, and the deviatoric strain, ∇ a W ba . (70) where G is the shear modulus of the equivalent isotropic continuum and the 2 is there because the shear modulus is defined in terms of the 'engineering' shear strain rather than the 'mathematical' shear strain. K is the bulk modulus and it appears as NK because K is defined using the volumetric strain, which is N̄ in N dimensions. In terms of Young's modulus E and Poisson's ratio , In 3D, for two-dimensional plane stress, and for two-dimensional plane strain, in which the strain perpendicular to the plane is zero, In 3D the elastic constant appearing in Eq. (77), which is zero when = 1 4 , as expected from the above. Similarly for plane stress which is zero when = 1 3 , again as expected. For plane strain, which is zero when = 1 4 , as was the case in 3D. Thus we can control the Poisson's ratio by including the effect of ̄ upon the arm tensions, in exactly the same way as described in Section 15 of Silling's 2000 paper [3]. .

Newtonian Fluid
If we replace the strain by the strain rate ̇ , the shear modulus G by the viscosity and the bulk modulus K by the bulk viscosity in the elastic Eq. (77), and include the thermodynamic pressure P then we obtain where N is again the number of dimensions. If is the velocity of particle a, then for the particle pair a and b is equal to the rate of strain in the arm. This appears in equation (4.2) of Monaghan [13].
Clearly one has to be careful if r ab happens to be very small, but W ab r ab should tend to zero as r ab → 0 , meaning that the force is finite. Even so Monaghan includes the factor 2 in the denominator in his equation (4.2). An alternative would be to use a kernel with a flat top so that it has zero curvature at r ab = 0.

An Elastic-Plastic Material
It is relatively easy to model an elastic solid or viscous fluid with arm tensions. However, it is not easily possible to model a material with a given yield criterion such as the well-known Tresca and von Mises criteria or the recently introduced Bigoni and Andrea Piccolroaz criterion [20]. On the other hand it is easy to specify a yield condition in terms of peridynamic arm tensions and compressions, and then see what this means in terms of stress and a yield criterion. Let us imagine that the tension in an arm depends on 1. the volumetric strain of the particles at its ends which produces an elastic isotropic stress via the bulk modulus and 2. an elastic-perfectly plastic elongation of the arm acting on its own. The arms will not all begin to yield at the same time so that the transition from elastic to plastic will be gradual. However, let us assume that all the arms connected to a particle are yielding, either in tension or compression, such that the principal plastic strain increments are d plastic I , d plastic II and d plastic III in 3D satisfy (85) S( ) = −P + (N + 2) ⋅̇ ⋅ + (N − (N + 2))̄̇, d plastic I + d plastic II + d plastic III = 0, so that the plastic volumetric strain is zero. We can impose this requirement simply via the elastic isotropic stress. If we align the principal strain increments with the coordinate axes, we have the tensor so that, in the direction of the unit vector the strain increment is which is positive if the arm is yielding in tension and negative if it is yielding in compression. The boundary between the tension and compression zones is where For simplicity, let us assume that all the arms have the same yield strength, either in tension or compression so that the force flux density is equal to ±S y per unit solid angle where the constant ±S y is the yield value of S. The principal stresses are aligned with the principal strain increments and therefore, ignoring the isotropic stress from the elastic bulk modulus, (92) 3 cos boundary − cos 3 boundary cos 2 d , 3 cos boundary − cos 3 boundary sin 2 d , cos 3 boundary d .
These equations will probably have to be integrated numerically and this has been done to plot the yield surface on the -plane [21] in Fig. 4. The arrows indicate the corresponding plastic strain increments, and it can be seen that they are normal to the surface and so obey the normality condition, as one would expect. Figure 5 shows the regions on the sphere yielding in tension and compression.
However, we can perform the integration for the two simple special cases. For the case of simple tension or compression, to give (95) 3 cos boundary − cos 3 boundary cos 2 d 3 cos boundary − cos 3 boundary sin 2 d cos boundary − cos 3 boundary d .

Fracture Mechanics
Dimensional analysis tells us that no amount of adjustment of stress/strain relationships will enable us to model brittle fracture because it is not possible to produce a dimensionless group containing the crack length. The Griffith theory [1,2,22] uses the concept of surface energy, that is the work required to increase the area of a crack with units of J/ m 2 . Surface energy is conventionally given the symbol , which should not be confused with the deviatoric strain or its components. Youngs modulus E and the failure stress away from the crack, failure have units of N/m 2 . Thus if c is the length of an edge crack in a material (not to be confused with the micromodulus in Section 1), the only independent non-dimensional groups are Griffith used interchange of elastic and surface energy to show that so that However we do not need to postulate the relationship between surface energy and failure stress, we only need to be able to model surface energy. If we assume arm tensions are controlled by yield, then we need to specify an elongation to failure to give the surface energy. The corresponding strain in an arm will depend upon its length, with shorter arms needing higher strains. This is familiar from necking in a tensile test on a ductile metal where the failure strain depends upon the gauge length since all the elongation is in the neck region. Thus effectively we are imagining that the arms fail by necking.
Let us write the surface energy as where is a constant length which is related to the amount of deformation associated with the surface energy. Then (103) (105) = yield tension , Eq. (105) should be compared with the concept of the modulus of cohesion in Barenblatt [22]. Writing where yield tension is the yield strain or proof strain, where the non-dimensional quantity Thus we can focus our attention on which will be less than 1 in any situation involving brittle fracture. If = 1 and yield ∕tension = 1 500 , then c = 1000 , which gives an indication of the sort of elongation we require of the inter-particle arms.

Numerical Simulation of Brittle Fracture
Let us replace the constitutive equation for linear elasticity in Eq. (77) by where ab is the total engineering strain calculated using Pythagoras' theorem to give the current arm length from which we subtract the initial length and then divide by the initial length.
plastic ab is the plastic strain. We will continue to calculate ̄a from ab so that the plastic volumetric strain is zero.
We now use the algorithm where yield is the yield strain which is taken as positive.
10 Numerical Implementation

The Kernel
The choice of kernel W r b h b , h b is arbitrary, provided that it satisfies Eq. (53) and Eq. (54), although some kernels will work better than others. We have chosen (107) yield ∕tension = E yield ∕tension, (108) failure ≈ yield ∕tension, to satisfy Eq. (54). The integer k and the non-dimensional constant are a matter of choice and are the same for all particles, regardless of their size. If k → ∞ and → ∞ we obtain the bell shaped Gaussian kernel. We have taken k = 2,

Model Setup
The numerical implementation that has been used to validate the theory is a simple elliptical disk in 2D with a crack spanning between the two focal points. The model is generated using elliptical coordinates for variables ∈ (− ∕2, ∕2) and ∈ (0, ) according to The disk is stressed by an evenly distributed prescribed displacement at the boundary under the assumption of two-dimensional plane strain, which means no displacement in the z direction perpendicular to the x − y plane and displacements in the x and y directions which are independent of z.
This means that the strains z , yz and zx are all zero, the stresses yz and zx are both zero, but z is nonzero due to Poisson's ratio effects in the elastic region. In the plastic region, stresses are larger than the yield stress can maintained since z can adjust itself to limit the deviatoric stress.
The geometrical parameters that control the setup of the base geometry are shown in Tables 1 and 2, where n refers to the number of divisions in the direction and n refers to the number of division in the direction. The grid shown in Fig. 6 is used to define zones  in which the particles are stored. The connectivity between the zones is computed and utilised every time when particle neighbour calculations are necessary. The setup of particles and arms follows a four-step process after the underlying grid has been generated. This is illustrated with the four quadrants a) -d) in Fig. 7. The first step is to populate all zones with constant number of particles, resulting in an increasing particle density closer to the focal point of the ellipse as a result of the parametrisation of Eq. (114). For the example that follows, each zone is populated with a pattern of 5x5 particles. A Monte Carlo inspired algorithm is then used to introduce noise in the particle distribution as shown in Fig. 7b and Fig. 7c, according to the procedure described in Fig. 8.
After the re-positioning of the particles, the particle-zone topology and the particle size need to be updated. The new size is calculated from, The setup for the numerical example. The stress on the plate is introduced through an incrementally imposed displacement field on the particles at the perimeter (in the grey zones). The direction of the displacement is illustrated with the arrows where h a is the particle size, a are the coordinates of current particle, b are the coordinates of a neighbouring particle, and n is the number of neighbours which is set to n = 5 . In a perfect even hexagonal-based circle packing, each circle has 6 adjacent neighbours. In a suboptimal case with noise and lower density packing it seemed reasonable to assume a lower number, thus the choice to set n = 5.
From the particle size and the choice of kernel parameter , the horizon for each particle a can be computed and the results for a few particles are illustrated in Fig. 7c. The horizon H a is the distance of connectivity and will be unique for each particle due to the irregularity of the distribution. Finally, the arms between the particles can be created by connecting each particle a with the neighbours withing its horizon H a = h a , as illustrated in Fig. 7d and Fig. 9. Note that the zone size should be chosen such that , Fig. 7 Showing the process of setting up the numerical model, starting from a regular particle distribution in a), in which noise is introduce using a Monte Carlo inspired algorithm as shown in b). Resulting particle distribution without the underlying grid of zones is shown in c) where some particles are highlighted for which the horizon ( h a ) is drawn as circles. The last quadrant d) shows the resulting arms from connecting each particle with all the other particles within its horizon the shortest side of the zone is larger than the horizon of the largest particle within that zone.
In Table 1 nP is the number of particles, nA is the number of arms, and k are constants that are used in the kernel Eq. (112) and Eq. (113), c is half the crack length which is specified in Eq. (114) and n is the number of neighbours used to calculate the particle size from Eq. (115). Furthermore, in Table 1 is the plastic elongation limit, E is the Young's modulus, y is the yield strain, y is the yield stress, is the Poission's ratio, is the material density, G and K are the shear and bulk modulus calculated from Eq. (78) and Eq. (81), respectively.
yes End no Calculate max neighbour overlap Generate random vector such that Start Move the particle Update particle size from Eq.(104) Fig. 8 Flowchart algorithm for the introduction of noise in an initially regular particle distribution, exemplified for a particle a with neighbours b. The particle size is updated only in each fifth iteration

Analytical Solution & Applied Load
The numerical implementation seeks to simulate progressive fracture and the results are compared to Griffith's theory formulated in Eq. (103). The load that is applied to trigger fracture is a prescribed incremental displacement of the particles in the boundary zone. The elastic solution to the analytical equation as outlined in Eq. (116) is used to impose this prescribed displacement filed. The displacements are then converted to a boundary stress for comparison with Eq. (103). This is done by summing over the reaction forces of the particles in the boundary zone which results from the imposed displacements, and dividing that force with the area of the disk perimeter, where the thickness is set to 1.
The analytical solution is also used to study convergence of the elastic solution with increasing particle count as presented in Fig. 10. A derivation of the analytical solution is out of the scope of the paper but the interested reader is referred to p.(187 -204) in [23]. For the case of plain strain, the analytical solution for the elliptical plate with an Fig. 9 All arms that connect neighbouring particles are shown in the left figure. The variation in density comes as a result of the varying particle size where particles are smaller closer to the two crack tips and larger towards the boundary. The largest particle has a radius which is more than 400 times greater that the smallest one. A zoomed in view of the sub region which is marked with a rectangle in the left figure and concentrated around the predefined crack as shown in the figure to the right Fig. 10 Showing a convergence study where the root mean square (RMS) for the displacement for the elastic solution of the theory presented here is compared with the analytical solution from Eq. (116). The numerical simulations are run for 5 different particle densities and compared to the analytical solution for the highest density (39200 particles) elliptical crack of zero width between the focal points can expressed in terms of a complex potential as , v is the displacement in y-direction, and referrers to curvilinear coordinates as described in [23] pg.192 and our Eq. (114), G is the shear modulus for plain strain and is the Poisson's ratio of the material.

Fracture
In order to validate the fracture model the relationship between stress and surface energy was examined numerically and compared with with Griffith's prediction, Eq. (103). It does not matter whether we have an edge crack or an internal crack, as we have studied, since we are primarily concerned with the fact that we expect the failure stress to be proportional to the square root of the surface energy divided by the crack length, and are less concerned with the constant of proportionality. Following the set up of the geometry as described above, the plastic elongation limit was varied from c/500 to c/3250 with incremental steps of c/250. The dimensionless proportional relationship, between the stress to failure over the yield stress and the square root of the elongation limit over half the crack width c, is expected to be linear. In order to calculate particle damage a damage number at the level of the arm is introduced such that The particle damage is then computed from where n a is the number of arms within the horizon H a for particle a.

Results
The results from the simulations are first shown for the case when = c∕1000 in Figs. 11 and 12, followed by a compilation of 11 simulations with varying . The load is increased by small increments until progressive fracture occurs such that the precise level of fracture stress could be obtained. The simulation is aborted when the crack propagation has reached a distance of 10c from the centre of the disk, effectively leaving the boundary zone intact. The elongation of the arms is shown in Fig. 11, where the total elongation in a) is computed L ab and the plastic elongation limit is shown in b). Figure 12 shows the particle damage in a) which is calculated from Eq. (119) and the total arm strain ab over the yield strain y in b).
The results of 11 simulations with a fixed irregular particle distribution and varying plastic elongation limit can be seen in Figs. 13 and 14. Figure 14 shows the a comparison with the linear relationship predicted by the Griffith criterion. Note that one can only establish the slope of the graph, that is the constant in the Griffith criterion, by doing numerical experiments. The high values for the fracture stress in Fig. 14  All arms including broken ones are drawn. A relative colour scheme is applied in a) where the most elongated arm is red and the least elongated arm in blue. The colour scheme in b) relates the plastic arm elongation to the elongation limit , where arms that exceed the elongation limit are coloured red Fig. 12 Drawing the damage of the particles in a) with a colour scale form green to red where ≤ 0.5 is red and = 1 is green (no damaged arms). In b) the arm yield strain is drawn in a colour scheme where ab ∕ y > 1 is red and ab ∕ y = 0 is blue Fig. 13 Showing ab ∕ y for the arms as a result of different values for . ab ∕ y > 1 is red and ab ∕ y = 0 is blue. In the case of = c∕500 the whole plate yields before fracture occurs are as a result of the plain strain assumption, as discussed above. The arm yield strain for the 11 simulations (including a 12th attempt for = c∕500 where fracture did not occur as a limiting case) are shown in Fig. 13.

Conclusion
The separation of the influence of the fibre orientation and the length of fibres per unit volume from the arm length between particles simplifies the peridynamic theory. It enables us to simulate a homogeneous material with varying particle sizes which can be useful in the study of stress concentrations and fracture propagation. It also enables us to produce a criterion for the yield of a peridynamic continuum which has the form of the Tresca hexagonal prism, but with rounded corners.
We introduce the concept of arm elongation to fracture in order to model surface energy in fracture mechanics. Numerical implementation demonstrates that the fracture stress is inversely proportional to the square root of the crack length as predicted by the Griffith theory. Table 2 Nomenclature and abbreviations table which is split into sections to follow the structure of the paper. Some of the variables in one section may be given another meaning in another section. This is most evident with the variables in Section 1 since the ambitions have been to use typical notations for each theoretical domain, which is causing some overlap  Unit tensor a, b Current and neighbouring particle notation T ab Tension in the arm between particle a and b w ab Non-dimensional weighting function for the arm ab The total weighted length associated with particle a and b a Total force applied to particle a [N] ∇ a Gradient at particle ā ,̄̇The mean strain and the mean strain rate ,̇

Nomenclature abbreviations and symbols
The Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.