Recommended Steps for Seismic Analyses

A seismic analysis of an earth structure can be quite complicated, especially if there is potential for dynamic liquefaction of the soil materials. The following steps are recommended to perform a seismic analysis using a nonlinear effective-stress material model. In this procedure, the nonlinear stress/strain behavior is represented by a “simple” elastic/plastic (Mohr-Coulomb) model, with formulations added to simulate liquefaction and hysteretic damping of the soils during cyclic loading. In general, these steps for performing a seismic analysis can also be applied for a more comprehensive liquefaction model.

1. Determine representative static and dynamic material characteristics — The effective-stress analysis approach in FLAC3D is performed as a fully coupled mechanical-fluid flow analysis, which requires the selection of model configure fluid, and assignment of drained stiffness and strength material properties, and the unsaturated (“dry”) density (or unit weight).

   In addition, the fully coupled effective-stress calculation requires that the water bulk modulus be specified explicitly. The water bulk modulus must be selected carefully for this type of analysis. The behavior of the model depends on the stiffness ratio, \(R_k\), as defined in Solving Flow-Only and Coupled-Flow Problems and repeated here:

   (1)\[ R_k = {{K_w/n} \over {K + 4 G/3}}\]

   where \(K_w\) is the water bulk modulus, \(n\) is the porosity, and \(K\) and \(G\) are the bulk and shear moduli of the unsaturated soil. There is a temptation to decrease the water bulk modulus arbitrarily in order to increase the timestep (and reduce the simulation time). However, there are two cases to consider when selecting \(K_w\):

   1. If the calculated value of \(R_k\) based upon the given soil properties and an assumed value of \(K_w\) = 4.18×10⁷ psf, or 2 GPa (pure water), is greater than 20, then the water bulk modulus can be reduced such that \(R_k\) = 20 without affecting the results significantly (and reducing the simulation time).
   2. If the calculated value for \(R_k\), using the actual water bulk modulus, is less than 20, then that value of water bulk modulus should be used.

   The water bulk modulus can be made to depend on the local elastic moduli, and can be specified differently for different soil units, provided that the two conditions listed above are satisfied.

   The material response to dynamic cyclic loading is typically quantified by a shear modulus degradation curve and a damping ratio curve. Representative curves should be selected for each material in order to provide an accurate representation of wave attenuation and energy dissipation during dynamic loading. In addition, it may be necessary to make the curves depth-dependent to make the energy dissipation more realistic (e.g., see Darendeli 2001).

   As discussed in Comprehensive Liquefaction Constitutive Models, a comprehensive material model to represent liquefaction behavior can be quite complex. The simplest forms of liquefaction models, as discussed in Simple Model Formulations, are based on a volume-change modification to the Mohr-Coulomb model. In this case, the only additional property is the blow count associated with the material. For example, see the Finn/Byrne formula defined in Simple Model Formulations.

2. Evaluate seismic motion characteristics — A design earthquake ground motion is often provided as an acceleration record. It is typically an outcrop motion that is often recorded at a rock outcrop. It is important to know the location of the recorded motion because the motion may need to be modified for input to the FLAC3D model. The deconvolution analysis that is performed to obtain the appropriate input motion is discussed below in Step 6.

   The characteristics of the seismic motion should always be checked, because these characteristics can influence the model conditions. As discussed in Wave Transmission, the frequency content of the input motion affects the selection of mesh size for accurate wave propagation. The criteria \(\Delta l \le {\lambda / 10}\) should be used to choose the appropriate maximum zone size for an accurate representation. If the highest frequencies associated with the input motion necessitate an extremely fine mesh (and a correspondingly small timestep), it may be possible to remove the high frequency components to permit a coarser mesh. If most of the power for the input is contained in the lower frequency components (say 80% to 90%), then the history can be filtered to remove the higher frequencies without significantly affecting the results. The FISH function “fft.fis” (in “\datafile\FISH\library") can be used to evaluate the frequency content of the wave. The FISH function “filter.fis” (also in “\datafile\FISH\library") can be used to perform the filtering.

   For assessing frequency content, it is suggested that the input motion be evaluated in the form of a particle velocity (not acceleration) history. (An acceleration record can be converted into a velocity record using “int.fis” described in “\datafile\FISH\library".) For a plane wave propagating through a continuous medium, it is easy to show (e.g., see textbooks on wave propagation) that stress \(\sigma = C \rho v\) (where \(C\) is the wave speed, \(\rho\) is the density and \(v\) is the particle velocity). Therefore, when considering the potential for yield or failure (which is determined by the level of induced stress), it is the particle velocity that is most relevant when evaluating the dominant frequencies, because the induced stress is directly proportional to velocity, not acceleration.

   The input record should also be checked for baseline drift (i.e., continuing residual displacement after the motion has finished — see Baseline Correction). The FISH function “int.fis,” is used to integrate the velocity record to produce the displacement waveform related to the input acceleration. If needed, a baseline correction can be performed by adding a low frequency sine wave to the velocity record; the sine wave parameters are adjusted so that the final displacement is zero. An example is given in the example of Dam Foundation.

3. Estimate material damping parameters to represent inelastic cyclic behavior — When using a simple elastic-perfectly plastic material model (such as the Mohr-Coulomb model) in a seismic analysis, it is necessary to incorporate additional material damping to account for cyclic energy dissipation during the elastic part of the response, as discussed here. A difficult aspect of this is how to determine the appropriate material damping input. Generally, two different schemes are used, either Rayleigh damping or hysteretic damping).

   The equivalent-linear method, as applied in program the SHAKE, for example, is one way to estimate material damping input for the Rayleigh damping scheme and the hysteretic damping scheme. An equivalent-linear analysis is performed on a soil column, representative of the site conditions, using the shear wave speeds and densities for the different soil layers in the column, the modulus reduction and damping ratio curves selected as representative of the materials, and the target earthquake design motion for the site. Elastic strain-compatible values are then determined for the shear-modulus reduction factors and damping ratios. Average modulus-reduction factors and damping ratios can be estimated for each soil layer; these are the input parameters for the model with Rayleigh damping.

   The anticipated range of cyclic shear-strain magnitudes for the given site conditions is needed to specify a best-fit range for the modulus reduction and damping ratio curves used with hysteretic damping (as discussed here). This range can be estimated from the range of equivalent uniform cyclic strains provided from a SHAKE simulation.

   Note that these estimates are derived from the equivalent-linear analysis that assumes a low level of nonlinearity. Some adjustment in parameters may be required, especially if the actual model exhibits strong nonlinearity. This is discussed further in Step 7.

