General Comments

As described in the first section of 3DEC Theory and Background, 3DEC uses a dynamic algorithm for problem solution. Natural dynamic systems contain some degree of damping of the vibrational energy within the system; otherwise, the system would oscillate indefinitely when subjected to driving forces. Damping is due, in part, to energy loss as a result of slippage along contacts of blocks within the system, internal friction loss in the intact material, and any resistance caused by air or fluids surrounding the structure. 3DEC is used to solve two general classes of mechanical problems: quasi-static and dynamic. Damping is used in the solution of both classes of problems, but quasi-static problems require more damping. Two types of damping (mass-proportional and stiffness-proportional) are available in 3DEC. Mass-proportional damping applies a force which is proportional to absolute velocity and mass, but in the direction opposite to the velocity. Stiffness-proportional damping applies a force, which is proportional to the incremental stiffness matrix multiplied by relative velocities or strain rates, to contacts or stresses in zones. In 3DEC, either form of damping may be used separately or in combination. The use of both forms of damping in combination is termed Rayleigh damping (Bathe and Wilson 1976). For solution of quasi-static problems using finite difference schemes, mass-proportional damping is generally used (Otter et al. 1966). 3DEC allows use of an automatic “adaptive” viscous damping scheme developed by Cundall (1982) for solution of quasi-static problems. These schemes are discussed in 3DEC Theory and Background. For dynamic analyses, either mass-proportional or stiffness-proportional, or both (i.e., Rayleigh), forms of damping may be used, as described in the next section.

Rayleigh Damping

In performing dynamic analysis with any code, it is usually necessary to account for energy losses in the physical system (e.g., heat, hysteresis) which are not accounted for in the numerical algorithm. In general, very little damping is used for highly elastic systems, and more damping is used for geomechanical materials, especially soils.

In the continuum analysis of structures, proportional Rayleigh damping is typically used to damp the natural oscillation modes of the system. In dynamic finite-element analysis, a damping matrix, \(C\), is formed with components proportional to the mass (\(M\)) and stiffness (\(K\)) matrices,

\[C = \alpha\ M + \beta\ K\]


\(\alpha\) = the mass-proportional damping constant; and

\(\beta\) = the stiffness-proportional damping constant.

For a multiple degrees-of-freedom system, the critical damping ratio, \(\xi_i\), at any angular frequency of the system, \(\omega_i\), can be found from (Bathe and Wilson 1976)

\[\alpha + \beta\ \omega_i^2 = 2\ \omega_i\ \xi_i\]


\[\xi_i = {{1} \over {2}}\ \bigl( {{\alpha} \over {\omega_i}} + \beta\ \omega_i \bigr)\]

The critical damping ratio, \(\xi_i\), is also known as the fraction of critical damping for mode \(i\) with angular frequency, \(\omega_i\).

Figure 1 shows the variation of the normalized critical damping ratio with angular frequency, \(\omega_i\). Three curves are given: mass or stiffness components only, and the sum of both components. As shown, mass-proportional damping is dominant at lower angular frequency ranges, while stiffness-proportional damping dominates at higher angular frequencies. The curve representing the sum of both components reaches a minimum at

\[\xi_{\min} = (\alpha\ \beta)^{1/2}\]
\[\omega_{\min} = (\alpha/\beta)^{1/2}\]


\[\alpha = \xi_{\min}\ \omega_{\min}\]
\[\beta = \xi_{\min}\thinspace /\thinspace \omega_{\min}\]

The fundamental frequency is then defined as

\[f_{\min} = \omega_{\min}\thinspace /\thinspace 2 \pi\]

Figure 1: Variation of normalized critical damping ratio with angular frequency.

The required input parameters to specify Rayleigh damping in 3DEC are \(f_{\min}\) (input parameter freq) and \(\xi_{\min}\) (input parameter fcrit).

For the case shown in Figure 1, \(\omega_{\min}\) = 10 radians per second, and \(\xi_{\min}\) = 1. Note that the damping ratio is almost constant over at least a 3:1 frequency range (e.g., from 5 to 15). Since damping in geologic media is known to be predominately hysteretic (Gemant and Jackson 1937; Wegel and Walther 1935), and hence independent of frequency, \(\omega_{\min}\) is usually chosen to lie in the center of the range of frequencies present in the numerical simulation. In this way, hysteretic damping is simulated approximately.

From the preceding equations and Figure 1, it can be seen that at frequency \(\omega_{\min}\) (or \(f_{\min}\)), and only at that frequency, mass damping and stiffness damping each supply half of the total damping.

Example of Different Damping Techniques

In order to demonstrate how Rayleigh damping works in 3DEC, the results of the following four damping cases involving a block sitting on a fixed block with gravity suddenly applied can be compared:

  • undamped;
  • Rayleigh damping (both mass and stiffness damping);
  • mass damping only; and
  • stiffness damping only.

