We will now introduce a new numerical method called peridynamics, which appears to be very promising for fracking problems. The main difference between the peridynamic theory and classical continuum mechanics is that the former is formulated using integral equations as opposed to derivatives of the displacement components. This feature allows damage initiation and propagation at multiple sites, with arbitrary paths inside the material, without resorting to special crack growth criteria. In the peridynamic theory, internal forces are expressed through non-local interactions between pairs of material points within a continuous body, and damage is a part of the constitutive model. Interfaces between dissimilar materials have their own properties, and damage can propagate when and where it is energetically favourable for it to do so.
Gabriel Hattori1, Jon Trevelyan1, Charles E. Augarde1, William M. Coombs1, Andrew C. Aplin2
1School of Engineering and Computing Sciences, Durham University, South Road, Durham DH1 3LE, UK. 2Science Labs, Department of Earth Sciences, Durham University, Durham DH1 3LE, UK
Received: 2 September 2015, Accepted: 6 January 2016, Published online: 29 January 2016
©The Author(s) 2016.
The peridynamics formulation was first developed by Silling , where he tried to overcome the limitation of current theories dealing with discontinuity, such as in fracture mechanics problems. The main argument was that the difficulty of existing theories was due to the presence of partial derivatives in the formulation to represent the displacement and forces, making necessary specific approaches to eliminate the singularities which would arise. Silling proposed a new formulation based on particular interactions as in molecular dynamics, but applied to continuum mechanics. The term peridynamics was adopted to describe this formulation, and it comes from the Greek roots for near and force. The pairwise interaction between two particles can be defined as 
where ρ is the mass density, f is the pairwise force function that the particle x′ exerts on the particle x,H is the neighbourhood of x, u is the displacement vector field, b is a prescribed force vector field (per unit volume). It is usual to adopt the relative position ξ of the two particles in the reference configuration as
Analogously, the relative displacement η is stated as
The current relative position can be easily given as η+ξ. The function f must satisfy two conditions
which represents Newton’s third law and enforces conservation of linear momentum, and
which assures conservation of angular momentum.
The interaction between particles is defined as a bond, which in continuum mechanics could also be considered as a spring connecting two particles. This definition is fundamentally the difference between the classical theory and peridynamics, where the main idea is the direct contact between two particles. The area of influence of a particle is defined as the horizon δ and is stated as
Fig. 18 Particle interaction in a peridynamics solid
Figure 18 illustrates the horizon δ in an arbitrary body. Outside the horizon δ, a particle has no influence on the other particles. For this reason, the peridynamics formulation is considered as a non-local model.
A material is microelastic if the pairwise function can be obtained through derivation of a scalar micropotential w such as
The micropotential w is the energy present in a single bond (in terms of energy per unit volume squared). Thus, the local strain energy density is defined as
where the factor 1/2 is present since each particle possesses half of the energy of the bond between them. If a material is microelastic, then every pair of particles x and x′ is connected by a spring. The force in the spring depends only on the distance between the particles in the deformed configuration. Hence, there is a scalar function w^ such that
From Eqs. (88) and (90), the pairwise function f is restated as
From Eqs. (82) and (91), the peridynamics model is fully defined for a non-linear microelastic material. However, a linearised theory of the peridynamics microelasticity can be defined as
where C is the material’s micromodulus function. It will be seen that the micromodulus has similar function to the material constitutive law. For more details, see reference .
Boundary conditions in peridynamics are not completely alike to the classical theory. Although the essential boundary condition is still present (displacements), there are no natural boundary conditions (tractions) in the peridynamics framework. Forces at the surface of a body must be applied as body forces b acting through the thickness of some layer under the surface. Usually, the thickness is taken to be the horizon δ. The displacement boundary conditions also have to be imposed as a volume rather than a surface. For more details see .
8.2 Constitutive Modelling
We assume that the bond force f depends only on the bond stretch s, defined as
As expected, s is positive only when the bond is under tension. Failure is introduced into the peridynamics model through breakage of the bonds connecting two particles over some stretching limit. Once a bond fails, it never becomes reconnected (i.e. no healing is considered). An example of a history dependent model is given by the prototype microelastic brittle (PMB) material, and is given by
where g(s)=cs,c is a constant and μ is a history-dependent scalar-valued function, assuming either the values 0 or 1 according to
In this case, s0 is the critical stretch for bond failure. The local damage at a point can be defined as where x has been included as a reminder that the history model also depends on the position in the body. One can see that 0≤φ≤1, 0 representing the undamaged state and 1 representing full break of all the bonds of a given particle to all other particles inside the horizon δ. The broken bonds will eventually lead to some softening material response, since failed bonds cannot sustain any load.
There are only two parameters that define the PMB material, the spring constant c and the critical stretch s0. Assuming η=sξ and substituting in Eq. (89), the local strain energy can be expressed as
This relation must be identical to its equivalent in the classical theory, W=9ks2/2, where k represents the material bulk modulus . The spring constant of the PMB material model is obtained as
Now we describe the bond breakage formulation. Let the work G0 necessary to break all the bonds per unit fracture area be given as
Fig. 19 Fracture energy evaluation
The elements of Eq. (100) are depicted in Fig. 19. Equation (100) is the energy to break all points A, where 0≤z≤δ from the points B. After evaluation of the integrals we obtain
8.3 Anisotropic Materials in Peridynamics
The peridynamics formulation was initially presented for isotropic materials, in order to make some simplifying assumptions. It is expected then that the spring stiffness of the bonds does not vary over the direction of ξ. It was demonstrated in detail in  that for isotropic materials, the Poisson’s ratio in the peridynamics formulations is constrained to take the constant value of 1/4. The constant Poisson’s ratio is a consequence of the Cauchy relation for a solid composed of a lattice of points that interact only through a central force potential .
Refinements of the peridynamics theory can allow the dependence of strain energy density on local volume change in addition to two-particle interactions .
A composite material is formed by different materials, commonly a brittle and stiff material (fibre) embedded into a ductile one (matrix). In , the micromodulus C is redefined in order to accommodate the new variables arising from the material’s anisotropy, including the fibre and matrix bonds for a laminate, and the shear and interlayer bonds present between two different laminates. However, in real composite materials, the fibre and matrix present properties vary significantly with the direction, which was not the case in this work. Instead, different isotropic materials were employed to form the composite fibre and matrix.
In , the fracture in fibre-reinforced composites is tackled with more attention to the material modelling, where the differences between the fibre and matrix bonds are specifically defined. Moreover, the effect of arbitrary fibre orientation in the peridynamic model is taken into account, and it is shown that for a given particle x, the number of fibre bond particles within the horizon δ can vary considerably, which leads to large variation of the strain energy density, the parameter which describes the bond stiffness. To consider this modelling issue, a semi-analytical model was deduced for fibre orientation of 45∘, and also for random fibre orientation.
Fig. 20 Direction of a peridynamic bond in the principal axes
A recent work  has deduced a peridynamic formulation for orthotropic media. The micromodulus C is defined in terms of the orientation of the angles φ, as illustrated in Fig. 20. The dependency on the θ orientation can be suppressed since the material properties do not change over θ for a transversely isotropic material. After some mathematical manipulation, the new definition of the micromodulus is given as
where An0 represents constant coefficients and Pmn are the associated Legendre functions of degree n and order m
Equation (103) can be further simplified into
Assuming c(0)=c1 and c(π/2)=c2, it can be shown that Eq. (104) is also equivalent to
where c1 and c2 are constants of the material model and are given by
where C11,C22,C16 and C66 are elements of the constitutive matrix given in the Voigt notation. Note that an orthotropic material has only 2 independent material constants in the peridynamic model instead of the normal 4 independent constants. This restriction is used linked to the fact that a point is only able to interact to another one individually, while in the classical theory this condition does not apply (a disturbance in a continuous point will automatically induce some disturbance on the points around the body). This restriction on peridynamics theory has been addressed by Silling et al.  and will be detailed in the next section.
The critical bond stretch also depends on the direction of ξ and is given by
where Bn0 are constants and are detailed in . The critical strain energy release rates for mode I crack propagation in the planes normal to the principal axes 1 (GIc1) and 2 (GIc2) can be obtained from the following relations
After integration of Eqs. (110) and (111), the critical stretches s01 and s02 are given by
The fracture behaviour of the material is fully defined by using the mode I energy release rates. Hence, mode II energies are not independent from mode I, which is another consequence of the bond-based peridynamic theory.
An important issue has been highlighted in [95, 201], concerning the use of “unbreakable” bonds near to the regions where a traction boundary condition is applied. The possible reason for this would be crack initiation and propagation close to these regions, due to the high stresses that could be present. It is important to understand the physics of the analysed problem properly in order to use this type of assumption during a peridynamic simulation.
8.4 State-based Formulation
The peridynamics formulation assumes that any pair of particles interacts only through a central potential which is independent of all the other particles surrounding it. This oversimplification has led to some restrictions of the material’s properties, such as the aforementioned fixed Poisson’s ratio of 1/4 for isotropic materials. Also, the pairwise force is responsible for modelling the constitutive behaviour of the material, which is originally dependent on the stress tensor. To overcome this limitation, Silling et al.  have extended the peridynamics formulation to include vector states. The vector states allow us to consider not only a particle, but a group of particles in the peridynamics framework. Moreover, the direction of the vector states would not be conditioned to be in the same direction of the bond, as in the bond-based theory. This property is fundamental to consider truly anisotropic materials.
Let A be a vector state. Then, for any ξ∈H, the value of A⟨ξ⟩ is a vector in R3, where brackets indicate the vector on which a state operates. The set of all vector states is denoted V. The dot product of two vector states A and B is defined by
The concept of a vector state is similar to a second order tensor in the classical theory, since both map vectors into vectors. Vector states may be neither linear nor continuous functions of ξ. The characteristics of the vector states are listed in , and they imply the vector states mapping of H may not be smooth as in the usual peridynamic model, including the possibility of having a discontinuous surface.
In the state theory, the equation of motion (82) is redefined as
with T as the force vector state field, and square brackets denote that the variables are taken in the state vector framework.
To ensure balance of linear momentum, T must satisfy the following relation for any bounded body B
The balance of angular momentum for a bounded body B is also required
The deformation vector state field is stated as
The non-local deformation gradient for each individual node is given by
where ⊗ denotes the dyadic product of two vectors, and ω(|ξ|) is a dimensionless weight function, used to increase the influence of the nodes closes to x. The use of this factor is still under study , but the assumption of ω(|ξ|)=1 has been seen to provide good results.
The discretisation of Eqs. (120) and (121) can be expressed as a Riemann sum as 
where m is the number of nodes with the horizon of node j. xj must be connected to at least three other nodes in the system to ensure that B(xj) will not be singular.
In state vector peridynamics, there are two ways to determine how the force state depends on the deformation near a given point. The first consists of formulating a constitutive model in terms of the force vector T and the deformation state Y [x,t]. In this case, the force state is defined as
where W is the strain energy density and ∇ indicates the Fréchet derivative, which is defined as any infinitesimal change in the deformation state dY resulting in a change of the strain energy density dW such as
with Hx being a sphere centred at the point x with radius equal to the horizon δ. Note that the Fréchet derivative can be seen as an equivalent of the tensor gradient in classical theory.
The second approach to relating the force and deformation in a state vector framework is to adopt a stress-strain model as an intermediate step [42, 280]. For a strain energy density W(F), the stress tensor can be expressed as
The force vector is redefined as 
After evaluation of the Fréchet derivative, the force vector can be defined explicitly as
The processing of mapping a stress tensor as a peridynamic force state is the inverse of the process of approximating the deformation state by a deformation gradient tensor. A peridynamic constitutive model that uses stress as an intermediate quantity results in general in bond forces which are not parallel to the deformed bonds. This type of modelling was called “non-ordinary” by Silling .
8.5 Numerical Discretisation
The discretisation of the peridynamics model is quite straightforward. Equation (82) can be rewritten as a finite sum
where n is the time step and subscripts denote the node number, i.e., uni=u(xi,tn),Vp is the volume of node p. Equation (129) is taken over all p nodes which satisfy |xp−ui|≤δ. The grid spacing Δx is also an important parameter in the peridynamics discretisation.
The discretised form of the linearised peridynamics model is given by
The displacements uni are obtained using an explicit central difference formulation,
with Δt as the time step. Some studies of the stability of the numerical discretisation were described in [154, 249]. It has been established that the time step must not exceed a certain value in order for the numerical discretisation to be stable. Moreover, the error associated with the discretisation depends on the time step with (O(Δt)) and the grid spacing with (O(Δx2)),
Convergence in peridynamics is affected by two parameters: the grid spacing Δx and the horizon δ. Reducing the horizon size for a fixed grid spacing will lead to the peridynamics solution approximating the solution using classical theory. However, fixing the horizon size while increasing the grid spacing will lead to the exact non-local solution for that particular horizon size . As for domain discretisation methods, it is important to balance the size of the horizon so the damage features in the analysed body are properly considered, and the grid spacing should be sufficiently small for the results to converge to the non-local solution. Usually, it ranges from 1/3 to 1/5 of the size of the horizon.
In recent works, the peridynamics formulation is used conjointly with other discretisation methods, such as meshless formulation  and finite element formulation .
In , peridynamics is used only to obtain the prediction of failure of the composite material, where an FEM code is employed to solve the global problem. This type of combined approach is often necessary since the peridynamic formulation can demand significant computational power, a common problem in molecular dynamics simulations as well.
9 Conclusions and Prospective Work
We have seen that the hydraulic fracture problem presents several characteristics which makes its study complicated: the shale is not a homogeneous material, it is not isotropic, the nanoporosity may retard crack propagation as the fluid penetrates the rock, and a large fracture network has to be considered in which cracks develop at multiple length scales, all of which can greatly increase the computational solving time. Moreover, most current analytical and numerical methods do not take into account crack branching, a key factor in order to obtain a correct estimation of the extended fracture network.
The current fracture models for brittle rocks and fracking have been useful as a first step in offering a more realistic fracking model. There are of course other limitations attached to each of the numerical models discussed earlier: for instance, in cohesive models, the cohesive zone model is not a parameter to be found, so the crack propagation path is already known a priori. Most works on X-FEM and BEM models consider that the crack propagation path is unique; only recently have some works appeared considering crack branching [184, 244, 285].
Fracking models developed so far have not considered the full complexity of shale rocks. Ulm and co-workers [268, 269, 270] have established that shales are likely to be transversely isotropic materials, with the direction perpendicular to the bedding planes taken as the symmetry axis. This is mainly due to the deposition process. It was also stated that the shale anisotropy is due more to the interaction between the particles than the elastic behaviour of the shale components.
It was seen in  that the fluid penetrating the crack may retard crack propagation, so the material’s porosity has to be taken into account in the numerical model.