Chapter 1 – Introduction

Overview

Fractures at some scale are pervasive in all rock masses exposed at or near the earth surfaces. Pervasive fracturing induces the formation of rock blocks, which are discrete bodies of rock within the rock mass that mechanically interact along the fracture planes subdividing the rock mass. Mathematical models of rock masses must account for the fracture geometry and discrete nature of blocky rock masses, and must be capable of accurately simulating physically viable rock mass failure. Finite element analyses, which may provide useful measures of in-situ stress, are limited to problems where a rock mass deforms slowly and has few material discontinuities. However, large deformation in fractured rock masses is controlled by the geometry of the fracturing (sitar:n1997a, maclaughlin:mm2001a, sitar:n2002). Because analytic solutions do not exist for systems with large numbers of discrete elements, numerical methods implemented using computers provide a relatively inexpensive alternative for conducting simulations of rock mass kinematics (cundall:pa1971, cundall:pa1988, hatzor:y1998, kim:yi1999a).

Discrete element methods

As defined by Cundall (cf. bardet:jp1997), discrete element methods require: 1. finite displacement and rotation of particles, including interparticle contact disengagement, and 2. new particle contacts are recognized algorithmically as the simulation proceeds. Discrete element methods are useful when the interaction of large numbers of particles is of more interest than detailed knowledge of the internal stress state of the particle. Internal stress states are less important in materials such as sand or rock, which are relatively stiff compared to translation and rotation. When necessary, mean measures of stress may be determined, often sufficient due to the propensity of stiff materials for brittle rather than ductile behavior.

In rock engineering, an early discrete element method was introduced as the distinct element method (DEM) by Cundall (cundall:pa1971). The distinct element method starts from Newton’s equations of motion, using a centered difference formulation for time integration and a numerical penalty applied to control penetration of contacting blocks. Penalties are applied through the force vector; contacts between blocks are not coupled through the stiffness matrix. An advantage of the distinct element method is that the formulation is conceptually simple to understand and relatively easy to implement, especially with circular or spherical elements. DEM or DEM-like techniques are probably the most popular discrete element methods among engineers, especially for soil mechanics and granular materials research. Disadvantages of DEM include requiring small time steps and/or numerical damping to ensure stability thus convergence for the central difference time integration.

Shi (shi:gh1988a) introduced the discontinuous deformation analysis (DDA), basing the equations of motion on a minimum potential energy formulation using an implicit method for time integration with numerical penalties applied to contacting blocks. DDA provides a numerical scheme for simultaneous translation, rotation and homogeneous deformation of a physical body, bridging rigid body (analytical) dynamics and continuum mechanics (\fig~\ref{fig:ddavenn}). Unlike Cundall’s DEM, the block-to-block contacts in DDA are coupled through a stiffness matrix, thus requiring an implicit integration scheme which also allows relatively large time steps compared to DEM. Numerical damping in DDA is implicit, induced by the time integration scheme, and by the block contact formulation. The primary disadvantage of DDA is that the coupled contact formulation and implicit solution algorithm require expensive factorization and iteration for solution, and is more difficult both to understand and to implement than DEM. The matrix form of the coupled equations is sparse (in larger problems) but unstructured due to contact conditions that may change at successive time steps, thus changing the structure of the stiffness matrix. DDA is most popular among rock mechanics practitioners (e.g., lin:ct1996, hatzor:y1998, bicanic:n2001), but has also found application in soil mechanics and granular materials studies~\cite{ke:tc1995, thomas:pa1999, rein:gsy2001, osullivan:c2001}. A recent overview by Ma~\citeyear{ma:my1999} lists many extensions to DDA that will not be discussed here. The mathematical structure and numerical accuracy of either method has not been satisfactorily resolved, in part because classical mechanics has resolved neither contact mechanics nor friction~\cite{chatterjee:a1998, stewart:de2000a}.

Currently, DDA is being used for engineering analysis on a number of high-profile projects, including the Three Gorges Dam project in the People’s Republic of China~\shortcite{dong:x1996, zhu:w1999}, Pueblo Dam, Colorado~\cite{kottenstette:jt1999}, the Yerba Buena tunnel portal, San Francisco California~\cite{law:hk1999}, Norway’s Gj\/ovik Olympic Cavern~\shortcite{scheldt:t2002} and Israel’s Masada National Monument~\shortcite{hatzor:yh2002}. Other, more routine, projects using DDA include: performance of tunnel excavations in gravels~\shortcite{chen:s96}, stability of a laminated voussoir beam in an archeaologically significant excavation~\cite{hatzor:y1998} and modeling fault motion in northern China~\cite{liu:l1996}. These projects have significant public interest, and some have life-safety considerations. Assessing the results of these analysis requires a thorough understanding of the structure and numerical accuracy of the DDA method.

\begin{figure}
\begin{center}
\includegraphics[width=3.0in]{figs/ddavenn.eps}
\caption{Venn diagram illustrating the DDA problem domain.}
\label{fig:ddavenn}
\end{center}
\end{figure}

The art and science of validation

As already mentioned, the original DDA formulation has been used for a wide range of problems, especially slope stability. Attempts to extend the original DDA formulation used for this thesis to include capabilities of general geotechnical interest, such as rock bolting, pore pressure, seismic loading, etc. were only partially successful. Any individual analysis using these extensions could be solved in a physically plausible manner for short times and small number of time steps, by adjusting any of several analysis parameters such as time step size or penalty value, or material properties such as Young’s modulus, viscous damping and joint (element interface) friction. However, the DDA formulation appears to be susceptible to numerical instability when applied to large-scale dynamic problems, especially when run for long times, or large number of time steps. Thus the motivation and purpose of this research are the issues of numerical stability and accuracy of DDA in the context of large dynamic systems. Since non-physical parameters such as time step sizes, penalty values, etc. are unavoidable for numerical simulation, developing a suite of simple cases with known analytic solutions (or solutions well-approximated by other means) appears to be a necessary step on the way to understanding the behavior of larger multiblock simulations.

