An analytic element model for highly fractured elastic media

We present an analytic formulation for an elastic medium under uniform stress and in plane strain, with a very large number of cracks that are allowed to intersect. The singularities at the tips of the cracks are represented exactly; the solution is written in terms of series expansions of negative powers of a complex variable defined outside the unit circle in a complex reference plane. We modeled an assembly of 10,000 cracks and show stress trajectories both for the whole model and for individual cracks, chosen at random from the assembly.


INTRODUCTION
The exact solution for a pressurized crack in an elastic medium under plane strain conditions was developed by Griffith 1 . Westergaard 2 used complex variables and a stress function introduced by Carothers 3 to study problems involving cracks. Plane strain problems in linear elasticity can be solved efficiently using the complex variable method, based on the Muskhelishvili-Kolosov displacement functions, Muskhelishvili 4 , and are reported in numerous modern textbooks, for example, Sadd 5 , and research articles, for example, Strack and Verruijt 6 and Verruijt 7 . We apply the displacement functions to develop an analytic element for a pressurized crack under plane strain; the derivation follows that given by Strack 8 .
Since the solution by Westergaard 2 , numerous new solutions were developed. Sneddon 9 and Sneddon and Elliott 10 extended the solution to include a variable pressure acting in the crack. Sih 11 used the displacement functions given in Muskhelishvili 4 to generalize the solution to that of a crack in an infinite plane under biaxial tension and pure shear, and Sanford 12 modified the solution to be valid for a crack near a boundary. Zhang et al. 13 extended the work by Sneddon 10 to a formulation capable of modeling nonuniform fluid pressure.
Gerolymatou 14 used conformal mapping to map a domain outside of a given shape onto the exterior of the unit circle and used a series expansion to obtain the solution; the use of conformal mapping and series expansions showed great agreement with the analytic solution. The solution which was presented by Gerolymatou 14 has similarities with the one presented in this paper.
There exist analytic solutions for multiple cracks in an elastic plane, for example, Benthem and Koiter 15 and Sneddon 16 . Several approaches are presented by Kachanov 17 , who suggested to superimpose cracks by finding the average tractions acting on each crack. Datsyshin and Savruk 18 suggested to use the Fredholm integral equation, as did Chen 19,20 , while Horii and Nemat-Nasser 21 suggested to use the method of pseudo traction.
There exist various other solutions for multiple cracks that are based on numerical methods, for example, the boundary collocation method, for example, Cheung et al. 22 and Wang et al. 23 , and the boundary element method, for example, Denda and Dong 24 and Guo et al. 25 . Helsing 26 presents a solution for a large number of cracks (10,000) in a periodic setting used to gain efficiency, based on solving the Fredholm integral equation with a high degree of accuracy using quadrature.
Further reading on crack solutions can be found in the work by Tada 27 , Kachanov 17 , Sih 28 and Sneddon 29 where a great collection of crack solutions are presented.
We base our formulation on the Muskhelisvili-Kolosov functions, written in a somewhat different form. We transform the Cartesian ( , )-space into the (̄, )-space where the bar indicates the complex conjugate, as suggested by Green and Zerna 30 . This formulation is derived for general cases by Strack 8 ; we merely summarize the formulation here. Our solution is based on the analytic element method; an overview of the method as applied to groundwater mechanics is given in Strack 31 , and an approach similar to that reported in this paper is found in Strack 32 . We express the exact solution for a single crack under plane strain conditions in an elastic medium in terms of the Muskhelishvili-Kolosov displacement functions written in terms of a complex variable , defined outside the unit circle in its plane. These functions are holomorphic outside the unit circle and thus may be expanded in terms of a Laurent series with exclusively negative powers (asymptotic expansions). We determine the coefficients of these expansions from the boundary conditions, which include the effects of all cracks in the system plus a given uniform stress. We achieve accuracy of up to machine precision for very large numbers of cracks, and show examples of assemblies of 10,000 cracks; the number of cracks is not fundamentally limited. As the formulation is "embarrassingly parallel," significant gains can be obtained by using multiprocessing.
We imagine our medium to be idealized rock and adopt the convention that compressive stresses are positive.

