Site Response Analysis of a Layered Soil Column (Total Stress Analysis)

From OpenSeesWiki
Revision as of 01:16, 16 September 2010 by Mcganncr (talk | contribs)
Jump to navigationJump to search

Example posted by: Christopher McGann, University of Washington

This article describes the OpenSees implementation of a site response analysis a soil deposit located above an elastic half-space using total stress analysis. A single soil column is modeled in two-dimensions and is subject to an earthquake ground motion in a manner which accounts for the finite rigidity of the underlying medium.

Provided with this article are three input files which execute this analysis in different ways.

Also included is the velocity time history of the selected ground motion, velocityHistory.out. This file is required for the analysis.

Several helpful files are included as well:

  • GilroyNo1EW.out, the acceleration time history for the selected ground motion, as downloaded from the Peer NGA strong motion database, with the header information removed.
  • processMotion.m, a Matlab script which reshapes the GilroyNo1EW.out acceleration history into a column vector, then integrates this using the trapezoidal method to obtain the velocity time history needed for the site response analysis.
  • accelPlot.m and respSpectra.m, Matlab scripts which produce the verification plots included in this article from the recorded results.

Download them all in a compressed file: siteResponseAnalysis.zip

To run this example, the user must download the one of the input files files () and the velocity time history file, velocityHistory.out, and place them in a single directory. Once this has been done, the user can then run the analysis. The Matlab scripts and the acceleration time history files are not essential to the analysis, however, they are provided to demonstrate how an alternative acceleration time history can be converted into the analysis and how certain plots can be obtained from the recorded output.

A large set of similar examples has been developed at the University of California at San Diego. They are available at http://cyclic.ucsd.edu/opensees. These examples utilize a wide variety of element and material formulations, as well as different boundary and loading conditions than those used in this example. The user is referred to these examples for further information on how to set up a site response model in OpenSees.

Model Description

Fig 1. Schematic representation of the site response model. Node numbers in blue, element numbers in red.

In this article, the site response analyses are performed for soil deposits which are underlain by an elastic half-space, which simulates the finite rigidity an underlying medium such as bedrock. It is assumed that there is no groundwater, therefore, total stress analysis is used in this example. The soil is modeled in two-dimensions with two degrees-of-freedom using the plane strain formulation of the quad element element. The provided input files consider only a single layer of soil, however, the files freeFieldIndepend.tcl and freeFieldDepend.tcl, consider a soil profile which has a parabolically-increasing shear wave velocity profile with depth. This is accomplished using sublayers, thus, the use of multiple soil layers in an OpenSees site response analysis is demonstrated. A general schematic of the model is shown in Fig. 1. In this example, the horizontal direction is the first degree-of-freedom and the vertical is the second. The soil node, element, and layer numbering schemes all begin at the bottom.

Note: The remainder of the main model description in this article, when not otherwise specified, refers to the simpler of the three provided input files, freeFieldSingle.tcl, with the intention of highlighting the main site response analysis details. The other input files are essentially the same, with the main difference being the complexity introduced by the inclusion of multiple layers. The more complex files are explained later in this article during the discussion of the comparison study.

To account for the finite rigidity of the underlying half-space, a Lysmer-Kuhlemeyer (1969) dashpot is incorporated at the base of the soil column using a zeroLength element and the viscous uniaxial material. The Lysmer-Kuhlemeyer (1969) dashpot is assigned a dashpot coefficient equal to the product of the mass density and shear wave velocity of the underlying layer with the area of the base of the soil column. In this example, the properties of a bedrock layer are used for the half-space. The soil column is excited at the base by a horizontal force time history which is proportional to the known velocity time history of the ground motion. Further information on this modeling approach can be found in Joyner and Chen (1975) and Lysmer (1978) among others.

The horizontal force time history is applied as a Path TimeSeries object using the file, velocityHistory.out. The force history is obtained by multiplying the known velocity time history by a constant factor set as the product of the area of the base of the soil column (width x thickness) with the mass density and shear wave velocity of the underlying layer. The area of the soil column is included to ensure proportional loading for any desired horizontal element size. The provided Matlab script, processMotion.m, shows how the velocity time history was obtained from the selected acceleration time history.

