- Email: [email protected]
Modification of the DDA method for rigid block Problems
PII: S0148-9062(97)00319-7
Int. J. Rock Mech. Min. Sci. Vol. 35, No. 6, pp. 683±693, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0148-9062/98 $19.00 + 0.00
Modi®cation of the DDA Method for Rigid Block Problems C. Y. KOO$ J. C. CHERN$ This paper presents the development of the rigid block version of discontinuous deformation analysis and its preliminary application in rockfall simulation. Modi®cations, including the reduction of the error due to large rigid body rotation by using post-correction at each time step of calculation and the incorporation of mass proportional damping to account for energy losses, are made to meet the requirement for large block translation as well as rotation. To improve the program eciency, rigid body displacement function and preconditioned conjugate gradient solver are adopted. The results obtained show that the newly developed code, RIG-DDA, can appropriately and eciently model various modes of rockfall movement and predict the trajectory of falling rock debris, which is helpful in the design of protection measures. # 1998 Elsevier Science Ltd.
INTRODUCTION
Discontinuous deformation analysis (DDA) pioneered by Shi [1] in the late 1980 s has received considerable attention from researchers and practising engineers in recent years. With its eectual use of block kinematics in dealing with the complicated interaction between discrete blocks under ensured equilibrium conditions at anytime by minimizing the total potential energy of the system, it is probably the most suitable method for the analysis of discrete block system involving large movements. Since DDA is a displacement method, it requires six unknown displacement variables to describe the ®rst order linear displacement ®eld in a block. Twelve unknowns are required for the second order displacement ®eld. Therefore, the number of simultaneous equations for the block system in an engineering problem become large. The time marching calculation process involves intricate block contact detection and the solving of a large number of simultaneous equations by iteration at each time step. It becomes prohibitively expensive to solve a real engineering problem involving a large number of blocks. This perhaps may be improved by using a more ecient equation solver. In this study, the preconditioned conjugate gradient method (PCG) is adopted to replace the original successive over-relaxation (SOR) solving technique for improving the program eciency. However, with the large computation task, there is a limit on the improvement of calculation eciency. {Geotechnical Research Center, Sinotech Engineering Consultants, Inc., No. 7, Lane 26, Yat-Sen Road, Taipei 10547, Taiwan, Republic of China.
Accordingly, developing methods of analysis for certain types of problem which require less variables to describe the displacement of the block is desirable. In many engineering problems, the block is either relatively rigid or the stress level is low. The movement of the block system is mainly due to the sliding or rotation of blocks and the deformation of the block itself is relatively small. Stability of rock slopes, rockfall and the behavior of rock®ll belong to this category of problem. The development of a rigid block model can reduce the number of unknown in the simulated problem and make the analysis possible. In this paper, the development of rigid block DDA is described and applications of this method in solving the rockfall problem are also presented.
DISPLACEMENT FUNCTION FOR RIGID BLOCK
The displacements of block due to rigid body translation and rotation can be expressed as: u u0 x ÿ x0 cos r0 ÿ 1 ÿ y ÿ y0 sin r0 v v0 x ÿ x0 sin r0 y ÿ y0 cos r0 ÿ 1
1
in which u0, v0 are the rigid body translations and r0 is the rigid body rotation. When the rigid body rotation is very small, the equation reduces to linear displacement function as follows. u u0 ÿ y ÿ y0 r0 v v0 x ÿ x0 r0
2
The displacement function may be expressed in matrix form:
683
684
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 1. Rotation error due to linear displacement function.
hui T23 D31 v 21 where 2
3 u0 6 7 D 4 v0 5 T
r0 1 0
ÿ y ÿ y0
0 1
x ÿ x0
are the displacement variable matrix and the transformation matrix respectively.
If a line P0P rotates to P0P1 in a rotating block as illustrated in Fig. 1, error due to approximation using linear displacement function in a problem involving large rigid body rotation will cause block expansion. This is also illustrated by the change in block size when it falls from a slope as shown in Fig. 2. Signi®cant expansion in block size can be observed. To solve the problem of block expansion, two approaches have been proposed. 1. the use of exact displacement function incorporating nonlinear terms [2]; 2. the use of linear displacement function and postcorrection [3, 4].
Fig. 2. Block expansion accumulated as a block falls from a slope.
KOO and CHERN: DDA AND RIGID BLOCKS
685
Fig. 3. Reasonable results obtained by using post-correction.
In the ®rst approach, the nonlinear terms of trigonometric function can be expressed in terms of in®nite series as follows. sin r0
r0 r30 r50 ÿ ÿ 1! 3! 5!
cos r0 1 ÿ
r20 r40 ÿ 2! 4!
Due to the exact displacement function shown in Equation (1) including trigonometric functions with in®nite series, the integration of potential energy becomes almost impossible to include many high order terms. However, if the high order terms are omitted, the error will occur. In this study, the second approach, involving simpli®ed linear displacement function and post-correction is adopted. The algorithm of this approach can be divided into two parts. The ®rst part is the adoption of simpli®ed linear displacement function, shown in Equation (2), to derive the potential energy. As a result of using the simpli®ed linear displacement function, the rotation error arising from high order terms omitted will occur. Therefore, to correct the error, post-correction is used in the second part of this approach. The locations of block vertices are re-calculated using the exact displacement function after each time step of calculation. As a result, the diculty of integration can be avoided and reasonable results can be obtained as illustrated in Fig. 3. FORMULATION OF RIGID BLOCK DDA
Based on the rigid body displacement function, the stiness matrices and force matrices due to point load,
body force, line load, inertia force, ®xed point, rock bolt, normal and shear contact forces and damping force are derived by minimizing the total potential energy. Details of the formulation can be found in Ref. [1]. In this paper, only inertia force matrix and damping force matrix are depicted.
(1) Simultaneous equilibrium equations for rigid block system For a system with N blocks, simultaneous equations of equilibrium may be obtained by assembling the stiness matrix and force matrix obtained by minimizing the total potential energy produced by the forces acting on the individual block, and can be expressed in matrix form as follows. 2 3 K11 K12 K13 K1N 6 K21 K22 K23 K2N 7 6 7 6 K31 K32 K33 K3N 7 6 7 6 . .. 7 .. .. .. 4 .. . 5 . . . KN1 KN2 KN3 KNN 3N 3N 2
3 2 3 F1 D1 6 D2 7 6 F2 7 6 7 6 7 6 D3 7 6 F 7 6 3 7 6 7 6 . 7 6 . 7 4 .. 5 4 .. 5 DN 3N1 FN 3N1 Due to rigid body displacement function used, each block has three degrees of freedom and the unknown submatrix [Di] is a 3 1 submatrix containing rigid body translation and rotation. The stiness submatrix
686
KOO and CHERN: DDA AND RIGID BLOCKS
[ Kij] and the force submatrix [Fi] for the given ith block have the matrix size of 3 3 and 3 1. (2) Inertia force Because the DDA method is based on dynamic approach, inertia force plays an essential role in rigid body motion. Newton's second law of motion is used to control the block movement. When a block moves in dynamic mode, its initial velocity is equal to the ®nal velocity of previous time step, whereas in static mode, its initial velocity is set to zero in each Q time step calculation. The total potential energy I of inertia force (Fx, Fy) can be written as follows. 2 2 3 @ u t Z Z Y 6 @t2 7 7dxdy Mu v6 4 2 5 @ v t I @t2 Z Z @Di t dxdy MDi T Ti T Ti @t2 where M mass per unit area of ith block. Omitting high order terms of Taylor series expansion, @2 Di t 2 2 @2 Di t D ÿ i @t2 D @t D2 in which D is the increment of time step. The total potential energy of ith block due to inertia force is: Z Z Y 2M 2M T T Ti Ti dxdy V0 Di Di ÿ D D2 I where [ V0] is the initial velocity of the block. By minimizing the total potential energy, the stiness matrix may be derived in the following term. Z Z 2MK Ti T32 Ti 23 dxdy; r; s 1; 2; 3 Krs D2 which will be added to the global stiness matrix of the block system, and the force matrix Z Z 2MK Ti T32 Ti 23 dxdy V0 ; r; s 1; 2; 3 fr D which will be added to the global force matrix of the block system. (3) Damping force Since DDA method has the capability of analyzing dynamic problems, it requires the damping procedure to make the blocks converge to a stable state. The movement of a block system is often retarded by the energy loss due to friction of blocks with the surrounding media. The eect of dragging force and collision between blocks was considered by Ohnishi et al. [3]. In the method proposed, the damping coecient due to dragging force and velocity reduction factor due to block collision have to be calibrated using ®eld test results. The energy loss due to block collision has been considered by the modi®cation of the stiness of the
contact spring in open±close iteration in the original code. No modi®cation was made in this aspect. In this study, mass proportional damping is used to account for the energy loss due to dragging force. The damping value C takes the form of C aM where a is the damping coecient and M is the mass per unit area. The total potential energy due to damping can be expressed as: Z Z Y Fx ÿ u v dxdy Fy D where damping force matrix can be expressed as: 2@u t3 Fx 6 7 ÿC 4 @t 5 @v t Fy @t The total potential energy of the ith block due to damping force is 2@u t3 Z Z Y 6 7 ÿ aMu v4 @t 5dxdy @v t D @t Z Z @Di t dxdy aMDi T Ti T Ti @t By minimizing the total potential energy, the stiness matrix of the ith block can be obtained as follows: Z Z aM Ti T32 Ti 23 dxdy; r; s 1; 2; 3 Krs D in which D is the time increment of calculation. This stiness matrix is added to the global stiness matrix of the block system. IMPROVED CALCULATION EFFICIENCY
As the block kinematics of DDA involves complicated interaction between discrete blocks, use of the dynamic approach, and a large number of simultaneous equations, computation becomes an enormous task, and it is prohibitively expensive to solve a real engineering problem containing large number of blocks by using the DDA method. In order to reduce the massive computing task, it is necessary to improve the eciency of calculation. In this study, calculation eciency is improved by reducing the number of simultaneous equations and adopting a more ecient solver. The number of simultaneous equations can be reduced by using the rigid body displacement function. Since only rigid body translation and rotation are considered, the matrix size of simultaneous equations is 3N 3N which is four times less than the original matrix size of 6N 6N. The result of adopting rigid body displacement function compared to the original DDA method is shown
KOO and CHERN: DDA AND RIGID BLOCKS
687
Fig. 4. Relationship of CPU time used and block number for one calculation time step.
in Fig. 4. Signi®cant reduction in CPU time may be achieved by adopting rigid block. According to the studies done by Meyer and Valocchi [5], the precondition conjugate gradient method (PCG) is highly ecient when solving a large number of simultaneous equations. Hence, the PCG solver is used to replace the original successive overrelaxation (SOR) solving technique for improving the
program eciency. The inverse of diagonal matrix is the choice for preconditioner in this study, because it is easy to implement and can be surprisingly eective. Figure 5 shows the comparison of eciency between SOR and PCG solver. It can be seen that the PCG solver converges faster than the SOR solver originally used in DDA method. From the results obtained, it is evident that after adopting PCG solver, program e-
Fig. 5. Comparison of SOR and PCG solver eciency for one calculation time step.
688
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 6. Geometry of validation example.
ciency has increased through the reduction of CPU time and also through the decrease in iteration number. Therefore, by adopting the rigid body formulation and a more ecient solver, the amount of calculation may be signi®cantly reduced and it is possible to analyze actual engineering problems dealing with large number of blocks. VALIDATION OF RIGID BLOCK DDA CODE
DDA code (version 1993) provided by Shi [1] is used as the basic program for developing the rigid block code named RIG-DDA. The RIG-DDA code is validated by using an idealized example as shown in Fig. 6. In this example, a rock slope of 92.5 m high is
cut on a 56.68 slope in a layered rock mass dipping at 608 into the hill. A regular system of 16 blocks which is considered as rigid resting on a base stepped at 1 m in every 5 m, and the properties of the rock block are assumed to be j = 38.158 and g = 25 KN/m3. Static mode is used in this analysis, and the normal and shear contact forces of each block are recorded during the calculation. The results obtained by RIG-DDA analysis are shown in Figs 7 and 8. In Fig. 7, it can be seen that good agreements in normal and shear forces at the points of block-slope contacts are obtained. The three predicted modes of movement, that is stable, toppling and sliding modes as shown in Fig. 8, are also the
Fig. 7. Comparison of normal and shear forces with exact solutions.
KOO and CHERN: DDA AND RIGID BLOCKS
689
Fig. 8. Results obtained from RIG-DDA.
Table 1. After Hoek and Bray [6] Block No. 16 15 14 13
Mode Stable Stable Stable Toppling
Block No. 12 11 10 9
Mode
Block No.
Mode
Block No.
Mode
Toppling Toppling Toppling Toppling
8 7 6 5
Toppling Toppling Toppling Toppling
4 3 2 1
Toppling Sliding Sliding Sliding
Fig. 9. The movement modes of rockfall observed by Ritchie [12].
690
KOO and CHERN: DDA AND RIGID BLOCKS
been developed. However, most of the computer program adopt simple models, such as spherical boulder, initial horizontal and vertical velocities, boulder obeying conventional law of re¯ection with the possibility of dierent normal and tangential coecients for bounce (Branwer [10]). Complicated interaction between adjacent blocks could not be considered. Discontinuous models such as distinct element method (DEM, Cundall [11]) and DDA can model the complicated interaction between falling block, and have been used to simulate the rockfall problems. Rigid block DDA with improved computation eciency described herein includes the kinematics of all modes of motion and the interaction between blocks. It is possible to simulate the movement of rock mass containing numerous distinct blocks along an irregular slope surface in an acceptable computing time. To investigate the capability of RIG-DDA in simulating the rockfall problem, dierent modes of movement of a single block on slopes with various inclination angles were calculated. The predominant modes of rockfall motion shown in Fig. 9, observed by Ritchie [12], are as follows: 1. rolling on slope at <458 2. bouncing on slope in the range of 45±758 3. falling on slope steeper than 758.
Fig. 10. The simulation of dierent movement modes on slopes with various inclinations by RIG-DDA (a) rolling mode, (b) bouncing mode and (c) falling mode.
same as those obtained by analytical method which are shown in Table 1. APPLICATION TO ROCKFALL ANALYSIS
The design of protection measures against rockfall requires knowledge of block moving trajectory as well as the energy produced. In general, three approaches, that is, ®eld test, numerical modeling and empirical method, are often used. In the numerical approach various computer models, such as Pfeier and Higgins [7], Hoek [8] and Kobayashi et al. [9] have
From the calculated trajectories of block movement on slopes with various inclinations, as shown in Fig. 10, the results are consistent with the observation made by Ritchie [12]. The second example used is taken from a paper by Chen et al. [13]. It is an actual rockfall site on ShuHwa highway in eastern Taiwan. The approximate distribution of fallen rock debris remained on roadway and near the toe of the slope has been described, which enables the calibration of energy dissipation parameter required in the calculation. It is a 130 m high rock slope with a highway passing near the toe of the slope. The geometry of the rock slope and the location of highway are illustrated in Fig. 11. Near the top of the slope, heavily jointed, weathered rock mass containing numerous rock blocks can be found. In this study, the block mesh was generated by the DDA preprocessor with the average block size ca 2 m2 similar to the actual condition in the site, and was illustrated as shown in Fig. 11. The joint strength was assumed to be c = 0 and f = 288. The mass proportional damping coecient a = 0.6 was used by a trial and error method to match the ®eld observation on the distribution of fallen rock blocks. The calculated falling mass during the falling process and the ®nal positions of rock blocks are illustrated in Fig. 12Fig. 13 respectively. The distribution of debris obtained is quite similar to those observed in the ®eld. To mitigate the rockfall hazard, a numerical investigation has been made. Based on the same boundary conditions and parameters used above, a retaining structure was assumed to be constructed at a
KOO and CHERN: DDA AND RIGID BLOCKS
691
Fig. 11. Geometry of rock slope and location of highway.
location near the toe of the high slope. Computed using the RIG-DDA code, the results shown in Fig. 14 demonstrate that the rock blocks can be retained by the structure. CONCLUSIONS
In this paper, the development of the rigid block version of RIG-DDA is presented. The program has
been validated and applied to simulate the rockfall problem. From the results of these studies, the following conclusions may be drawn. 1. The new code adopts the rigid body displacement function and incorporates the damping eect in the formulation. With improved calculation eciency, RIG-DDA has expanded the capability of DDA in modeling in the low stress, small deformation, and
Fig. 12. Blocks in the falling process.
692
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 13. Final position of rock blocks.
Fig. 14. Final position of falling rock blocks stopped by retaining wall.
large displacement engineering problems involving large number of blocks, especially in dynamic analysis. 2. The adoption of linear displacement function with post-correction can eectively eliminate the error due to large rigid body rotation and avoid the diculty of integration of potential energy in using exact displacement function. 3. From the results of simulating the rockfall problem, various modes of block movement as those
observed in the ®eld may be correctly predicted. Therefore, it makes the analysis of an extremely complicated rockfall problem possible.
Accepted for publication 17 September 1997
REFERENCES 1. Shi, G. H. and Goodman, R. E., Int. J. Numerical Analyt. Methods Geomech., 1989, 13, 359±380.
KOO and CHERN: DDA AND RIGID BLOCKS 2. Maclaughlin, M. M. and Sitar, N., Rigid body rotations in DDAProceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, 620. Berkeley, CA, 1996. 3. Ohnishi, Y., Yamamukai, K. and Chen, G. Q., Application of DDA in rockfall analysisProceedings of the 2nd North American Rock Mechanics Symposium. Montreal, Quebec, Canada, 1996, pp. 2031±2037. 4. Ke, T. C., The issue of rigid-body rotation in DDAProceedings of the ®rst International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media. Berkeley, CA, U.S.A., 1996, pp. 318±325. 5. Meyer, P. D. and Valocchi, A. J. A., Water Resources Res., 1989, 25, 1440±1446. 6. Hoek, E. and Bray, J. W., Rock Slope Engineering. Institute of Mining and Metallurgy, London, 1981.
693
7. Pfeier, T. J. and Higgins, J. A., Rockfall hazard analysis using the Colorado rockfall simulation program. TRB Transportation Research Board, 1990. 8. Hoek, E., RockfallÐA Program in Basic for the Analysis of Rockfalls from Slopes. Golder Associates, Vancouver, B.C., 1987. 9. Kobayashi, Y., Harp, E. L. and Kagawa, T., Rock Mech. Rock Eng., 1990, 23, 1±20. 10. Branwer, C. O., Participant Workbook For Rockfall Hazard Mitigation Methods. National Highway Institute, U.S.A, 1994. 11. Cundall, P. A., Rock Fracture, Proceedings of International Symposium Rock Fracture, 1971, Nancy, 11±18. 12. Ritchie, A. M., Highway Record, 1963, 17, 13±28. 13. Chen, H. G., Chen, R. H. and Huang, T. H., Journal of Civil and Hydraulic Engineering, 1993, 20, 19±30.
Int. J. Rock Mech. Min. Sci. Vol. 35, No. 6, pp. 683±693, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0148-9062/98 $19.00 + 0.00
Modi®cation of the DDA Method for Rigid Block Problems C. Y. KOO$ J. C. CHERN$ This paper presents the development of the rigid block version of discontinuous deformation analysis and its preliminary application in rockfall simulation. Modi®cations, including the reduction of the error due to large rigid body rotation by using post-correction at each time step of calculation and the incorporation of mass proportional damping to account for energy losses, are made to meet the requirement for large block translation as well as rotation. To improve the program eciency, rigid body displacement function and preconditioned conjugate gradient solver are adopted. The results obtained show that the newly developed code, RIG-DDA, can appropriately and eciently model various modes of rockfall movement and predict the trajectory of falling rock debris, which is helpful in the design of protection measures. # 1998 Elsevier Science Ltd.
INTRODUCTION
Discontinuous deformation analysis (DDA) pioneered by Shi [1] in the late 1980 s has received considerable attention from researchers and practising engineers in recent years. With its eectual use of block kinematics in dealing with the complicated interaction between discrete blocks under ensured equilibrium conditions at anytime by minimizing the total potential energy of the system, it is probably the most suitable method for the analysis of discrete block system involving large movements. Since DDA is a displacement method, it requires six unknown displacement variables to describe the ®rst order linear displacement ®eld in a block. Twelve unknowns are required for the second order displacement ®eld. Therefore, the number of simultaneous equations for the block system in an engineering problem become large. The time marching calculation process involves intricate block contact detection and the solving of a large number of simultaneous equations by iteration at each time step. It becomes prohibitively expensive to solve a real engineering problem involving a large number of blocks. This perhaps may be improved by using a more ecient equation solver. In this study, the preconditioned conjugate gradient method (PCG) is adopted to replace the original successive over-relaxation (SOR) solving technique for improving the program eciency. However, with the large computation task, there is a limit on the improvement of calculation eciency. {Geotechnical Research Center, Sinotech Engineering Consultants, Inc., No. 7, Lane 26, Yat-Sen Road, Taipei 10547, Taiwan, Republic of China.
Accordingly, developing methods of analysis for certain types of problem which require less variables to describe the displacement of the block is desirable. In many engineering problems, the block is either relatively rigid or the stress level is low. The movement of the block system is mainly due to the sliding or rotation of blocks and the deformation of the block itself is relatively small. Stability of rock slopes, rockfall and the behavior of rock®ll belong to this category of problem. The development of a rigid block model can reduce the number of unknown in the simulated problem and make the analysis possible. In this paper, the development of rigid block DDA is described and applications of this method in solving the rockfall problem are also presented.
DISPLACEMENT FUNCTION FOR RIGID BLOCK
The displacements of block due to rigid body translation and rotation can be expressed as: u u0 x ÿ x0 cos r0 ÿ 1 ÿ y ÿ y0 sin r0 v v0 x ÿ x0 sin r0 y ÿ y0 cos r0 ÿ 1
1
in which u0, v0 are the rigid body translations and r0 is the rigid body rotation. When the rigid body rotation is very small, the equation reduces to linear displacement function as follows. u u0 ÿ y ÿ y0 r0 v v0 x ÿ x0 r0
2
The displacement function may be expressed in matrix form:
683
684
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 1. Rotation error due to linear displacement function.
hui T23 D31 v 21 where 2
3 u0 6 7 D 4 v0 5 T
r0 1 0
ÿ y ÿ y0
0 1
x ÿ x0
are the displacement variable matrix and the transformation matrix respectively.
If a line P0P rotates to P0P1 in a rotating block as illustrated in Fig. 1, error due to approximation using linear displacement function in a problem involving large rigid body rotation will cause block expansion. This is also illustrated by the change in block size when it falls from a slope as shown in Fig. 2. Signi®cant expansion in block size can be observed. To solve the problem of block expansion, two approaches have been proposed. 1. the use of exact displacement function incorporating nonlinear terms [2]; 2. the use of linear displacement function and postcorrection [3, 4].
Fig. 2. Block expansion accumulated as a block falls from a slope.
KOO and CHERN: DDA AND RIGID BLOCKS
685
Fig. 3. Reasonable results obtained by using post-correction.
In the ®rst approach, the nonlinear terms of trigonometric function can be expressed in terms of in®nite series as follows. sin r0
r0 r30 r50 ÿ ÿ 1! 3! 5!
cos r0 1 ÿ
r20 r40 ÿ 2! 4!
Due to the exact displacement function shown in Equation (1) including trigonometric functions with in®nite series, the integration of potential energy becomes almost impossible to include many high order terms. However, if the high order terms are omitted, the error will occur. In this study, the second approach, involving simpli®ed linear displacement function and post-correction is adopted. The algorithm of this approach can be divided into two parts. The ®rst part is the adoption of simpli®ed linear displacement function, shown in Equation (2), to derive the potential energy. As a result of using the simpli®ed linear displacement function, the rotation error arising from high order terms omitted will occur. Therefore, to correct the error, post-correction is used in the second part of this approach. The locations of block vertices are re-calculated using the exact displacement function after each time step of calculation. As a result, the diculty of integration can be avoided and reasonable results can be obtained as illustrated in Fig. 3. FORMULATION OF RIGID BLOCK DDA
Based on the rigid body displacement function, the stiness matrices and force matrices due to point load,
body force, line load, inertia force, ®xed point, rock bolt, normal and shear contact forces and damping force are derived by minimizing the total potential energy. Details of the formulation can be found in Ref. [1]. In this paper, only inertia force matrix and damping force matrix are depicted.
(1) Simultaneous equilibrium equations for rigid block system For a system with N blocks, simultaneous equations of equilibrium may be obtained by assembling the stiness matrix and force matrix obtained by minimizing the total potential energy produced by the forces acting on the individual block, and can be expressed in matrix form as follows. 2 3 K11 K12 K13 K1N 6 K21 K22 K23 K2N 7 6 7 6 K31 K32 K33 K3N 7 6 7 6 . .. 7 .. .. .. 4 .. . 5 . . . KN1 KN2 KN3 KNN 3N 3N 2
3 2 3 F1 D1 6 D2 7 6 F2 7 6 7 6 7 6 D3 7 6 F 7 6 3 7 6 7 6 . 7 6 . 7 4 .. 5 4 .. 5 DN 3N1 FN 3N1 Due to rigid body displacement function used, each block has three degrees of freedom and the unknown submatrix [Di] is a 3 1 submatrix containing rigid body translation and rotation. The stiness submatrix
686
KOO and CHERN: DDA AND RIGID BLOCKS
[ Kij] and the force submatrix [Fi] for the given ith block have the matrix size of 3 3 and 3 1. (2) Inertia force Because the DDA method is based on dynamic approach, inertia force plays an essential role in rigid body motion. Newton's second law of motion is used to control the block movement. When a block moves in dynamic mode, its initial velocity is equal to the ®nal velocity of previous time step, whereas in static mode, its initial velocity is set to zero in each Q time step calculation. The total potential energy I of inertia force (Fx, Fy) can be written as follows. 2 2 3 @ u t Z Z Y 6 @t2 7 7dxdy Mu v6 4 2 5 @ v t I @t2 Z Z @Di t dxdy MDi T Ti T Ti @t2 where M mass per unit area of ith block. Omitting high order terms of Taylor series expansion, @2 Di t 2 2 @2 Di t D ÿ i @t2 D @t D2 in which D is the increment of time step. The total potential energy of ith block due to inertia force is: Z Z Y 2M 2M T T Ti Ti dxdy V0 Di Di ÿ D D2 I where [ V0] is the initial velocity of the block. By minimizing the total potential energy, the stiness matrix may be derived in the following term. Z Z 2MK Ti T32 Ti 23 dxdy; r; s 1; 2; 3 Krs D2 which will be added to the global stiness matrix of the block system, and the force matrix Z Z 2MK Ti T32 Ti 23 dxdy V0 ; r; s 1; 2; 3 fr D which will be added to the global force matrix of the block system. (3) Damping force Since DDA method has the capability of analyzing dynamic problems, it requires the damping procedure to make the blocks converge to a stable state. The movement of a block system is often retarded by the energy loss due to friction of blocks with the surrounding media. The eect of dragging force and collision between blocks was considered by Ohnishi et al. [3]. In the method proposed, the damping coecient due to dragging force and velocity reduction factor due to block collision have to be calibrated using ®eld test results. The energy loss due to block collision has been considered by the modi®cation of the stiness of the
contact spring in open±close iteration in the original code. No modi®cation was made in this aspect. In this study, mass proportional damping is used to account for the energy loss due to dragging force. The damping value C takes the form of C aM where a is the damping coecient and M is the mass per unit area. The total potential energy due to damping can be expressed as: Z Z Y Fx ÿ u v dxdy Fy D where damping force matrix can be expressed as: 2@u t3 Fx 6 7 ÿC 4 @t 5 @v t Fy @t The total potential energy of the ith block due to damping force is 2@u t3 Z Z Y 6 7 ÿ aMu v4 @t 5dxdy @v t D @t Z Z @Di t dxdy aMDi T Ti T Ti @t By minimizing the total potential energy, the stiness matrix of the ith block can be obtained as follows: Z Z aM Ti T32 Ti 23 dxdy; r; s 1; 2; 3 Krs D in which D is the time increment of calculation. This stiness matrix is added to the global stiness matrix of the block system. IMPROVED CALCULATION EFFICIENCY
As the block kinematics of DDA involves complicated interaction between discrete blocks, use of the dynamic approach, and a large number of simultaneous equations, computation becomes an enormous task, and it is prohibitively expensive to solve a real engineering problem containing large number of blocks by using the DDA method. In order to reduce the massive computing task, it is necessary to improve the eciency of calculation. In this study, calculation eciency is improved by reducing the number of simultaneous equations and adopting a more ecient solver. The number of simultaneous equations can be reduced by using the rigid body displacement function. Since only rigid body translation and rotation are considered, the matrix size of simultaneous equations is 3N 3N which is four times less than the original matrix size of 6N 6N. The result of adopting rigid body displacement function compared to the original DDA method is shown
KOO and CHERN: DDA AND RIGID BLOCKS
687
Fig. 4. Relationship of CPU time used and block number for one calculation time step.
in Fig. 4. Signi®cant reduction in CPU time may be achieved by adopting rigid block. According to the studies done by Meyer and Valocchi [5], the precondition conjugate gradient method (PCG) is highly ecient when solving a large number of simultaneous equations. Hence, the PCG solver is used to replace the original successive overrelaxation (SOR) solving technique for improving the
program eciency. The inverse of diagonal matrix is the choice for preconditioner in this study, because it is easy to implement and can be surprisingly eective. Figure 5 shows the comparison of eciency between SOR and PCG solver. It can be seen that the PCG solver converges faster than the SOR solver originally used in DDA method. From the results obtained, it is evident that after adopting PCG solver, program e-
Fig. 5. Comparison of SOR and PCG solver eciency for one calculation time step.
688
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 6. Geometry of validation example.
ciency has increased through the reduction of CPU time and also through the decrease in iteration number. Therefore, by adopting the rigid body formulation and a more ecient solver, the amount of calculation may be signi®cantly reduced and it is possible to analyze actual engineering problems dealing with large number of blocks. VALIDATION OF RIGID BLOCK DDA CODE
DDA code (version 1993) provided by Shi [1] is used as the basic program for developing the rigid block code named RIG-DDA. The RIG-DDA code is validated by using an idealized example as shown in Fig. 6. In this example, a rock slope of 92.5 m high is
cut on a 56.68 slope in a layered rock mass dipping at 608 into the hill. A regular system of 16 blocks which is considered as rigid resting on a base stepped at 1 m in every 5 m, and the properties of the rock block are assumed to be j = 38.158 and g = 25 KN/m3. Static mode is used in this analysis, and the normal and shear contact forces of each block are recorded during the calculation. The results obtained by RIG-DDA analysis are shown in Figs 7 and 8. In Fig. 7, it can be seen that good agreements in normal and shear forces at the points of block-slope contacts are obtained. The three predicted modes of movement, that is stable, toppling and sliding modes as shown in Fig. 8, are also the
Fig. 7. Comparison of normal and shear forces with exact solutions.
KOO and CHERN: DDA AND RIGID BLOCKS
689
Fig. 8. Results obtained from RIG-DDA.
Table 1. After Hoek and Bray [6] Block No. 16 15 14 13
Mode Stable Stable Stable Toppling
Block No. 12 11 10 9
Mode
Block No.
Mode
Block No.
Mode
Toppling Toppling Toppling Toppling
8 7 6 5
Toppling Toppling Toppling Toppling
4 3 2 1
Toppling Sliding Sliding Sliding
Fig. 9. The movement modes of rockfall observed by Ritchie [12].
690
KOO and CHERN: DDA AND RIGID BLOCKS
been developed. However, most of the computer program adopt simple models, such as spherical boulder, initial horizontal and vertical velocities, boulder obeying conventional law of re¯ection with the possibility of dierent normal and tangential coecients for bounce (Branwer [10]). Complicated interaction between adjacent blocks could not be considered. Discontinuous models such as distinct element method (DEM, Cundall [11]) and DDA can model the complicated interaction between falling block, and have been used to simulate the rockfall problems. Rigid block DDA with improved computation eciency described herein includes the kinematics of all modes of motion and the interaction between blocks. It is possible to simulate the movement of rock mass containing numerous distinct blocks along an irregular slope surface in an acceptable computing time. To investigate the capability of RIG-DDA in simulating the rockfall problem, dierent modes of movement of a single block on slopes with various inclination angles were calculated. The predominant modes of rockfall motion shown in Fig. 9, observed by Ritchie [12], are as follows: 1. rolling on slope at <458 2. bouncing on slope in the range of 45±758 3. falling on slope steeper than 758.
Fig. 10. The simulation of dierent movement modes on slopes with various inclinations by RIG-DDA (a) rolling mode, (b) bouncing mode and (c) falling mode.
same as those obtained by analytical method which are shown in Table 1. APPLICATION TO ROCKFALL ANALYSIS
The design of protection measures against rockfall requires knowledge of block moving trajectory as well as the energy produced. In general, three approaches, that is, ®eld test, numerical modeling and empirical method, are often used. In the numerical approach various computer models, such as Pfeier and Higgins [7], Hoek [8] and Kobayashi et al. [9] have
From the calculated trajectories of block movement on slopes with various inclinations, as shown in Fig. 10, the results are consistent with the observation made by Ritchie [12]. The second example used is taken from a paper by Chen et al. [13]. It is an actual rockfall site on ShuHwa highway in eastern Taiwan. The approximate distribution of fallen rock debris remained on roadway and near the toe of the slope has been described, which enables the calibration of energy dissipation parameter required in the calculation. It is a 130 m high rock slope with a highway passing near the toe of the slope. The geometry of the rock slope and the location of highway are illustrated in Fig. 11. Near the top of the slope, heavily jointed, weathered rock mass containing numerous rock blocks can be found. In this study, the block mesh was generated by the DDA preprocessor with the average block size ca 2 m2 similar to the actual condition in the site, and was illustrated as shown in Fig. 11. The joint strength was assumed to be c = 0 and f = 288. The mass proportional damping coecient a = 0.6 was used by a trial and error method to match the ®eld observation on the distribution of fallen rock blocks. The calculated falling mass during the falling process and the ®nal positions of rock blocks are illustrated in Fig. 12Fig. 13 respectively. The distribution of debris obtained is quite similar to those observed in the ®eld. To mitigate the rockfall hazard, a numerical investigation has been made. Based on the same boundary conditions and parameters used above, a retaining structure was assumed to be constructed at a
KOO and CHERN: DDA AND RIGID BLOCKS
691
Fig. 11. Geometry of rock slope and location of highway.
location near the toe of the high slope. Computed using the RIG-DDA code, the results shown in Fig. 14 demonstrate that the rock blocks can be retained by the structure. CONCLUSIONS
In this paper, the development of the rigid block version of RIG-DDA is presented. The program has
been validated and applied to simulate the rockfall problem. From the results of these studies, the following conclusions may be drawn. 1. The new code adopts the rigid body displacement function and incorporates the damping eect in the formulation. With improved calculation eciency, RIG-DDA has expanded the capability of DDA in modeling in the low stress, small deformation, and
Fig. 12. Blocks in the falling process.
692
KOO and CHERN: DDA AND RIGID BLOCKS
Fig. 13. Final position of rock blocks.
Fig. 14. Final position of falling rock blocks stopped by retaining wall.
large displacement engineering problems involving large number of blocks, especially in dynamic analysis. 2. The adoption of linear displacement function with post-correction can eectively eliminate the error due to large rigid body rotation and avoid the diculty of integration of potential energy in using exact displacement function. 3. From the results of simulating the rockfall problem, various modes of block movement as those
observed in the ®eld may be correctly predicted. Therefore, it makes the analysis of an extremely complicated rockfall problem possible.
Accepted for publication 17 September 1997
REFERENCES 1. Shi, G. H. and Goodman, R. E., Int. J. Numerical Analyt. Methods Geomech., 1989, 13, 359±380.
KOO and CHERN: DDA AND RIGID BLOCKS 2. Maclaughlin, M. M. and Sitar, N., Rigid body rotations in DDAProceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, 620. Berkeley, CA, 1996. 3. Ohnishi, Y., Yamamukai, K. and Chen, G. Q., Application of DDA in rockfall analysisProceedings of the 2nd North American Rock Mechanics Symposium. Montreal, Quebec, Canada, 1996, pp. 2031±2037. 4. Ke, T. C., The issue of rigid-body rotation in DDAProceedings of the ®rst International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media. Berkeley, CA, U.S.A., 1996, pp. 318±325. 5. Meyer, P. D. and Valocchi, A. J. A., Water Resources Res., 1989, 25, 1440±1446. 6. Hoek, E. and Bray, J. W., Rock Slope Engineering. Institute of Mining and Metallurgy, London, 1981.
693
7. Pfeier, T. J. and Higgins, J. A., Rockfall hazard analysis using the Colorado rockfall simulation program. TRB Transportation Research Board, 1990. 8. Hoek, E., RockfallÐA Program in Basic for the Analysis of Rockfalls from Slopes. Golder Associates, Vancouver, B.C., 1987. 9. Kobayashi, Y., Harp, E. L. and Kagawa, T., Rock Mech. Rock Eng., 1990, 23, 1±20. 10. Branwer, C. O., Participant Workbook For Rockfall Hazard Mitigation Methods. National Highway Institute, U.S.A, 1994. 11. Cundall, P. A., Rock Fracture, Proceedings of International Symposium Rock Fracture, 1971, Nancy, 11±18. 12. Ritchie, A. M., Highway Record, 1963, 17, 13±28. 13. Chen, H. G., Chen, R. H. and Huang, T. H., Journal of Civil and Hydraulic Engineering, 1993, 20, 19±30.