format: text

text:
# Abstract

## Abstract

We review microstructural fracture growth models suitable for the study of hydraulic fracture processes in disordered porous materials and present some basic results. It is shown that microstructural models exhibit certain similarities to corresponding theories of continua. These similarities are most easily demonstrated for simple crack geometries, i.e., straight cracks (finite size scalings). However, there exist even scaling relations which are completely independent of the particular employed crack structure. Furthermore it is demonstrated that disorder in cohesional/flow properties can influence the crack growth and the resulting fracture geometry in an essential way.

## 1 Introduction

The classification and prediction of breaking processes within solid materials forms an extremely difficult problem in engineering science as well as in solid state physics. This is even true for brittle materials exhibiting the simplest rheological behavior, i.e. no remanent deformations. Continuum mechanical treatments of brittle solids based on Griffith’s theory$`^{\text{References}}`$ are quite widespread and accepeted in scientific papers, and have been relatively successful for practical and engineering purposes.$`^{\text{References}}`$ There are, however, a couple of problems being associated with classical fracture mechanics limiting its range of application. Certain restrictions arise from the assumption that the deformations can be conducted in a thermodynamical reversible manner. Connected to this is the basic assumption that the involved fracture processes are thermodynamically reversible too. This is usually expressed by minimizing thermodynamic potentials with respect to crack/fracture growth, which is a perturbational scheme. The typical modelling situation is such that some governing field equations are describing the global state of the deformation (i.e. the free energy density) according to some boundary/loading conditions while the free surface energy field (intrinsic strength) is classically assumed to be a constant property throughout the solid material. The assumption of a constant intrinsic strength constitutes a quite strong idealization for most experimental situations because it eliminates the field nature of cohesional properties. The assumption of a constant theoretical strength is most easily justified for homogeneous materials, i.e., perfect single crystals. However, even for such classes of materials there exist certain limitations: a) the crystal’s cohesion is in general non-isotropic (variation in crystal strength according crystallographic directions) and b) the cohesion is subject of thermal influences (lattice vibrations) introducing thermal fluctuations into the problem.

The survey becomes more difficult for most solid materials of practical interest as they are frequently composite materials. Variations of cohesional properties occur already on macroscopic length scales comparable to those of the microstructural elements. There is no simple argument to replace the theoretical strengths of the isolated solid phases by some average effective strength in connection with fracture problems. This is because the fracture and its associated natural geometry advances on a microstructural level where the heterogeneity in grain orientation and strength is present. We mention this in order to illustrate that the assumption of constant cohesion/binding properties throughout a solid and therefore the representation employing a homogeneous continuum needs careful consideration.

The investigation of geometrical properties of irreversibly grown cracks in the presence of disorder has recently received considerable attention amoung physicists. Basic ideas go back to statistical physics, theory of phase transitions, scaling theory, selfsimilarity and self-organized criticality. These concepts are described in good standard books of condensed matter physics, i.e. Chaikin & Lubensky.$`^{\text{References}}`$
Perhaps one of the most well studied phenomena exhibiting self-organized criticality is Diffusion Limited Aggregation (DLA), theoretically first implemented by Witten & Sander.$`^{\text{References}}`$ Its outline is quite simple and it is well suited to exemplify geometrical selfsimilarity.
Particles being small enough to be subject to Brownian motion diffuse in a large gas or fluid reservoir. The particle concentration is assumed to be very low (infinite dilution), i.e., there are at each time step precisely one diffusing particle and one cluster where the diffusing particle is allowed to be attached to. If the diffusing particle touches the particle-cluster it sticks to it forever and a new diffusing particle is started far away from the cluster. The process is called diffusion limited because the particle diffusion is the rate limiting step in this reaction scheme. Figure 1 shows such a DLA cluster.

The cluster has a dendritic structure exhibiting branches on all length scales between particle and cluster length scale. This supposes that the average number $`N`$ of particles belonging to the cluster follows a power law in cluster radius $`R`$, i.e.,

$$NR^{d_f}.$$
(1)
The exponent $`d_f`$ is called the fractal dimension of the cluster. The value for $`d_f`$ is in general always smaller equal the embedding dimension; in the above example $`d_f=1.71<d=2`$. The physical meaning of a power law relation Eq.(1) is that the structure does not contain any natural length scale and that it looks similar to itself (selfsimilar) under a change of scale. The reader interested into concepts and applications of Fractals is referred to References References and References. Nowadays it is widely believed that the value of the fractal dimension does only depend on the symmetries of the governing equations but not on particularities such as the particle shape or the employed discretisation scheme in numerical simulations. Deviations from scaling behavior are anticipated for systems of finite size.

The main interest in the DLA problem results from the fact that the governing equation is the Laplace equation being one of the most widespread and simplest equations in physics. A diffusing particle sticks with probability proportional to $`P/n`$ to the clustersurface with $`P`$ being the solution of the diffusion equation. The sticking-probability is highest near the clusters extremities, so that the cluster is growing at a higher rate at the tips than in its interior. The tips screen the cluster interior to a certain amount leading to a diluted dendritic structure. Other examples belonging to the same universality class as DLA are the dielectric breakdown (DBM) (see Fig. 2), the viscous fingering (VF) (see Fig. 3a) and the Stefan problem of solidification.$`^{\text{References},\text{References}}`$ All this problems lead experimentally and numerically to similar geometrical patterns.

There has been considerable discussion whether fracture processes in presence of fluctuations would lead to DLA-like breakdown pattern too. The conducted investigations do, however, not support this idea in an unambiguous way. Rather the structure of the basic elastic equations (Lamé or Cosserat equations) cannot be reduced to a Laplace type equation even in limiting cases for the elastic constants. The experiments of Van Damme and collaborators on twodimensionally confined clay pastes display a sharp transition between DLA morphologies (viscous fingering, compare Fig. 3a) and the much more spiky and diluted hydraulic fractures, e.g. see Fig. 3b.

Further below we will present numerically obtained results mainly concerned with the problem of hydraulic fracturing of disordered materials in two dimensions.
We will first consider the case in which the model material contains only elastic interactions. We will investigate the cases of a constant injection (borhole-) pressure and of a constant fluid influx.
In a later chapter we consider a liquid saturated porous elastic material being fractured in a similar way. However, in such a case is the pressure distribution originating from the pore fluid flow directly coupled to the elastic equations.

## 2 The Model

In the following we will outline the employed model. We give a brief description of the basic elastic and flow equations. Thereafter we explain in detail the employed boundary conditions. After this we demonstrate how heterogeneous cohesive properties are taken into account. Finally we present the breaking rules which contain the physics of the considered breaking process.

### 2.1 Pore space

There are at least two different ways of implementing porosity numerically. In both cases certain microstructural elements (beams, plates or shells) will form a stress carrying elastic backbone. One can think, however, of two different realizations for the microstructural elements. The elements could be ’made’ of a material which already follows the constitutive laws for porous materials, e.g. Biot’s theory.$`^{\text{References}}`$ This would be most satisfying from a continuum mechanical point of view as long as no fracturing occurs. However, if fracturing at the pore level is to be taken into account the above mentioned approach does not suffice (or even becomes inconsistent) and a micromechanical approach must be employed. The pore space/volume has to be modelled explicitely in that case. We have focussed on the latter implementation: the microstructural elements follow the classical constitutive equations of linear elasticity, and the pore space is viewed as the remaining vacant volume. The breaking of microstructural elements is associated with a coallescence of pore volume under question. Because the pore volume is modelled explicitely the couplings between elastic and flow interactions ought to be derived explicitely from the microstructural configuration.

In this paper we consider a two-dimensional model-microstructure laying on a square lattice with properties of the pore space (hydrostatic pressure field) defined on its dual lattice.$`^{\text{References}}`$

### 2.2 Elastic equations

We consider the beam model (as defined in p. 232 of Ref. References) on a two dimensional square lattice of linear size $`L`$. This model is preferable compared to Hookean spring models from a microstructural point of view because elastic beams due to their finite crossection can be physically interpreted as fibers, while Hookean springs (having zero cross section) can not. In the two-dimensional beam model each of the lattice sites $`i`$ carries three real variables: the two translational displacements $`x_i`$ and $`y_i`$ and a rotational angle $`\phi _i`$ in the xy-plane. Neighbouring sites are rigidly connected by elastic beams of length $`l`$. The beams are assumed to have the same quadratic cross section, $`A=d^2`$, and the same elastic behavior, governed by three material and geometry dependent constants $`a=l/(Ed^2)`$, $`b=l/(Gd^2)`$ and $`c=12l^3/(Ed^4)`$ with $`E`$ and $`G`$ being the Young and shear moduli. With characteristic values for rock materials $`E=10^{11}Pa`$, $`G=410^{10}Pa`$, $`l=10^5m`$ and $`d=10^6m`$ the above three constants become $`a=10^4(Pam)^1`$, $`b=2.510^4(Pam)^1`$ and $`c=0.12(Pam)^1`$. When a site is rotated ($`\phi _i0`$) the beams bent accordingly always forming tangentially $`90^{}`$ angles with each other. In this way local momenta are taken into account. For a horizontal beam between sites $`i`$ and $`j`$ one has the longitudinal force acting at site $`j`$,