Soil Profile Geometry

For the simple soil profile considered in the file freeFieldSingle.tcl, the soil profile geometry is controlled entirely by the thickness of the soil deposit. In this example, this value is set at 40 meters.

In the multi-layered soil profiles in the other input files, the soil profile geometry is controlled by a specified number of layers and the thicknesses of each individual soil layer. In this case the overall thickness of the soil profile is the sum of the layer thicknesses.

Mesh Geometry

The geometry of the mesh is based upon the concept of resolving the propagation of the shear waves at or below a particular frequency by ensuring that an adequate number of elements fit within the wavelength of the chosen shear wave. This ensures that the mesh is refined enough such that the desired aspects of the propagating waves are well captured in the analysis.

This is accomplished in the model by specifying the highest frequency which the user wishes to be well resolved and the number of elements that the user desires to be in one wavelength of a shear wave propagating at this frequency. The wavelength used to define the mesh geometry is determined as the product of the specified cutoff frequency and the minimum shear wave velocity of the soil profile. In the provided examples there are either 8 or 10 elements per wavelength. The largest value possible value for the specified frequency is the sampling frequency of the ground motion. The provided examples use 100 Hz for the 'cutoff' frequency.

The horizontal size of the elements is set to be the minimum vertical element size in the soil column. The number of nodes and the total number of elements are computed automatically. The number of elements is based on the computed element size and the thickness of the soil deposit. For n total elements, there are 2n + 1 total nodes.

Soil Nodes

The soil nodes are created automatically from the input geometry and meshing information. As shown in Fig. 1, the node numbering scheme is left-to-right, top-to-bottom. Nodes with even numbers fall on the y-axis, and the odd-numbered nodes are spaced horizontally by the computed horizontal element size.

Dashpot Nodes

A single zeroLength element is used to define the Lysmer-Kuhlemeyer (1969) dashpot, therefore, only two nodes are required. These nodes are arbitrarily assigned numbers 2000 and 2001. If the user has modified the meshing information such that there are more than 2000 nodes in the soil column, the dashpot node numbers will need to be changed.

Boundary Conditions and Equal Degrees-of-Freedom

The nodes at the base of the column are fixed against displacements in the y-direction in accordance with the assumption that the soil layers are underlain by bedrock. The remaining soil nodes are then tied together using the equalDOF command in order to achieve a simple shear deformation pattern. This is done by declaring equalDOF for every pair of nodes which share the same y-coordinate.

One of the dashpot nodes is fully fixed (node 2000), while the other is fixed only against displacements in the y-direction (node 2001). To incorporate the dashpot element into the total model, equalDOF is again used, this time linking the horizontal degrees-of-freedom of the partially fixed dashpot node and one of the nodes at the base of the soil column.

Soil Material Properties and Objects

A series of material properties are required to define the constitutive behavior of the soil and the underlying elastic half-space. The main soil properties are the mass density, the shear wave velocity, and Poisson's ratio. From these, the elastic, shear, and bulk moduli are computed. Poisson's ratio is set as zero in these examples for the purposes of emulating a one-dimensional analysis. With Poisson's ratio at zero, no vertical accelerations are generated.

The remaining soil properties correspond to the particular constitutive model selected in each analysis, either PressureDependMultiYield or PressureIndependMultiYield. For these examples, these properties are largely based on the recommendations detailed on the PressureDependMultiYield or PressureIndependMultiYield pages of the OpenSees command manual.

For the case of multiple layers, one material object is created for each layer using the material properties defined in the input file. With the exception of a few global properties, each layer is given a separate set of properties.

