Damping
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: quasistatic and dynamic. Damping is used in the solution of both classes of problems, but quasistatic problems require more damping. Two types of damping (massproportional and stiffnessproportional) are available in 3DEC. Massproportional damping applies a force which is proportional to absolute velocity and mass, but in the direction opposite to the velocity. Stiffnessproportional 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 quasistatic problems using finite difference schemes, massproportional 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 quasistatic problems. These schemes are discussed in 3DEC Theory and Background. For dynamic analyses, either massproportional or stiffnessproportional, 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 finiteelement analysis, a damping matrix, \(C\), is formed with components proportional to the mass (\(M\)) and stiffness (\(K\)) matrices,
where:
\(\alpha\) = the massproportional damping constant; and
\(\beta\) = the stiffnessproportional damping constant.
For a multiple degreesoffreedom system, the critical damping ratio, \(\xi_i\), at any angular frequency of the system, \(\omega_i\), can be found from (Bathe and Wilson 1976)
or
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, massproportional damping is dominant at lower angular frequency ranges, while stiffnessproportional damping dominates at higher angular frequencies. The curve representing the sum of both components reaches a minimum at
or
The fundamental frequency is then defined as
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 largestrain on
block create brick 0 10 0 1 0 10
block cut jset or 0 0 5
block fix range posz 0 5
model grav 0 0 10
block prop dens 1000
block contact gensub
block contact jmodel assign mohr
block contact property stiffnessnormal 5e7 stiffnessshear 5e7
block contact materialtable default property stiffnessnormal 5e7 ...
stiffnessshear 5e7
block mechanical damp rayleigh 0 0
hist interval 1
block his velz pos 0 0 10
block his dispz pos 0 0 10
model his dynamic timetotal
model save "damp"
;no damping
model cycle 250
model title "Vertical displacement v time (undamped)
model save "dampundamped"
;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 "dampmassstiff"
;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 "dampmass"
;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 "dampstiff"
;
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
where:
\(a\) = joint area (10 m^{2}, 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).
Guidelines for Selecting Rayleigh Damping Parameters for Dynamic Analysis
What is normally attempted in a dynamic analysis is the reproduction of the frequencyindependent damping of materials at the correct level (e.g., 25% for geological materials, and 210% for structural systems (Biggs 1964)). Rayleigh damping is frequencydependent, 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.
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 stiffnessproportional damping.
The stiffnessproportional 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.
Maxwell Damping
Maxwell damping provides an alternative to Rayleigh damping for timedomain seismic deformation analyses. The numerical calculation time using Maxwell damping can be a small fraction of that required for Rayleigh damping. Maxwell damping makes practical large 3D site response and soilstructure interaction analyses.
Using Rayleigh damping, a target damping ratio can obtain at two specific frequencies; the resulting damping is slightly smaller than the target damping ratio when the frequency is within the range bounded by these two specific frequencies, while damping becomes dramatically larger outside of this range. While Rayleigh damping is helpful to reduce artificial numerical high frequency noise, it can also over damp high frequency content because the damping ratio of the Rayleigh scheme can be more than 20% at high frequencies, even if it is low in the frequency range of interest. As pointed out by Hall (2006), the massproportional component of Rayleigh damping is essentially a nonphysical mechanism and the stiffnessproportional component for softening materials can result in unrealistically large damping forces under certain conditions. Hall (2006), thus, suggested eliminating the massproportional damping contribution and bounding the stiffnessproportional damping contribution for Rayleigh damping. Another serious drawback of Rayleigh damping is the drastically reduced numerical time step. The stiffness damping component of Rayleigh damping leads to increasingly stiff response at higher frequencies, so that a smaller numerical time step is needed for stable numerical integration (Belytschko 1983). For large 3D models with many degrees of freedom, the requirement of significantly longer machine calculation time makes the applications of the stiffnessproportional Rayleigh damping less practical.
Bielak et al. (2011) and others have shown that an improved form of damping for timedomain analyses can be developed by replacing the dashpot of stiffness damping, with one or more Maxwell components (a spring in series with a dashpot). While each Maxwell component produces frequency dependent damping, peaking at a particular frequency, a number of Maxwell components acting in parallel can be used to produce relatively constant damping over a desired frequency range. Dawson and Cheng (2021) implemented a practical Maxwell damping scheme with three components to provide small strain damping for timedomain seismic deformation analyses. It has been showed (Dawson and Cheng, 2021; Cheng and Dawson, 2022) that this scheme with three Maxwell components can produce relatively constant damping in the frequency range of interest for most seismic deformation analyses. Damping drops off gradually outside of this range.
Formulation
Standard Linear Solid
The basic building block of the Maxwell damping scheme is what is known as the Standard Linear Solid (SLS) as shown in Figure 4. This consists of a primary spring acting in parallel with a single Maxwell component – the Maxwell component itself being a spring in series with a viscous dashpot. Note that the stiffness of the Maxwell component spring is specified here as a fraction of the stiffness K of the primary spring.
We first derive an expression for the complex (frequency domain) stiffness of the Maxwell component only. The differential equation describing the stressstrain response of the Maxwell component is:
where \(\alpha K\) is the spring stiffness and η is the viscosity. This differential equation can be solved by Fourier transform, or equivalently, by assuming harmonic motions of the form:
where \(\omega\) is the circular frequency and \(\varepsilon_0\) and \(\sigma_0\) are complex amplitudes, accounting for phase as well as magnitude. The complex stiffness, or modulus, of the Maxwell component is then
where
is the relaxation time, the ratio of viscosity to elastic stiffness. Dividing relation (12) by the stiffness K of the primary spring and using (13), results in the normalized, nondimensional complex modulus for the Maxwell component
Next, consider the complex modulus of the Standard Linear Solid (Figure 4). Since the primary spring and the Maxwell component act in parallel, the stiffness of this system is simply the sum of the individual stiffnesses. The normalized complex modulus is then
or
The damping ratio for a viscoelastic material can be computed from the ratio of the imaginary and real parts of the complex stiffness (Bland, 1960), so that for the Standard Linear Solid, the damping ratio is:
It can be shown that maximum damping occurs at a circular frequency of
with a maximum damping value of
The damping versus frequency response of the SLS is illustrated in Figure 5 which shows the damping ratio for an SLS with 5% damping with a peak frequency of 3.5 Hz. The damping of the SLS has a narrow peak and broad tails.
Full Maxwell Damping Scheme
Next, we consider the full Maxwell damping scheme with three Maxwell components placed in parallel with the primary spring, as shown in Figure 6. Here the primary spring represents the elastic stiffness of the base constitutive model to which damping is to be added.
According to Equation (14), the normalized complex stiffness of an individual Maxwell component with index k is
where
The normalized complex stiffness of the system as a whole, with a primary spring and three Maxwell components acting in parallel is then
and the damping ratio of the system is
With careful selection of the properties for the three Maxwell components, damping of the system can be made almost constant over a wide frequency range, as shown in Figure 7. Using three Maxwell components, relatively constant damping can be obtained between 0.5 Hz and 25 Hz, the relevant frequency range for most seismic deformation analyses. Outside this range, damping drops off gradually.
It can be seen in Figure 7 that damping of the system as a whole is approximately the superposition of the damping of the individual Maxwell components. This rough approximation provides a helpful starting point for selecting frequencies and damping values for the components to match a target damping value. Damping of the components can then be adjusted to obtain a closer fit. Figure 7 illustrates that this superposition is only approximate, as the peak damping for the 25 Hz component exceeds the damping of the combined system. In fact, there is interaction between the three components, as would be expected from Equation (16).
The lower damping at high frequencies shown in Figure 7 might be seen as a drawback of this technique, as timedomain numerical analyses can be subject to highfrequency numerical noise. However, high frequency vibrations, by definition, experience more cycles per unit time, with each cycle dissipating energy according to the damping ratio. Thus, high frequency noise can be damped out rapidly (in time) with even a very small damping ratio.
User Inputs and Maxwell Component Properties
For specifying input properties for the damping scheme, it is convenient to think of the system as composed of three SLS, each consisting of a Maxwell component in parallel with the primary spring. Damping of the complete system is then very roughly the superposition of the damping of each of these SLS system. Convenient input parameters for the damping scheme are the center frequencies \(f_k\) (in hertz) of each SLS and the damping ratios \(\xi_k\) of each SLS. The \(\alpha_k\) constants can then be computed from the damping ratio with Equation (15), using the quadratic formula
The relaxation time for the individual Maxwell components is then
and the dashpot viscosities are
If pairs of \((f_k,\xi_k)\) are known, the closedform solution for the target damping ratio can be straightforwardly obtained. However, calibration requires working in the opposite direction  finding pairs of \((f_k,\xi_k)\) that produce a suitable target damping ratio. A practical guideline is suggested here:
Select a target damping ratio, e.g., 5%.
Estimate a target frequency range over which the approximate target damping ratio is desired. The lower and upper bounds of this range is \(f_1\) and \(f_3\). Then \(f_2\) can be a value so that \(f_2 \approx \sqrt{f_1 f_3}\). For example, if we wish a target damping ratio of 0.05 and constant damping in the range of 0.5 and 25 Hz, then \(f_1\), \(f_2\), \(f_3\) can be 0.25, 3.5, 25 Hz, respectively. Generally, the frequency range of interest should cover the first primary natural frequency of the soil profile, the third primary natural frequency of the soil profile, and the dominate frequency of the input motion.
Using a plotting tool (e.g., Excel), adjust \((\xi_1,\xi_2,\xi_3)\) to approximate the target damping ratio over the target range. Experience has shown that with only a few trails, a satisfactory target damping ratio can be obtained.
Table 1 lists Maxwell component properties for various damping values, all fits using the same set of center frequencies. The resulting damping versus frequency curves are shown in Figure 8. Note that other values can be assigned to the center frequencies \(f_k\) in order to change the position or shape of the damping versus frequency curve.
Maxwell Component Properties 