$$F_j=\alpha (x_ix_j),$$
(2)
the shear force

$$S_j=\beta (y_iy_j)+\frac{\beta }{2}l(\phi _i+\phi _j),$$
(3)
and the flexural torque at site $`j`$,

$$M_j=\frac{\beta }{2}l(y_iy_j+l\phi _j)+\delta l^2(\phi _i\phi _j),$$
(4)
using $`\alpha =1/a`$, $`\beta =1/(b+c/12)`$ and $`\delta =\beta (b/c+1/3)`$. The corresponding equations for vertical beams are similar. The outline of the beam network is displayed in Fig.4. We will not consider inertial or gravity forces, however, the flow of a pore fluid amounts in general to non-vanishing body forces.

If there is a pressure drop $`p_2p_1`$ across a beam (pore wall) the resulting line force density can be replaced by a shear force dipole, $`S=(p_1p_2)ld/2`$, acting at both ends of that beam. This follows directly from the equilibrium conditions for beams. Hence the dipoles couple the pressure field in form of body forces to the mechanical equilibrium conditions.

In mechanical equilibrium the sum over all elastic and body forces (torques) acting on site $`j`$ must vanish. Thus if the pressure field is known one can calculate the displacements.

To determine the hydrostatic pressure field one needs to consider flow equations.

### 2.3 Flow equations

The beams exhibit a typical length $`l`$ and thickness $`d`$ and the liquid saturated equilibrium pore volume is simply $`(ld)d^2`$. This states explicitly the considered material to be a binary system, i.e. any volume element is either occupied by elastic material or by pore fluid. By varying the ratio of $`l/d`$ different values for the volumetric porosity can be achieved (ranging from nearly zero to nearly one). In order to obtain most simplified flow equations, we make the following assumtions: a) the pore fluid is incompressible (which appears well justified for water) and b) the fluid velocity relative to the bulk is proportional to the gradient of the pressure field (Darcy’s relation). The continuity equations for solid and liquid mass density lead to the required pressure equation (in the continuum limit);

$$\text{div}\left(\frac{k}{\mu }\text{grad}p\right)=\frac{}{t}\text{div}𝐮,$$
(5)
with $`p`$ being the pressure, $`𝐮`$ the displacement vector, $`k`$ the mechanical permeability and $`\mu `$ the fluid viscosity. This equation is a Poisson equation for the pressure field for constant $`k`$ and $`\mu `$. The right hand side of Eq.(5) represents in general a time-dependent source term for the pressure field. It characterizes the local rate of relative volume change due to elastic deformations, and hence the net fluid flux itself. In general these net fluid fluxes will control the pressure field’s time-evolution. However, in this work we will only consider the stationary case. The physical motivation for this is the observation, that the characteristic internal relaxation time $`\tau =\mu /E`$ for water-rock systems is typically of order $`10^6s`$. It appears justified to assume that the time periods between sucessive breaking events are much larger than the internal relaxation time $`\tau `$, and hence that breaking events are triggered by stationary pressure and displacement distributions. In connection with Eq.(5) this implies to balance the net fluid flux for a given cell $`i`$ to zero (Laplacian);

$$\frac{1}{\mu }\underset{j}{}k_{ij}(p_ip_j)+\text{boundary conditions}=0.$$
(6)
Here $`k_{ij}`$ denotes the local permeabilities across the beams (pore walls).

The pressure equation Eq.(6) is solved without referring to any mechanical properties, which is directly a consequence of explicitely dropping all dynamic terms. Its solution for given boundary conditions represents the body force under which the elastic equations, among their own boundary conditions, are solved. Formally the corresponding continuum problem is analogous to the standard problem of thermoelasticity$`^{\text{References}}`$ which was already noted by Biot$`^{\text{References}}`$. From the resulting stresses/strains and breaking thresholds (see Sec.2.5) we then determine which element(s) are broken next.

Here we would just like to mention what happens to the pressure equations, when a beam is broken. Because the beams represent the pore walls, and because the fluid is assumed to be incompressible the pressure within a connected crack has to be constant. In this sense a crack is a macro-pore. Numerically this can be accomplished in different ways. We have chosen to set the permeabilities for broken beams very high, i.e. $`k_{ij}=10^9`$. Hence the pressure drop along a crack is negligible small.

### 2.4 Boundary conditions

The results presented in this review are mainly for hydraulic fracturing conditions in impermeable$`^{\text{References},\text{References}}`$ and permeable$`^{\text{References}}`$ materials. The only exception is Sec. 3 where we compare results for tensile and hydraulic loading configurations.$`^{\text{References}}`$ The corresponding tensile boundary conditions are therefore discussed in Sec. 3.

Let us make few remarks about terminology. A boundary will be called ‘free’ if its normal stress component vanishes and will be called (hydrostatically) ‘loaded’ if the normal stress component follows some known non-zero stress distribution (von Neuman condition). A boundary is called ‘fixed’ or ‘clamped’ if the displacement values are assigned to the boundary (Dirichlet condition). A particularity in calculations are the sometimes employed ‘periodic’ boundary conditions for displacements/stresses. They enforce periodic displacement/stress field solutions. Periodic boundary conditions are especially useful in numerical calculations when the asymptotic behavior for infinite systems is of interest. We will encounter them in Sec. 3.

For most physical fracture problems the boundary conditions for the displacements are von Neuman conditions, or mixed boundary conditions. In presence of an additional flow equation boundary conditions for the fluid pressure field need to be specified. In our case homogeneous Dirichlet conditions are sufficient, i.e., the hydraulic crack pressure, $`p`$, and an equilibrium pressure, $`p_0`$, for the external boundary.

In Sec. 3.2 and Sec. 4 we describe simulations where the hydraulic crack is ‘loaded’ and propagated by a constant crack pressure. The crack pressure does neither vary in space nor in time. The considered external boundaries (lattice edges) are either free or periodic.

A, for practical purposes, somewhat more realistic boundary condition propagating a hydraulic crack is a constant fluid injection rate, $`q_{in}`$. To keep things tractable we only consider an incompressible hydraulic fluid. The boundary condition for an impermeable crack, Sec. 5.1, is a particular case of the more general permeable crack boundary condition, see below. A constant fluid injection rate $`q_{in}`$ corresponds to a certain driving pressure difference, $`pp_0`$, where $`p`$ denotes the hydraulic pressure within the crack and $`p_0`$ corresponds to an equilibrium pressure. We have used $`p_0=0`$. The hydraulic pressure does depend on $`q_{in}`$ and the crack opening volume. For permeable cracks/materials $`p`$ also depends on the hydraulic losses. In the following we will outline how to determine the hydraulic crack pressure $`p`$ for given injection rate $`q_{in}`$ and crack structure (contour) $`V`$.

Because for an incompressible fluid the conservation of fluid mass directly implies a conservation of fluid volume, it is possible to formulate a continuity equation for the crack opening volume $`V`$,

$$\frac{V}{t}+_V𝐯𝑑𝐀=q_{in}.$$
(7)
Here $`𝐯`$ denotes the Darcy field integrated over the crack surface $`V`$ with surface element $`d𝐀`$. As source term appears the rate of injected (fluid) volume, $`q_{in}`$. Because the Darcy vectors are linear funtionals of the pressure, $`𝐯=k/\mu \text{grad}p`$, and the crack opening volume is a linear fuctional of displacements, $`V=_V𝐮𝑑𝐀`$, and because all considered equations are linear, a scaling of the form $`p=\lambda p^{}`$ implies $`𝐯=\lambda 𝐯^{}`$, $`𝐮=\lambda 𝐮^{}`$ and $`V=\lambda V^{}`$. Note that such a scaling is in general impossible if one considers Eq. (5) instead of the Laplacian, $`\mathrm{\Delta }p=0`$, as flow equation.

There exists an interesting application of Eq.(7) for undercritical, non-growing cracks. With $`V^{}`$ and $`𝐯^{}`$ depending only on the constant crack contour $`V`$ and crack pressure $`p^{}`$ one obtains a first order differential equation for the scaling factor $`\lambda (t)`$,

$$\frac{d\lambda }{dt}_V𝐮^{}𝑑𝐀+\lambda _V𝐯^{}𝑑𝐀=q_{in},$$
(8)
with general solution,

