Comparison of FLAC3D to SHAKE for a Layered, Nonlinear-Elastic Soil Deposit

Note

To view this project in FLAC3D, use the menu command Help ‣ Examples…. The main data files used are shown at the end of this example. The remaining data files can be found in the project.

In this section, we compare a simulation of wave propagation through a geological profile using FLAC3D with hysteretic damping to SHAKE-91. The FLAC3D and SHAKE models simulate the following problem: a horizontally layered soil deposit overlying a rigid bedrock is subjected to a horizontal acceleration base motion. The soil deposit is 150 feet (45.72 m) deep and contains 10 different soil types. The soils are treated as nonlinear elastic materials by assuming that shear moduli and damping are strain-dependent. Table 1 summarizes the properties and locations of the soil layers in the deposit. The dynamic characteristics of these soils are governed by two sets of modulus reduction factor (\(G/G_{\rm max}\)) and damping ratio (\(\lambda\)) curves: the first set (set 1) for clay, and the second set (set 2) for sand. [1] These curves are shown in Figure 1 through 4, and are denoted by “SHAKE91” in the legend.

Table 1: Soil Deposit Profile and Properties

Soil	1	2	3	4	5	6	7	8	9	10
Shear Modulus (MPa)	186	150	168	186	225	327	379	435	495	627
Density (Kg/m3)	2000	2000	2000	2000	2000	2082	2082	2082	2082	2082
Dynamic Property (set)	2	2	2	1	1	2	2	2	2	2
Location (feet)	1-5	5-20	20-30	30-50	50-70	70-90	90-110	110-130	130-140	140-150

Figure 1: Modulus reduction curve for FLAC3D default model versus dynamic property set 2 (Seed & Idriss data — sand — “upper range”).

Figure 2: Damping ratio curve for FLAC3D default model versus dynamic property set 2 (Seed & Idriss data — sand — “upper range”).

Figure 3: Modulus reduction curve for FLAC3D default model versus dynamic property set 1 (Seed & Sun data — clay — “upper range”).

Figure 4: Damping ratio curve for FLAC3D default model versus dynamic property set 1 (Seed & Sun data — clay — “upper range”).