During the creation of these examples is was noted that using too few yield surfaces in these nested surface models, or using poorly spaced user-generated surfaces, can cause high frequency shaking to develop in the model. This appears to be due to the way in which these constitutive models work. There is not one consistent tangent over the length of the backbone curve, instead, a smooth backbone curve is approximated using a series of linear segments. If there are too few, or poorly spaced, yield surfaces, a dramatic change in stiffness can occur as the stress path crosses a particular yield surface, resulting in a shock which is manifested as the observed high frequency shaking. This behavior can be alleviated by using a larger number of yield surfaces, resulting in a smoother approximation to the backbone curve of the soil. Care should be taken, however, depending upon the intentions of the analysis. As the number of surfaces increases, the time needed to complete the analysis also increases, and in the end, the inclusion of this high frequency shaking does not significantly alter the results of the simulation. It is recommended that the user keep these ideas in mind when selecting the number of surfaces to use in a particular analysis.

Soil Elements

Four-node quad elements are used to model the soil using the plane strain formulation of the quad element. The element connectivity uses a counterclockwise pattern for the previously-described node numbering scheme (see Fig. 1). The soil elements in each layer are assigned the material tag of the material object corresponding to that layer. A unit thickness is used in all examples. The self-weight of the soil is considered as a body force acting on each element. The body force is set as the unit weight of the soil in each layer, which is determined from the respective mass density input value.

Dashpot Material and Element

The viscous uniaxial material is used to define the Lysmer-Kuhlemeyer (1969) dashpot. This material model requires a single input, the dashpot coefficient, c. Following the method of Joyner and Chen (1975), the dashpot coefficient is defined as the product of the mass density and shear wave velocity of the underlying medium, which, in this example, is assumed to have the properties of bedrock. The dashpot coefficient must also include the area of the base of the soil column to maintain proportional results for any horizontal element size. Since the elements are created with unit thickness, this area is simply the computed horizontal element size.

A zeroLength element is used for the Lysmer-Kuhlemeyer (1969) dashpot. This element connects the two previously-defined dashpot nodes and is assigned the material tag of the Viscous uniaxial material object in the first degree-of-freedom (horizontal direction).

Recorders

The recorders defined for this example document the following items:

  • the nodal displacements, velocities, and accelerations in both degrees-of-freedom
  • the stress and strain response at each Gauss point in each element

The recorded nodal values are the true values for each parameter, not relative values as with the uniformExcitation command.

Two sets of recorders are included to capture the response of the soil during the gravity analysis (which is sometimes useful to ensure proper model generation) and the response of the soil during the application of the ground motion. These files are differentiated by a naming scheme in which the records corresponding to the gravity analysis are named 'Gacceleration', 'Gstress', etc ...

Gravity Loading and Analysis

The gravity analysis in this example is conducted as a transient analysis with very large time steps, thus simulating a static analysis while avoiding the conflicts which may occur when mixing static and transient analyses. The self-weight of the soil elements provides the loads, therefore, no loading object is required. Gravity is applied for ten steps with entirely elastic constitutive behavior. This allows the material objects to update various parameters to account for confining pressure. Once these steps have converged, the material objects are updated using the updateMaterialStage command to consider elastoplastic behavior, and the gravity analysis is repeated for 40 steps.

Horizontal Loading and Analysis

Using the method of Joyner and Chen (1975), dynamic excitation is applied as a force time history to the base of the soil column, at the node which shares equal degrees-of-freedom with the Lysmer-Kuhlemeyer (1969) dashpot. This force history is obtained by multiplying the velocity time history of the recorded ground motion by the mass density and shear wave velocity of the underlying bedrock layer and the area of the base of the soil column (the horizontal element size in these examples). This technique considers the finite rigidity of the underlying layer by allowing energy to be radiated back into the underlying material.

The force history is applied to the model as a Path TimeSeries object using a Plain load pattern object. The actual force applied to the node in each time step is the product of the load factor indicated in the pattern object (1.0 in this example), the additional load factor included in the timeSeries object, and the value found in the file, velocityHistory.out, at that time step. The load factor included in the timeSeries object is used to create a force history from the velocity history found in the namesake file.

The transient analysis is conducted with the Newmark integrator using the gamma and beta coefficients defined near the top of the input file. These values are set at 0.5 and 0.25, respectively to ensure there is no numerical damping in the analysis.