$$\lambda (t)=\lambda _{\mathrm{}}+(\lambda _0\lambda _{\mathrm{}})e^{t/\tau },$$
(9)
using

$$\lambda _{\mathrm{}}=\frac{q_{in}}{_V𝐯^{}𝑑𝐀},\lambda _0=\frac{p_0}{p^{}},\tau =\frac{_V𝐮^{}𝑑𝐀}{_V𝐯^{}𝑑𝐀}.$$
(10)
The important point is that one can choose an arbitrary boundary value $`p^{}`$ within a calculation; the resulting $`p=\lambda (t)p^{}`$ will be unique. However, one needs to calculate the contour integrals in Eq. (10) in order to do the rescaling.

The magnitude $`\lambda _{\mathrm{}}`$ does depend only on the injected fluid flux $`q_{in}`$ and the total fluid flux through the crack contour. It does not depend on elastic properties. The typical time scale $`\tau `$, above which the crack pressure and the crack opening volume tend to approach their asymptotic values, is given by the ratio of crack volume and fluid volume flux through the crack surface. It does not depend on the injected fluid flux. Again we would like to stress that the crack contour is assumed to be constant. However, Eq.(9) is quite general in the sense that all dependencies on the actual crack contour have been implicitly taken into account within the constants $`\lambda _{\mathrm{}}`$ and $`\tau `$. Equation (9) should therefore hold for arbitrarily shaped cracks.

There are two limiting cases for the practical important situation $`\lambda _00`$. For large times, $`t\tau `$, the scaling factor $`\lambda (t)\lambda _{\mathrm{}}`$ becomes constant and so do the crack pressure $`p_{\mathrm{}}=q_{in}\tau (p^{}/V^{})`$ and the crack volume $`V_{\mathrm{}}=q_{in}\tau `$. The other limiting case is for small times, $`t\tau `$, where the scaling factor becomes linear in time, $`\lambda (t)=\frac{q_{in}}{V^{}}t`$ with corresponding pressure $`p=q_{in}\frac{p^{}}{V^{}}t`$ and volume $`V=q_{in}t`$. The latter case also holds for impermeable materials because $`lim_{k0}\lambda (t)=q_{in}/V^{}t`$.

Equations (7) and (9) hold for arbitrary crack contours, i.e., branched or ramified cracks.

Equation 7 has a corresponding discrete update scheme. One just needs to be careful with the time derivative in Eq. (7) because the scaling factors $`\lambda `$ are in general different before and after the time increment, $`\mathrm{\Delta }t`$,

$$\lambda (t+\mathrm{\Delta }t)=\lambda (t)\left(1\frac{D^{}}{V^{}}\mathrm{\Delta }t\right)+\frac{q_{in}}{V^{}}\mathrm{\Delta }t.$$
(11)
Here $`D^{}=_V𝐯^{}𝑑𝐀`$ and $`V^{}=_V𝐮^{}𝑑𝐀`$ denote the crack contour integrals/sums for the Darcy flows and displacements respectively, as obtained from the solution of the boundary value problem at time $`t`$ with boundary value $`p^{}0`$. As soon as the crack growths and its boundary changes one needs to recaculate the flow and the elastic equations in order to obtain the corresponding $`D^{}`$ and $`V^{}`$. The ratio $`D^{}/V^{}`$ being the inverse of a typical time will therefore change during the crack growth process.

We summarize:

1. Choose an initial crack boundary.
2. Set the initial value $`\lambda (0)`$, i.e., $`\lambda (0)=0`$ (closed crack).
3. Solve the pressure equation Eq. (6) for an arbitrary hydraulic pressure $`p^{}`$ (boundary value for the current crack boundary).
4. Insert the corresponding pressure gradients into the elastic equations Eqs. (2)-(4) as inhomogeneities and solve the elastic equations.
5. Calculate from the obtained elastic solution the stress/force distribution $`\sigma ^{}`$.
6. Calculate from the elastic solution the crack opening volume $`V^{}`$ and from the pressure solution the total Darcy flux $`D^{}`$.
7. Update $`\lambda `$ according to Eq.(11).
8. Use the updated stresses/forces, $`\sigma =\lambda \sigma ^{}`$, in the breaking rule (see Sec. 2.6) and break the corresponding element(s).
9. Repeat steps 3-8 until a stopping criterion is fulfilled.

### 2.5 Heterogeneities

The problem of fracture processes in elastic media has undergone much attention in the last years, in particular the role of disorder in such processes.$`^{\text{References}\text{References}}`$ An important question of interest is how a variation in the local properties affects the global elastic properties in comparison to ideal homogeneous media. It was recognized that elastic properties could considerably be changed and therefore interest was attracted how failure mechanisms like fracture, rupture or damage processes are influenced, changed or possibly controlled by disorder.

In the literature are two essentially different ways to take the disorder into account.$`^{\text{References}}`$

The first possibility is to use exclusively the strain/stress field in order to calculate the probability for breaking a local constitutive element. This is known as stochastic disorder. Prominent scalar examples for this kind of disorder are known in the literature as Diffusion Limited Aggregation (DLA), Dielectric Breakdown or Viscous Fingering.$`^{\text{References},\text{References}}`$ We will present results for stochastic disorder in Sec. 4.

A second different kind of disorder is to initially assign a statically set of random values to certain material specific constants (structural disorder). Once these constants like strength, bending thresholds, elastic properties and others have been chosen the system behaves in a fully deterministic way. Corresponding results are presented in Sec. 5.

### 2.6 Breaking rules

Lattice fracture models like the beam model are spatially discrete. This is because one would like to have direct access to each constitutive elastic element. The medium is represented by a lattice and nearest neighbouring lattice points are connected by constitutive elastic elements called bonds. The bonds represent interactions not on an atomistic scale but rather on a mesoscopic scale, for example the elastic coupling between two adjacent grains where the grains themselves are linear elastic bodies. New features are the definitions of physically motivated ”selection” and ”breaking” rules for lattice bonds. Either a bond is broken according to some breaking rule and stays broken for the future or it remains intact and is further susceptible for the breaking process. The simulation ends when the elastic lattice fully breaks apart. Moreover, one has to define in the selection rule those bonds which could break in principle at a given stage of rupture. In single crack models only the surface bonds of the crack are allowed to break. However, in the many crack version in general all bonds are eligible. This allows in principle the nucleation of several microcracks in the system. We have used for our simulations the first kind of selection rule. In a recent published work$`^{\text{References}}`$ it was shown that if the disorder for the breaking thresholds is not to weak there are two different regimes in the macroscopic breaking characteristic. For small macroscopic strains first the weakest bonds are broken without leading to catastrophic failure. The breaking process is stable in this region and one has to increase the external applied stress in order to proceed breaking bonds. In finite systems however, this disorder controlled regime must end at the breaking point. The breaking process localizes, becomes unstable and a single large crack controlls the entire breaking. It was shown that for different probability distributions of breaking thresholds and different loading conditions two nontrivial possibly universal exponents are sufficient to rescale the macroscopic breakdown characteristics with respect to the system size. Breaking rules which determine the breaking conditions are crucial in lattice models because one has to put additional information about the fracture mechanism into the model which could not be obtained from elasticity theory. The employed breaking rules are described in the corresponding sections Sec. 3, Sec. 4, and Sec. 5.

## 3 Fracturing without disorder

We investigate both the case of uniaxial loading due to remote tensile forces and the loading due to a hydrostatic pressure acting from the interior on the crack surface. The former situation is typically encountered when a material is pulled at its outer faces, while the latter one can be found in engineering applications of hydraulic fracturing.$`^{\text{References}}`$ We will neglect inertial forces, assuming a sufficiently slow fracture process. The cohesive properties of the beams are considered to be homogeneous, i.e. the beams behave linear-elastic up to their theoretical strength $`F_{th}`$ (cohesion force), above which they irreversibly break (zero elastic constants). At the beginning of each simulation we break one vertical beam located at the center of the lattice. This is the initial crack. Corresponding to the employed boundary conditions described below we calculate the internal force distribution $`F_{ij}`$ of beams connecting sites $`i`$ and $`j`$ using a conjugate gradient method.$`^{\text{References}}`$ From this longitudinal forces we determine and break the most over stressed beam, i.e. the beam carrying the highest force, $`F_{i_0j_0}=\mathrm{max}_{\{ij\}}F_{ij}`$. The maximum force allows us to calculate a scaling factor $`\lambda =F_{th}/F_{i_0j_0}=\sigma _{th}/\sigma _{i_0j_0}`$ which in turn is used to determine the macroscopic breaking stress $`\sigma _c`$ or breaking pressure $`P_c`$ (measured in units of $`\sigma _{th}`$) necessary to fulfill the local fracture criterion $`\sigma _{i_0j_0}=\sigma _{th}`$. Such scaling is possible due to the linearity of the employed elastic equations. Breaking the beam connecting sites $`i_0`$ and $`j_0`$ destroys the balance of forces at those sites and one has to relax the system to its new equilibrium configuration. Then the above steps are repeated until the lattice breaks apart.

