MDOF Dynamic Analysis

mdofDynamicAnalysis.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 form.
# Please see the Programming note for how to get started, and notice
# that you must copy certain functions into the file terjesfunctions.py
# Also notice that for this particular code you need the file
# ElCentroGroundMotion.txt, which is posted nearby.
# ------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
import terjesfunctions as fcns


# ------------------------------------------------------------------------
# INPUT [N, m, s, kg]
# ------------------------------------------------------------------------

# Constants
E = 2.0e+11
I = 8.09e-5
A = 0.0091
M = 400.0
storeyHeight = 4.0
bayWidth = 6.0

# Nodes (x, y)
nodes = np.array([(0.0,          0.0),
                  (bayWidth,     0.0),
                  (2.0*bayWidth, 0.0),
                  (0.0,          storeyHeight),
                  (bayWidth,     storeyHeight),
                  (2.0*bayWidth, storeyHeight),
                  (0.0,          2.0*storeyHeight),
                  (bayWidth,     2.0*storeyHeight),
                  (2.0*bayWidth, 2.0*storeyHeight)])

# Frame elements (node1, node2, E, I, A)
frameElements = np.array([(1, 4, E, I, A),
                          (2, 5, E, I, A),
                          (3, 6, E, I, A),
                          (4, 5, E, I, A),
                          (5, 6, E, I, A),
                          (4, 7, E, I, A),
                          (5, 8, E, I, A),
                          (6, 9, E, I, A),
                          (7, 8, E, I, A),
                          (8, 9, E, I, A)])

# Mass (node, value, dof)
mass = np.array([(4, M, 1),
                 (5, M, 1),
                 (6, M, 1),
                 (7, M, 1),
                 (8, M, 1),
                 (9, M, 1)])

# Boundary conditions (node, 0=free, 1=fixed)
constraints = np.array([(1, 1, 1, 1),
                        (2, 1, 1, 1),
                        (3, 1, 1, 1)])


# ------------------------------------------------------------------------
# ASSEMBLE THE STIFFNESS MATRIX Ka
# ------------------------------------------------------------------------

# Number of DOFs
numDOFsPerNodeInStructure = 3
numNodes = nodes.shape[0]
totalNumDOFs = numDOFsPerNodeInStructure * numNodes

# Count elements
numFrameElements = frameElements.shape[0]

# Initialize the big stiffness matrix
Ka = np.zeros((totalNumDOFs, totalNumDOFs))

# Assemble frame elements
KlStorageFrame = []
TlgStorageFrame = []
for i in range(numFrameElements):

    # Information about element
    E = frameElements[i, 2]
    I = frameElements[i, 3]
    A = frameElements[i, 4]

    # Information about nodes
    node1 = int(frameElements[i, 0])
    node2 = int(frameElements[i, 1])
    nodeList = np.array([node1, node2])
    nodalcoords = np.array([nodes[node1-1, 0], nodes[node1-1, 1], nodes[node2-1, 0], nodes[node2-1, 1]])

    # Element length and orientation
    delta_x = nodalcoords[2] - nodalcoords[0]
    delta_y = nodalcoords[3] - nodalcoords[1]
    L = np.sqrt(delta_x ** 2 + delta_y ** 2)

    # Get Kl and Tlg
    nu = 0
    beta = 0
    [Kl, Kgeometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, 0)
    KlStorageFrame.append(Kl)
    TlgStorageFrame.append(Tlg)

    # Calculate Kg = Tlg^T * Kl * Tlg
    Kg = (np.transpose(Tlg).dot(Kl)).dot(Tlg)

    # Assemble
    numDOFsPerNodeInElement = 3
    Ka = fcns.assembleStiffnessMatrix(Ka, Kg, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)


# ------------------------------------------------------------------------
# ASSEMBLE MASS MATRIX
# ------------------------------------------------------------------------

# Add the mass to DOF 1 and 2 of the node, not the rotational DOF
Ma = np.zeros((totalNumDOFs, totalNumDOFs))

numMasses = mass.shape[0]

for i in range(numMasses):
    dof = int((mass[i, 0]-1)*numDOFsPerNodeInStructure + mass[i, 2] - 1)
    Ma[dof, dof] = mass[i, 1]


# ------------------------------------------------------------------------
# INTRODUCE BOUNDARY CONDITIONS & LOCK ZERO-STIFFNESS DOFs
# ------------------------------------------------------------------------

# DOFs locked by input
dofStatus = np.zeros(totalNumDOFs)
nfix = constraints.shape[0]
for k in range(nfix):
    fixedNode = int(constraints[k, 0])
    for m in range(numDOFsPerNodeInStructure):
        if constraints[k, m + 1] == 1:
            dofStatus[(fixedNode - 1) * numDOFsPerNodeInStructure + m] = 1

# DOFs without stiffness
for i in range(totalNumDOFs):
    if Ka[i, i] == 0:
        dofStatus[i] = 1

numFreeDOFs = int(totalNumDOFs - np.sum(dofStatus))