4. Create appropriate model grid for accurate wave propagation — The characteristics of the input motion are used to help select the appropriate mesh size and adjust the input wave record in order to provide an accurate solution in the seismic analysis. (See Wave Transmission for further information on the relation between wave propagation characteristics and mesh size.)

   In most seismic simulations, especially if excessive motion and strains are anticipated, a large-strain simulation should be performed in order to provide a more accurate deformation solution. When significant deformation and distortion of the grid is anticipated, it is important to minimize the number of tetrahedral-shaped zones in the mesh and, in particular, along slope faces. Tetrahedral zones along slope faces are prone to become badly distorted during large-strain calculations, because tetrahedral zones do not contain overlaid sets of subzones. Hexahedral-shaped zones contain two overlaid sets of subzones, which provide a more accurate calculation for materials undergoing plastic yield. (See Mixed Discretization for a Hexahedral Grid.) [1]

5. Calculate static equilibrium state for site — It is important to model the construction sequence of the earth structure as closely as possible in order to provide a reasonable representation of the initial, static shear stresses in the structure. This is important, particularly in a liquefaction analysis, because the initial static shear stresses can affect the triggering of liquefaction.

   Simple analyses typically assume that the initial shear modulus (\(G_{max}\)) is uniform throughout each material unit. It may be more appropriate to vary \(G_{max}\) for soils as a function of the in-situ effective stress (e.g., see Kramer 1996). For example, it may be considered that the maximum shear modulus varies as a function of effective stress as defined by the Seed et al. (1986) expression,

   (2)\[G_{max} = 21.7 \times P_a \times K_{2,max} \left( {\sigma_m^{\prime} \over P_a} \right)^{0.5}\]

   where \(P_a\) is the atmospheric pressure, \(K_{2,max}\) is a constant determined from the relative density, and \(\sigma_m^{\prime}\) is the effective mean stress. This initial modulus variation can be implemented via FISH during the static loading stage.

   It is important, if hysteretic damping is applied in the seismic simulation, to check the initial shear stress at the static equilibrium state. The hysteretic damping formulation is assumed to initiate hysteresis from an initial value of zero shear strain. If the initial shear stresses are high, then the shear stress and shear strain state for hysteretic damping may not be compatible. If shear stresses are high at the static equilibrium state of the model, then in order to ensure that shear stresses and strains are consistent during the dynamic phase, hysteretic damping should be invoked before the model is brought to the initial equilibrium state. See here for a recommended procedure to incorporate hysteretic damping during the static solution stage.