It is well known that the highest stress enhancement factors occur at the crack tips for the pressure as well as for the tensile problem. Because we consider only homogeneous cohesion forces the terms ‘most over stressed’ and ‘highest stress’ are equivalent and the cracks grow linear in direction of the highest local tensile force. In our case the direction of crack growth is the x-direction because of the uniaxiality of imposed boundary conditions, see below.

### 3.1 Tensile fracturing

To impose an external strain we attach at the bottom of the lattice a zeroth line on which for all sites $`x_i=y_i=\phi _i=0`$ are fixed, and on the top we attach a $`(L+1)`$st line on which all sites have the same fixed values $`x_i=\phi _i=0`$ and $`y_i=1`$. With this boundary conditions the external displacement in y-direction is fixed to unity (Dirichlet boundary conditions) and one can imagine them as being represented by rigid bars attached at the top and the bottom of the lattice. For the boundary conditions in x-direction we consider periodic as well as free lattice boundaries. The externally applied stress, $`\sigma _0`$, necessary to maintain a unity displacement between the two rigid bars can be easily calculated by summing up all $`y_i`$ displacements over the first line, $`\sigma _0=\frac{\alpha }{AL}y_i`$, where $`\alpha `$ is the longitudinal force constant, $`A`$ the cross section of the beam and $`L`$ the (dimensionless) lattice size. Employing the above defined scaling factor $`\lambda `$ we express the externally applied breaking stress, $`\sigma _c`$, in the form,

$$\sigma _c=\lambda \sigma _0,$$
(12)
or in dimensionless form,

$$\frac{\sigma _c}{\sigma _{th}}=\frac{1}{F_{i_0j_0}}\frac{\alpha }{L}y_i.$$
(13)
With this notation the external breaking stress $`\sigma _c`$ is just the stress necessary to yield the local breaking condition $`\sigma _{i_0j_0}=\sigma _{th}`$.

We will use Eq.(13) for our finite-size scaling purposes. It expresses the breaking stress, measured in units of cohesion stress, in terms of the beam force $`F_{i_0j_0}`$ at the ‘crack tip’ and the average force $`\frac{\alpha }{L}y_i`$ on a beam due to the global unit displacement. In general both terms depend on the number of broken beams $`N`$ (crack size) and on the lattice size $`L`$, see Sec. 3.3.1.

### 3.2 Hydraulic fracturing

Experimentally hydraulic fracturing is encountered when an incompressible fluid is injected under high pressure into an existing crack. The characteristics of this particular loading mode is that the loading of the crack only happens on the crack surface itself, but not on remote surfaces like in the case of tensile fracturing. Recently hydraulic fracturing has been numerically addressed by introducing for each broken beam a force dipole of strength $`F_0`$ simulating an inner pressure $`P_0=F_0/A`$.$`^{\text{References},\text{References}}`$ In the present section we will follow this implementation. The basic assumptions are a) the pore fluid is highly compressible (gas), b) the permeability is negligible and c) the injected fluid is incompressible. Starting from one vertical broken beam (zero elastic constants, one vertical force dipole) we calculate the internal stress distribution $`F_{ij}`$ using the conjugate gradient. Calculating the scaling factor $`\lambda `$ as described above we determine and break the beam carrying the highest stress, $`F_{i_0j_0}`$. Thereafter the lattice has to be relaxed again (now two broken beams, two force dipoles) and one repeats the above mentioned steps until the lattice breaks apart. The crack growth is single connected. As always the two vertical beams in front of both crack tips are most overstressed and broken. Hence the existence of a connected path for the liquid is garantued.

We have performed calculations for free external boundaries in x- and y-direction as well as for periodic boundaries in x-direction and free boundaries in y-direction. As we only impose stresses the elastic solution (displacement field) is unique up to a general translation and rotation. Without restriction we therefore fix the displacement of an arbitrary lattice site to zero. Analogous to Eqs. (12) and (13) one has for the breaking pressure,

$$P_c=\lambda P_0,$$
(14)
and

$$\frac{P_c}{\sigma _{th}}=\frac{F_0}{F_{i_0j_0}},$$
(15)
respectively. It should be noted that $`\sigma _{th}`$ is the tensile strength and not the compressional strength. In fact the stresses at the tips are tensile stresses as the crack is hydraulically opened.

We will use Eq. (15) for the finite-size scaling of the hydraulic problem.

### 3.3 Results

As explained in Sec. 3.1 and 3.2 we obtain from our simulations the macroscopic breaking stresses and pressures, $`\sigma _c(N,L)`$ and $`P_c(N,L)`$, measured in units of the theoretical strength $`\sigma _{th}`$, in terms of the number of broken beams $`N`$ and the employed lattice sizes $`L`$. As we are interested in the finite-size scaling properties of the breaking stresses ($`\mathrm{\Sigma }_c=\sigma _c`$) and pressures ($`\mathrm{\Sigma }_c=P_c`$) we introduce two new functions, $`\mathrm{\Psi }(N)`$ and $`\mathrm{\Phi }(N/L)`$,

$$\mathrm{\Sigma }_c(N,L)=\sigma _{th}\mathrm{\Psi }(N)\mathrm{\Phi }(N/L).$$
(16)
The function $`\mathrm{\Phi }(N/L)`$ describes the finite-size corrections to the scaling $`\mathrm{\Psi }(N)`$ between the breaking stress/pressure and the crack length for given boundary conditions in the infinite system. The quantity $`N/L`$ is for our purposes the most suitable scaling variable as it represents the length of a linear crack measured in units of the employed lattice size. In general both functions, $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ will depend on the employed boundary conditions. The dependence of the asymptotic scaling $`\mathrm{\Psi }`$ on employed boundary conditions appears to be contrary to what is expected from the theory of critical phenomena. In fact the asymptotic scalings $`\mathrm{\Psi }`$ for breaking stresses and pressures of continua without internal length scale are identical. However, as we will see below this is no longer necessarily true for systems with internal length scale. One obtains rather finite-size corrections originating from this scale which, depending on the employed boundary conditions, can result in very broad, non universal transient regions.

In order to show numerically that the finite-size scaling function $`\mathrm{\Phi }(N/L)`$ exists, one needs to know (or to conjecture) the explicit form of $`\mathrm{\Psi }(N)`$, i.e. the asymptotic behavior of $`\mathrm{\Sigma }_c`$ as $`L\mathrm{}`$.

We will access the latter information from simple continuum mechanical considerations of symmetric elasticity.

#### 3.3.1 Tensile fracturing

Previous continuum mechanical scaling considerations have successfully been applied to fracture models. It has been shown$`^{\text{References}}`$ that the breaking stress $`\sigma _c`$ for tensile fracturing of a finite, periodic beam lattice follows very well the known ‘Tangent-Formula’ of fracture mechanics $`^{\text{References},\text{References}}`$,

$$\frac{\sigma _c}{\sigma _{th}}N^{1/2}\sqrt{\frac{N}{L}\mathrm{cot}\frac{\pi }{2}\frac{N}{L}}.$$
(17)
By comparison with Eq. (16) we see that the asymptotic scaling for the breaking stress $`\sigma _c`$ is given by $`\mathrm{\Psi }N^{1/2}`$, a result that originally was found by Griffith. In Fig. 5 we show the finite-size scalings for the numerically obtained tensile breaking stresses (a) for free boundary conditions and (b) for periodic boundary conditions. For both curves the data collapse is quite acceptable. While for the periodic problem the finite-size correction $`\mathrm{\Phi }(N/L)`$ follows Eq. (17)$`^{\text{References}}`$, no analytical expression for $`\mathrm{\Phi }`$ is known in the case of free boundary conditions. However, the differences for the breaking stresses at given values of $`N/L`$ are less than $`0.1\sigma _{th}`$, see Fig. 5, and Eq. (17) represents a reasonably good approximation for finite solids for many practical purposes.

#### 3.3.2 Hydraulic fracturing

In general analytical solutions of crack problems are difficult to obtain and even for two-dimensional problems in infinite, semi-infinite or periodic domains one is often faced with complicated integral equations. The reader interested in this is referred to Muskhelishvili$`^{\text{References}}`$ and Sneddon$`^{\text{References}}`$.

It has been shown analytically that the breaking pressure $`P_c`$ of a semi-infinite two-dimensional continuum with periodic boundary conditions is given by Eq. (17) where one has to replace $`\sigma _c`$ by $`P_c`$.$`^{\text{References}}`$ We show in Fig. 6 the analogous plot to Fig. 5 for the breaking pressure. One does not obtain a data collapse, indicating that $`\mathrm{\Psi }N^{1/2}`$, see Eqs. (16) and (17), is not the scaling of the breaking pressure in the asymptotic limit of infinite lattices.