TRANSFORMED TENSORS AND TRACTIONS
We formally transform the stress tensor components , (where or equal to one represents 1 ≡ and , = 2 represents 2 ≡ ) into the contravariant components in the = ( 1 , 2 ) ≡ (̄, )-space, We consider only two components in complex space because of the following relations: The expressions for the components of the transformed strain tensor, ( = 1, 2; = 1, 2) in terms of the Cartesian strain tensor components are 11 = 11 − 22 − 2i 12 , These expressions differ from those used by Muskelishvili; the components of the transformed stress and strain tensors are the complex conjugates of the ones used by Muskelishvili. This is the consequence of applying the formal transformations that follow from the characteristics of the governing partial differential equations for divergence and curl, see Strack 8 .
where [F/L 2 ] is the shear modulus, and [.] is Poisson's ratio. We write the displacement components in terms of local Cartesian coordinates ( , ), that are tangential ( ) and normal ( ) to a plane at an angle to the -direction. We represent these displacements as and and combine them into a complex displacement We transform the Cartesian tangential, , and normal, , components of the traction in terms of the normal, = (cos , sin ), and transform the result in terms of the contravariant components of the traction, see Strack 8 We apply the transformed equilibrium conditions and after some manipulation obtain the following expression for the complex displacement in the absence of body forces, see Strack 8 The expressions for the transformed stress tensor components are and These equations are equivalent to the Muskhelishvili-Kolosov functions, see Muskhelishvili 4 . The functions Φ and Ψ are holomorphic functions in terms of a single complex variable.

A PRESSURIZED CRACK
The analytic element for a pressurized crack is based on an alternative expression for a single pressurized crack in an infinite domain. We write the classical analytic solution for a single pressurized crack in an infinite domain in a different, but equivalent, form; it simplifies the formulation for the analytic element that permits superposition.

Exact solution
We transform the crack shown in Figure 1(A) into the unit circle in a reference complex plane, the -plane in Figure 1(B), in such a way that one side of the crack, indicated with a + and shown in blue, corresponds F I G U R E 1 Transformation of the crack in the -plane onto a circle in the -plane to the upper half of the circle in the -plane and the other side, marked with a − and shown in red, corresponds to the lower half. This property of the transformation is useful when applying and verifying the boundary conditions. We first transform the crack in the physical plane, extending from 1 to 2 , into a line of length 2 in the dimensionless -plane where is the angle between the crack and the -axis, and is the length of the crack. The transformation ( ) maps the end points of the crack onto = −1 and = 1, and 0 is the complex coordinate of the center of the crack The corresponding expression for in terms of is The function that transforms the exterior of a slot from the -plane onto the -plane is a standard conformal mapping A point on the unit circle, = e i , corresponds to = cos , which is real and varies between −1 and 1. The inverse of ( ) is we choose the positive sign to lead the square root, so that infinity in the -plane corresponds to infinity in the -plane.
We compute the square root as the product of two square roots to obtain the proper mapping, see Strack 8 The displacement components in terms of local coordinates

F I G U R E 3 The traction components in terms of local coordinates
We express̄in terms of̄in the expressions for , 11 and 12 , Equation (9) through Equation (11) and The complex displacement in the plane of the element is given by Equation (6), with Equation (18) where and are defined in Figure 2.
The expression for the traction 1 in terms of 11 and 12 is given by Equation (8) where and are defined in Figure 3.
We introduce a function = ( ), defined by the equation and use it in Equation (21) We introduce a function differentiate Ψ with respect to and use this expression for Ψ ′ in Equation (19) for 11 We obtain the expression for 12 in terms of ′ from Equation (20) with Equation (25) The analytic element for a pressurized crack is based on an alternative expression for a single pressurized crack in an infinite domain, written in terms of . The following choice for and renders a solution that meets the boundary conditions along a pressurized crack The boundary conditions are that the normal and shear components of the tractions along the inner crack faces are equal to − and 0, respectively, taking the normal to each crack face to point out of the domain and into the crack. This implies that the complex traction 1 is We use Equations (27) and (20) with Equation (25) in expression (22) for the traction We chose ′ ( ) = −̄′( ) and equals̄along the element This expression is purely imaginary, so that = 0, as required. We differentiate with respect to and obtain an expression for ∕ from so that Equation (33) becomes We use this in the expression for − i , Equation (32) and use that̄= 1∕ This solution indeed meets the boundary condition that = 0 and = − .

ADDING A UNIFORM STRESS FIELD
We subject the crack to a uniform stress field with its principal axis in the -direction, and with zero minor principal stress. We modify the solution for a single crack by replacing the constant in the solution by a complex constant̃. We represent the constant stress as We integrate this, and require that the displacement be zero at the origin, which gives The contribution of the uniform far-field to the complex traction along a crack at an angle to the -axis, is given by We use Equation (41) or, writing cos(2 ) − 1 as −2 sin 2 ( ) and sin(2 ) as 2 sin( ) cos( ) We add the contributions of the far-field and the crack and solve for the coefficients.