Kf = np.zeros((numFreeDOFs, numFreeDOFs))
Mf = np.zeros((numFreeDOFs, numFreeDOFs))
iota = np.zeros(numFreeDOFs)
iCounter = 0
jCounter = 0
for i in range(totalNumDOFs):

    if dofStatus[i] == 0:

        jCounter = 0

        for j in range(i + 1):

            if dofStatus[j] == 0:
                Kf[iCounter, jCounter] = Ka[i, j]
                Kf[jCounter, iCounter] = Ka[i, j]
                jCounter += 1

        Mf[iCounter, iCounter] = Ma[i, i]
        iota[iCounter] = Ma[i, i]
        iCounter += 1


# ------------------------------------------------------------------------
# EIGENVALUE ANALYSIS
# ------------------------------------------------------------------------

# Here we need scipy, to solve the *generalized* eigenvalue problem [K-omega^2*M]{u} = 0
import scipy

# Determine eigenvalues, which are gamma = omega^2
allEigenvalues = scipy.linalg.eigvals(Kf, Mf)

# Sort them
sortedEigenvalues = np.sort(allEigenvalues)

# Check them
eigenvalues = []
for i in range(len(sortedEigenvalues)):

    if sortedEigenvalues[i].real == np.inf:

        # Remove infinite eigenvalues (zero mass DOFs)
        pass

    else:

        eigenvalues.append(sortedEigenvalues[i].real)

# Calculate corresponding frequencies and periods
numEigenvalues = len(eigenvalues)
naturalFrequencies = []
naturalPeriods = []
print('\n'"Found the following", numEigenvalues, "eigenvalues:"'\n')
for i in range(numEigenvalues):
    frequency = np.sqrt(eigenvalues[i])
    naturalFrequencies.append(frequency)
    period = 2*np.pi/frequency
    naturalPeriods.append(period)
    print("   Frequency: %.2frad, Period: %.3fsec" % (frequency, period))


# ------------------------------------------------------------------------
# SET RAYLEIGH DAMPING MATRIX
# ------------------------------------------------------------------------

# Set target damping ratio
targetDampingRatio = 0.05

# Enforce that damping on the first two modes
omega1 = naturalFrequencies[0]
omega2 = naturalFrequencies[1]

# Calculate alpha and beta in C = alpha*M + beta*K
factor = 2 * targetDampingRatio / (omega1 + omega2)
alpha = omega1 * omega2 * factor
beta = factor

# Create the damping matrix
Cf = np.multiply(Mf, alpha) + np.multiply(Kf, beta)

# Calculate the damping ratio at the different periods
dampingRatios = []
for i in range(numEigenvalues):
    dampingRatios.append(alpha / (2 * naturalFrequencies[i]) + beta * naturalFrequencies[i] / 2)

# Plot those damping ratios
plt.ion()
fig = plt.figure(1)
plt.title("Damping Ratio at Different Natural Periods")
plt.plot(naturalPeriods, dampingRatios, 'ks-', label='Actual damping')
plt.plot([naturalPeriods[0], naturalPeriods[-1]], [targetDampingRatio, targetDampingRatio], 'rs-', label='Target damping')
plt.legend()
plt.xlabel("Period [sec]")
plt.ylabel("Damping Ratio")
print('\n'"Press any key to continue...")
plt.waitforbuttonpress()


# ------------------------------------------------------------------------
# READ GROUND MOTION [Unit of g]
# ------------------------------------------------------------------------

# Load the El Centro ground motion record
groundAcceleration = []
f = open("ElCentroGroundMotion.txt", "r")
lines = f.readlines()
for oneline in lines:
    splitline = oneline.split()
    for j in range(len(splitline)):
        value = 9.81 * float(splitline[j])
        groundAcceleration.append(value)

# Output
print('\n'"Max ground acceleration in units of g:", np.max(np.abs(groundAcceleration)) / 9.81)

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

# Plot the ground motion
plt.ion()
fig = plt.figure(2)
plt.title("El Centro Ground Motion")
plt.plot(t, groundAcceleration, 'k-')
plt.xlabel("Time [sec]")
plt.ylabel("Ground Acceleration [m/s^2]")
print('\n'"Press any key to continue...")
plt.waitforbuttonpress()


# ------------------------------------------------------------------------
# INITIALIZE PLOT
# ------------------------------------------------------------------------

# Initial declarations
Pplus1 = np.zeros(numFreeDOFs)
displacement = np.zeros(numFreeDOFs)
velocity = np.zeros(numFreeDOFs)
acceleration = np.zeros(numFreeDOFs)

# Initialize plot
plt.ion()
fig = plt.figure(3)
plt.autoscale(True)

# Find outer plot limits
xmin = 0.0
xmax = 0.0
ymin = 0.0
ymax = 0.0

for i in range(numNodes):

    if nodes[i, 0] < xmin:
        xmin = nodes[i, 0]

    elif nodes[i, 0] > xmax:
        xmax = nodes[i, 0]

    elif nodes[i, 1] < ymin:
        ymin = nodes[i, 1]

    elif nodes[i, 1] > ymax:
        ymax = nodes[i, 1]

