A rheological model of a fractured solid

187 Tectonophys& 226 (1993) 187-198 Elsevier Science Publishers B.V., Amsterdam A rheological model of a fractured solid Vladimir Lyakhovsky a, Yury Podladchikov b and Alexei Poliakov ’ ’ Department of Geo~ysi~ and P~aneta~ Sciences, Raymo~ and Beverly Sackier PacUrryof zyxwvutsrqponmlkjihgfedcbaZYXWVUT Exact Sciences, Tel Aviv Vniue~~~, 69978 Tel Aviv, Israel b institute zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA ofZ.%perimentalMineralogy, Chernogolovka, Moscow District 142432, Russia ’ Hans Ramberg Tectonic Laboratory, Institute of Geology, Uppsala University, Box 555, 751 22 Uppsala, Sweden (Received April 23, 1993; revised version accepted May 26, 1993) ABSTRACT Experiments to study 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 elasticity of a fractured material will be discussed for a physically non-linear model of an elastic solid. The strain-dependent moduli model of material, presented in this paper, makes it possible to describe this feature of a solid. ft also permits to simulate a dilatancy of rocks. A damage parameter, introduced into the model using a the~o~amical approach, allows to describe a rheologicai transition from the ductile regime to the brittle one, and to simulate the rock’s memory, narrow fracture zone creation and strain rate localization. Additionally, the model enables the investigation of the final geometry of fracture zones, and also to simulate their creation process, taking into account pre-existing fracture zones. The process of narrow fracture zone creation and strain rate localization was simulated numerically for single axis compression and shear flow. 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~~-1951/93/$~.~ 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, 0 1993 - Elsevier Science Publishers B.V. All rights reserved 18X zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 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- V. LYAKHOVSKY b-1’ 41 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 Lockner et al. (1991) 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 (1958, 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; Lyakhovsky and Myasnikov, 1987, 1988) and investigate a faulting process in the northern Dead Sea rift (Ben-Avraham and Lyakhovsky, 1992). RHEOLOGICAL MODEL OF A FRACTURED 189 SOLID Deformational properties of rocks Following Lyakhovsky and Myasnikov (1984), the dependence of the elastic potential Ue on the invariants I, = eii and I, = lijeij of the strain tensor lij may be written as: U”= Jji~z:+pz2-vzI~ 1 (1) 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: qj = ( Azr - “~)6ij + (2Z-6- VZ*/g)Eij (2) Fig. 1. Average pressure versus volume strain of Westerly granite under different type ~+/a, of 3-D load (dash lines with markers). Simulated average pressure versus volume strain (heavy lines). 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: ~[~-(3A+2~)~+3v=O and has the following solution: 5=3A;v2p _/m 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: where Sij = 1, if zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA i = j, and 0, if i Zj. t= 3v Using the function 6 = I,/ fi to describe the 3A+2~ form of stress state, effective elastic moduli can The obtained expression for 5 is strictly positive. be introduced by: This means that a shear stress in the solid causes lV=A-V/t; ~e=~---_v~ extensional strain (dilatation). The effect of dilatancy in rocks was discovered by Bridgeman in To find the value ,$ in terms of vii, the following 1949 and has been studied by Reynolds (1885) equation should be solved: and has been studied by many researchers. Figure vii (3A + 2p)e - 3v - v[’ 1 shows the average pressure p = (2~7, + a,)/3 of -= (3) Westerly granite under 3-D load versus volume \is5 (2P - US){_ strain of sample @chock, 1977); a, is the maximal tension stress; ui = a2 are compressive where: stresses. The results of analogous simulations uss, = uijuij - +ui; ing the above described scheme are also shown 190 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA *wo, Fig. 2. Stress-strain relation for marble (0) and diabase (A ). 3-D compression - heavy lines; 1-D compression, (T= ~~,a, = cA = 0, E, - compression (dash lines), Ed - tension (dotted lines). 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/ \ t.YAKtlOLSh> I.1 ,\i 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 ‘.@(a) 1 b Fig. 3. Elastic moduli A (a) and p (b) vs. type of deformation of granite from Kola Peninsula. RHEOLOGICAL MODEL OF A FRACTURED SOLID 191 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFED entropy can be written (e.g., Malvern, 1969). In state variables, including “hidden variables” (e.g., the case of lack of damage variation and viscous Malvern, 1969) not available for macroscopic obflow, the process of deformation must be reservation. Coleman (1964) and Coleman and versible. Mizel (1964) considered the class of materials in Using this condition we arrive at the ordinary which thermodynamical state functions depend equation for the stress tensor (Myasnikov et al., not only on the instantaneous values of the vari1990): zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPON ables, but also on their gradients. Coleman’s approach promises to be useful since it deals with a aF limited numbers of explicitly enumerated macro(5) Oij= PG scopic state variables and not a vague collection It is usually assumed in irreversible thermodyof unspecified substate variables. Edelen and Laws (1971) and Edelen et al. (1971) demonnamics (e.g., Fitts, 1962; De Groot and Mazur, strated the dependence of thermodynamical state 1962) according to a principle of maximum rate functions on gradient of “hidden variables” reof entropy production or maximum dissipation sults in nonlocal properties of continuum. This power, that the constitutive equations give the approach was used by Truskinovsky (1991) to fluxes as functions of the forces, and that at least describe the process of phase transition. “in the neighborhood of equilibrium” the constiIn order to simulate a fracturing process in tutive equations (called phenomenological equaterms of fractal analysis an non-dimensional damtions) are equations of the form: age parameter (Ywas incorporated. This parameter lies in the interval from 0 to 1 to describe the evolution of the medium from the ideal undewhere C is a positive constant describing the stroyed elastic solid to the absolutely destroyed temporal scale of a damage process. material. This equation is well-known as the GinzburgFollowing this approach we represent the free Landau equation, widely used in the theory of energy of a solid zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA F as: phase transitions (e.g., Honenberg and Halperin, F=F,+U”+/ 3Z,(T-T,,) (4) 1977; Umantsev and Roitburd, 1988; Bronsard and Kohn, 1991; Truskinovsky, 1991). Analogous where F. = F,, CT, (~1; U” is the elastic potential to the J-integral in fracture mechanics, the ther(1) with the elastic moduli A, p, u to be functions of a; 0 is the coefficient of temperature expanmodynamic force aF/ aa can be interpreted sivity, and T the temperature. (Lemaitre, 1985) as the energy release rate resultWe express the elastic strain tensor eij in ing from the increase of damage. terms of a metric tensor describing the current To describe a viscous flow we use a common state of medium gij and a metric tensor gly. rheological law for a incompressible Maxwell describing the irreversible viscous deformation: body: 1 aF - -aij zyxwvutsrqponmlkjihgfedcbaZYX -=- 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): 1 avi aVj I eii= 5 -+G. i axj I Since the free energy is a function of T, lij, (Y, taking into account the irreversible changes of T, (Y and g& the balance equations of energy and 3 a%ck To close the set of eqns. (5-7), the equation of motion have to be added <fi forces): dvi aaij z +.fi”Pd, I 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 192 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 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 lij 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: - a 1 a 2 a F of elastic + A,a 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: F, = K( U I_ I’..\1 k-a-l Parameters and main features of the model A = A,( 1 -a) LYAKHOVSKY FI were discussed in Myasnikov et al. (1990) and Lyakhovsky et al. (1991). Using this assumption the dependence moduli on a may be written as: V Ul)2(U - a~)’ Fig. 4. Free energy as a function of damage parameter a for different stresses. where: Equation (10) has a number of stationary soluwhere zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA K is a constant and (Ye, a2 are levels of tions. The dependence of the free energy F on a damage corresponding to the ductile and brittle in the case of an unloaded sample (aij = 0) is regimes, respectively. Using eqns. (6) and (9) the shown in Figure 4a. There are two equivalent damage evolution is: stable solutions corresponding to the ductile (a = a,) and brittle (a = a,> regimes and one unstable (u,e - u;) + zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 4 K( (Y - a,)( (Y - (Y2) [a = zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPO (a, + a*)/21 solution. Thus, the damage of unloaded homogeneous solid comes to one of a1 +a2 these two solutions. This process may be very x a-(10) 2 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA slow and the solid can “remember” its initial ( )I RHEOLOGICAL MODEL OF A FRACNRED 193 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA SOLI D to describe the process of creation of narrow high fractured zones with the strain rate localization. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA Since the modeled solid has “memory” it is possiIOble 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 0 -1 these old zones or provoke some new fracturing. ll- x -2 p-I A -3 Fig. 5. 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~ 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 GaussLegendre 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 two pieces. into Shear flow simulation 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: q =A exp( -Ba) Fig. 6. Distribution of damage parameter LYfor the simulated evolution of the narrow fracture zone in rectangular area under I-D compression. (a) Geometry of the pre-existing damage zone to initiate the fracture process. (b) intermediate stage of the damage zone propagation. Cc) Final stage of simulation. 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- 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 equa- strain 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. RHEOLOGICAL MODEL OF A FRACTURED 195 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 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA L, (k = 1,3) are defined as: Lk=ak+xbk+yck, k=1,3 and (x,y) are grid where ak, b,, k are instants 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: $)( t + At) = u!“‘(t) I + c Fig. 8. Computer simulation of the shear deformation of initially homogeneous rectangle area. (a) Initial shape of the area. (b) The beginning of the deformational process. Cc) The final stage demonstrates the three rigid blocks which moves separately. ,@)E P The FLAC method has advantage over implicit methods because it is computationally inexpenIf a body is at mechanical equilibrium, the net sive for each time step and it is memory-efficient forces on each node are zero; otherwise, the because matrices storing the system of equations nodes are accelerated. New c~rdinates of the are not required. grid nodes are computed by: A computer simulation of the shear deformax$“)( t + At) =x;“)(t) + u,‘“‘ht tion of an initially homogeneous rectangular body and then calculations are repeated for the new and the damage distribution is shown in Figures 8 configuration. Together with the deformation, a and 9. The boundary conditions of this example new damage distribution is calculated using were as follow: shear movement at the top and known stresses. Thus, visco-elastic deformation of the bottom and free surface at the left and right. a damage material is simulated. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA The first step of deformation (Fig. 8b) results in 196 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 a zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 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 C and ductile regimes. The present study of the Fig. 9. Same as in Fig. 8 showing the damage increase (isorheological model of rocks is based on the therlines correspond to a level of the damage). 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, homogeneous deformation similar to a shear flow and the process of rheological transition (bifurcain a linear viscous body. The deformation causes tion). The simulated state of stress at the a change of the shape of the area and creates ductile-brittle transition (Fig. 5) is close to the conditions for stress concentration at the two empirical Coulomb-Mohr failure criterion and is corners (Fig. 9b). The damage propagates from in a good agreement with laboratory observathe corners into the area (Fig. 9c, d) and finally tions. The description of slow flow of material in yields two parallel zones separating the sample a ductile regime (Fig. 7) is analogous to the into three rigid blocks (Fig. 9e). This result of power law viscous model. three rigid blocks separated by fractured zones The results of some numerical experiments are can be also seen in Figure 8c. presented in order to study mechanisms involved in strain localization and fracturing processes. Conclusions These numerical examples show one of the main features of the model, i.e., the possibility to simuExperiments for studying the behavior of varilate the process of fracture zone propagation. In ous materials point to the relation that exists contrast to the model of single crack propagation between elastic properties and the type of stress. in elastic media, there is no singularity at the tip The influence of the state of stress on the elastic- RHEOLOGKAL MODEL OF A FRACTURED SOLID 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 Lyakhovsky, 1992). References Aim, O., Jaktlund, L.L. and Kou, S., 1985. The influence of microcrack density on the elastic and fracture mechauical properties of Stripa granite. Phys. Earth Planet. Inter., 40: 161-179. Ambartsumyan, S.A., 1982. Raznomodulnaya Teoria Uprugosti (Strain-Dependence Elastic Theory). Nauka, Moscow, 320 pp. (in Russian). Aviles, CA., Scholz, C.H. and Boatwright, J., 1987. Fractal analysis applied to characteristic segments of the SanAndreas fault. J. Geophys. Res., 92: 331-344. Ben-Avraham, Z. and Lyakhovsky, V., 1992. Faulting process along the northern Dead Sea transform and the Levant margin. Geology, 20: 1143-1146. Bolshoi, A.N., Cundall, P.A., Podladchikov, Y.Y. and Lyakhovsky, V.A., 1992. An explicit inertial method for the simulation of viscoelastic flow: an evaluation of elastic effects on diapiric flow in two and three layers models. UMSI Rep. 92/110, 20 pp. Bronsard, L. and Kohn, R.V., 1991. Motion by mean curvature as a singular limit of Ginzburg-~ndau Dynamics. 3. Diff. Eqn., 90: 211-237. Bruner, W.M., 1976. Comments on “Seismic velocities in dry and saturated cracked solids” by O’Connel and Budiansky. J. Geophys. Res., 81: 2573-2576. Budiansky, B. and Q’Connel, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct., 12: 81-97. Byerlee, J.D., 1968. Brittle-ductile transition in rocks. J. Geophys. Res., 73: 4741-4750. Christensen, R.M., 1979. Mechanics of Composite Materials. Wiley, New York, NY, 330 pp. Coleman, B.D., 1964. Thermodynamics of materials with memory. Arch. Ration. Mech. Anal., 17: l-46. Coleman, B.D. and Mizef, V.J., 1964. Existence of caloric equations of state in thermodynamics. J. Chem. Phys., 40: 1116-1125. 197 Collins, J.A., 1981. Failure of Materials in Mechanical Design. Wiley, New York, NY, 620 pp. Cundall, P.A., 1989. Numerical experiments on localization in frictional materials. tng. Arch., 59: 148-159. Cundall, PA. and Board M., 1988. A microcomputer program for modeling large-strain plasticity problems. In: C. Swoboda (Editor), Numerical Methods in Geomechanics. Proc. 6th Int. Conf. on Numerical Methods in Geomechanics, Innsbruck. Balkema, Rotterdam, pp. 2101-2108. De Groot, S.R. and Mazur, P., 1962. Nonequilibrium Thermodynamics. North”Holland Publishing Co., Amsterdam. Edelen, D.G.B. and Laws, N., 1971. On the thermodynamics of systems with nonlocality. Arch. Ration. Mech. Anal., 43: 24-35. Edelen, D.G.B., Green, A.E. and Laws, N., 1971. Nonlocal continuum mechanics. Arch. Ration. Mech. Anal., 43: 36-44. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. R. Sot. London, Ser. A 241: 376-396. Fitts, D.D., 1962. Nonequilibrium Thermodynamics. McGrawHill, New York, NY. Griggs, D.T., 1936. Defo~ation of rocks under high confining pressure, J. Geol., 44: 541-577. Goetz, C. and Evans, B., 1979. Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics. Geophys. J.R. Astron. Sot., 59: 436-478. Hale, J. and Kocak, H., 1991. Dynamics and Bifurcations. Springer-Verlag, New York, NY, 568 pp. Henyey, F.S. and Pomphrey, N., 1982. Self-consistent moduli of a cracked solid. Geophys. Res. Lett., 9: 903-906. Herrmann, H.J. and Roux, S. (Editors), 1990. Statistical Models for the Fracture of Disordered Media. Elsevier, Amsterdam, 353 pp. Hoff, N.J., 1953. The necking and rupture of rods subjected to constant tensile loads. J. Appl. Mech., 20: 105-108. Hohenberg, P.C. and Halperin, B.I., 1977. Theory of dynamic critical phenomena. Rev. Mod. Phys., 49: 435-480. Jones, R.M., 1977. Stress-strain relation for materials with different moduli in tension and compression. AIAA J., 15: 62-73. Kaehanov, L.M., 1958. Time of the rupture process under creep conditions. Izv. Acad. Nauk USSR, Otd. Tech., 8: 26-31 (in Russian). Kachanov. L.M., 1986. Introduction to Continuum Damage Mechanics. Martinus Nijhoff, Dordrecht, 135 pp. Kachanov, L.M., 1987. Elastic solids with many cracks: a simple method of analysis. Int. J. Solids Struct., 23: 23-43. Kirby, S.H., 1980. Tectonic stress in the lithosphere: Constrains provided by the experimental deformation of rocks. J. Geophys. Res., 85: 6353-6363. Lemaitre, J., 1985. A continuous damage mechanics model for ductile fracture. J. Eng. Mater. Technol., 107: 83-89. Lockner, D.A., Byerlee, J.D., Kuksenko, V., Ponomarev, A. and Sidorin, A., 1991. Quasi-static fault growth and shear fracture energy in granite. Nature, 350: 39-42. 19x Lockner, D.A. model and Madden, of brittle T.R., fracture. 1991. A multiple-crack 1. Non-time-dependent tion; 2. Time-dependent simulation. Reynolds. simula- J. Geophys. Res.. 96: 19,623-19,654. Lomakin, E.V. and Rabotnov, Yu.N.. dla isotropnogo 1978. Sootnosheniya raznomodulnogo Nauk SSSR. Mekh. Tverd. tela. lzv. Tela. 6: 29-34 V.A. and Myasnikov, of elastic cracked (in Rus- Lyakhovsky, V.P., 1984. On the behavior solid. Phys. Solid Earth, V.A. and Myasnikov, of visco-elastic Lyakhovsky, cracked V.A. and V.P., 1985. On the behavior Myasnikov, V.P., 4: 28-35. 1987. Relation V.A. rheologically Inter., and solids. V.P., beJ. 1988. Acoustics Int. J. Phys. Earth of Planet, 50: 60-64. Lyakhovsky, Rheological 91/296, Malvern. model 30 of Y. and a fractured Poliakov, L.E., 1969. Introduction Medium. solid. to the Prentice-Hall, Kolpakov, Lyakhovsky, problem V.A. N.I., and of the nature (interpretation Lanev, V.S., Myasnikov, V.P., of subhorizontal of Kola super-deep V.P., Lyakhovky, 1990. Non-local media. Rep. Acad. O’Connel, of a Cliffs, NJ, MS., On the dry and saturated cracked Sci. Yu.Yu., Springer, C.B. and Paterson, mation Geophys. A.N. and Nedra, Stavrogin, A.N. Rocks. of serpentinite solids Int. J. Solids Struct.. of rocks to large stresses, and R.B. Merrill Cratering. Protosenya, Moscow, and Nedra, (Editors), Pergamon Press. New Turcotte, A.G., 1979. Plasticity of Strength of 301 pp. (in Russian). Protosenya, Moscow. A.G., 1985. pp. (in Russian). 450 L., lY9l. Kinks versus shocks. In: R. Fosdick, and M. Slemrod (Editors), Shock in General D.L., 1986. Fractals Induced Media. E. Transi- IMA Vol., and fragmentation. J. Geophys. Res., 91: 1921-1926. A.P. relaxation and Roitburd, in nonlocal A.L.. media. 1988. Nonisothermal Sov. Phys. Solid (FIT), 4: 613. B., Dubois, Fractal patterns Lett., 104: 25-35. Von Karman, Druck. visco-elastic Volarovich, 312: 302-305. solids. fault system. J. Geophys. Appl. Mech., Raleigh, Stavrogin, J., Moore, D. and Touchard, of fractures in granites. T., 1911. Festigkeits Verh. Dtsch. G.. Earth versuche 1991. Planet. Sci. unter all seitigem Ingr., 55: 1749-1757. M.P. (Editor), J. Geophys. Res., 79: Yukutake, H., 1988. Physical 1989. Fracturing from measurements Y.N., 1969. Creep fat of cracked Properties of Rocks. Handbook, Nedra, Moscow. 255 pp. (in Russian). B., 1974. Seismic velocities zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIH in P.G. and Aki, K., 1987. Fractal Andreas for constants York, NY, pp. 657-688. Velde. seismic boundaries hole). Rep. Acad. of strain-dependent Sci. USSR, Explosion Berlin. V.A. and Podladchikov, model R.J. and Budiansky, Rabotnov, and Structures 1987. M., 19YI. A simple technique Pepin Springer-Verlag, 5412-5426. Okubo. Impact R.O. tions and Phase Rusanov, Trans. on the Sec. 1977. The response Rep. Mechanics variation of steels. statistics. UMSI Englewood ot Am. orientation Dumn 296: 71-76. Myasnikov, R.N., Umantsev. M.V., USSR, crack 1991. 713 pp. Mints, arbitrary elastic A., pp. Continuum effective Truskinovsky, V., Podladchikov, of temperature strength C.M. and Kachanov, finding Rocks. Myasnikov, nonlinear composed 174: 777-781. In: D.L. Roddy, Sot., 2: 429-437. Lyakhovsky, Sayers, Schock, tween seismic wave velocity and state of stress. Geophys. Astron. E.L., 1952. Effect rupture of media Philos. May. Ser. 5, 20: 469-481. 27: 67 I-680. 10: 71-75. solid. Phys. Solid Earth, in contact. Mech. Eng.. sian). Lyakhovsky, Robinson, long-term teorii uprugosti Akad. O., 1X85. On the dilatancy rigid particles rupture. geometry in the San 12th Int. Congr. MS., Zienkiewics, neering Berlin. and 1965. Experimental its tectonic Res., 70(16): 3965-3985. during triaxial process of granite and temporal deformation. inferred variations J. Geophys. in Res., 94: 15,639-15,651. Res., 92: 345-355. Proc. velocity of spatial defor- implications. J. Zobak, O., 1971. The Finite Science. McGraw M.D. and Byerlee, rack dilatancy Geophys. Element Method in Engi- Hill, New York, NY, 541 pp. J.D., 1975. The effect on the permeability Res.. 80-B: 7522755. of Westerly of microcgranite. J.