The data file for each of these cases follows.

Data File

; -------------------------------------------
;  Block under gravity - 3 critically damped cases
; -------------------------------------------
model new
model config dynamic
model large-strain on

block create brick 0 10 0 1 0 10
block cut j-set or 0 0 5

block fix range pos-z 0 5
model grav 0 0 -10

block prop dens 1000

block contact gen-sub
block contact jmodel assign mohr
block contact property stiffness-normal 5e7 stiffness-shear 5e7 
block contact material-table default property stiffness-normal 5e7 ...
                                              stiffness-shear 5e7 

block mechanical damp rayleigh 0 0

hist interval 1
block his vel-z pos 0 0 10
block his disp-z pos 0 0 10
model his dynamic time-total

model save "damp"

;no damping
model cycle 250
model title "Vertical displacement v time (undamped)
model save "damp-undamped"

;mass and stiffness
model restore "damp"
hist interval 20
block mech damp rayleigh 1 16
model cycle 4321

model title "Vertical displacement v time (MASS & STIFFNESS DAMPED; damp 1 16)"
model save "damp-mass-stiff"

;mass only
model rest "damp"
block mech damp rayleigh 2 16 mass
model cycle 250
model title "Vertical displacement v time (MASS DAMPED; damp 2 16 mass)"
model save "damp-mass"

;stiffness only
model rest "damp"
hist interval 50
block mech damp rayleigh 2 16 stiff
model cycle 8623
model title "Vertical displacement v time (STIFFNESS DAMPED; damp 2 16 stiff)"
model save "damp-stiff"
program return

In the first case, with no damping, a natural frequency of oscillation of approximately 16 cycles/sec can be observed (see Figure 2). The theoretical period of oscillation is given by

frequency = \({{1} \over {2 \pi}}\ \bigl( {{ka} \over {m}}\bigr)^{1/2}\) = 15.9 cycles/second


\(a\) = joint area (10 m2, in this case);

\(k\) = joint stiffness (5e7 Pa/m); and

\(m\) = mass of upper block (5e4 kg).

The problem is critically damped if: a fraction of critical damping = 1 is specified; the natural frequency of oscillation = 16 is specified; and both mass and stiffness damping are used.

The results in Figure 3 show that the problem is critically damped using these parameters. If only mass (Figure 4) or stiffness (Figure 5) damping is used, then the damping must be doubled (since each contributes 1/2 to the Rayleigh damping) to obtain critical damping. Figures 4 and 5 show, again, that the problem is critically damped using mass and stiffness damping with twice the damping specified (see the data file above).


Figure 2: Plot of vertical displacement versus time, for a single block contacting on a rigid base with gravity suddenly applied (no damping).


Figure 3: Plot of vertical displacement versus time, for a single block contacting on a rigid base with gravity suddenly applied (mass and stiffness damping).


Figure 4: Plot of vertical displacement versus time, for a single block contacting on a rigid base with gravity suddenly applied (mass damping only)


Figure 5: Plot of vertical displacement versus time, for a single block contacting on a rigid base with gravity suddenly applied (stiffness damping only).

Guidelines for Selecting Rayleigh Damping Parameters for Dynamic Analysis

What is normally attempted in a dynamic analysis is the reproduction of the frequency-independent damping of materials at the correct level (e.g., 2-5% for geological materials, and 2-10% for structural systems (Biggs 1964)). Rayleigh damping is frequency-dependent, but has a “flat” region that spans about a 3:1 frequency range, as shown in Figure 1. For any particular problem, a Fourier analysis of typical velocity records might produce a response such as that shown in Figure 6.


Figure 6: Plot of velocity spectrum versus frequency.

If the highest predominant frequency is three times greater than the lowest predominant frequency, then there is a 3:1 span or range that contains most of the dynamic energy in the spectrum. The idea in dynamic analysis is to adjust \(\omega_{\min}\) of the Rayleigh damping so that its 3:1 range coincides with the range of predominant frequencies in the problem. \(\xi_{\min}\) is adjusted to coincide with the correct physical damping ratio. The predominant frequencies are neither the input frequencies nor the natural modes of the system, but a combination of both. The idea is to try to get the right damping for the important frequencies in the problem.

For some problems involving large movements of blocks, it is improper to use any mass damping because the block motion might be artificially restricted. Examples of such problems include any problems involving free flow or fall of blocks under gravity, and impulsive loading of blocks due to explosions. In such cases it may be appropriate to use only stiffness-proportional damping.

The stiffness-proportional component of Rayleigh damping does affect the critical timestep for the explicit solution scheme in 3DEC. The controlling timestep, therefore, will be reduced as the stiffness damping component is increased (see Belytschko 1983). For problems involving free fall and “bounce” of blocks from a fixed base, the coefficient of restitution is required for accurate modeling. Onishi et al. (1985) provide a method for estimating stiffness damping parameters based on the coefficient of restitution.