Damping 
\(f_1=0.5\) Hz 
\(f_2=3.5\) Hz 
\(f_3=25\) Hz 
\(\xi_1\) 
\(\xi_2\) 
\(\xi_2\) 

0.5% 
0.0039 
0.0029 
0.0041 
1.0% 
0.0080 
0.0060 
0.0085 
2.0% 
0.0160 
0.0120 
0.0180 
3.0% 
0.0235 
0.0190 
0.0284 
4.0% 
0.0311 
0.0260 
0.0398 
5.0% 
0.0385 
0.0335 
0.0520 
To provide low strain damping for linear elastic, as well as elastoplastic constitutive models. Maxwell damping has also been implemented for structural elements such as beams, shells, cables and piles. The incremental equation describing the stressstrain response of a Maxwell component is:
where \(\Delta \varepsilon\) and \(\Delta t\) are known variables based on the FLAC3D dynamic implementation scheme. \(\sigma^n\), which is the stress of the Maxwell component calculated in the previous step, n, is also known. In order to calculate the stress of the Maxwell components in the current step, n+1, explicit or implicit algorithms can be employed.
Explicit Algorithm
If \(\sigma=\sigma^n\) in Equation (28), it can be shown that
The time step required for numerical stability is \(\Delta t < 2 \eta / \alpha K\).
Implicit Algorithm
If \(\sigma=\sigma^{n+1}\) in Equation (28), it can be shown that
This algorithm is unconditionally stable. Even if \(\Delta t \to +\infty\), \(\sigma^{n+1}\) is still convergent to a bounded value \(\eta \dot{\varepsilon}\). Thus, the time step for integration of the Maxwell component does not have to be smaller than the time step used for the element dynamic integration. The implicit algorithm, Equation (30), is used for the FLAC3D implementation (Dawson and Cheng, 2021). For solid elements, the actual implementation considers as a stress tensor, as a strain tensor and K as a constitutive stiffness matrix. For structural elements K is the element stiffness matrix.
Verification Problem
The following example demonstrates how Maxwell damping works in 3DEC.
Was this helpful? ...  Itasca Software © 2022, Itasca  Updated: Mar 09, 2023 