This result is somewhat unexpected, because the continuum limit predicts identical scaling relations for breaking stresses and pressures and on the other hand we find good agreement between lattice and continuum scalings for the tensile problem.

However, the proposed model for hydraulic fracturing is defined on a lattice and not in a continuum. As lattices always show additional structure, higher order gradient terms of the continuum displacement field appear in the elastic solution, which have no counterparts in simple continua.$`^{\text{References},\text{References},\text{References}}`$

In the following we present an argument why the lattice finite-size scalings for the pressure and tensile problem are different.

In the tensile problem the loading forces are acting at the lattice boundaries remote to the stress free crack surfaces. Therefore local stresses close to the crack tips are carried through nearly the whole elastic volume. Contrary to this in the hydraulic problem loading only happens due to force dipoles acting at the crack surfaces close to the crack tips. One therefore might expect that the breaking pressures $`P_c`$ show a more pronounced ‘(micro structure) lattice behavior’ than the breaking stresses $`\sigma _c`$. For very large cracks these differences should diminish as the discrete loading approaches a force density.

In order to extent a linear crack of length $`N`$ in an infinite continuum one has to impose a breaking pressure $`P_cN^{1/2}`$. However, employing a single double force of magnitude $`F_c`$ one finds, $`F_cN^{1/2}`$ for the critical loading force (related Boussinesq problem).$`^{\text{References}}`$

We argue now that for a single crack in an infinite disk loaded by equidistant force dipoles the asymptotic scaling for $`\mathrm{\Psi }`$ is equivalent to that of an infinite lattice. Knowing $`\mathrm{\Psi }`$ we plot the finite-size scalings for the investigated lattices.

The asymptotic continuum scaling for the breaking pressure (as calculated from equidistant force dipoles) is most easily obtained using the complex stress function of Westergaard.$`^{\text{References},\text{References}}`$ We find for the asymptotic scaling in an infinite disk,

$$\mathrm{\Psi }(N)\frac{\pi }{2}\frac{N^{1/2}}{1+Nf(N)},$$
(18)
with

$$f(N)=\underset{k=1}{\overset{(N1)/2}{}}\frac{1}{\sqrt{(\frac{N+1}{2})^2k^2}},N\text{ odd}.$$
(19)
The sum $`f(N)`$ converges rapidly towards $`\mathrm{arcsin}(\frac{N1}{N+1})`$ and for large $`N`$ towards $`\pi /2`$. Hence in the limit of a large number of force dipoles we obtain the aforementioned asymptotic Griffith scaling for the breaking pressure, $`P_c/\sigma _{th}=N^{1/2}`$.

In Fig. 7 we show the resulting finite-size scalings for the breaking pressure, based on Eqs. (16) and (18). In this plot a reasonable scaling is obtained, which proves our earlier statement that the asymptotic limits ($`L\mathrm{}`$) of the breaking pressures for continua and lattices are different.

Comparing Fig. 7 with Fig. 5 one also notes a significant deviation from the scaling form Eq. (17), demonstrating a modified correction to scaling behavior.

## 4 Hydraulic Fracturing with stochastic disorder

We briefly discussed in the introduction, Sec. 1, Diffusion Limited Aggregation and its associated fractal aggregates. Essentially DLA clusters growth by adding one particle per “simulation step” to the clusters surface with a probability, $`p_{ij}`$, proportional to the local field gradient $`|P|`$ obtained from Laplace’s equation, $`\mathrm{\Delta }P=0`$. Similar considerations apply for the quasistatic growth of electrical discharge patterns, where an elementary discharge event is assumed to happen with a probability proportional to the local electrical field.$`^{\text{References}}`$ Typically the boundary conditions are Dirichlet conditions for the external and the cluster surface.

There have been attempts in literature to generalize the DLA model in several ways. A main proposition was the introduction of a growth exponent, $`\eta `$, describing how strong the growth probabilities are correlated to the gradients, $`p_{ij}|P|^\eta `$. For $`\eta \mathrm{}`$ the probability distribution becomes arbitrarily narrow (a delta peak) and the model is completely deterministic, i.e., the cluster grows always at the surface location of highest field gradient (current). The resulting structures are usually linear in two dimensions. The case $`\eta =1`$ corresponds to fractal DLA cluster growth, as already mentioned above. For $`\eta =0`$ the cluster growth is completely controlled by fluctuations and independent of the field solution. The corresponding structures are usually spatially compact.

It is important to note that the growth rule and in particular the employed growth exponent $`\eta `$ is completely independent of the investigated field equation. Rather the growth rule bases on plausible considerations.

There has been subsequently work been done to apply stochastic growth rules onto elastic systems, especially applying them to crack growth problems. The “Gradient” in the growth rule is then typically replaced by the hydrostatic strain (eigenvalue) component. For certain classes of tensile problems fractal crack growth was found on two dimensional Hookean spring networks, exhibiting a fractal dimension $`d_f1.6`$ for $`\eta =1`$.$`^{\text{References},\text{References},\text{References}}`$ This value for $`d_f`$ is relatively close to the one observed in a DLA problem. However, the first author of this review conducted simulations on two dimensional beam lattices (see Sec.2.2) where a small preexisting crack was loaded at constant pressure from inside, and found non-fractal rather compact crack growth, already for $`\eta =1`$.$`^{\text{References}}`$ It was concluded that compressed microstructural elements (beams) should not be considered for the breaking process because of the very asymmetry of interatomic potentials. With this modified breaking rule also non-fractal cracks were found for $`\eta =1`$. The cracks were, however, now linear in shape. In order to obtain in a hydraulic fracture simulation persisting branching cracks, much stronger fluctuations need to occur than obtainable with $`\eta =1`$. This is not so surprising as almost the whole elastic bulk is in a compressed state, while just around the crack tips there exist very narrow cones of elements being in a stretched state. In Fig. 8 we show the crack of a single, very large simulation for the case of maximum fluctuations, i.e., $`\eta =0`$ using the above mentioned asymmetric breaking criterion. Note that the breaking probability does not depend on the strain magnitude, however, all compressed beams are assigned a breaking probability identically to zero. The hydraulic pressure inside the crack was kept constant and the external boundaries are periodic in displacements, i.e., spanned on a torus. Furthermore the system is assumed as impermeable and there are no losses of the hydraulic fluid. The resulting crack structures are scaleinvariant (fractal) as can be seen in Fig.9. The measured fractal dimension $`d_f=1.56\pm 0.05`$ is slightly below the numerical value for tensile problems, and slightly above the experimental value $`d_f=1.41`$ obtained from hydraulic fracturing of clay pastes in two dimensions, compare Fig.3.

The total crack opening volume of the ramified cracks has been calculated and monitored during the growth, and plotted versa the average crack radius. This is shown in Fig.10. The crack opening volume for fractal cracks shows the same dependency on (average) cracks size in comparison to straight cracks, namely $`VR^2`$.

We note that the foregoing considerations and results are for impermeable elastic materials, i.e., the influence of possible effects due to hydraulic fluid flow (hydraulic losses) have been neglected.

## 5 Hydraulic Fracturing with structural disorder

The concept of a local cohesion strength has been used in a number of papers.$`^{\text{References}}`$ One assumes that a deformed elastic element connecting sites $`i`$ and $`j`$ breaks above a certain material specific threshold force $`F_{th}^{ij}`$ (’cohesion strength‘). This kind of disorder is also called deterministic or quenched disorder in the literature, because once the breaking thresholds have been chosen the model behaves in a completely deterministic way, contrary to the model described in the foregoing section. However, we prefer the term ‘structural disorder’ as a spatial variation in material properties is frequently associated with certain micro-structural particularities within a material. If the inner stresses are above their thresholds the beams are broken and removed, i.e. their elastic moduli are set to zero. Since the cohesional strength for compression is much higher than for tension we assume that compressed beams cannot break.$`^{\text{References}}`$ The breaking thresholds are initially distributed randomly according to some probability density function, $`\rho (F_{th})`$. At this point it would be ideal to employ empirical strength distributions of the considered material for $`\rho (F_{th})`$. However, at the present stage of research we are basically interested in generic features of such crack growth models. The width of the strength distribution is most easily controlled employing power-law distributions,

$$\rho (F_{th})F_{th}^r,$$
(20)
with $`F_{th}[0,F_{max}]`$ and $`r>1`$. Negative exponents $`r`$ are used to describe strong cohesive disorder (broad distributions) while large positive exponents correspond to weak disorder (narrow distributions). It is convenient to express the normalization factor and $`F_{max}`$ by the distribution’s expectation value $`F_{th}`$ and the exponent $`r`$. This allows to investigate disorder effects due to deviations around $`F_{th}`$ for constant $`F_{th}`$.

