Modelling of Granular Fracture in Polycrystalline Materials Using Ordinary State-Based Peridynamics

An ordinary state-based peridynamic formulation is developed to analyse cubic polycrystalline materials for the first time in the literature. This new approach has the advantage that no constraint condition is imposed on material constants as opposed to bond-based peridynamic theory. The formulation is validated by first considering static analyses and comparing the displacement fields obtained from the finite element method and ordinary state-based peridynamics. Then, dynamic analysis is performed to investigate the effect of grain boundary strength, crystal size, and discretization size on fracture behaviour and fracture morphology.


Introduction
Polycrystalline materials are widely used in many different industrial applications. Amongst the various existing polycrystalline materials, metals and ceramics are common examples. Polycrystalline materials are composed of individual crystals that have a particular crystal orientation and are separated from neighbouring crystals via grain boundaries. Microscopic features of polycrystalline materials such as crystal orientation, grain boundary strength, etc. may have a significant effect on the overall macroscopic behaviour of the material, especially on the fracturing behaviour of these materials. Therefore, it is essential to analyse this type of material at the microscopic scale. Although experimental approaches can be useful for this purpose [1][2][3], it is not always possible to perform such experiments and they can also be very expensive. On the other hand, numerical approaches can be a very good alternative. There are various numerical studies available in the literature. For the majority of these studies, cohesive zone elements (CZM) [4][5][6], extended finite element methodology (XFEM) [7,8], and boundary element method (BEM) [9,10] were utilized. Although these techniques are widely used and powerful, they may have limitations for some specific cases and conditions. Instead, a new continuum mechanics formulation, peridynamics (PD) [11], can be considered.
As opposed to partial differential equations that traditional approaches are based on, peridynamics utilizes integro-differential equations without containing any spatial derivatives. Hence, these equations are always applicable regardless of discontinuities such as cracks. Peridynamics has been used for the fracture analysis of many different types of materials and material behaviours [12][13][14][15][16][17][18][19]. It has also been applied for the analysis of polycrystalline materials [20][21][22]. However, these studies used either original bond-based formulation (BB) [11] or non-ordinary state-based (NOSB) [23] formulations. Bond-based formulation has limitations on material constants whereas non-ordinary state-based formulations may encounter the zero-energy mode problem. In order to overcome all these issues, an ordinary state-based (OSB) peridynamic formulation [23,24] can be utilized. The numerical solution of peridynamics is generally obtained by using a meshless scheme. Therefore, the formulation does not where ρ(x) represents the density and ..
u(x, t) is the acceleration of material point x at time t. Moreover, t(u − u, x − x, t) and t (u − u, x − x, t) denote the force density vectors of the material points x and x , and u − u represents the difference of displacements of the material points x and x at time t. In Equation (1), H represents the peridynamic horizon that defines the range of interaction of a particular material point, as shown in Figure 1. In the peridynamic literature, the size of the horizon is usually represented by the symbol δ.
Materials 2016, 9,977 2 of 24 numerical solution of peridynamics is generally obtained by using a meshless scheme. Therefore, the formulation does not suffer from issues such as mesh distortion, especially for problems involving large deformations and fractures. There are also other mesh-free methods available in the literature and used for modelling shear bands in metals [25][26][27][28]: concrete fragmentation [29][30][31], dynamic fracture in thin shells [26,32], and fluid-structure interaction [33]. A good example of a meshless approach used for fracture is the cracking particles method (CPM) [29,34,35]. In CPM, the crack path is represented by a set of cracked particles. CPM is especially useful for complex fracture patterns such as crack branching and coalescence. Other promising techniques used for fracture modelling include smoothed particle hydrodynamics [36][37][38], molecular dynamics [39], the discrete element method [40], and the force potential-based particle method [41]. Although the aforementioned approaches may have certain advantages for particular conditions, in this study peridynamics was chosen for modelling granular fracture in polycrystalline materials. Moreover, the ordinary state-based formulation was utilized for the first time in the literature in order to overcome the limitations of bond-based formulation and zero-energy mode problem of non-ordinary state-based theory. After validating the formulation, several demonstration cases are considered to investigate the effect of grain boundary strength, crystal (grain) orientation, and grain size.