THE ANALYTIC ELEMENT
We generalize the special form of the exact solution for a pressurized crack to create the analytic element, which can be combined with other analytic elements to create a new solution. The objective is to generalize the exact solution such that the boundary conditions can be met precisely, while maintaining the singular form of the element at the tips. The element provides the degrees of freedom necessary to account for the effect of other elements, such as other cracks in the medium. The functions and are holomorphic outside the unit circle; the analytic element only represents the crack itself. The functions and can thus be expanded in terms of asymptotic expansions, valid outside the unit circle We demonstrate that these expressions can be used to model a crack with discontinuous displacements and with normal and shear stresses along the crack expressed as series expansions in terms of powers of along the crack boundary, that is, in terms of a Fourier series. We determine the coefficients and from the conditions that the sum of all analytic elements in the solution plus the element itself meet the boundary conditions along each crack We apply the expressions for , 11 , and 12 as given by Equations (24), (27), and (28) with expressed in terms of .
We differentiate Equation (46) with respect to where ∕ is given by Equation (34), and we obtain ∕ from Equation (14) = We obtain the expression for ∕ in the same way We differentiate ∕ with respect to to obtain ′′ ( ) or, with ∕ , Equation (34), and ∕ , see Equation (14) 2 and simplify

Tractions
Expression (31) for the tractions along the element, where =̄, becomes We use Equations (50) and (51) to write this in terms of We use that̄= 1 along the element, so that and rewrite Equation (56)

Continuity of the tractions
The tractions must be continuous across the element; the same value must be obtained for as for̄= −1 where + and − labels in the exponent position indicate the + and − sides of the crack. We multiply by and combine terms This condition must be fulfilled for all values of on the unit circle, excluding the singular points at the crack tips ( = ±1) so that expression (58) for the traction along the element becomes We write = e i along the unit circle and separate real and imaginary parts The complex coefficients must be chosen such that the boundary conditions are met. This implies for the case of a single crack that = 0, or ℑ = 0, = 1, … , ∞ and that = − , or ℜ = 0, = 2, 3, … , ∞ and which matches the exact solution for a single crack along the real axis (i.e., the angle = 0), compare Equation (29).

The stresses
We find the components of the stress tensor due to a crack at any point of the domain from Equations (27) and (28), with the expressions for and , and their first derivatives. The second derivative of is given by Equation (54). Tractions on any plane in the domain are found from Equation (8).

The displacements
We obtain the expressions for and from Equations (46) and (47) with Equation (62) We find the displacements along the crack from Equation (24) We set =̄to obtain the displacement along the crack or, with Equations (68) and (69), and̄= 1 There is a jump in the displacements or For the special case of a single crack oriented along the real axis, where 1 = ∕2, and = 0, = 2, … , ∞, this becomes The jump in displacement is purely imaginary for this case, and equals The normal displacement component on the positive side minus that on the negative side is positive; the crack opens. This solution matches that obtained from the exact solution written in terms of .

IMPLEMENTATION
We superimpose the analytic elements to obtain the expressions for both the displacements and the stress tensor components for a set of fractures. We use the degrees of freedom of each element, the complex coefficients , to meet the boundary conditions; the accuracy with which these conditions are met depends only on the precision of the machine used. We use an iterative scheme introduced by Janković and Barnes 33 , and an approach similar to that proposed by Strack 32 , which leads to a solution that converges to the exact one with up to machine accuracy. We represent the number of elements by , and truncate the infinite series at . We number the elements as = 1, 2, … , and denote their end points as 1 and 2 . We introduce functions ( ), one for each crack, using definition (12) and define according to Equation (16) ( ) = ( ( )).
We introduce functions ( ) and ( ), one for each element, and label their coefficients as

Global displacements
The contribution of element to the global displacement is , obtained from the local displacement , (70), by multiplication by e −i , where is the angle between element and the -axis where is the length of element . We add all contributions to obtain the global displacement

Global stress tensors
We find the contributions of element to the global expressions for the transformed stress tensors from Equations (27) and (28) We obtain the expressions for the stress tensor components due to all element by superposition and The tractions along element due to element are given by Equations (65) and (66), with the angle defining a point on the unit circle in the -plane of the element. The sum of the contributions of all elements, including element , must satisfy the condition that the shear traction is zero and the normal traction is equal to minus the pressure specified for the element, − where is a point of the element, and 11 and 12 represent the stress tensor components obtained by adding the contributions of all elements, excluding element and We apply the complex condition (87) for the complex constants of element at a minimum of 2 points of the unit circle of the -plane of element . This leads to 2 equations for the 2 unknowns of element . In case over specification is used, that is, if we apply the condition at more points than the minimum, the number of equations exceeds the number of unknowns and we minimize the sum of the squares of the errors. For each of these points, we compute the value of the global coordinate that corresponds to the point of the circle, and use it in the global functions inside the brackets of Equation (87).