In the following we will discuss the modelling of crack growth in impermeable and permeable materials separately because the employed breaking criteria slightly differ for both cases.

### 5.1 The impermeable medium

The hydraulic fracturing of an impermeable solid represents a limiting case due to the vanishing hydraulic losses. For a porous material saturated by an incompressible fluid this would imply that only pure shear deformations would be possible. This follows directly from the time-dependent flow equation, i.e., for $`k0`$ or $`\mu \mathrm{}`$ Eq. (5) reduces to the condition $`\text{div}𝐮=\text{const.}`$ preserving the volume of a deformation. This is not a very typical experimental situation. However, if the pores are filled by very compressible gases one might even neglect the deformation dependent volume forces for impermeable solids.

In the following we consider the case where an incompressible fluid hydraulically propagates a crack through an impermeable material. We require that the pores are filled by a very compressible substance compared to the solids compressibility. Such situation is approximately encountered for magma driven cracks.$`^{\text{References}}`$

We consider a constant fluid injection rate, $`q_{in}=510^2l^2d`$, into the hydraulic crack and a vanishing pore pressure, $`p_0=0`$. The boundary conditions have been discussed in Sec. 2.4, i.e., see Eq. (11) for vanishing Darcy losses, $`D^{}=0`$. The simulation procedure follows the mentioned steps 1-9 in Sec. 2.4 with the exception that in the present case no pressure equations need to be solved. Instead the ‘test pressure’ $`p^{}`$ does appear directly as inhomogeneity for the elastic equations along the crack contour. Initially one small crack (one broken vertical beam) is prepared on the grid center. For the distribution of breaking thresholds, Eq. (20), the values $`F_{th}=10^2`$ (arbitrary units), $`r=0.5`$ and $`r=0.7`$ (strong disorder), and $`r=+\mathrm{}`$ (no disorder) are used. After each time increment, $`\mathrm{\Delta }t=1`$, all beams carrying tensile forces larger than the associated breaking thresholds are broken, i.e., their elastic moduli are set to zero.

For $`r=+\mathrm{}`$ (no disorder) straight cracks are observed and the hydraulic crack pressure drops in form of a power-law in time, $`P(t)t^{1/3}`$.$`^{\text{References}}`$

For strong disorder the situation becomes more complicated. In Fig. 11 we show a calculated sample crack. The obtained crack is a branching structure exhibiting many small scale branches but only a very few branches on the larger length scale. A comparison with the crack structure shown in Fig. 8 reveals a much lower crack density in the present case. This is directly measured in Fig. 12 by plotting the number of broken beams $`N`$ in dependence of the average crack size $`R`$,$`^{\text{References}}`$ yielding the cracks fractal dimension $`d_f`$ using the relationship $`NR^{d_f}`$. The fractal dimensions for the investigated threshold distributions $`r=0.7`$ $`()`$ and $`r=0.5`$ $`(+)`$ take the values $`d_f=1.44\pm 0.10`$ and $`d_f=1.39\pm 0.10`$ respectively. These values for $`d_f`$ are consistent with the fractal dimension found in the two dimensional experiments on hydraulic fracturing of viscoelastic clays.$`^{\text{References}}`$

The time record of the hydraulic pressure at the injection point is of considerable interest because on the one hand it is often experimentally measurable and on the other hand does it contain information about rheological and breaking properties. If the considered crack is at rest the hydraulic pressure $`P`$ will generally build up in time (loading). For an impermeable medium the pressure build up will correspond to the materials elastic response on the injected fluid (and possible fluid flow effects within the crack which are neglect in the present case). We consider the fluid as incompressible implying a linear relationship between crack volume and time, and due to the linearity and time independence of employed elastic equations a linear increase of pressure in time is the consequence. The associated slope is related to the materials mechanical stiffness (for the present crack). However, if the crack starts to move the stiffness decreases and so does the hydraulic pressure. There exists a cutoff for this pressure drop because the materials average cohesion is non-zero. Therefore the crack comes to rest again after some time period. Figure 13 displays such an ‘erratic’ pressure record. The strongest pressure fluctuations occur for small cracks (early times) at high frequency whereas lateron the loading periods become larger and larger. This is due to the rapidly decreasing stiffness as soon as the crack advances. The loading periods are periods of ‘quiescence’ (no breaking); the unloading periods correspond to ‘bursts’ (temporally localized groups of breaking events). This is exemplified in Fig. 14. The temporal correlations of such bursts have been analyzed in detail in Ref. References.

### 5.2 The permeable medium

This section considers hydro-fracturing of a permeable porous medium, where an incompressible fluid is injected at the center of a quadratic grid at a constant rate, $`q_{in}`$. The boundary conditions have been discussed in Sec. 2.4.

The dimensionless formulation of the pressure equation is first explained before we show some examples. The considered pressure equation Eq. (6) is made dimensionless by considering the forces acting on the beams surrounding the injection cell at the center of the grid. Darcy’s law yields the Darcy velocity $`v=(k/\mu )dp/dx`$ in the $`x`$-direction. The Darcy velocity is also $`v`$ in the $`y`$-direction when assuming symmetry. The forces acting on the beams surrounding the injection cell are $`F=l^2ddp/dx`$, because the pressure difference between two neighbor cells is $`dp=(dp/dx)l`$. The length of a beam is $`l`$ and the width of a beam is $`d`$. The beam force $`F`$ can be written $`F=\mu vl^2d/k`$, when Darcy’s law is used to replace the pressure gradient $`dp/dx`$ with the Darcy velocity $`v`$. The volume flux into the center cell is initially $`q_{in}=vl^2`$. The force $`F`$ is written in terms of $`q_{in}`$ rather than $`v`$, and we have $`F=\mu q_{in}d/k`$. This force is scaled with the unit force $`F_0=Ed^2`$, (where $`E`$ is Young’s modulus), which can be written $`F/F_0=q_{in}/q_0`$. The unit injection rate becomes $`q_0=kEd/\mu `$.

This scaling shows two things. Firstly, it shows that a dimensionless injection rate of order one or larger implies strong enough forces on the center cell for deformations on grid to be noticeable. Secondly, the (dimensionless) average beam strength must be less than $`q_{in}/q_0`$ for a crack to nucleate.

The characteristic time becomes $`t_0=\mu l^2/(Ek)`$, which will be as low as nano seconds for a grid with beam length $`l=10^5`$ m, permeability $`k10^{15}`$ m<sup>2</sup>, Young’s modulus $`E10^{11}`$ Pa, and fluid viscosity $`\mu 10^3`$ Pa s. We consider a constant permeability throughout the material.

The full pressure equation for an incompressible fluid Eq. (5) is not solved. The pressure equation is therefore not fully coupled to the volume changes in both the bulk and the fracture. However, it is still possible to couple the hydraulic pressure to the elastic deformations of the crack by considering conservation of fluid, as shown in Sec. 2.4. The (Laplacian) pressure equation is first solved for a unit pressure $`p^{}`$ on the crack contour. The total Darcy flux $`D^{}`$ is determined from this solution; the corresponding crack volume $`V^{}`$ is computed from the associated elastic solution. The scaling parameter $`\lambda (t)`$ yielding the physical pressures and strains is computed each time step according Eq. (11), see Sec. 2.4 steps 1-9.

The beams on the crack surface are checked after each time step to see if any of them are stretched beyond their strength. If there are any over-stretched beams, then the most over-stretched beam will be broken. We showed in the foregoing section that the advance of fracture is localized in time in form of bursts, separated by large time periods of quiescence (loading). The time-step used during the bursts is therefore “small”. On the other hand, the quiet periods between two bursts can be very long compared to the short time span of a burst. The time-step is therefore multiplied by a factor larger than $`1`$ after time-steps where no beams were broken. This factor could be chosen so that ten such steps together yield a combined factor equal to for instance $`2`$ or $`10`$. Increasing each time-step between the bursts by a factor makes it possible to obtain a step length sufficiently long to get through these relatively long periods of pressure build-up. However, the time-step is immediately set to its minimum value once a new beam breaks.

Three numerical examples are studied, which differ in the average beam strength $`F_{th}`$. The average beam strength in these cases and the injection rate $`q_{in}`$ are given in the Table 1.

All three cases have the grid size $`100\times 100`$ and a beam strength distribution characterized by the exponent $`r=0.9`$, compare Eq. (20).

A value for $`r`$ close to $`1`$ indicates that the beams have very broadly distributed strengths. The employed material and geometry dependent constants are $`a=1`$, $`b=2.5`$ and $`c=1200`$, when expressed dimensionless. (The scaling of these constants is such that $`a`$ becomes $`1`$.)