Ordinary State-Based Peridynamic Formulation for a Cubic Crystal
The equation of motion of Ordinary State-Based (OSB) peridynamics (PD) can be written as: , (1) where ( )   [22], the ordinary state-based model for a cubic crystal can be represented using two types of interactions (bonds), as shown in Figure 2.  [22], the ordinary state-based model for a cubic crystal can be represented using two types of interactions (bonds), as shown in Figure 2. These are: where  represents the angle between the orientation of the bond and the crystal (grain) orientation. As an example, bonds within the horizon of a particular material point for a grain orientation of = are shown in Figure 2.
According to OSB PD theory, the strain energy density of a material point can be written as [24]: where J is the total number of material points within the family of material point   k x .
By using the strain energy density expression given in Equation (2), the peridynamic force densities t and ' t can be obtained as: where with and These are: 1. Type 1 bonds (green dashed lines)-interactions along all directions ( φ = 0 ∼ 2π), 2.
Type 2 bonds (red solid lines)-interactions along the directions of φ = 1 4 π, 3 4 π, 5 4 where φ represents the angle between the orientation of the bond and the crystal (grain) orientation.
As an example, bonds within the horizon of a particular material point for a grain orientation of ϕ = π 4 are shown in Figure 2. According to OSB PD theory, the strain energy density of a material point can be written as [24]: where J is the total number of material points within the family of material point x (k) . By using the strain energy density expression given in Equation (2), the peridynamic force densities t and t can be obtained as: where with µ T2 = 1 Type-2 bonds 0 otherwise (5) and where Materials 2016, 9, 977 4 of 23 In Equations (3) and (6), y and y represent the location of material points x and x after deformation, i.e., y = x + u and y = x + u (see Figure 1). The PD dilatation, θ, for a crystal can be expressed as: and the parameter, Λ, is defined as: The stretch parameter s can be expressed as: The PD material parameter a is associated with the deformation involving dilatation, θ (k) . The remaining material parameters, b T1 and b T2 , are associated with deformation of the bonds along the Type 1 and Type 2 bond directions, respectively, as shown in Figure 2. All PD material constants can be expressed in terms of material constants of a cubic crystal, Q ij , from classical theory. The procedure for obtaining these relationships is presented in Section 3.
When the stretch, s (k)(j) between material points k and j exceeds a critical stretch value, s c , the interaction breaks and damage occurs. Hence, there will no longer be any interaction between these two particles. The critical stretch parameter (2D) can be expressed as [24]: For this loading condition, the length of the relative position of material points   j y and   k y in the deformed state becomes: The PD dilatation term can be evaluated as: in which   x j    x k   . As expected, there is no dilatation for this loading condition. Therefore, the strain energy density can be calculated as: For this loading condition, the length of the relative position of material points y (j) and y (k) in the deformed state becomes: The PD dilatation term can be evaluated as: in which ξ = x (j) − x (k) . As expected, there is no dilatation for this loading condition. Therefore, the strain energy density can be calculated as: or Equating expressions of strain energy density from CCM and OSB PD formulations, i.e., Equations (15) and (19), results in: 3.2. Second Loading Condition (Uniaxial Stretch in Crystal Orientation Direction: ε 11 = ζ, ε 22 = 0) In the second loading condition, a constant strain is applied along the direction of crystal orientation (Figure 4), and the corresponding dilatation and strain energy density from CCM can be expressed as: and W CCM 3.2. Second Loading Condition (Uniaxial Stretch in Crystal Orientation Direction: 11 22 , 0 In the second loading condition, a constant strain is applied along the direction of crystal orientation (Figure 4), and the corresponding dilatation and strain energy density from CCM can be expressed as: and Due to this deformation, the dilatation of PD can be evaluated as: The length of the relative position of material points y (j) and y (k) in the deformed state becomes: Due to this deformation, the dilatation of PD can be evaluated as: By equating expressions of the dilatation term from CCM and PD formulations, i.e., Equations (21) and (25), results in: The PD strain energy density for this loading condition can be calculated as: or Hence, by equating expressions of strain energy density from Equations (22) and (28), the following relationship can be obtained: 3.3. Third Loading Condition (Biaxial Stretch: In the third loading condition, a constant strain is applied in all directions ( Figure 5), and the corresponding dilatation and strain energy density from CCM can be expressed as: and The length of the relative position under this loading condition can be evaluated as: Hence, the dilatation term in PD formulation can be expressed as: By equating the expressions of dilatation from both CCM and PD, i.e., Equations (30) and (34), the same expression given in Equation (26) can be obtained. Moreover, the PD strain energy density under this loading condition can be evaluated as: By equating Equations (31) and (35), a new relationship can be obtained, as: The length of the relative position under this loading condition can be evaluated as: Hence, the dilatation term in PD formulation can be expressed as: By equating the expressions of dilatation from both CCM and PD, i.e., Equations (30) and (34), the same expression given in Equation (26) can be obtained. Moreover, the PD strain energy density under this loading condition can be evaluated as: By equating Equations (31) and (35), a new relationship can be obtained, as: Hence, the OSB PD material parameters can be expressed in terms of engineering constants of CCM by utilizing the three relationships given in Equation (20), (29), and (36) as: For BB PD, the parameter a associated with dilatation should vanish, leading to the constraint equation Q 12 = Q 44 , which is a limitation of BB PD in cubic polycrystal analysis.

Numerical Results and Discussion
In this section, the results generated from static and dynamic PD analyses are presented, and comparisons with finite element method (FEM) results are also provided. For static analysis (Section 4.1), a single cubic Niobium (Nb) crystal model is considered first (Section 4.2.1) and displacement fields of PD and FEM are compared. Then, a cubic Molybdenum (Mo) polycrystal model with 18 Voronoi grains is analysed (Section 4.2.2), and the PD and FEM displacement fields are compared. For the dynamic analysis, the influence of the discretization size and the interface strength coefficient (β) on the results is considered first (Section 4.3.1). Then, the influence of crystal size on fracture behaviour is investigated (Section 4.3.2).

Material Data
Two different materials are considered in this study: niobium (Nb) for single crystal static analysis, and molybdenum (Mo) for polycrystal static and dynamic analysis. According to [43], the local stiffness matrix of each individual crystal can be written as: However, for plane stress configuration, the material matrix given in Equation (38) can be written by using reduced stiffness matrix following the procedure given in [44] as where Q ij are the reduced stiffness coefficients and can be expressed in terms of the elements of the stiffness matrix, C ij as Therefore, the material properties of Nb and Mo are given in Table 1 as shown below [44]: The fracture toughness of Mo is specified as K IC = 24.2 MPa √ m [45], which corresponds to a critical stretch value of 0.008127 for plane stress configuration.

Static Analysis
A cubic crystal model with a length of 152.4 µm and a width of 76.2 µm is considered and the number of particles along the horizontal and vertical directions is 240 and 120, respectively. The values of grid spacing and horizon are ∆x = 0.635 µm and δ = 1.915 µm, respectively. A uniform discretization scheme is used throughout this study. However, non-uniform discretization can also be possible by using small grid sizes at critical regions such as interfaces and utilizing the "Dual-horizon peridynamics" concept, as introduced by Ren et al. (2016) [46,47]. A horizontal tension loading of, P, is applied on the right edge of the model, and the left edge is fully fixed as shown in Figure 6. The tension loading is specified as a body load and applied to a single layer of material points at the right edge of the model. The displacement constraint condition at the left edge is also imposed to a single layer of material points. To reach the steady-state condition and perform static analysis, an adaptive dynamic relaxation technique was utilized as described in Madenci and Oterkus (2014) [24].
Therefore, the material properties of Nb and Mo are given in Table 1 as shown below [44]:  , respectively. A uniform discretization scheme is used throughout this study. However, non-uniform discretization can also be possible by using small grid sizes at critical regions such as interfaces and utilizing the "Dual-horizon peridynamics" concept, as introduced by Ren et al. (2016) [46,47]. A horizontal tension loading of, P, is applied on the right edge of the model, and the left edge is fully fixed as shown in Figure 6. The tension loading is specified as a body load and applied to a single layer of material points at the right edge of the model. The displacement constraint condition at the left edge is also imposed to a single layer of material points. To reach the steady-state condition and perform static analysis, an adaptive dynamic relaxation technique was utilized as described in Madenci and Oterkus (2014) [24].  Based on the results presented in Figures 7 and 8, good agreement is obtained between PD and FEM analyses. Therefore, it can be concluded that the OSB PD crystal model presented in this study can produce accurate results for different crystal orientations for a single crystal.

Static Analysis of Mo Polycrystal
In the second case, a polycrystal model with 18 randomly orientated grains is generated by using Voronoi tessellation. A uniform discretization is utilized. Depending on the location of the material point, corresponding grain orientation is determined. Hence, Type 2 bonds will exist in many different directions according to the random orientation of the crystals. The average crystal size is 645.16 μm 2 , and the amount of tension loading applied on the right edge is 374.5 MPa P  . The layout of the polycrystal model is shown in Figure 9.  Based on the results presented in Figures 7 and 8, good agreement is obtained between PD and FEM analyses. Therefore, it can be concluded that the OSB PD crystal model presented in this study can produce accurate results for different crystal orientations for a single crystal.

Static Analysis of Mo Polycrystal
In the second case, a polycrystal model with 18 randomly orientated grains is generated by using Voronoi tessellation. A uniform discretization is utilized. Depending on the location of the material point, corresponding grain orientation is determined. Hence, Type 2 bonds will exist in many different directions according to the random orientation of the crystals. The average crystal size is 645.16 µm 2 , and the amount of tension loading applied on the right edge is P = 374.5 MPa. The layout of the polycrystal model is shown in Figure 9. Based on the results presented in Figures 7 and 8, good agreement is obtained between PD and FEM analyses. Therefore, it can be concluded that the OSB PD crystal model presented in this study can produce accurate results for different crystal orientations for a single crystal.

Static Analysis of Mo Polycrystal
In the second case, a polycrystal model with 18 randomly orientated grains is generated by using Voronoi tessellation. A uniform discretization is utilized. Depending on the location of the material point, corresponding grain orientation is determined. Hence, Type 2 bonds will exist in many different directions according to the random orientation of the crystals. The average crystal size is 645.16 μm 2 , and the amount of tension loading applied on the right edge is 374.5 MPa P  . The layout of the polycrystal model is shown in Figure 9.  Figures 10 and 11 and a good agreement is obtained between two solutions, which confirms that the present model can also accurately represent polycrystal material behaviour.

Dynamic Analysis of Mo Polycrystals
For the dynamic analysis, a 5 mm by 5 mm square plate with randomly oriented grains is considered as shown in Figure 12. A horizontal velocity boundary condition of 5 m/s V  is applied on both the left and right edges of the model. Three layers of virtual particles are placed along the left and right edges to impose this condition, as suggested in [24]. A no-fail zone is also imposed on virtual particles and their neighbouring particles in order to allow the load to be accurately transferred inside the plate. Two pre-existing cracks with a length of 0.4 mm are applied at the centre of the bottom and top edges, as shown in Figure 13

Dynamic Analysis of Mo Polycrystals
For the dynamic analysis, a 5 mm by 5 mm square plate with randomly oriented grains is considered as shown in Figure 12. A horizontal velocity boundary condition of V = 5 m/s is applied on both the left and right edges of the model. Three layers of virtual particles are placed along the left and right edges to impose this condition, as suggested in [24]. A no-fail zone is also imposed on virtual particles and their neighbouring particles in order to allow the load to be accurately transferred inside the plate. Two pre-existing cracks with a length of 0.4 mm are applied at the centre of the bottom and top edges, as shown in Figure 13. The time step size is specified as dt = 0.05 ns and the total number of time steps is 100,000, i.e., the total simulation time of 5.0 µs. The study considers three different interface strength coefficient, β, values (0.5, 1.0, and 2.0), three different mesh sizes (74 × 74, 150 × 150, and 300 × 300) and three different total numbers of grains (25, 100, and 400; i.e., different crystal sizes).

Dynamic Analysis of Mo Polycrystals
For the dynamic analysis, a 5 mm by 5 mm square plate with randomly oriented grains is considered as shown in Figure 12. A horizontal velocity boundary condition of 5 m/s V  is applied on both the left and right edges of the model. Three layers of virtual particles are placed along the left and right edges to impose this condition, as suggested in [24]. A no-fail zone is also imposed on virtual particles and their neighbouring particles in order to allow the load to be accurately transferred inside the plate. Two pre-existing cracks with a length of 0.4 mm are applied at the centre of the bottom and top edges, as shown in Figure 13

Effect of PD Discretization Size and Interface Strength Coefficient (β)
The aim of this analysis is to investigate the effect of the peridynamic discretization size on the crack pattern predicted by PD model and the morphology of intergranular and transgranular fracture modes when changing the value of the interface strength coefficient, β. The horizon is specified as 3.015 x    , which means that it is controlled by changing the PD discretization (74 × 74 particles, 150 × 150 particles, and 300 × 300 particles). Moreover, three different interface strength coefficient, β, values are considered to investigate the intergranular and transgranular fracture modes of the polycrystal.

Effect of PD Discretization Size and Interface Strength Coefficient (β)
The aim of this analysis is to investigate the effect of the peridynamic discretization size on the crack pattern predicted by PD model and the morphology of intergranular and transgranular fracture modes when changing the value of the interface strength coefficient, β. The horizon is specified as δ = 3.015 · ∆x, which means that it is controlled by changing the PD discretization (74 × 74 particles, 150 × 150 particles, and 300 × 300 particles). Moreover, three different interface strength coefficient, β, values are considered to investigate the intergranular and transgranular fracture modes of the polycrystal. Figures 14-16 show the fracture pattern of the polycrystal under plane stress configuration at five different times (1.5 µs, 2.0 µs, 2.5 µs, 3.0 µs, and 3.5 µs) for β = 0.5 with 74 × 74 particles, 150 × 150 particles, and 300 × 300 particles, respectively.

Effect of PD Discretization Size and Interface Strength Coefficient (β)
The aim of this analysis is to investigate the effect of the peridynamic discretization size on the crack pattern predicted by PD model and the morphology of intergranular and transgranular fracture modes when changing the value of the interface strength coefficient, β. The horizon is specified as 3.015 x    , which means that it is controlled by changing the PD discretization (74 × 74 particles, 150 × 150 particles, and 300 × 300 particles). Moreover, three different interface strength coefficient, β, values are considered to investigate the intergranular and transgranular fracture modes of the polycrystal.     The results show that with an increasing total number of particles, the intergranular crack pattern can be predicted more accurately and in more detail. However, the simulation time will increase rapidly as well. Therefore, it is important to find a good balance between accuracy and time. In this study, 150 × 150 particles can provide appropriate results, which is the reason why most of the simulations in this paper are chosen by using this number of particles.    The results show that with an increasing total number of particles, the intergranular crack pattern can be predicted more accurately and in more detail. However, the simulation time will increase rapidly as well. Therefore, it is important to find a good balance between accuracy and time. In this study, 150 × 150 particles can provide appropriate results, which is the reason why most of the simulations in this paper are chosen by using this number of particles.   The results show that with an increasing total number of particles, the intergranular crack pattern can be predicted more accurately and in more detail. However, the simulation time will increase rapidly as well. Therefore, it is important to find a good balance between accuracy and time. In this study, 150 × 150 particles can provide appropriate results, which is the reason why most of the simulations in this paper are chosen by using this number of particles. Figures 17-22 show the fracture patterns of the polycrystal at five different times (1.5 µs, 2.0 µs, 2.5 µs, 3.0 µs, and 3.5 µs) for β = 1.0 and β = 2.0 with 74 × 74 particles, 150 × 150 particles, and 300 × 300 particles, respectively.  The results show that with an increasing total number of particles, the intergranular crack pattern can be predicted more accurately and in more detail. However, the simulation time will increase rapidly as well. Therefore, it is important to find a good balance between accuracy and time. In this study, 150 × 150 particles can provide appropriate results, which is the reason why most of the simulations in this paper are chosen by using this number of particles.                     As described above, similar conclusions can be found in these simulations. For instance, branching of cracks can be obtained more clearly by increasing the total number of particles, but the simulations become more time-consuming. Moreover, the transgranular fracture mode becomes more dominant as the interface strength coefficient increases.

Effect of the Crystal Size
The aim of this section is to investigate the effect of the crystal size on fracture pattern. The plate is discretized by 150 × 150 particles, containing three different numbers of randomly orientated grains (25 grains, 100 grains, and 400 grains). Three different grain boundary strength coefficients, 0.5    As described above, similar conclusions can be found in these simulations. For instance, branching of cracks can be obtained more clearly by increasing the total number of particles, but the simulations become more time-consuming. Moreover, the transgranular fracture mode becomes more dominant as the interface strength coefficient increases.

Effect of the Crystal Size
The aim of this section is to investigate the effect of the crystal size on fracture pattern. The plate is discretized by 150 × 150 particles, containing three different numbers of randomly orientated grains (25 grains, 100 grains, and 400 grains). Three different grain boundary strength coefficients, β = 0.5, β = 1.0 and β = 2.0, are considered to investigate the effect of crystal size for different fracture modes. As described above, similar conclusions can be found in these simulations. For instance, branching of cracks can be obtained more clearly by increasing the total number of particles, but the simulations become more time-consuming. Moreover, the transgranular fracture mode becomes more dominant as the interface strength coefficient increases.

Effect of the Crystal Size
The aim of this section is to investigate the effect of the crystal size on fracture pattern. The plate is discretized by 150 × 150 particles, containing three different numbers of randomly orientated grains (25 grains, 100 grains, and 400 grains). Three different grain boundary strength coefficients, 0.5    As described above, similar conclusions can be found in these simulations. For instance, branching of cracks can be obtained more clearly by increasing the total number of particles, but the simulations become more time-consuming. Moreover, the transgranular fracture mode becomes more dominant as the interface strength coefficient increases.

Effect of the Crystal Size
The aim of this section is to investigate the effect of the crystal size on fracture pattern. The plate is discretized by 150 × 150 particles, containing three different numbers of randomly orientated grains (25 grains, 100 grains, and 400 grains). Three different grain boundary strength coefficients, 0.5        According to the damage plots shown in Figure 23 with 25 grains at 2.0 µs, the propagation does not always occur from pre-existing cracks. Only the top pre-existing crack propagates in the 100-grain model and both pre-existing cracks propagate in the 400-grain model. This is because with an increase in the total number of grains, the probability of the pre-existing cracks being located on a grain boundary increases. In other words, since the grain boundary strength β = 0.5 promotes intergranular fracture mode, the crack can more easily propagate if it is located on the grain boundary. However, for the grain boundary strength values of β = 1.0 and β = 2.0, there is no such difference observed in fracture behaviour and both pre-existing cracks propagate as shown in Figures 26-31. According to the damage plots shown in Figure 23 with 25 grains at 2.0 μs, the propagation does not always occur from pre-existing cracks. Only the top pre-existing crack propagates in the 100-grain model and both pre-existing cracks propagate in the 400-grain model. This is because with an increase in the total number of grains, the probability of the pre-existing cracks being located on a grain boundary increases. In other words, since the grain boundary strength    According to the damage plots shown in Figure 23 with 25 grains at 2.0 μs, the propagation does not always occur from pre-existing cracks. Only the top pre-existing crack propagates in the 100-grain model and both pre-existing cracks propagate in the 400-grain model. This is because with an increase in the total number of grains, the probability of the pre-existing cracks being located on a grain boundary increases. In other words, since the grain boundary strength

Conclusions
In this paper, a new ordinary state-based peridynamic formulation is presented and related derivations are provided. The current model does not have any limitations on material constants as in the bond-based peridynamic theory. Static analyses were carried out for validation purposes and a comparison of results between PD and FEM shows that the proposed PD model can accurately capture the deformation behaviour of cubic polycrystalline materials. Then, dynamic analyses were carried out by considering different configurations to investigate the effect of interface strength coefficient, discretization size, and crystal size. The observations based on the evaluated results can be summarized as:

1.
Intergranular and transgranular fracture modes can be captured by changing the interface strength coefficient. As a future study, by comparing the experimental and PD results of crack morphology, actual interface strength coefficients can be estimated.

2.
The accuracy of simulation can be improved by increasing the total number of particles for intergranular fracture. However, the difference is not significant for transgranular fracture.
In order to prevent the simulation from being time-consuming, a good balance should be considered between accuracy and simulation time.

3.
Pre-existing cracks can propagate more easily with decreasing crystal size for inter-granular fracture mode, since there is a higher probability of a pre-existing crack interacting with a grain boundary.
As a future study, experimental studies can be used to validate and refine the damage predictions of the proposed PD model. Moreover, as the current study is mainly focused on a 2D model, the formulation can be extended to a 3D model. Acknowledgments: Results were obtained using an EPSRC-funded ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk). EPSRC grant no. EP/K000586/1. critical stretch x vector defining the position of a generic particle x x vector defining the position of a generic neighbour of particle x y vector defining the position of particle x in the deformed configuration y vector defining the position of particle x in the deformed configuration b (x, t) body