SDOF Dynamic Analysis

sdofDynamicAnalysis.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.
# Also notice that for this particular code you need these three files:
#    GroundAccelerationX.txt
#    GroundAccelerationY.txt
#    GroundAccelerationZ.txt
# ------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
from SDOFnewmark import *


# ------------------------------------------------------------------------
# FUNCTION TO INTEGRATE RECORDS
# ------------------------------------------------------------------------
def trapezoidalRecordIntegration(delta_t, record):

    integrated = [0]

    for i in range(1, len(record)):

        integrated.append(integrated[i-1] + delta_t * 0.5 * (record[i-1]+record[i]))

    return integrated

# ------------------------------------------------------------------------
# READ AND PLOT GROUND MOTIONS [Records are in unit of g]
# ------------------------------------------------------------------------

# Load x-direction acceleration
ax = []
f = open("GroundAccelerationX.txt", "r")
lines = f.readlines()
for oneline in lines:
    splitline = oneline.split()
    for j in range(len(splitline)):
        value = 9.81 * 100 * float(splitline[j])
        ax.append(value)

# Load y-direction acceleration
ay = []
f = open("GroundAccelerationY.txt", "r")
lines = f.readlines()
for oneline in lines:
    splitline = oneline.split()
    for j in range(len(splitline)):
        value = 9.81 * 100 * float(splitline[j])
        ay.append(value)

# Load z-direction acceleration
az = []
f = open("GroundAccelerationZ.txt", "r")
lines = f.readlines()
for oneline in lines:
    splitline = oneline.split()
    for j in range(len(splitline)):
        value = 9.81 * 100 * float(splitline[j])
        az.append(value)

# Create time axis
numTimePoints = len(ax)
delta_t = 0.01
t = np.linspace(0, delta_t*numTimePoints, numTimePoints)

# Integrate acceleration to obtain velocity
vx = trapezoidalRecordIntegration(delta_t, ax)
vy = trapezoidalRecordIntegration(delta_t, ay)
vz = trapezoidalRecordIntegration(delta_t, az)

# Integrate velocity to obtain displacement
ux = trapezoidalRecordIntegration(delta_t, vx)
uy = trapezoidalRecordIntegration(delta_t, vy)
uz = trapezoidalRecordIntegration(delta_t, vz)

# Plot
plt.ion()
fig1 = plt.figure(1)
plt.axis('off')
plt.title('Ground Motion')

axes = fig1.add_subplot(331)
axes.get_xaxis().set_visible(False)
plt.plot(t, ax, 'k-', linewidth=1.0)
plt.ylabel('$\ddot{u}_g [cm/s^2]$')

axes = fig1.add_subplot(332)
axes.get_xaxis().set_visible(False)
plt.plot(t, ay, 'k-', linewidth=1.0)

axes = fig1.add_subplot(333)
axes.get_xaxis().set_visible(False)
plt.plot(t, az, 'k-', linewidth=1.0)

axes = fig1.add_subplot(334)
axes.get_xaxis().set_visible(False)
plt.plot(t, vx, 'k-', linewidth=1.0)
plt.ylabel('$\ddot{u}_g [cm/s]$')

axes = fig1.add_subplot(335)
axes.get_xaxis().set_visible(False)
plt.plot(t, vy, 'k-', linewidth=1.0)

axes = fig1.add_subplot(336)
axes.get_xaxis().set_visible(False)
plt.plot(t, vz, 'k-', linewidth=1.0)

axes = fig1.add_subplot(337)
plt.plot(t, ux, 'k-', linewidth=1.0)
plt.ylabel('${u}_g [cm]$')
plt.xlabel('Time [sec]')

axes = fig1.add_subplot(338)
plt.plot(t, uy, 'k-', linewidth=1.0)
plt.xlabel('Time [sec]')

axes = fig1.add_subplot(339)
plt.plot(t, uz, 'k-', linewidth=1.0)
plt.xlabel('Time [sec]')

print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()

# ------------------------------------------------------------------------
# CALCULATE AND PLOT SDOF RESPONSE
# ------------------------------------------------------------------------
m = 100
k = 10000
c = 2 * np.sqrt(k * m) * 0.03

[structuralXdisplacement, structuralXvelocity, structuralXacceleration] = sdofNewmark(m, c, k, delta_t, ax)
[structuralYdisplacement, structuralYvelocity, structuralYacceleration] = sdofNewmark(m, c, k, delta_t, ay)
[structuralZdisplacement, structuralZvelocity, structuralZacceleration] = sdofNewmark(m, c, k, delta_t, az)

plt.ion()
fig2 = plt.figure(2)
plt.axis('off')
plt.title('Structural Responses')

axes = fig2.add_subplot(331)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralXacceleration, 'k-', linewidth=1.0)
plt.ylabel('$\ddot{u} [cm/s^2]$')

axes = fig2.add_subplot(332)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralYacceleration, 'k-', linewidth=1.0)

axes = fig2.add_subplot(333)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralZacceleration, 'k-', linewidth=1.0)