The crack in case I, where $`q_{in}/q_0F_{th}/F_0`$, is shown in Fig. 15a. This figure shows the crack after it stopped growing. The hydraulic pressure is shown in Fig. 15b, where it is seen that the pressure reaches a stationary value. A small crack nucleated in this case, even though the average beam strength is $`10`$ times larger than the injection rate. This is due to the strong heterogeneities ($`r=0.9`$). The pressure follows an exponential rise in time according to the expression for the scaling parameter $`\lambda `$ in Eq. (9). The first smooth rise of the pressure ends by breaking beams until no more beams are overstressed, and a second rise of pressure follows. This second rise leads to a stationary pressure. Notice that the first rise of pressure points towards a higher stationary pressure than the second rise. This is expected because the maximum possible pressure in a constant crack geometry, given a constant injection rate, is decreasing with increasing crack size. A large crack requires a lower stationary hydraulic pressure compared to a small crack in order to balance the injection rate to the hydraulic losses. The small crack seen in Fig. 15a is almost entirely due to the first burst. After the first burst, the pressure cannot reach a sufficiently high level any more to break more beams. The time $`\tau `$ (see Eq. (10)) is plotted in Fig. 18a. Because $`\tau `$ measures the typical elastic response versus the typical flow response, it defines the characteristic time scale below which the (re)loading process appears to happen within an impermeable medium. It is seen that $`\tau `$ is increasing in crack size, which is consistent with the pressure plot Fig. 15b. The characteristic time increases with the crack size even though the stationary pressure is decreasing.

The case II shows crack propagation for an average beam strength equal to the injection rate, see Table 1. The forces from the fluid gradients on the beam structure are then comparable with the strength of the beams. The pressure during the cracking process is shown in Fig. 16b, and it is seen that the pressure is rising almost linearly between the bursts of beam breaking. The small initial bursts were sufficient to generate a crack with a time constant $`\tau `$ much larger than the time span shown in the pressure plot. This explains why the pressure rise appears almost linear rather than exponential in time as given by Eq. (9). The characteristic time $`\tau `$ is shown in Fig. 18b. The crack structure at the end of the time period in the pressure plot (Fig. 16b) is shown in Fig. 16a. The crack in this case of broadly distributed beam strengths shows the growth of a branching crack structure. Note that within the crack ‘fjords’ the Darcy flow almost vanishes.

Case III has an average beam strengths $`F_{th}/F_0`$ one order of magnitude lower than the injection rate $`q_{in}/q_0`$. The injection rate is therefore high compared to strength of the structure. This is also seen from the pressure plot in Fig. 17b, which is characterized by more frequent bursts compared to case II. Only a slight pressure increase is enough to trigger a new burst. The hydraulic pressure level of beam breaking in Fig. 17b is about one order of magnitude lower than the corresponding pressure level in figure Fig. 16b. This is as expected from the average beam strengths of these cases. The characteristic time $`\tau `$ for case III is shown in Fig. 18c, which also shows a large number of bursts. The crack in its final state is shown in Fig. 17a. A branching crack structure is seen, where the branch close to the boundary seems to be growing faster than the other branches. This seems to be a boundary effect because the pressure gradients close to a boundary are larger. This effect favourizes growth of cracks close to the boundaries.

The presented cases demonstrate several features of hydraulic fracturing in a permeable medium with broadly distributed cohesion, when a fluid is injected at a constant rate at the center of the test sample. They clearly show how the hydraulic pressure is increasing during “quiet” periods. The pressure increases according to Eq. (9). The characteristic time $`\tau `$ in Eq. (10) is computed, and it is seen that $`\tau `$ becomes large compared to the time period between successive bursts. The pressure therefore appears to rise linearly between the bursts. The characteristic time $`\tau `$ is also increasing with the crack size, which implies that the time periods between successive bursts get longer as the crack growths. This can also be seen from the slopes in Figs.15b-17b which are decreasing with increasing crack size. The pressure drops during the bursts. Furthermore, the bursts also become more frequent when the average beam strength $`F_{th}/F_0`$ gets low compared to the injection rate $`q_{in}/q_0`$.

## Conclusion

We have reviewed some numerical results for hydraulic fracturing in porous materials employing microstructural fracture growth models. The considered models are not finite difference schemes of certain continuum mechanical equations. Rather they implement an approach to crack growth problems in the spirit of Monte Carlo methods of statistical mechanics. This allows to track quasistatic crack growth under very general conditions, i.e., for disordered systems on mesoscopic length scales.

## Acknowledgments

F. T. would like to acknowledge financial support from CEC under grant number ERBFMBICT 950009.


text:
# Test Beams and Polarized Fixed Target Beams at the NLC

## I Introduction

During the last 10 years, colliding beam physics with the SLC/SLD was SLAC’s main program. During most of the SLC colliding beam running, SLAC was also engaged in a highly successful fixed target program. The main emphasis was on the spin content of the nucleon (E-142 slace142  and E-143 slace143  at 30 GeV, and E-154 slace154  and E-155 slace155  at 50 GeV), but other fundamental experiments like E-144 also ran slace144 . Typically two months per year were set aside from the main program. Test beamscavalli93 ) were used parasitically for smaller experiments (E-146, LPM ) and for detector development, in particular GLAST glast96 .

This report discusses the possibilities for using the NLC nlczdr  primary electron beam for fixed target physics or to produce secondary hadron test beams. While dedicated operation of a fixed target program is obviously the most efficient, we will assume that test beam or fixed target running is done in parallel with colliding beam physics. Other design conditions are that there must be permitted access to the fixed target and test beam experimental building during colliding beam operation, there should be useful hadron and electron yields up to about 80% of the primary beam energy, and a single beam transport line and experimental area can be used for both tests and fixed target running.

Figure 1 shows possible locations of such test beam and fixed target beam facilities: they could be at the end of the electron linac (dedicated running only), or behind either IP following the NLC extraction line. The scheme shown is based on the current NLC design of a low energy IR (up to 500 GeV center of mass (cms), i.e., 250 GeV per beam), and a high energy IR (up to 1 TeV cms) expandable to higher energies. To keep the option for multi-TeV beams open, the beam path to the high energy IR has minimal bends; this requires that the linacs have a small angle of $`20`$ mrad between them to enable the outgoing bunches to be separated from the incoming beam.

Figure 2 shows a schematic layout of the NLC low energy IR extraction line, including the 5 MW water beam dump, and the proposed test beam and fixed target take-off. For colliding beam physics the beam dump intercepts both the disrupted, primary electrons and the beamstrahlung $`\gamma `$’s. For test beams, a 0.3 radiation length Be target is placed in the disrupted beam for production of secondary e-, $`\pi ^{}`$, K-, p̄ (or their anti-particles). To extract the primary beam for fixed target physics, a dipole is placed at the location of the Be target to bend the beam into the secondary transport line. The dipole aperture must be large enough so as not to intercept the beamstrahlung photons. The 5 MW dump (a larger dump will be required for the high energy IR) does double duty as a collimator of the disrupted beam energy tail, see Section II-B. Some design requirements of the transport beam line to the test beam/fixed target experimental area are: 1) able to transport the full primary beam energy, 2) able to collimate the disrupted beam energy tail within 1% of the undisrupted energy, 3) no energy dispersion in the experimental area, 4) magnet apertures sized to allow useful secondary yields, 5) zero net bend so that the longitudinal polarization is preserved no matter what beam energy is being used in NLC, 6) the experimental area has a sufficient transverse displacement from the straight-ahead NLC dump line so that the muon flux from the 5 MW dump is not a problem.

## II Test Particle Production and Tail Collimation

### II.1 Test Particle Production

Figure 3 shows simulations fluka  of pion and proton yields at production angles of 0.25, 0.5 and 1.0 for a 250 GeV electron beam of intensity 1012/pulse train on a 0.3 rl Be target in a secondary beam line of acceptance 4 $`\mu `$sr and $`\mathrm{\Delta }P/P`$ = 4%. The curves show that the yield at 1 is large enough for lower energies, but that the maximum energy of particles is a function of angle, so smaller angles give more latitude in choosing the right yield at higher secondary energies. Kaon yields are somewhere in between. Not enough test particles were used in these simulations to determine if enough yield is produced at an angle of 1; this will be investigated at a later time.

### II.2 Fixed Target Experiments

Figure 4 shows the concept of a combined test beam and fixed target facility. The beams after the dump are bent and sent through a system of slits. These slits will be designed to define the momenta of the particles for the test beam application and to clean up the remainder of the tails of the disrupted beams for fixed target experiments.