Since these models consider elastoplastic soil behavior, there is inherent hysteretic damping which occurs, however, a small amount of Rayleigh damping is used so there is still some damping at low strain values. The level of Rayleigh damping is controlled by the damping ratio and the circular frequencies of two modes of vibration. This type of approach is discussed in most dynamic analysis texts, such as Chopra (2007).

The time step used in the analysis is selected to meet stability considerations using the Courant-Friedrichs-Levy condition (discussed well by LeVeque, 2007). Meeting the requirements of this condition ensures that the time step is small enough for stability given the maximum shear wave velocity and the minimum vertical element size in the model. The example files are set up such that the analysis time step and the corresponding number of steps are automatically generated to satisfy this condition. If the required time step is larger than the time step increment found in the recorded ground motion, the value specified in the ground motion is used. If the user alters the shear wave velocity profile such that the lowest layer does not have the highest shear wave velocity, then the value of $vsMax used in this routine will need to be updated.

Representative Results

Fig 2. Acceleration time history response of at the ground surface of the soil column. Figure generated using accelPlots.m


Several plots displaying results obtained from running this example site response analysis are provided for use as a means of verification. These plots were generated using the Matlab scripts getAccel.m, accelPlots.m, and depthPlots.m, which are available for download with the rest of the files in a compressed folder here.

Fig. 2 shows the acceleration response at the ground surface (top of the soil column) due to the applied force history. This figure was created by simply plotting the recorded response at each time increment. The m-file, accelPlots.m was used to create the figure.

Fig. 3 shows the computed acceleration, velocity, and displacement response spectra at the ground surface. The figure was created using the m-file, accelPlots.m.

Fig. 4 shows the variation of several recorded parameters with depth. The dashed lines represent the soil layer boundaries. The m-file, depthPlots.m was used to generate this figure.

Fig 3. Acceleration, velocity, and displacement response spectra at the ground surface of the soil column. Figure generated using accelPlots.m
Fig 4. Variation of maximum nodal acceleration and displacement and maximum shear stress and strain with depth. Figure generated using depthPlots.m

Comparison of Results with Shake

A comparison study was accomplished for the example soil profile using the one-dimensional, equivalent linear site response analysis program EduShake (http://www.proshake.com/). The results of this study to not correlate well with the OpenSees site response analysis results due to several factors such as the way the programs handle modulus reduction and damping.

In the example problem, the default modulus reduction (G/Gmax) curves were used. These curves do not match those used by EduShake. User-defined modulus reduction curves can be used in the PressureDependMultiYield material, however, this can create difficulties with the friction angle. A good approach is to define PressureIndependMultiYield materials for many sublayers within each of the soil layers, each with different parameters, and utilize user-defined G/Gmax curves.


The example contained in this article is for the purposes of familiarizing users with how to conduct a site response analysis using OpenSees. A forthcoming example will describe some techniques which can aid in achieving results which are more similar to an EduShake (or ProShake) analysis.

References

  1. Chopra, A.K. (2007). Dynamics of Structures, Third Edition, Pearson Prentice Hall.
  2. Joyner, W.B. and Chen, A.T.F. (1975). "Calculation of nonlinear ground response in earthquakes," Bulletin of the Seismological Society of America, Vol. 65, No. 5, pp. 1315-1336, October 1975.
  3. LeVeque, R.J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia.
  4. Lysmer, J. (1978). "Analytical procedures in soil dynamics," Report No. UCB/EERC-78/29, University of California at Berkeley, Earthquake Engineering Research Center, Richmond, CA.
  5. Lysmer, J. and Kuhlemeyer, A.M. (1969). "Finite dynamic model for infinite media," Journal of the Engineering Mechanics Division, ASCE, 95, 859-877.
Retrieved from "https://opensees.berkeley.edu/wiki/index.php?title=Site_Response_Analysis_of_a_Layered_Soil_Column_(Total_Stress_Analysis)&oldid=5125"