Previous validation studies

A growing suite of validation studies for DDA has emerged over the last decade. Yeung~\citeyear{yeung:mcr1991} studied frictional sliding and toppling, a three-hinged beam problem, and the behavior of rock bolts in a mine roof, finding that DDA appeared to capture physically viable behavior. MacLaughlin~\citeyear{maclaughlin:mm1997} studied sliding friction for short displacements using straight and doubly sloped inclines, finding that DDA computed displacements consistently within a few percent of analytic solutions. Hatzor and Feintuch~\citeyear{hatzor:yh2001} studied short (meter scale) displacements of a sliding block subjected to cyclic loading to investigate the applicability of DDA to seismic analysis. Their results indicate that the computed solution ultimately diverges, but in many cases the rate of divergence is slow enough to allow DDA some use for seismic site investigation in blocky rock masses. Bi{\’c}ani{\’c} and Stirling~\citeyear{bicanic:n2001} proposed the Couplet$/$Heyman arch for a DDA benchmark problem. Their analysis found that DDA appeared to reproduce the behavior expected from the analytic solution for arch stability.

However, most of the existing validation for discrete element methods is more qualitative than quantitative, typically plots comparing results to analytic solutions (e.g.,~\shortciteNP{bangash:t2002}). Moreover, most of these exercises were conducted for small time steps, for short run times, short displacements, or short number of time steps (on the order of 100s to 1000s of time steps). The behavior of the limited validation is then assumed to hold for more demanding conditions.

While the behavior in the first few time steps provides information on how the actual simulation is perturbed from the mathematical model, running simulations for very long times helps to expose instabilities due to cumulative error, if present. In some cases, initial perturbations may “damp out” as the system evolves in time, providing for covergence. In cases of divergence, if the rate of divergence can be bounded over the time scale of interest, or can be shown to diverge sufficiently slowly, the simulation may well be of practical use (e.g.,~\shortciteNP{tsesarsky:m2002}).

Two components to validation are algorithmic approximation to a physical model, and validation of the model to the data. In the first case, the algorithm may approximate the mathematical model of some phenomenon. In discrete element methods, this is usually the case. For example, a mathematical block sliding on a mathematically inclined plane with frictional resistance turns out to have a number of issues associated with algorithmic implementation, in part due to enforcing a frictional, no-penetration contact condition between the plane and the sliding block. One popular method, used in many discrete element methods, is to allow a small amount of non-physical penetration between elements, then impose a penalty force, usually modeled as a linearly elastic spring, to push the elements apart. As a result, the algorithm approximates the mathematical model.

One of the best techniques for validating an algorithmic representation of a mathematical model is to compare an analytic solution to the numerical output. This is usually possible only for systems of unrealistic simplicity. When analytical solutions are not available, other quantities, such as kinetic and potential energy, may be required to be conserved. The ability of an algorithm to conserve such quantities may be used (with caution, {\em cf.}~\citeNP{ortiz:m1986}) as an indicator of efficacy. For problems amenable to finite element analysis, solution accuracy may be investigated by increasing the amount of discretization.

Once an algorithmic representation is validated, model validation may be performed. Typically, this involves comparing the simulation results to laboratory test data, or to field data such as landslide runout. Insofar as the numerical simulation results are satisfactory, the simulation can be said to have predictive capability. If the software implementation is verifiably correct, poor comparisons between the software and the data suggest that the mathematical model is incorrect.

Pr{\’e}cis

The goal of this dissertation is twofold. The first goal is to explore in detail the mathematical structure of DDA in order to understand the inherent mathematical simplifications of the physical problem. The second goal is to experimentally examine the numerical behavior of the DDA formulation applied to simple problems for which analytic or semi-analytic solutions exist. Both of these investigations are necessary for a complete understanding of the DDA method. Simulation behavior of the model problems is examined by varying both material properties (e.g., Young’s modulus, friction angle) and analysis parameters (e.g., timestep, penalty values). The analyses are conducted systematically to uncover trends in the data. These trends in turn form an empirical basis for future development of analytical tools to estimate the accuracy of DDA simulations. All of the analyses are described in detail sufficient for independently verifying the results.

The overall contributions of this dissertation include:

• Developing a framework applicable for both DDA and the distinct element method (DEM), in two- and three-dimensional form with either rigid body or first order deformation, as a result of:
  • establishing DDA as a numerical method for the theory of pseudo-rigid bodies.
  • expressing DDA in terms of a Lagrangian functional rendered stationary in an integral variational statement.
  • explicitly incorporating initial conditions into the integral variational statement.
  • deriving an algorithm for open-close iteration which provides no-penetration constraints between blocks.
  • constructing an explicit expression for kinetic energy.
  • verifying that the mixed dimensional structure (displacements and strains) of the DDA method is consistent.
• Verifying that DDA can be an effective method for simulating large displacement frictional sliding.
• Verifying that DDA captures contact forces at equilibrium conditions, at least in an average sense.
• Developing techniques for assessing the accuracy of DDA results.
• Exposing limitations on collisional and impact capability.
• Developing a method for analyzing displacement dependent properties.

Chapter references

[bibtex file=dem.bib key=cundall:pa1971]
[bibtex file=dem.bib key=cundall:pa1988]
[bibtex file=dda.bib key=hatzor:yh1998]
[bibtex file=dda.bib key=kim:yi1999a]

bardet:jp1997
lin:ct1996,
bicanic:n2001
ke:tc1995,
thomas:pa1999,
rein:gsy2001,
osullivan:c2001