Figure 5 shows the energy distribution of the 250 GeV disrupted beam and the disrupted power in 1% bins, calculated using the beam-beam program GUINEAPIGguineapig . In order to achieve a $`\mathrm{\Delta }E/E`$ of 1%, the quality required for fixed target experiments (see below), 1.6 MW have to be collimated in the transport line to the experimental area. As an existence proof, the energy defining slits (SL-10) in the SLAC A-line are designed for such power. In the NLC case, not all power has to be absorbed by the energy defining slits because some will be absorbed by the 5 MW dump.

At higher cms energy (1 TeV), the beams at collision are smaller, and the beamstrahlung and coherent pair production effects more severe. The use of disrupted beams at higher energy also seems possible, but needs more study.

The disrupted beam depolarization due to the beam-beam interaction at the IP was calculated using the program CAINcain . Figure 6 shows that the depolarization in the highest energy bins is small (0.7%), even for 500 GeV beams (1 TeV cms). The depolarization for 250 GeV beams would be 0.3%.

## III Requirements on the Fixed Target Beams and Possible Experiments

### III.1 Requirements on the Beam Parameters

The demands on beams for fixed target physics are very modest compared to the demands of luminosity production with colliding beams. For colliding beams both beams must have a very small emittance, and both the electron and positron system must be working at the highest level.

These requirements do not apply for fixed target experiments. Here the luminosity is determined by the number of beam particles and the density of the target, and not the emittance of the beams. Luminosities of $``$ 1039 cm-2 sec-1 can be achieved; this is 4-5 orders of magnitude larger than in colliding beam experiments. Furthermore, only the system for producing polarized electrons is needed.

Traditionally, in preparing collider experiments a “Snowmass year” = 107 sec, corresponding to 32% efficiency, is used for run time calculations. Because of the relaxed conditions we will use $`\sqrt{2}`$\*32% = 45% efficiency for fixed target experiments at the NLC. This number compares well with the experience at SLAC.

| Parameters | Value Required |
| --- | --- |
| Energy spread $`\mathrm{\Delta }E`$ | $``$ 1% E |
| Polarization | $``$ 80% |
| Beam size | $`<`$$``$ 1 mm |
| Charge jitter | $`<`$ 2% |
| Position jitter on target | $`<`$ 100 $`\mu `$m |
| \*Beam asymmetries $`\mathrm{\Delta }Q`$ | $``$ 10-9 Q |
| \*Position asymmetries $`\mathrm{\Delta }x`$ | $`<`$10 nm |

Requirements of fixed target electron beams: These requirements are similar to present day spin physics experimental needs, properties already achieved in End Station A at SLAC. The two last quantities (\*) are the uncertainties with respect to the polarization state integrated over an entire run.

As shown above, the specific requirements of energy spread and polarization for fixed target experiments can be achieved, even with disrupted beams. This opens up parasitic use of otherwise unused beams, thus enhancing the physics value of operations at moderate cost.

### III.2 Fixed Target Area Design Parameters

The fixed target area, including the targets, must be accessible during colliding beam operations to make installation and operation of the fixed target experiments efficient and independent from the colliding beam operation. The main colliding beam dump at $``$ 200 m is assumed to be hermetic, except for $`\mu `$’s. Muons range out with 0.4 m per GeV in earth and a transverse spread of about 10 m. At 250 GeV they will reach $``$600 m in the earth beyond the main colliding beam dump, so the transverse separation of the experimental area from the primary line-of-flight of the $`\mu `$’s has to be at least 10 m.

This optimal scattering angles at 250 GeV are 1-4 mrad for Møller experiments and $``$ 20 mrad for spin structure experiments. This in turn requires large drift sections to have enough transverse space available for the spectrometers. The target and spectrometer locations can be separated by a divergent tunnel to minimize the excavation costs of an experimental hall.

### III.3 Possible Experiments

Below is a list of possible experiments which should continue to be interesting ten years from now, thereby opening a window of opportunity for NLC. The availability of high energy, high flux and high polarization $`\gamma `$-beams will spawn new and exciting proposals because the figure of merit of $`\gamma `$-absorption experiments increases with energy.

Experiments with (e,e):

* ALR by Møller scattering (sin2$`\mathrm{\Theta }`$W)
* Spin structure functions at very low $`x`$

Experiments with $`\gamma `$-absorption :

* Polarized gluon density $`\mathrm{\Delta }`$G
* Gerasimov-Drell-Hearn sum rule (GDH)
* General charm physics

With a backscattered laser beam, quasi-monochromatic $`\gamma `$’s of 200 GeV energy can also be produced, but at much lower intensities.

As listed above, with $`\gamma `$-absorption, the gluon spin contribution $`\mathrm{\Delta }`$G in the nucleon can be determined. This is an interesting possibility via open charm production from $`\gamma `$g fusion, because the cross section for $`\gamma `$N -$`>`$ cc̄X at 200 GeV is 700 nb, so that the event rate with a primary $`\gamma `$-beam will be very high. An experiment using muons, but similar in kinematics with a much lower event rate, is COMPASS at CERN compass . A new proposal to do such measurements at SLAC with a 45 GeV coherent bremsstrahlung beam (SLAC-Proposal E-161, slace161 ) has just been approved.

Another experiment which could have a very large impact on the physics of charm mixing, is $`\gamma `$N -$`>`$ cc̄X. Here the cross section is 3 orders of magnitude larger than at B-factories. The kinematic boost is large compared even to asymmetric B-factories, opening up a new regime for lifetime determinations, D0-mixing and CP-violation. This experiment would improve on the work of the FOCUS collaboration (FermiLab E831, focus ) which used a 250 GeV endpoint broad band bremsstrahlung beam from secondary electrons and positrons produced from the 800 GeV Tevatron proton beam.

### III.4 Specific Example: The Running of sin2$`\mathrm{\Theta }`$W by Polarized Møller Scattering

SLAC is currently finishing construction of a 50 GeV fixed target experiment to measure the weak mixing angle Kumar95  at low momentum transfer (off the Z) with Møller scattering with great precision (E-158, slace158 ). The first runs are scheduled for early 2001.

Figure 7 shows a theoretical prediction and several experimental measurements of sin2$`\mathrm{\Theta }`$W as well as the predicted errors for two potential NLC fixed target experiments. Like other Weinberg angle measurements SLD , they are based on measuring ALR, the left-right asymmetry.

The table below shows the effect of beam energy on the Møller measurement error.

| Experiment | E-158 | NLC-I | NLC-II |
| --- | --- | --- | --- |
| E/GeV | 46.4 | 250 | 500 |
| ALR/10-7 | 3.2 | 16.1 | 32.2 |
| relative $`\delta `$ from ALR alone | 1 | 1/5.4 | 1/10.8 |

TABLE: Impact of beam energy on the Møller measurement error: as explained in the text the statistical error decreases with increasing beam energy. 100% polarization was assumed for the table.

The small errors result mainly from the fact that the left-right asymmetry is proportional to the beam energy: A$`{}_{LR}{}^{}`$E. This has a unique and beneficial consequence for Møller scattering. Generally in electron scattering the cross section is proportional to the inverse of the beam energy squared, $`\sigma E^2`$, but in the case of Møller scattering the cross section is $`E^1`$ (see Ref. slace158 ). The figure of merit for this experiment, however, is $`A^2\sigma E`$. The consequence is that the statistical error decreases with increasing beam energy. The point in Figure 7 at NLC-250 assumes a total charge of 170 Coul on target. Assuming the nominal NLC beam of 1012 electrons per bunch train, this corresponds to eight months running at 45% efficiency using an undisrupted beam. For the disrupted beam, 57% of the beam is within $`\mathrm{\Delta }`$E/E = 1%. In this case 120 Coul can be collected in nine months running, or 170 Coul in 13 months.

In summary for Møller scattering, the use of disrupted beams is possible and the beam depolarization due to colliding beams is small. Within about one year of running at NLC energies, the compositeness of the electron can be probed at the 60 TeV scale and the existence of an extra Z-boson, Z, can be verified or excluded up to 2.7 GeV. An important cross check on the Higgs mass can be performed at the 10% level kumar96 .

## IV Summary

Single beam use at the NLC offers possibilities for a wide variety of fixed target experiments at (relatively) modest cost. We have investigated the possibilities of attaching a single beam facility to the low energy IR at 250 GeV maximum energy for polarized electrons, 200 GeV for polarized quasi-monochromatic photons. Other locations are possible. NLC detector development and testing, in particular for a detector in the HE IR, is possible with the early availability of a 250 GeV test beam.

Acknowledgments
We would like to thank: Peter Bosted (University of Massachusetts, Amherst), David Burke (SLAC), Don Crabb (University of Virginia), Krishna Kumar (University of Massachusetts, Amherst), Mike Olson (St. Norbert College), Tom Markiewicz (SLAC), Makis Petratos (Kent State University) and Terry Toole (American University), for much and diverse help and contributions.