6. Match dynamic loading to target seismic record for site — The input seismic record to the FLAC3D model should produce a calculated motion that can be matched to the target design earthquake motion. If the target motion is at a different location than the model input motion, then a deconvolution analysis should be performed to apply an input motion that will produce a motion that is comparable to the target motion at the target’s location.

   Note that a different form of the deconvoluted motion is applied for a rigid-base model than for a compliant-base model in order to reproduce the target motion. (See here for further explanation.) For the rigid-base boundary, the input motion to FLAC3D is the within motion calculated by SHAKE-91 at the specified depth. For the compliant-base boundary, the input motion is the upward-propagating motion, which is half of the outcrop motion output by SHAKE at the specified depth. Figure 1 illustrates the different base input motions:

   \(x\)-velocities at model base (GP 189) and top of wall (GP 4715)

   
   
   

   Figure 1: Deconvolution analysis to obtain base input motion.

7. Apply dynamic loading and boundary conditions — The dynamic boundary and loading conditions should be specified in a specific order for a seismic analysis. This is because the free-field boundaries that are used to represent the lateral extent of the far-field are assigned properties and initial conditions that are transferred from the main-grid zones adjacent to the free-field boundaries. Any changes to these zones or gridpoints after the free field is assigned are not seen by the free field.

   The following sequence is recommended to assign conditions for a dynamic calculation to ensure that these conditions are applied properly:

   1. Turn on the dynamic calculation mode (model dynamic active on) and the large-strain calculation mode (model large-strain on).
   2. Input the filtered and baseline-corrected input velocity that corresponds to the acceleration at the base of the model derived from the deconvolution analysis of the target acceleration. The input velocity can be called into FLAC3D via the table import command.
   3. Using the zone history command, select various parameters to monitor during the dynamic simulation, such as gridpoint accelerations and zone pore pressures and stresses.
   4. Initialize the displacements and velocities in the grid (using the zone gridpoint initialize command), and specify the type of mechanical damping and parameters (using the zone dynamic damping command).
   5. Apply the free-field boundary conditions on the sides of the model using the zone dynamic free-field on command. Then, assign the dynamic boundaries for the seismic loading using the zone face apply command with the table or fish multiplier.

   In order to apply quiet boundary conditions along the same boundary as the dynamic input, the dynamic input must be applied as a stress boundary, because the effect of the quiet boundary will be nullified if the input is applied as an acceleration (or velocity) wave. A velocity record is converted into a shear-stress boundary condition using a two-step procedure:

   1. Convert the velocity wave into a shear stress wave and/or a normal stress wave, depending on the direction of the applied wave relative to the model base, using the formulae on \(\sigma_n\) and \(\sigma_s\) given in Dynamic Loading and Boundary Conditions.
   2. Monitor the acceleration at the foundation surface during the dynamic run to compare this acceleration to the target acceleration. If material failure occurs within the model, this will affect the attenuation of the wave. Also, if the model is shallow, the free surface can cause an increase in the velocity and acceleration that could extend to the base. [2] Some adjustment to the input stress wave may be required in order to produce an acceleration that is comparable to the target acceleration.

8. Perform undamped, elastic-material simulation — Before running a dynamic model with actual material strength and damping properties, an elastic simulation should be made without damping to estimate the maximum levels of cyclic strain and frequency response ranges of the model system. If the cyclic strains are large enough to cause excessive reductions in shear modulus, then the use of additional damping is questionable. In such a case, the damping will be performing outside of its intended range of application. The model properties and input amplitude should be checked if excessive strains are calculated throughout the model.

   Velocity histories should be monitored throughout the model to provide an estimate of the dominant frequencies produced in the model system. Also, shear-strain histories should be recorded to estimate maximum cyclic shear-strain levels when no material damping is provided.

9. Perform simulations with damping and actual strength properties — The results of the undamped run can be used to help select appropriate damping characteristics for the materials in the model. Additional damping may be prescribed for the model in order to damp the oscillation modes identified from the undamped simulation.

   Acceleration histories should be checked first in the initial run with damping and actual strength properties and compared to the target motion for the site. Some adjustment to the model input motion may be required to improve the comparison.

It is important to monitor several different variables during the seismic shaking phase. For example, shear stress/shear strain plots can illustrate the effect of hysteretic damping on the model response. Excess pore-pressure histories can help quantify the potential for liquefaction, and contours of cyclic pore-pressure ratio can delineate regions of liquefaction in a model.

[1]

A badly distorted zone can cause a calculation to stop prematurely in a mesh containing a few tetrahedral zones. It may be possible to prevent this by increasing the strength of the individual zone. This should not significantly affect the model results, provided that the strengths of only a few zones are changed.

[2]

This is a result of the velocity-doubling effect of the free surface. Note that the effect of a free surface on a simple sinusoidal velocity wave extends beneath the surface to, approximately, a depth of one-fourth of the wavelength of the wave transmitted through the medium. (The extent of velocity doubling can be shown simply by applying a wave to the base of a column of zones with a free surface and monitoring the maximum amplitude experienced by each zone as the wave travels through the column.)