A rheological model of a fractured solid

Introduction

The early experiments of Adams and his colleagues Von Karman (1911) and Griggs (1936) demonstrated that rheological properties of rocks depend on the applied load and temperature. Two basic regimes of rock behavior have been revealed in these and later experiments (Kirby, 1980):

(1) The brittle regime, in which specimen deformation results from localized displacements or shear fractures or faults.

(2) The ductile regime, which is characterized by a very small effect of pressure on strength, no significant dilatation and negligable evidence of intercrystalline glide or other mechanisms of plastic deformation.

The transition between brittle and ductile behavior involves a broad region of semibrittle behavior in which both stable microfracturing and ductile processes occur.

Empirical relationships between the strength and the critical least principal stress of the brit-tle-ductile transition were investigated by Byerlee (19681, Kirby (1980) and many others. Goetz and Evans (1979) paid attention to the fact that the application of a nnmber of well-known constitutive relationships, such as perfectly elastic, visco-elastic, elasto-plastic to plastic behavior of the lithosphere, fails to describe a number of basic phenomena. Among them the dependence of the deformational parameters of the material on the type of loading, the ability of the material to remember the history of a fracture process, etc. The main difficulty lies in describing the brittle regime. There are different ways of expressing the constitutive relationships for non-homogeneous or fractured media. One of them is based on the assumption that a body consists of a matrix and a set of cracks or inclusions. All physical characteristics of the matrix and the inclusions are supposed to be known. This approach was put forward by Eshelby (1957) and has been successfully applied to composite materials (Christensen, 1979). For higher crack densities when crack interactions cannot be neglected, O'Connel and Budiansky (1974) and Budiansky and O'Connel (1976) proposed a self-consistent method for the calculation the elastic moduli for random crack orientation statistics. In this method the effect of crack interactions is included by assuming that each crack is embedded in a medium with the effective stiffness of the cracked body. Bruner (1976) and Henyey and Pomphrey (1982) have pointed out that this scheme may overestimate the crack interactions and have proposed an alternative, differential scheme. Another method has been proposed by Kachanov (19871, Sayers and Kachanov (1991) for finding the effective properties for solids having interacting cracks with arbitrary crack interactions. However, these results are obtained for each given arrangement of cracks rather than in statistical terms. Thus, the only schemes available at present for ~al~lating the effective elastic constants at finite crack densities are the self-consistent and differential schemes. Both of them are not applicable for the description of a fracture distribution and fragmentation of rocks which have been identified as power law phenomena and indicates a fractal relation (e.g., Turcotte, 1986;Herrmann and Roux, 1990;Velde et al., 1991). Fractal analysis has been used to characterize the nature of the pattern of fractures in natural rocks and particularly the pattern in major fault systems (e.g., Okubo and Aki, 1987;Aviles et al., 1987). In this case, there is another method to simulate effective rheological parameters based on a phenomenological approach.

Experimental studies of the behavior of various granular materials show noticeable dependence of elastic properties on the type of loading. This is particularly manifest in the difference between strain characteristics of the material under tension and under compression. Opening and closing of cracks under tension and under compression produce variations of elastic moduli and including abrupt changes when the strain reverses. The dependence of elastic properties on the type of loading has been observed in experiments with a variety of design materials (Collins, 1981) and rocks (e.g., Zoback and Byerlee, 1975;Stavrogin and Protosenya, 1979;Alm et al., 1985).

Several mechanical models with strain-depen-dent moduli have been suggested; the one of Ambartsumyan-Khachatryan (Ambartsumyan, 1982) assumes that in the stress-strain law each of the compliances (Poisson's ratio divided by Young's modulus) changes whenever the stress for which this compliance acts as a factor reverses. Jones (1977) suggested a model of a material with a matrix of weighted compliances; the weights depend on absolute values of the principal stresses. For materials with a weakly non-linear response it is possible to assume that the elastic moduli depend only on the type of stress. This assumption was the basis for the model of Lomakin and Rabotnov (1978).

A long-term action of stress and temperature will cause growth or healing of microcracks resulting in a variation of the rheological properties of a solid. The recent investigations by Yukutake (1989) and of a fracturing process in granite show that this process may hardly be described in terms of the classical Griffits model of single crack propagation. Lockner and Madden (1991) put forward a multiple-crack model of brittle fracture. Many attempts have been made to define and estimate suitable variables to describe the extent of fracturing. Such variables should be able to characterize the processes of fracturing quantitatively and should be suitable for modeling from the point of view of continuum mechanics. Among such approaches are Robinson's (1952) linear cumulative creep damage law, Hoff's (1953) ductile creep rupture theory, Kachanov's (1958Kachanov's ( , 1986 brittle rupture theory, Rabotnov's (1969) coupled damage creep theory and many modifications of these theories. Following this approach and using the strain-dependent moduli model, a description of the damage evolution of a fractured medium under longterm stress has been put forward (Lyakhovsky and Myasnikov, 1985;Myasnikov et al., 1990;Lyakhovsky et al., 1991). This model has made it possible to simulate the evolution of fracturing under tectonic stress, resulting in the appearance of seismic boundaries of rheological nature (Mints et al., 1987;Myasnikov, 1987, 1988) and investigate a faulting process in the northern Dead Sea rift (Ben-Avraham and .

MODEL OF A FRACTURED SOLID

tions. An application of this method to the simulation of a visco-elastic flow is described in Bolshoi et al. (1992). The computational mesh consists of quadrilateral elements, which are subdivided into pairs of constant-strain triangles, with different diagonals. This overlay scheme ensures symmetry of solution by averaging the results obtained on two meshes (Cundall and Board, 1988). Linear, triangularelement shape functions L, (k = 1,3) are defined as:

where ak, b,, k are instants and (x,y) are grid coordinates. These shape functions are used to interpolate the nodal velocities ZJ!~' within each triangular element (e): k=t This formula enables the calculation of the strain increments in each triangle and the element stresses using eqn. 6). When the stresses in each triangle are known, the forces at the node number n @!"I) are calculated by projecting the stresses from all elements surrounding the node:

where nj is the jth component of the unit vector normal to each side of the two element adjacent to the node n. The length of each side is denoted by Al. The minus sign is a consequence of Newton's Third Law. New velocities are computed by integrating eqn. (8) over a given time step At:

If a body is at mechanical equilibrium, the forces on each node are zero; otherwise, nodes are accelerated. New c~rdinates of grid nodes are computed by: The FLAC method has advantage over implicit methods because it is computationally inexpensive for each time step and it is memory-efficient because matrices storing the system of equations are not required.

A computer simulation of the shear deformation of an initially homogeneous rectangular body and the damage distribution is shown in Figures 8 and 9. The boundary conditions of this example were as follow: shear movement at the top and the bottom and free surface at the left and right. The first step of deformation (Fig. 8b) results in a C Fig. 9. Same as in Fig. 8 showing the damage increase (isolines correspond to a level of the damage).

homogeneous deformation similar to a shear flow in a linear viscous body. The deformation causes a change of the shape of the area and creates conditions for stress concentration at the two corners (Fig. 9b). The damage propagates from the corners into the area (Fig. 9c, d) and finally yields two parallel zones separating the sample into three rigid blocks (Fig. 9e). This result of three rigid blocks separated by fractured zones can be also seen in Figure 8c.

Deformational properties of rocks

Following Lyakhovsky and Myasnikov (1984), the dependence of the elastic potential Ue on the invariants I, = eii and I, = l ijeij of the strain tensor l ij may be written as:

where i, j = 1, 2, 3; p is the density, A, /.L are Lame constants; the coefficient v is an additional elastic modulus; here and later, repeating indexes imply summation.

Apart from the quadratic terms containing Z: and Z2 common to elastic solids, the potential (1) includes the non-analytical second-order term vZ,& which accounts for the effects of fracturing in the material. Such an addition to the elastic potential is the simplest one for which the requirement of being analytical is rejected. It allows for differing properties of a solid under tension and under compression. For models in which non-linearity is taken into account by means of higher-order terms, non-linear effects vanish at small strains. The suggested model preserves its properties under arbitrarily small deformations since the principal terms and the additional term are of the same order.

In a solid described by potential (11, the relation between the stress tensor aij and the strain tensor is:

where Sij = 1, if i = j, and 0, if i Zj.

Using the function 6 = I,/ fi to describe the form of stress state, effective elastic moduli can be introduced by:

To find the value ,$ in terms of equation should be solved:

where:

vii, the following

(3) Constructing a numerical method to find the strain tensor in terms of a known stress tensor to solve eqn. (31, we have to search for roots of a polynomial with real coefficients. In the case of pure shear stress (aii = 01, eqn. (3) has the form:

and has the following solution:

The sign in front of the radical is chosen to give regularity of the solution as v + 0. For small v we get the solution by neglecting higher terms:

The obtained expression for 5 is strictly positive. This means that a shear stress in the solid causes extensional strain (dilatation). The effect of dilatancy in rocks was discovered by Bridgeman in 1949 and has been studied by Reynolds (1885) and has been studied by many researchers. Figure 1 shows the average pressure p = (2~7, + a,)/3 of Westerly granite under 3-D load versus volume strain of sample @chock, 1977); a, is the maximal tension stress; ui = a2 are compressive stresses. The results of analogous simulations using the above described scheme are also shown on Figure 1. The elastic moduli were assumed to be h = CL, Y = 0.2/t. Comparison of the slopes of the average pressure shows a good agreement between the simulated description of the dilatancy and the experiment. Figure 2 shows the results of sample testing with marble and diabase under 1-D and 3-D compression (Stavrogin and Protosenia, 1979). The connection between stress and strain for small deformation is linear, but the value of elastic moduli depends on the type of loading. According to this experiment bulk modulus of marble under 3-D compression is 3.3 X lo5 kg/cm', but under 1-D compression it is 2.5 X lo5 kg/ '.@ (a) 1 b

cmL. This means a necessity to take into account the dependence of elastic moduli on the type of loading and makes it possible to estimate the value of additional modulus incorporated in eqn.

(1). Using the results of experiment shown in Figure 2, the estimation of v/(A + 2~) for marble is about 10% and for diabase is about 2.5%. This difference is very natural because the strength of diabase is much higher than that of marble. Figure 3 shows results of estimations of effective elastic moduli A, p for different types of loading 5 and the approximation of a strain-dependent moduli model. Experimental points were calculated using results of granite sample testing collected by Volarovich (1988). Each point represents an independent measurement of the load and the resulting strain tensor that allows an easy calculation of the value of elastic moduli. The approximation represented by Figure 3 is much better then the one by Hook which provide constant (independent on type of loading) estimation of elastic moduli. Evaluation of the three elastic moduli as a best fit to the experimental data gives: A = 3.95 X lo4 MPa; p = 0.11 X lo4 MPa; v = 0.12 x lo4 MPa.

Model of medium with damage

Many writers on continuum thermodynamics have postulated it as a state function of all the state variables, including "hidden variables" (e.g., Malvern, 1969) not available for macroscopic observation. Coleman (1964) and Coleman and Mizel (1964) considered the class of materials in which thermodynamical state functions depend not only on the instantaneous values of the variables, but also on their gradients. Coleman's approach promises to be useful since it deals with a limited numbers of explicitly enumerated macroscopic state variables and not a vague collection of unspecified substate variables. Edelen and Laws (1971) and Edelen et al. (1971) demonstrated the dependence of thermodynamical state functions on gradient of "hidden variables" results in nonlocal properties of continuum. This approach was used by Truskinovsky (1991) to describe the process of phase transition.

In order to simulate a fracturing process in terms of fractal analysis an non-dimensional damage parameter (Y was incorporated. This parameter lies in the interval from 0 to 1 to describe the evolution of the medium from the ideal undestroyed elastic solid to the absolutely destroyed material.

Following this approach we represent the free energy of a solid F as:

where F. = F,, CT, (~1; U" is the elastic potential (1) with the elastic moduli A, p, u to be functions of a; 0 is the coefficient of temperature expansivity, and T the temperature. We express the elastic strain tensor eij in terms of a metric tensor describing the current state of medium gij and a metric tensor gly. describing the irreversible viscous deformation:

The expression for the whole strain rate tensor eij = i(dgij/d t) is similar to the common used form, usually formulated in terms of velocities (vi):

Since the free energy is a function of T, l ij, (Y, taking into account the irreversible changes of T, (Y and g& the balance equations of energy and entropy can be written (e.g., Malvern, 1969). In the case of lack of damage variation and viscous flow, the process of deformation must be reversible.

Using this condition we arrive at the ordinary equation for the stress tensor (Myasnikov et al., 1990):

It is usually assumed in irreversible thermodynamics (e.g., Fitts, 1962;De Groot and Mazur, 1962) according to a principle of maximum rate of entropy production or maximum dissipation power, that the constitutive equations give the fluxes as functions of the forces, and that at least "in the neighborhood of equilibrium" the constitutive equations (called phenomenological equations) are equations of the form:

where C is a positive constant describing the temporal scale of a damage process.

This equation is well-known as the Ginzburg-Landau equation, widely used in the theory of phase transitions (e.g., Honenberg and Halperin, 1977;Umantsev and Roitburd, 1988;Bronsard and Kohn, 1991;Truskinovsky, 1991). Analogous to the J-integral in fracture mechanics, the thermodynamic force aF/aa can be interpreted (Lemaitre, 1985) as the energy release rate resulting from the increase of damage.

To describe a viscous flow we use a common rheological law for a incompressible Maxwell body:

a%ck

To close the set of eqns. (5-7), the equation of motion have to be added <fi forces):

The final system of eqns. (5-8) describes elastic deformation, viscous flow, and damage evolution of the media. Detals of the application of the thermodynamics equations to the damage process were discussed in Myasnikov et al. (1990) and Lyakhovsky et al. (1991).

Parameters and main features of the model

The simplest formulation of the model is to assume a linear dependence of elastic moduli A, w, v on the damage parameter a. For a = 0 the material is expected to be an ideal undestroyed solid, i.e., a linear elastic Hook's model. So, for (Y = 0 the third elastic modulus v in eqn. (1) have to be equal to zero. For a = 1 the material is absolutely destroyed, i.e., it is impossible to create stresses resulting in expansion. Mathematically this means that elastic potential U' (eqn. 1) is not convex as a function of strain tensor l ij for I, 2 0. This condition gives us the maximum value for Y = v1 as a function of A =hl, p =pI for a= 1:

Using this assumption the dependence of elastic moduli on a may be written as:

The next step is to describe the existence of two different stable regimes which may be interpreted as ductile and brittle, and the process of rheological transition (bifurcation). From the theory of catastrophe it is known that to describe the bifurcation process we have to assume the energy to be a function of the fourth order of its parameter (e.g., Hale and Kocak, 1991). Following this recommendation, we represent F, as:

where K is a constant and (Ye, a2 are levels of damage corresponding to the ductile and brittle regimes, respectively. Using eqns. (6) and (9) the damage evolution is: where:

Equation 10has a number of stationary solutions. The dependence of the free energy F on a in the case of an unloaded sample (aij = 0) is shown in Figure 4a. There are two equivalent stable solutions corresponding to the ductile (a = a,) and brittle (a = a,> regimes and one unstable [a = (a, + a*)/21 solution. Thus, the damage of unloaded homogeneous solid comes to one of these two solutions. This process may be very slow and the solid can "remember" its initial ll-

. Differentiat stress versus confining pressure at the transition to brittle behavior; I= granite Kola Peninsula Wavrogin and Protosenia, 1985); 2 = Spruce Pine Dunite; 3 = Cabramurra Serpentinite; 4 = Nahant Gabbro (Raleigh and Paterson, 1965). Solid line is the boundary between the brittle and ductile regions determined from the model. level of damage for a very long time. Loading of a sample results in a variation of its elastic energy causing an increase or decrease of damage. A relatively high shear stress promotes the increase of damage. The energy (Fig. 4b) corresponding to the ductile regime becomes much higher than the energy of the brittle one. The subsequent increase of the stress leads to a disappearance of the first solution ( Fig. 4c>. It is impossible for the sample to keep its initial ductile regime and it yields an abrupt change of the damage. At this moment a rheological transition from the ductile to the brittle regime occurs. Under relatively high pressures this process may reverse and will result in a healing of damage.

The level of stresses at the ductile-brittle transition depends on the type of loading, Figure 5 represents the comparison of the simulated differential stress versus confining pressure at the ductile-brittle transition with results of experiments for different crystalline rocks.

The problem of space-temporal damage evolution is of special interest for various geological problems. One of the main features of the model, as it will be shown nume~cally, is the possibili~ to describe the process of creation of narrow high fractured zones with the strain rate localization. Since the modeled solid has "memory" it is possible to introduce some pre-existing high-and lowdamage zones and to simulate their temporal evolution. After the unloading, in contrast to the models of plasticity, the simulated material "remembers" information about the created damage zones. The next stage of loading may reactivate these old zones or provoke some new fracturing.

Numerical simulation of damage evolution under single axis compression

First, we start the simulation with the most simple example: single axis compression of the rectangle sample neglecting the process of viscous relaxation of stresses. The process of damage evolution is assumed to be quasi-static, i.e., the stress distribution satisfies the equilibrium equation of an elastic solid for each time step. The new pattern of damage is calculated using eqn. (10) for the known stress distribution and the new distribution of cr or new elastic moduli results in a new stress field at the next step. To simulate elastic stresses at each time step of the simulation, the unite-element method based on an algorithm by Zienkiewicz (1971) was used. Application of the finite-element method to the solution of different variations of the elastic problem was widely discussed in literature.

In this case, quadratic interpolation of displacements within two-dimensional triangular elements was used. The area integrals were evaluated by Gauss-Legendre quadrature. To take into account the dependence of effective elastic moduli on the type of loading simple iteration method was used. Figure 6 shows different time steps of the damage evolution in a rectangle area under one axis compression. To initiate the beginning of the fracture, a small zone of high damage value was located at the boundary. Stress distribution and evolution of damage at the tip of the fracture zone is similar to the problem of crack growth in elastic media. In contrast to the problem of crack evolution in elastic media, due to the finite size of a gradient zone between destroyed and undestroyed area, there is no sin~lari~. Stress distri-possible. The initially continuous area breaks into two pieces. bution has a regular solution at every point and eqn. (10) provides the fracture zone with a finite rate of growth. At the intermediate stage (Fig. 6bf, similar to the process of crack propagation in elastic media, there are two different possible directions for the future damage increase. The selection of the direction is very sensitive to the stress distribution and the accuracy of the numerical scheme strongly influences this process. During the final stage of the simulation (Fig. 6~1, the damage zone reaches the opposite boundary of the area and the ~ntinuat~on of the simulation in the frame of elastic displacements becomes im-To simulate any flow it is necessary to define an effective viscosity as a function of the parameters of a solid. Power law viscosity is usually assumed for rocks in the ductile regime. For the brittle regime hardly anything is known about the viscosity, except its very strong dependence on the level of fracturing. To take into account stronger viscosity variations in comparison with elastic moduli, we assume the viscosity 9 to be an exponential function of the damage parameter LU:

Shear flow simulation

where A, B are constants. Figure 7 shows a comparison of a stress-strain rate relation obtained by this model and by power law for the order n = 3 and n = 10.

Numerical simulation has been done using "FLAC" algorithm (Cundall and Board, 1988;CundaH, 1989). The formulation is explicit-intime, using an updated Lagrangian scheme to provide the capability for large strains. FLAC is believed to offer advantages over conventional finite-element schemes in cases where flow and material instabili~ occur. Physical instability can be modeled without numerical instability if internal terms are included in the equilibrium equastrain rate Fig. 7. Stress versus strain rate for ductile regime. Heavy fine suggested model; dotted and dash lines -power law viscosity for n = 3 and n = 10, respectively.

Conclusions

Experiments for studying the behavior of various materials point to the relation that exists between elastic properties and the type of stress. The influence of the state of stress on the elastic-ity in a fractured material can be taken into account in a physically non-linear model of an elastic solid. The additional non~analytica~ second-order term included in the elastic potential makes it possible to simulate the dependence of effective elastic moduli on the type of loading, particularly under tension and compression. Such an addition to the elastic potential is the simplest one for which the requirement of being analytical is rejected. For models in which non-linearity is taken into account by means of higher-order terms, non-linear effects become negligable at small strain. The suggested model preserves its properties under arbitrarily small deformations since the principal terms and the additional term are of the same order. It is also shown that a shear stress in the simulated solid causes extensional strain (dilatation). A comparison of the slopes of the average pressure for the sample testing ( Fig. la) and numerical experiment (Fig. lb) proves a good agreement between the model and experiment.

A long-term action of stress and temperature will cause growth or healing of microcracks resulting in a change in rheological properties of a solid and rheological transition between brittle and ductile regimes. The present study of the rheological model of rocks is based on the ther-mod~amical approach in continuum mechanics. The crucial point of the model is the description of the existence of two different steady regimes which may be interpreted as ductile and brittle, and the process of rheological transition (bifurcation). The simulated state of stress at the ductile-brittle transition (Fig. 5) is close to the empirical Coulomb-Mohr failure criterion and is in a good agreement with laboratory observations. The description of slow flow of material in a ductile regime (Fig. 7) is analogous to the power law viscous model. The results of some numerical experiments are presented in order to study mechanisms involved in strain localization and fracturing processes. These numerical examples show one of the main features of the model, i.e., the possibility to simulate the process of fracture zone propagation. In contrast to the model of single crack propagation in elastic media, there is no singularity at the tip of a propagating fracture zone and a finite rate of its growth. The final stage of the evolution of the fracture zones is similar to the same numerical experiments for visco-plastic models, but the present model makes it possible to investigate all intermediate stages of the process too. This approach makes it also possible to take into account the effect of the rock's memory, or the ability of the material to remember a history of fracture process. Thus, some pre-existing fracture zone might be introduced into the simulation. This feature of the model has been used for the investigation of a faulting process in the northern Dead Sea rift (Ben-Avraham and .