This example is derived from the data file provided with the downloaded SHAKE-91 code (see https://nisee.berkeley.edu/elibrary/software.html). The *SHAKE-91* data file is modified such that the deepest layer is given a large wave speed to correspond to rigid bedrock. Note that to simulate a rigid base, we increase the stiffness (shear modulus) of the bedrock from 3.33 GPa to 2080 GPa in the SHAKE-91 data file for the purpose of comparison.

The base acceleration input is a set of seismic data recorded in the Loma Prieta Earthquake, which is also downloadable (see https://nisee.berkeley.edu/elibrary/getpkg?id=SHAKE91) with the SHAKE-91 code (i.e., “DIAM.ACC”). The input accelerogram is shown in Figure 5. The record has a peak acceleration of approximately 0.11 g, duration of 40 seconds, and dominant frequency of 2.47 Hz.

Figure 5: Input accelerogram.

SHAKE model — The free-field response of the layered soil deposit to the base-motion acceleration is calculated using SHAKE-91. An iterative solution is performed until the secant shear moduli and the damping ratios throughout the column converge to within a given tolerance of the values corresponding to 50% of the maximum shear strain. Calculations are made for different input accelerations, with the record scaled to a maximum base acceleration. [2]

FLAC3D model with hysteretic damping — A 1D column FLAC3D model composed of 30 cubic zones, each with a depth of 5 feet (1.524 m), is assigned the soil properties listed in Table 1. Model conditions are prescribed to simulate the free-field motion of a layered soil deposit with a rigid base. The horizontal acceleration record is applied as an input boundary condition at the base of the model. FLAC3D takes the scaled base motion from SHAKE-91. For example, note that the input table “Diam-flac-0001.acc” is scaled with a maximum value of 0.0001 g. Vertical movement is prevented in the model.

Hysteretic damping is added to this model by fitting a secant-modulus function curve to the (\(G/G_{\rm max}\)) versus strain data, as shown in Figures 1 and 3. Different fitting functions are available in FLAC3D. In this example, the default (two-parameter) function is used. The parameters of the secant-modulus functions are determined by best-fitting (e.g., least squares regression analysis) the laboratory modulus-reduction curves, assuming that the modulus reduction factor versus shear strain follows the default relation. The best-fit values for dynamic property set 1 — sand are listed here, and for dynamic property set 2 — clay here CS: this is not cool, but understandable. Preceding x-refs want better labeling. A comparison between the secant-modulus function curves and the laboratory curves is shown in Figure 1 and 3.

Damping ratio curves that correspond to the secant-modulus function curves are then calculated. [3] Both shear modulus reduction factor (\(G/G_{\rm max}\)) versus strain and damping ratio (\(\lambda\)) versus strain should be monitored, and parameters adjusted, if necessary, to provide a reasonable fit to both curves over a specified strain range. As shown in Figures 1 through 4, fairly reasonable fits for both modulus-reduction factor and damping-ratio curves are obtained with the default damping function and the parameters listed in Numerical Fits to Seed and Idris over a strain range of 0.0001% to 0.1% for both the sand and the clay.

The fitting curves of the FLAC3D default model (Figures 1 through 4) show that, for small strain (corresponding to small input acceleration), the behavior is approximately linear (i.e., both shear modulus and damping ratio are constants for both FLAC3D and SHAKE-91). Both codes are thus expected to give similar results in this circumstance. Here we compare accelerograms and response spectra at the top of the model for very low input acceleration. Figure 6 shows the horizontal acceleration at the top of the model (gridpoint 121 in FLAC3D and sub-layer 1 in SHAKE-91) as a function of time with maximum input acceleration amplitude of 0.0001 g. Both records are very similar; maximum accelerations calculated by FLAC3D and SHAKE-91 have less than 0.1% difference.

Figure 6: Accelerograms at the top of the model with small input.

Figures 7 and 8 provide pseudo-acceleration and pseudo-velocity spectra, calibrated in SHAKE-91 and FLAC3D when the maximum input acceleration amplitude is small (0.0001 g). In SHAKE-91, the response spectra are calculated using option 9. In FLAC3D, response spectra are computed using a FISH function, spectra, implemented in file “spectra.dat”. Here the damping ratio, minimum period and maximum period of interest are 5%, 0.01 and 10, respectively. From these plots it can clearly be seen that the FLAC3D and SHAKE-91 results correspond quite closely for small-amplitude input motion.

Figure 7: Pseudo-acceleration spectrum at the top of the model.

Figure 8: Pseudo-velocity spectrum at the top of the model.

For validation purposes, acceleration amplifications at the top of the model, defined as the ratio of the maximum acceleration value at the top to the maximum base input acceleration value, are compared between FLAC3D and SHAKE-91. We measure the ratios when applying scaled maximum input acceleration value in “DIAM.ACC” of 0.0001 g, 0.0005 g, 0.001 g, 0.005 g, 0.01 g, 0.05 g, 0.1 g, 0.5 g, and 1 g, where g is the acceleration of gravity (9.81 m/s²). FLAC3D simply takes the scaled base motion from SHAKE-91. According to the calculation results shown in Figure 9, it is confirmed that FLAC3D and SHAKE-91 predict similar acceleration amplification.

Figure 9: Acceleration amplifications comparison at the top of the model.

References

Hardin, B. O., and V. P. Drnevich. “Shear Modulus and Damping in Soils: I. Measurement and Parameter Effects, II. Design Equations and Curves,” Technical Reports UKY 27-70-CE 2 and 3, College of Engineering, University of Kentucky, Lexington, Kentucky. [These reports were later published in the Journal of Soil Mechanics and Foundation Division, ASCE, Vol. 98, No. 6, pp. 603-624 and No. 7, pp. 667-691, in June and July 1972.]

Idriss, I. M., and Joseph I. Sun. User’s Manual for SHAKE-91. University of California, Davis, Center for Geotechnical Modeling, Department of Civil & Environmental Engineering (November 1992).

Seed, H. Bolton, and I. M. Idriss. “Soil Moduli and Damping Factors for Dynamic Response Analysis,” Earthquake Engineering Research Center, University of California, Berkeley, Report No. UCB/EERC-70/10, p. 48 (December 1970).

Seed, H. Bolton, et al. “Moduli and Damping Factors for Dynamic Analyses of Cohesionless Soils,” Journal of the Geotechnical Engineering Division, ASCE, Vol. 112, No. GT11, pp. 1016-1032 (November 1986).

Sun, J. I., R. Golesorkhi and H. Bolton Seed. “Dynamic Moduli and Damping Ratios for Cohesive Soils,” Earthquake Engineering Research Center, University of California, Berkeley, Report No. UCB/EERC-88/15, p. 42 (1988).

Vucetic, M., and R. Dobry. “Effect of Soil Plasticity on Cyclic Response,” in J. Geotech. Eng. Div., ASCE, Vol. 111, No. 1, pp. 89-107 (January 1991).

Data Files

F3ShakeNonLinearElasticCase.dat

;---------------------------------------------------------------------
;              Comparison of Hysteretic Damping with Shake91
;              This file uses SHAKE91-comparison-flac3d.f3dat
;              to create the acceleration
;              amplification ratio at the top of layer 1. It changes
;              the acceleration amplifications for different runs.
;              It also creates the graphs and checks for errors.
;---------------------------------------------------------------------
model new
fish automatic-create off
model title 'Comparison of Hysteretic Damping with Shake91'
model configure dynamic

;----------------------------
; Load applied acceleration
;----------------------------
table 'iacc' import 'Diam-flac-0001.acc'

; Run 9 different cases of acceleration multiplier,
;   storing the acceleration amplification result for each
[global all_mult = map(1,1,2,5,3,10,4,50,5,100,6,500,7,1000,8,5000,9,10000)]
             ; Map of acceleration multiplier cases, 9 cases from 1 to 10000
fish define run_all_cases
    global acc_case
    loop foreach acc_case map.keys(all_mult) ; All case numbers
        global acc_mult = all_mult(acc_case) ; Acceleration multiplier
                                             ; for that case 
        command              ; program call the data file that runs it
            program call 'comparison'
        end_command
    end_loop
end
[run_all_cases]

model save 'final'

comparison.dat

model large-strain off
;---------------------------------------------------------------------
; Comparison of Hysteretic Damping with Shake91
;---------------------------------------------------------------------

; First, reset model without clearing everything
zone delete ; Remove all zones
history delete ; Remove all histories
model dynamic time-total 0 ; Reset dynamic time to 0

;-------------------------------------
;Grid generation and model properties
;-------------------------------------
zone create brick size 1 1 30 
zone cmodel assign elastic
zone property bulk 300e6 shear 186e6 density 2000 ...
              range position-z 29 30   
zone property bulk 200e6 shear 150e6 density 2000 ...
              range position-z 26 39 ;  5- 20 ft
zone property bulk 200e6 shear 168e6 density 2000 ...
              range position-z 24 26 ; 20- 30 ft
zone property bulk 270e6 shear 186e6 density 2000 ...
              range position-z 20 24 ; 30- 50 ft
zone property bulk 350e6 shear 225e6 density 2000 ...
              range position-z 16 20 ; 50- 70 ft
zone property bulk 480e6 shear 327e6 density 2082 ...
              range position-z 12 16 ; 70- 90 ft
zone property bulk 550e6 shear 379e6 density 2082 ...
              range position-z  8 12 ; 90-110 ft
zone property bulk 600e6 shear 435e6 density 2082 ...
              range position-z  4  8 ;110-130 ft
zone property bulk 750e6 shear 495e6 density 2082 ...
              range position-z  2  4 ;130-140 ft
zone property bulk 900e6 shear 627e6 density 2082 ...
              range position-z  0  2 ;140-150 ft

;------------
; Histories
;------------
model history name 'time' dynamic time-total ; 2
zone history name 'atop' acceleration-x position (0,0,30) ;top accel 0'

;--------------------
;Boundary Conditions
;--------------------
zone gridpoint fix velocity-y
zone gridpoint fix velocity-z
zone face apply acceleration-x [acc_mult*9.81] table 'iacc' ...
                time dynamic range position-z = 0

;--------------------
; Hysteretic damping
;--------------------
zone dynamic damping hysteretic default -3.325 0.823 range position-z  0 16
zone dynamic damping hysteretic default -3.156 1.904 range position-z 16 24
zone dynamic damping hysteretic default -3.325 0.823 range position-z 24 30

;----------------------------------------------
; Increase size by 1.524 to match shake
;----------------------------------------------
zone gridpoint initialize position (1.524,1.524,1.524) multiply

;----------------------------------------------
; Solve to 40.48s, taking histories every 0.02s.
;----------------------------------------------
model dynamic timestep fix 0.0002
history interval 100
model solve time-total 40.48

;top accel hist to table in order to generate response spectra
history export 'atop' vs 'time' table ['tacc'+string(acc_case)]

NonlinearElasticCase-SHAKE.dat

option 1 -- dynamic soil properties -- (max is thirteen):
    1
    3
   11    #1 modulus for clay (Seed & Sun 1989) upper range
0.0001   0.0003   0.001    0.003    0.01     0.03     0.1      0.3
1.           3.           10.
1.000    1.000    1.000    0.981    0.941    0.847    0.656    0.438
0.238    0.144    0.110
   11      damping for clay (Idriss 1990) --
0.0001   0.0003   0.001    0.003    0.01     0.03     0.1      0.3
1.           3.16       10.
0.24     0.42     0.8      1.4      2.8       5.1      9.8     15.5
21.       25.      28.
   11    #2 modulus for sand (Seed & Idriss 1970) -- upper range
0.0001   0.0003   0.001    0.003    0.01     0.03     0.1      0.3
1.           3.           10.
1.000    1.000    0.990    0.960    0.850    0.640    0.370    0.180
0.080    0.050    0.035
   11      damping for sand (Idriss 1990) -- (about LRng from SI 1970)
0.0001   0.0003   0.001    0.003    0.01     0.03     0.1      0.3
1.           3.          10.
0.24     0.42     0.8      1.4      2.8       5.1      9.8     15.5
21.       25.       28.
    8     #3 ATTENUATION OF  ROCK  AVERAGE
 .0001    0.0003   0.001    0.003    0.01     0.03      0.1    1.0
 1.000    1.000    0.9875   0.9525   0.900    0.810    0.725    0.550
    5     DAMPING IN ROCK
 .0001    0.001    0.01       0.1       1.
  0.4       0.8        1.5         3.0       4.6
    3   1    2    3
Option 2 -- Soil Profile
    2
    1   17      Example -- 150-ft layer; input:Diam [] .1g
    1    2           5.00                .050      .125     1000.
    2    2           5.00                .050      .125      900.
    3    2          10.00                .050      .125      900.
    4    2          10.00                .050      .125      950.
    5    1          10.00                .050      .125     1000.
    6    1          10.00                .050      .125     1000.
    7    1          10.00                .050      .125     1100.
    8    1          10.00                .050      .125     1100.
    9    2          10.00                .050      .130     1300.
   10    2          10.00                .050      .130     1300.
   11    2          10.00                .050      .130     1400.
   12    2          10.00                .050      .130     1400.
   13    2          10.00                .050      .130     1500.
   14    2          10.00                .050      .130     1500.
   15    2          10.00                .050      .130     1600.
   16    2          10.00                .050      .130     1800.
   17    3                               .010      .140   100000.
Option 3 -- input motion:
    3
 1900 4096   .02     diam.acc                     (8f10.6)
           .0001      25.        3    8
Option 4 -- sublayer for input motion within (1) or outcropping (0):
    4
   17    0
Option 5 -- number of iterations & ratio of avg strain to max strain
    5
    0    8    0.50
Option 6 -- sublayers for which accn time histories are computed & saved:
    6
    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
    0    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    0    0    0    0    1    0    0    0    1    0    0    0    1    0
Option 6 -- sublayers for which accn time histories are computed & saved:
    6
   16   17   17
    1    1     0
    0    1     0
option 7 -- sublayer for which shear stress or strain are computed & saved:
    7
    4    1    1    0 1800           -- stress in level 4
    4    0    1    0 1800           -- strain in level 4
option 7 -- sublayer for which shear stress or strain are computed & saved:
    7
    8    1    1    0 1800           -- stress in level 8
    8    0    1    0 1800           -- strain in level 8
option 9 -- compute & save response spectrum:
    9
    1    0
    1    0     981.0
   0.05
option 10 -- compute & save amplification spectrum:
   10
   17    0    1    0     0.125      -- surface/rock outcrop
execution will stop when program encounters 0
    0

Endnotes

[1]

More variations of \(G/G_{\rm max}\) and \(\lambda\) for soils are available in the literature mentioned in the SHAKE-91 manual (Idriss and Sun 1992) and duplicated here for easy reference (e.g., Hardin and Drnevich 1972, Seed and Idriss 1970, Seed et al. 1986, Sun et al. 1988, and Vucetic and Dobry 1991).

[2]

The scaling is accomplished through the built-in functionality in SHAKE-91, using either of two scaling parameters in OPTION 3.

[3]

The data file in the project Several Cyclic Strain Levels can be used to perform a series of simple shear tests on a single-zone FLAC3D model in order to develop damping ratio curves that correspond to the secant-modulus functions.