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

Example posted by: <span style="color:blue"> Christopher McGann, University of Washington</span>

----

This article describes the OpenSees implementation of a site response analysis for a layered soil profile 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 is the main input file needed to execute this analysis in OpenSees, [[FreeFieldDamp.tcl|freeFieldDamp.tcl]], along with several additional necessary and/or helpful files:
* ELCENTRO.EQ, the acceleration time history for the considered ground motion
* forceHistory.out, the path timeSeries information needed to run the analysis
* getForceHistory.m, a Matlab script which computes the timeSeries information from the acceleration time history
* getAccel.m, accelPlots.m, and depthPlots.m, Matlab scripts which produce the plots included in this article from the recorded results

Download them all in a compressed file: [[Media:siteResponseTotal.zip|siteResponseTotal.zip]]

To run this example, the user must download the files freeFieldDamp.tcl and forceHistory.out and place them in a single directory.  Once this has been done, the user can then type "source freeFieldDamp.tcl" into the interpreter of the OpenSees.exe application to run the analysis.  The [http://www.mathworks.com 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, and utilize 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 to incorporate effective stress analysis and/or cohesive materials.

== Model Description ==

[[File:siteRespModel.png|frame|Fig 1.  Schematic representation of the site response model.  Node numbers in blue, element numbers in red.]]

The site response analysis is performed for a soil profile with three layers of cohesionless soil.  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.  The nDMaterial model, [http://opensees.berkeley.edu/OpenSees/manuals/usermanual/1550.htm PressureDependMultiYield], is used as a constitutive model for the cohesionless soil, and each layer is assigned separate material properties.  A 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. 

To account for the finite rigidity of the underlying medium (assumed to be bedrock in this example), a Lysmer-Kuhlemeyer (1969) dashpot is incorporated at the base of the soil column using a [[zeroLength Element|zeroLength element]] and the [[Viscous Material| 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 bedrock layer.  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, forceHistory.out.  This force time history has been precomputed from the acceleration time history detailed in the file ELCENTRO.EQ.  The provided Matlab script, getForceHistory.m, shows how this computation was accomplished, and can be modified to produce the required force time history for any given acceleration time history.

=== Soil Profile Geometry ===

There are three layers of cohesionless soil in this example.  The user can specify the thickness of each layer.  The default layer thickness values are:
* uppermost layer, 2 m thick
* middle layer, 8 m thick
* lower layer, 40 m thick

=== Mesh Geometry ===

In this example, each of the soil elements is a one meter square.  This is controlled by defining the element size in the x-direction and the number of elements in each layer.  These are the only required inputs in this section.  The number of nodes and the total number of elements are computed automatically.  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 input horizontal element size (1 m in this example).

=== Dashpot Nodes ===

A single [[zeroLength Element|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 ===

The [http://opensees.berkeley.edu/OpenSees/manuals/usermanual/1550.htm PressureDependMultiYield] nDMaterial is used to define the constitutive behavior of the soil.  One material object is created for each of the three layers in this example 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.  For this example, these properties are largely based on the recommendations detailed on the [http://opensees.berkeley.edu/OpenSees/manuals/usermanual/1550.htm PressureDependMultiYield] page of the OpenSees command manual.  

=== Soil Elements ===

Four-node quad elements are used to model the soil using the plane strain formulation of the [[Quad Element |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.  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 Material | 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 be bedrock having a mass density of 2.4 Mg/m^3 and a shear wave velocity of 914.4 m/s.

A [[zeroLength Element |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 [[Recorder Command |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.  A real time display recorder is included for visualization of the response of the soil column during the analysis.

=== Gravity Loading and Analysis ===

The gravity analysis in this example is conducted as a [[analysis Command |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 [http://opensees.berkeley.edu/OpenSees/manuals/usermanual/1559.htm 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.  This technique considers the finite rigidity of the bedrock 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 Pattern |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 (also 1.0 in this example), and the value found in the file, forceHistory.out, at that time step.  The loading could be scaled by modifying the load factors to produce the desired effect.

The transient analysis is conducted with the [[Newmark Method | Newmark integrator]] using the same gamma and beta coefficients defined in the gravity analysis.  The time step increment and number of steps used in the analysis match those of the original acceleration time history. 

== Representative Results ==

[[File:surfaceAccel.png|frame|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 [[Media:siteResponseTotal.zip|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.

[[File:surfaceSpectra.png|frame|left|Fig 3.  Acceleration, velocity, and displacement response spectra at the ground surface of the soil column.  Figure generated using accelPlots.m]]

[[File:depthSummary.png|frame|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.

It is possible to define the material behavior in the OpenSees analysis such that the results are comparable to those from a site response analysis program such as EduShake, however, it much simple to do so using the PressureIndependentMultiYield material model.  

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 ==

# 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.
# 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.
# Lysmer, J. and Kuhlemeyer, A.M. (1969). "Finite dynamic model for infinite media," ''Journal of the Engineering Mechanics Division, ASCE'', '''95''', 859-877.