INTERSECTING CRACKS
An assembly of cracks generally includes intersections. Each crack element creates a displacement discontinuity, and the displacement discontinuity at the point of intersection is the algebraic sum of the two displacement discontinuities of the intersecting cracks. Moving through the point of intersection along either one of the cracks, the displacement discontinuity jumps. We apply the conditions that are imposed at any point of the cracks for each of the two cracks also at the point of intersection. Since the conditions are valid for the tractions along the planes of the two cracks, the stress state at the point of intersection is isotropic; normal stresses are equal at two different planes, and the shear stresses along these planes are zero; Mohr's circle reduces to a point.

Results
We applied the approach described in this paper to a set of 10,000 cracks in plane strain in an elastic medium subjected to a uniform stress in the -direction of magnitude 10 units, and pressures in the cracks that vary between zero and 1 units. A plot of the major and minor stress trajectories are shown at progressively smaller scales in Figures 4 and 5.

F I G U R E 6 Stress discontinuity across a crack
The transformed stress tensor component 11 offers a direct way to obtain the values of both the stress deviator and the directions of the major and minor principal stresses. We introduce a function so that where is the stress deviator, and the principal direction. There may be a visible inconsistency as stress trajectories intersect cracks. Although the tractions are continuous across the cracks, the stress tensor components are not. The principal stress that is tangential to the crack can be, and often is, different at the two sides of the crack; such a discontinuity does not violate the conditions of static equilibrium as illustrated in Figure 6. The major and minor principal stresses switch as a result of the discontinuity of the stress component̃̃, wherẽand̃are local coordinates as indicated in the figure. This causes a sign change in the arctangent in Equation (90).
Errors reported here are relative; they are the computed minus the desired shear traction (zero) and the difference between computed and specified normal traction. All of the errors are divided by the stress tensor component 12 , which is 2 along the crack. The relative mean error as a function of the number of coefficients is shown in Figure 7. It illustrates that the solution preforms well even for a small number of coefficients; in this example, the relative mean error is 2.6E−5 when using 100 coefficients. This case applies to two intersecting cracks, as displayed in Figure 8.

CONCLUSION
We constructed analytic elements capable of modeling very large assemblies of pressurized cracks in an elastic medium under plane strain conditions. We demonstrated that the solution is both very accurate and potentially extremely efficient, with no theoretical limitation on the number of cracks, thanks to an iterative scheme. Computational time depends linearly on the number of processors, a property known as "embarrassingly parallel". Computation of the coefficients also benefits greatly from parallel processing, but not quite as much as for the computation of results. We demonstrated that the analytic element method, developed originally for groundwater flow modeling by the lead author, also applies to rock mechanics and other problems governed by the equations for linear elasticity. This is encouraging, because for the case of groundwater flow, the method applies also to three-dimensional flow through large numbers F I G U R E 7 Relative mean error in (red curves) and (blue curves) as a function of the number of coefficients F I G U R E 8 Displacement vector field for a pair of intersecting cracks of elliptical inhomogeneities, as demonstrated by Janković and Barnes 34 . That this approach, consisting of using series expansions in ellipsoidal coordinates in three dimensions, can be applied to problems of three-dimensional cracks is by no means certain, but the present paper suggests that it is at least worth a serious attempt.
The objective of creating analytical models of fissured rock and other fractured elastic media is the ability to study the behavior of highly fractured rock with detailed data regarding stresses and displacements available to a high degree of accuracy. Numerical laboratory experiments could be used, for example, to define meaningful definitions of average stresses in sub domains, and define aggregate material properties. This goal is very similar to that of the distinct element method presented in a paper by Cundall and Strack 35 , which was originally also restricted to two-dimensional assemblies, and is the first step in developing highly accurate and sophisticated three-dimensional distinct element models (the method was later renamed to "discrete element method"). The advantage of such models over both laboratory studies and discrete numerical models is that all local values are available with high precision at any point in the medium, and that singularities and their effect on the material behavior are represented with a very high degree of precision. We hope that the work presented in this paper is a first step in a lengthy process of developing increasingly realistic analytic element models for elastic media.

A C K N O W L E D G M E N T S
The corresponding author wishes to thank Auli Niemi and Ramin Moghadasi for their comments on this work. The author Erik A.L. Toller was supported by the Swedish Transport Administration [Grant Number 7017].

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in AE_LE_crack at https://github.com/eriktoller/ AE_LE_crack.