xmin = xmin - (xmax - xmin) * 0.3
ymin = ymin - (ymax - ymin) * 0.3
xmax = xmax + (xmax - xmin) * 0.3
ymax = ymax + (ymax - ymin) * 0.3

# Make the plotting area square
xrange = xmax - xmin
yrange = ymax - ymin

if xrange > yrange:
    ymin -= 0.5 * (xrange - yrange)
    ymax += 0.5 * (xrange - yrange)
else:
    xmin -= 0.5 * (yrange - xrange)
    xmax += 0.5 * (yrange - xrange)

# ------------------------------------------------------------------------
# MDOF NEWMARK ALGORITHM
# ------------------------------------------------------------------------

# Newmark parameters (0.167<beta<0.25) (gamma=0.5)
newmarkBeta = 0.25
newmarkGamma = 0.5

# Newmark constants
a1 = 1.0 / (newmarkBeta * delta_t**2)
a2 = -1.0 / (newmarkBeta * delta_t**2)
a3 = -1.0 / (newmarkBeta * delta_t)
a4 = (1.0 - 1.0 / (2.0 * newmarkBeta))
a5 = newmarkGamma / (newmarkBeta * delta_t)
a6 = -newmarkGamma / (newmarkBeta * delta_t)
a7 = (1.0 - newmarkGamma / newmarkBeta)
a8 = delta_t * (1.0 - newmarkGamma / (2.0 * newmarkBeta))

# Track max displacement
uMax = 0

# Start the loop
for step in range(numTimePoints):

    # Keep track of time
    time = step * delta_t

    # Set the external force vector
    Pplus1 = np.multiply(iota, groundAcceleration[step])

    # Establish the system of equations using the Newmark constants
    newDispFactor = np.multiply(a1, Mf) + np.multiply(a5, Cf)
    accelFactor = np.multiply(a4, Mf) + np.multiply(a8, Cf)
    velocFactor = np.multiply(a3, Mf) + np.multiply(a7, Cf)
    dispFactor = np.multiply(a2, Mf) + np.multiply(a6, Cf)
    Kefficient = newDispFactor + Kf
    Fplus1 = Pplus1 - accelFactor.dot(acceleration) - velocFactor.dot(velocity) - dispFactor.dot(displacement)

    # Solve the system of equilibrium equations
    uPlus1 = np.linalg.solve(Kefficient, Fplus1)

    # Maximum displacement
    uaMax = np.max(np.abs(uPlus1))
    if uaMax > uMax:
        uMax = uaMax

    # Determine increment in displacement
    delu = uPlus1 - displacement

    # Determine acceleration at time step i+1
    accelerationPlus1 = np.multiply(a1, delu) + np.multiply(a3, velocity) + np.multiply(a4, acceleration)

    # Determine velocity at time step i+1
    velocityPlus1 = np.multiply(a5, delu) + np.multiply(a7, velocity) + np.multiply(a8, acceleration)

    # Set response values at time step i equal to response at time step i+1
    displacement = uPlus1
    velocity = velocityPlus1
    acceleration = accelerationPlus1

    # Get ua by adding zeros to the array
    ua = np.zeros(totalNumDOFs)
    iCounter = 0
    for i in range(totalNumDOFs):
        if dofStatus[i] == 0:
            ua[i] = displacement[iCounter]
            iCounter += 1

    # Clear previous lines from plot
    plt.cla()

    # Set scale factor
    scale_factor = 500

    # Get coordinates of displaced nodes
    xcoord = np.zeros(2)
    ycoord = np.zeros(2)
    for i in range(numFrameElements):
        node1 = int(frameElements[i, 0])
        node2 = int(frameElements[i, 1])
        xcoord[0] = nodes[node1 - 1, 0] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
        ycoord[0] = nodes[node1 - 1, 1] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
        xcoord[1] = nodes[node2 - 1, 0] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
        ycoord[1] = nodes[node2 - 1, 1] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor

        # Plot each line segment
        plt.plot(xcoord, ycoord, 'ks-')

    # Polish and hold plot
    plt.title("MDOF Dynamics (Time: %.1f sec, Max displacement: %.4f mm)" % (time, 1000*uMax))
    plt.axis([xmin, xmax, ymin, ymax])
    plt.pause(0.0001)

# Hold the plot a few seconds at the very end
print('\n'"Press any key to continue...")
plt.waitforbuttonpress()

Hello, thank you for this code, i want to help me please.
NameError: name 'fcns' is not defined
How can i fix it @terjehaukaas

Dear Professor Haukaas

Thank you very much for your code!
I'm facing a problem when I try to run the code.
The third figure doesn't plot a displaced shape of the frame and also glitches a long time while plotting.
The title also prints a total duration of the ground motion instead of the corresponding time of maximum displacement (as in your screenshot).
Can you suggest where can be a mistake?