axes = fig2.add_subplot(334)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralXvelocity, 'k-', linewidth=1.0)
plt.ylabel('$\dot{u} [cm/s]$')

axes = fig2.add_subplot(335)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralYvelocity, 'k-', linewidth=1.0)

axes = fig2.add_subplot(336)
axes.get_xaxis().set_visible(False)
plt.plot(t, structuralZvelocity, 'k-', linewidth=1.0)

axes = fig2.add_subplot(337)
plt.plot(t, structuralXdisplacement, 'k-', linewidth=1.0)
plt.ylabel('${u} [cm]$')
plt.xlabel('Time [sec]')

axes = fig2.add_subplot(338)
plt.plot(t, structuralYdisplacement, 'k-', linewidth=1.0)
plt.xlabel('Time [sec]')

axes = fig2.add_subplot(339)
plt.plot(t, structuralZdisplacement, 'k-', linewidth=1.0)
plt.xlabel('Time [sec]')

print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()

plt.ion()
fig3 = plt.figure(3)
axes = fig3.add_subplot(111, projection='3d')
axes.plot(structuralXdisplacement, structuralYdisplacement, structuralZdisplacement, label='Structural displacement', color='red', linewidth=1.0)
axes.plot(ux, uy, uz, label='Ground displacement', color='black', linewidth=1.0)

axes.set_xlabel('x-displacement [cm]')
axes.set_ylabel('y-displacement [cm]')
axes.set_zlabel('y-displacement [cm]')

axes.legend()
axes.view_init(60, 20)

print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()

# ------------------------------------------------------------------------
# CALCULATE AND PLOT RESPONSE SPECTRA
# ------------------------------------------------------------------------
dampingRatio = 0.03
minPeriod = 0.05
maxPeriod = 10
numSpectrumPoints = 200
xDisplacementSpectrum = []
xVelocitySpectrum = []
xAccelerationSpectrum = []
yDisplacementSpectrum = []
yVelocitySpectrum = []
yAccelerationSpectrum = []
zDisplacementSpectrum = []
zVelocitySpectrum = []
zAccelerationSpectrum = []
periodAxis = []

for i in range(numSpectrumPoints):

    period = minPeriod + i * maxPeriod/numSpectrumPoints

    periodAxis.append(period)

    k = 4 * m * np.pi**2 / period**2

    c = 2 * np.sqrt(k * m) * dampingRatio

    [disp, vel, accel] = sdofNewmark(m, c, k, delta_t, ax)

    xDisplacementSpectrum.append(np.max(np.abs(disp)))
    xVelocitySpectrum.append(np.max(np.abs(vel)))
    xAccelerationSpectrum.append(np.max(np.abs(accel))/981) # ... in units of g

    [disp, vel, accel] = sdofNewmark(m, c, k, delta_t, ay)

    yDisplacementSpectrum.append(np.max(np.abs(disp)))
    yVelocitySpectrum.append(np.max(np.abs(vel)))
    yAccelerationSpectrum.append(np.max(np.abs(accel)))

    [disp, vel, accel] = sdofNewmark(m, c, k, delta_t, az)

    zDisplacementSpectrum.append(np.max(np.abs(disp)))
    zVelocitySpectrum.append(np.max(np.abs(vel)))
    zAccelerationSpectrum.append(np.max(np.abs(accel)))

plt.ion()
fig4 = plt.figure(4)
plt.axis('off')
plt.title('Response spectra')

axes = fig4.add_subplot(331)
plt.loglog(periodAxis, xAccelerationSpectrum, 'k-', linewidth=1.0)
plt.ylabel('$S_a [cm/s^2]$')

axes = fig4.add_subplot(332)
plt.loglog(periodAxis, yAccelerationSpectrum, 'k-', linewidth=1.0)

axes = fig4.add_subplot(333)
plt.loglog(periodAxis, zAccelerationSpectrum, 'k-', linewidth=1.0)

axes = fig4.add_subplot(334)
plt.plot(periodAxis, xVelocitySpectrum, 'k-', linewidth=1.0)
plt.ylabel('$S_v [cm/s]$')

axes = fig4.add_subplot(335)
plt.plot(periodAxis, yVelocitySpectrum, 'k-', linewidth=1.0)

axes = fig4.add_subplot(336)
plt.plot(periodAxis, zVelocitySpectrum, 'k-', linewidth=1.0)

axes = fig4.add_subplot(337)
plt.plot(periodAxis, xDisplacementSpectrum, 'k-', linewidth=1.0)
plt.ylabel('$S_d [cm]$')
plt.xlabel('Period [sec]')

axes = fig4.add_subplot(338)
plt.plot(periodAxis, yDisplacementSpectrum, 'k-', linewidth=1.0)
plt.xlabel('Period [sec]')

axes = fig4.add_subplot(339)
plt.plot(periodAxis, zDisplacementSpectrum, 'k-', linewidth=1.0)
plt.xlabel('Period [sec]')

print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()