Skip to content

Instantly share code, notes, and snippets.

@geosharma
Last active November 8, 2019 02:12
Show Gist options
  • Save geosharma/de1cede130bb39c8b29d51563f452b4b to your computer and use it in GitHub Desktop.
Save geosharma/de1cede130bb39c8b29d51563f452b4b to your computer and use it in GitHub Desktop.
Code snippets
#!/usr/bin/env python
# -*-coding:utf-8 -*-
import numpy as np
import openseespy.opensees as ops
from pathlib import Path
# Example problem:
# Description: 1D Total Stress Site Response Analysis, Single soil layer
# Original authors of this example: Christopher McGann and Pedro Arduino, University of Washington
# ref: https://opensees.berkeley.edu/wiki/index.php/Site_Response_Analysis_of_a_Layered_Soil_Column_(Total_Stress_Analysis)
# Units: kN, m, sec
def roundup(num):
"""Check if a float is a whole number if not round up"""
return int(num) if num.is_integer() else int(np.ceil(num))
def data_recorders(fpath, fname_prefix, number_of_nodes, number_of_elements, motiondt):
# this is a deviation from the original solution, since only ground
# accelerations, velocity and displacements are considered, for now only ground
# node values will be recorded
# fmt: off
# recorders for nodal displacements, velocity and acceleration at each time step
restypes = ["disp", "vel", "accel"]
for restype in restypes:
fname = fname_prefix + restype + ".out"
print(f"Filename: {fname}")
fout = str(fpath.joinpath(fname))
ops.recorder(
"Node",
"-file", fout,
"-time",
"-dT", motiondt,
"-nodeRange", 1, number_of_nodes,
"-dof", *[1, 2],
restype,
)
# # ref: http://opensees.berkeley.edu/wiki/index.php/PressureIndependMultiYield_Material
# # For 2D problems, the stress output follows this order: σxx, σyy, σzz, σxy, ηr, where ηr
# # is the ratio between the shear (deviatoric) stress and peak shear strength at the current
# # confinement (0<=ηr<=1.0). The strain output follows this order: εxx, εyy, γxy.
# stresses and strains at each Gauss point in the soil elements
gausspoints = ["1", "2", "3", "4"]
for gpt in gausspoints:
fname = fname_prefix + "_stress" + gpt + ".out"
fout = str(fpath.joinpath(fname))
ops.recorder(
"Element",
"-file", fout,
"-time",
"-dT", motiondt,
"-eleRange", 1, number_of_elements,
"material", gpt,
"stress",
)
fname = fname_prefix + "_strain" + gpt + ".out"
fout = str(fpath.joinpath(fname))
ops.recorder(
"Element",
"-file", fout,
"-time",
"-dT", motiondt,
"-eleRange", 1, number_of_elements,
"material", gpt,
"strain",
)
# fmt: on
def main():
# file path and filenames
path_root = Path(r"./")
path_data = path_root / "data"
# general filename for recorder outputs
fname_prefix = "ex_g_02_siteresponse_totalstress_"
# velocity time history file
velfile = "velocityhistory.out"
fvel = str(path_data.joinpath(velfile))
print(f"Ground motion file path: {fvel}")
fvel = "data/velocityHistory.out"
# 1. Define analysis parameters
# soil profile, thickness of the soil layer [m]
depth_soillayer = 40.0
# acceleration due to gravity [m/s2]
accg = 9.81
# material properties
# # soil saturated density [Mg/m3], or soil mass density
soil_satdensity = 1.7
# soil shear wave velocity [m/s]
soil_vs = 250.0
# Poisson's ratio of the soil
nu = 0.0
# the undrained shear strength, c of Mohr-Coulomb failure envelop [kPa]
cohesion = 95.0
# soil friction angle, phi of Mohr-Coulomb failure envelop [deg]
friction_angle = 0.0
# peak shear strain
peak_shear_strain = 0.05
# reference pressure [kPa]
reference_pressure = 80.0
# pressure dependency coefficient
pressure_dependence = 0.0
# bedrock shear wave velocity [m/s]
rock_vs = 760.0
# bedrock mass density [Mg/m^3]
rock_density = 2.4
# soil shear modulus [kPa]
modG = soil_satdensity * soil_vs ** 2
# soil elastic modulus [kPa]
modE = 2 * modG * (1 + nu)
# soil bulk modulus [kPa]
modK = modE / (3 * (1 - 2 * nu))
# number of yield surfaces
nsurf = 25
print(f"Soil elastic properites:")
print(f"G: {modG/1000:6.1f} MPa, E: {modE/1000:6.1f} MPa, K: {modK/1000:6.1f} MPa")
# ground motion paramters
# time step in ground motion record [s]
gm_dt = 0.005
# number of steps in ground motion record
# this depends on the number of points in the ground motion file
gm_nsteps = 7990
# wavelength parameters
# highest frequency desired to be well resolved [Hz]
fmax = 100.0
# wavelength of the highest resolved frequency
lamdamin = soil_vs / fmax
# number of elements per one wavelength
nele_wave = 10
# rayleigh damping parameters (read on this)
# damping ratio
damping_ratio = 0.02
# lower frequency [Hz]
f_lower = 0.2
omega_lower = 2 * np.pi * f_lower
# upper frequency [Hz]
f_upper = 20
omega_upper = 2 * np.pi * f_upper
# rayleigh damping coefficients
alpham = 2.0 * damping_ratio * omega_lower * omega_upper / (omega_lower + omega_upper)
betak = 2.0 * damping_ratio / (omega_lower + omega_upper)
betakinit = 0.0
betakcomm = 0.0
print(f"Damping coefficients: alpham = {alpham:6.3f}, betak = {betak:6.5f}")
# Newmark parameters
newmark_gamma = 0.5
newmark_beta = 0.25
# acceleration due to gravity in the x and y direction
accg_x = 0.0
accg_y = -accg
wgt_x = 0.0
wgt_y = accg_y * soil_satdensity
# 2. define mesh geometry
# number of element in x and y direction
nx = 1
# trial for vertical element size
# initial element size
h_trial = lamdamin / nele_wave
# if the number of elements is not an integer round up
ny = roundup(depth_soillayer / h_trial)
print(f"Number of elements: nx = {nx}, ny = {ny}")
# size of elements in x and y direction of the nodes
ly = float(depth_soillayer / ny)
lx = ly
print(f"Element size: lx = {lx:6.3f} m, ly = {ly:6.3f} m")
# soil element thickness
thickness = 1.0
# analysis type, "PlaneStrain" or "PlaneStress"
analysis_type = "PlaneStrain"
# surface pressure on elements
surface_pressure = 0.0
# start node number
start_node_number = 1
# number of nodes in the vertical direction
num_node_y = 2 * (ny + 1)
# wipe previous analyses
ops.wipe()
# corner nodes are 2D with 2 DOF (displacements in x1 and x2)
ops.model("basic", "-ndm", 2, "-ndf", 2)
ycoord = 0.0
count = 0
# loop over nodes
a = range(1, num_node_y + 1, 2)
for j in range(1, num_node_y, 2):
ycrd = ycoord + count * ly
ops.node(j, 0.0, ycrd)
ops.node(j+1, lx, ycrd)
print(f"Node: {j}, x: 0.0, y: {ycrd}")
print(f"Node: {j + 1}, x: {lx}, y: {ycrd}")
count += 1
print("Finished creating all soil nodes.....")
# 4. Define dashpot nodes
ops.node(2000, 0.0, 0.0)
ops.node(2001, 0.0, 0.0)
print("Finished creating dashpot nodes....")
# 5. Define boundary conditions and equalDOF
# definie fixity of the base nodes
ops.fix(1, 0, 1)
ops.fix(2, 0, 1)
# define fixity of the dashpot nodes
ops.fix(2000, 1, 1)
ops.fix(2001, 0, 1)
# define equalDOF for simple shear deformation of soil elements
for k in range(3, num_node_y, 2):
ops.equalDOF(k, k + 1, 1, 2)
print(f"equalDOF: rnode: {k}, cnode: {k + 1}, dofs: 1 2")
# define equalDOF for dashpot and base soil nodes
ops.equalDOF(1, 2, 1)
ops.equalDOF(1, 2001, 1)
print(f"Created all boundary conditions and equalDOF:")
# 6. Define soil materials
soil_mattag = 1
# nd = 2 for plane strain and 3 for 3-D
nd = 2
# number of dimensions (nd) 2 for plain-strain 3 for 3D
ops.nDMaterial(
"PressureIndependMultiYield",
soil_mattag,
nd,
soil_satdensity,
modG,
modK,
cohesion,
peak_shear_strain,
friction_angle,
reference_pressure,
pressure_dependence,
nsurf,
)
print(f"Soil material created")
element_mass_density = 0.0
print(
"{0:5s}, {1:6s}, {2:6s}, {3:6s}, {4:6s}".format(
"Ele", "Node i", "Node j", "Node k", "Node l"
))
for j in range(1, ny + 1):
nodei = 2 * j - 1
nodej = nodei + 1
nodek = nodei + 3
nodel = nodei + 2
ops.element("quad", j, nodei, nodej, nodek, nodel, 1.0, "PlaneStrain", 1.0, 0.0, 0.0, wgt_x, wgt_y)
print("{0:5d}, {1:6d}, {2:6d}, {3:6d}, {4:6d}".format(j, nodei, nodej, nodek, nodel))
print(f"All soil elements created")
# material and elements for viscous damping
# dashpot coefficients, again figure out how this coefficient comes about
dashpot_c = lx * rock_density * rock_vs
# for linear damping, power factor = 1
dashpot_alpha = 1.0
# dashpot material
dashpot_mattag = 4000
ops.uniaxialMaterial("Viscous", 4000, dashpot_c, 1.0)
# element, connecting two dashpot nodes
dashpot_elenum = 5000
ops.element(
"zeroLength",
5000,
*[2000, 2001],
"-mat",
4000,
"-dir",
*[1],
)
print("Finished creating dashpot material and element...")
print(f"Gravity analysis:")
# update material stage to 0 for elastic behavior
# ops.updateMaterialStage("-material", soil_mattag, "-stage", 0)
print(f"Material updated for elastic behavior during gravity loading")
# gravity loading
ops.constraints("Transformation")
ops.test("NormDispIncr", 1e-5, 30, 1)
ops.algorithm("Newton")
ops.numberer("RCM")
ops.system("ProfileSPD")
ops.integrator("Newmark", newmark_gamma, newmark_beta)
ops.analysis("Transient")
ops.analyze(10, 500.0)
print(f"Elastic gravity analysis over")
# update material stage to include plastic analysis
ops.updateMaterialStage("-material", soil_mattag, "-stage", 1)
ops.analyze(40, 500.0)
print(f"Plastic gravity analysis over")
ops.setTime(0.0)
ops.wipeAnalysis()
print(f"Create data recorders")
# specify the node and element at which the acceleration, velocity, displacement,
# stress and strain time histories are required, for surface node
last_node_number = 322
last_element_number = 160
data_recorders(
path_data, fname_prefix, last_node_number, last_element_number, gm_dt
)
print(f"Finished creating data recorders")
# destory previous analysis object and start consolidation
print(
f"Wiping gravity analysis, resetting time to 0.0, and starting dynamic analysis:"
)
# define constant factor for applied velocity, find out how this comes about?
# this is what I found, the link might not work the next time
# http://www.bgu.ac.il/geol/hatzor/articles/Bao_Hatzor_Huang_2012.pdf
cfactor = dashpot_c
# load pattern and timeseries for applied force history
print(f"Ground motion file: {fvel}, dt: {gm_dt}, cfactor: {cfactor}")
# load timeseries tag
tstag_eq = 20
# load pattern tag
pattag_eq = 10
# node at which the excitation is applied
nodetag_eq = 1
loadvals = [1.0, 0.0]
ops.timeSeries(
"Path", tstag_eq, "-dt", gm_dt, "-filePath", fvel, "-factor", cfactor
)
ops.pattern("Plain", pattag_eq, tstag_eq)
ops.load(nodetag_eq, *loadvals)
# determine the analysis time step using CFL condition
# CFL stands for Courant-Friedrichs-Levy
# ref: http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf
# duration of groundmotion [s]
duration = gm_dt * gm_nsteps
# trial analysis step, I do not know how this came about
ktrial = lx / np.sqrt(rock_vs)
if gm_dt < ktrial:
nsteps = gm_nsteps
dt = gm_dt
else:
nsteps = int(floor(duration / ktrial) + 1)
dt = duration / nsteps
print(f"Number of steps: {nsteps}, time step, dt: {dt}")
print(f"Newmark: gamma: {newmark_gamma}, beta: {newmark_beta}")
print(f"Rayleigh: alpham: {alpham}, betak: {betak}, betakinit: {betakinit}, betakcomm: {betakcomm}")
ops.constraints("Transformation")
ops.test("NormDispIncr", 1.0e-3, 15, 5)
ops.algorithm("Newton")
ops.numberer("RCM")
ops.system("ProfileSPD")
ops.integrator("Newmark", newmark_gamma, newmark_beta)
ops.rayleigh(alpham, betak, betakinit, betakcomm)
ops.analysis("Transient")
ops.analyze(nsteps, dt)
print(f"Dynamic analysis over")
ops.wipe()
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment