Building response to earthquake

Solving ordinary differential equation of dynamic motion involving mass, damping and spring in response to earthquake

A building, considered here as a rigid mass, is subject to an earthquake. The foundation of the building, interacting with the ground, is represented by a spring and a dashpot (damper).

The dynamic motion of the building is defined by the following differential equation:

\[ k \left(u_b(t)-u_g(t) \right) + c \left(\frac{d u_b}{d t} - v_g(t) \right) + m \left(\frac{d^2 u_b}{d t^2} \right) = F \]

where:

• \(u_b\) = horizontal displacement of the building

• \(u_g\) = horizontal displacement of the ground (= earthquake outcrop motion)

• \(v_g\) = \(du_g / dt\) = horizontal velocity of the ground

• \(t\) = time

• \(k\) = (shear) stiffness between the building foundation and the ground

• \(c\) = damping between the building foundation and the ground

• \(m\) = mass of the building

• \(F\) = external force acting on the building

In the remainder of this assignment we consider \(F\) = 0

Some initialisation:

import numpy as np
import matplotlib.pyplot as plt

k = 10  # Stiffness
c = 10  # Damping
m = 10  # Mass

Part 1: Reading and interpretion of the earthquake record

The earthquake outcrop motion is given in terms of acceleration \(a_g\), velocity \(v_g\) and displacement \(u_g\) as time series at regular time intervals (0.005 seconds). In the code cell below the earthquake data is read from the respective files and stored in corresponding arrays. Note that the unit of length in earthquake records is usually centimeter, so \(a_g\) is given in \(cm/s^2\), \(v_g\) in \(cm/s\) and \(u_g\) in \(cm\).

import os
from urllib.request import urlretrieve


def findfile(fname):
    if not os.path.isfile(fname):
        print(f"Downloading {fname}...")
        urlretrieve(
            "https://github.com/TUDelft-MUDE/source-files/raw/main/file/" + fname, fname
        )

# Download the data files (this may take some time if the server is busy)
findfile("Earthquake_acceleration.txt")
findfile("Earthquake_velocity.txt")
findfile("Earthquake_displacement.txt")

Downloading Earthquake_acceleration.txt...

Downloading Earthquake_velocity.txt...

Downloading Earthquake_displacement.txt...

from itertools import islice


def readfile(filename):
    """Read data from file and return the data as a 1-D numpy array."""
    values = []
    with open(filename) as f:
        for line in islice(f, 3, None):  # Skip first 3 rows
            values.append(np.fromstring(line, sep=" "))
    return np.concatenate(values)


ag = readfile("Earthquake_acceleration.txt")
vg = readfile("Earthquake_velocity.txt")
ug = readfile("Earthquake_displacement.txt")

num = vg.shape[0]  # number of datapoints in the earthquake record
dt = 0.005  # time interval between subsequent data points
t = np.linspace(start=0, stop=dt * num, num=num)

Plotting the earthquake ground velocity

Task 1.1

Write the Python code to plot the ground velocity \(v_g\) as a function of time \(t\) using the standard plotting facilities in Matplotlib.

Solution 1.1

plt.figure(figsize=(12, 4))
plt.plot(t, vg, linewidth=0.5)
plt.plot([0, t[-1]], [0, 0], color="black", linewidth=0.2, linestyle=":")
plt.xlim(0, t[-1])
plt.xlabel("time [s]")
plt.ylabel("velocity [cm/s]")
plt.show()

Numerical differentiation and integration

As an exercise, we consider the earthquake ground velocity \(v_g\) as the only given data, from which we calculate the ground acceleration \(a_g\) by numerical differentiation and the ground displacement \(u_g\) by numerical integration.

For the numerical differentiation we use the Forward Difference scheme:

Task 1.2

Formulate \(a_g(t)\) = \(\cfrac{d v_g(t)}{d t}\) by numerical differentiation according to the Forward Difference scheme, and write the necessary Python code.

Solution 1.2

Forward Difference scheme: $\( a_g(t) = \frac{v_g(t+\Delta t) - v_g(t)}{\Delta t} \)\( or \)\( a_g^i = \frac{v_g^{i+1} - v_g^i}{\Delta t} \)$

# Forward Difference scheme to differentiate ground acceleration from given velocity

ag_FD = np.zeros([num])

for i in range(0, num - 1):
    ag_FD[i] = (vg[i + 1] - vg[i]) / dt

For the numerical integration we use the Trapezoidal Rule:

Task 1.3

Formulate \(u_g(t)\) = \(\int{v_g(t) d t}\) by numerical integration according to the Trapezoidal Rule, and write the necessary Python code.

Solution 1.3

Trapezoidal Rule: $\( u_g(t + \Delta t) = u_g(t) + \frac{\Delta t}{2} \left[v_g(t+\Delta t) + v_g(t) \right] \)\( or \)\( u_g^{i+1} = u_g^i + \frac{\Delta t}{2} \left(v_g^{i+1} + v_g^i \right) \)$

# Trapezoidal Rule to integrate ground displacement from given velocity

ug_TR = np.zeros([num])

for i in range(0, num - 1):
    ug_TR[i + 1] = ug_TR[i] + dt / 2 * (vg[i + 1] + vg[i])

Comparing the calculated data with the original data

Task 1.4

Plot the calculated and original data in one plot to enable visual comparison. Do this for the acceleration \(a_g\) as well as the displacement \(u_g\) (in separate plots).

Solution 1.4

# Plotting ground acceleration

plt.figure(figsize=(12, 4))
plt.plot(t, ag_FD, linewidth=0.3, label="calculated")  # Calculated ground acceleration
plt.plot(t, ag, linewidth=0.3, label="original")  # Original ground acceleration
plt.plot([0, t[-1]], [0, 0], color="black", linewidth=0.2, linestyle=":")
plt.xlim(0, t[-1])
plt.xlabel("time [s]")
plt.ylabel("acceleration [cm/s2]")
plt.legend()
plt.show()

# Plotting ground displacement

plt.figure(figsize=(12, 4))
plt.plot(t, ug_TR, linewidth=0.3, label="calculated")  # Calculated ground displacement
plt.plot(t, ug, linewidth=0.3, label="original")  # Original ground displacement
plt.plot([0, t[-1]], [0, 0], color="black", linewidth=0.2, linestyle=":")
plt.xlim(0, t[-1])
plt.xlabel("time [s]")
plt.ylabel("displacement [cm]")
plt.legend()
plt.show()

Interpretation of the earthquake data

Task 1.5

Conclude on the accuracy of the calculated results.

Are the results according to your expectation? For example, comment on the acceleration, velocity and displacement at the end of the recorded time series.

Solution 1.5

Calculated results are visually the same as the original data, which means: they seem accurate.

Acceleration and velocity are expected to be zero at the end of the recorded time series. Indeed, they seem to be very small.

Displacement at the end of the time seris is not zero. This may be due to the following causes:

1. The ground surface may have moved permanently as a result of the earthquake.

2. The earthquake records were not corrected for drift.

Part 2: Numerical solution of the dynamic problem

Now, we proceed with the numerical solution of the above differential equation for the dynamic motion of the building \(u_b\). Since the second derivative already requires two subsequent timesteps, we use the Central Difference scheme both for the first and second derivative to obtain maximum accuracy. Note that the ground acceleration is not used in the differential equation and the numerical solution; only the ground displacement and velocity! (Why?)

Task 2.1

Formulate the differential equation in a numerical scheme using the Central Difference scheme for the first and second derivative and by considering \(u_b(t)\) = \(u_b^i\) as the central point.

Elaborate the numerical scheme such that the solution for the building displacement at the next time step \(u_b(t + \Delta t)\) = \(u_b^{i+1}\) is expressed explicitly in the other (known) components.

Hint: Refer to next week’s lecture or Section 1.6 in the book for an application involving the second derivative and verify that: $\( \frac{d^2 u_b}{d t^2} \approx \frac{u_b^{i+1} - 2 u_b^i + u_b^{i-1}}{\Delta t^2} \)$

Solution 2.1

\[ k u_b^i + c \frac{u_b^{i+1} - u_b^{i-1}}{2 \Delta t} + m \frac{u_b^{i+1} - 2 u_b^i + u_b^{i-1}}{\Delta t^2} = k u_g^i + c v_g^i \]

\[ \left( \frac{c}{2 \Delta t} + \frac{m}{\Delta t^2} \right) u_b^{i+1} = \left( -k + \frac{2 m}{\Delta t^2} \right) u_b^i + \left( \frac{c}{2 \Delta t} - \frac{m}{\Delta t^2} \right) u_b^{i-1} + k u_g^i + c v_g^i \]

\[ u_b^{i+1} = \left[ \left( -k + \frac{2 m}{\Delta t^2} \right) u_b^i + \left( \frac{c}{2 \Delta t} - \frac{m}{\Delta t^2} \right) u_b^{i-1} + k u_g^i + c v_g^i \right] / \left( \frac{c}{2 \Delta t} + \frac{m}{\Delta t^2} \right) \]

Task 2.2

How many initial conditions are required to solve the differential equation?

Formulate these initial conditions for the numerical scheme.

Solution 2.2

Since the differential equation is second order, \(two\) initial conditions are required:

1. Initial building displacement is zero

2. Initial building velocity is zero

\[ u_b^0 = 0 \]

and $\( v_b^0 \approx \frac{u_b^1 - u_b^0} {\Delta t} = 0; u_b^1 = u_b^0 = 0 \)$

Implementation of the numerical solution

Task 2.3

Based on your elaboration above, write the necessary Python code to implement the numerical solution for the building displacement \(u_b\) in time. Store the results in a list or (Numpy) array.

Use a time step \(\Delta t\) equal to the time interval of the earthquake data points.

Solution 2.3

def solution(k, c, m, num):
    """Solve the ODE for the building displacement in time
        
        Parameters
        ----------
        k : float
            Stiffness of the building
        c : float
            Damping of the building
        m : float
            Mass of the building
        num : int
            Number of time steps to solve for
    """
    u = np.zeros([num])  # Initial conditions at t=0 automatically satisfied

    for i in range(1, num - 1):
        u[i + 1] = (
            (-k + 2 * m / (dt * dt)) * u[i]
            + (c / (2 * dt) - m / (dt * dt)) * u[i - 1]
            + k * ug[i]
            + c * vg[i]
        ) / (c / (2 * dt) + m / (dt * dt))
    return u


ub = solution(k, c, m, num)

As the implementation of this numerical solution is a bit trickier, you can use the code below to check if your output is correct:

# checks your output against the expected output

assert len(ub) == 7827, "number of timesteps is not correct."
assert ub[0] == 0, "Initial condition not satisfied."

assert np.allclose(ub[100], -0.09866151298075838), "Incorrect value at index 100."
assert np.allclose(ub[1000], 6.55694268229219), "Incorrect value at index 1000."
assert np.allclose(ub[-1], 1.9276415440472443), "Incorrect value at last index."

print("All tests passed!")

All tests passed!

Plotting the results

Task 2.4

Plot the results (building displacement \(u_b\)) as a function of time, together with the earthquake motion (ground displacement \(u_g\)).

Solution 2.4

plt.figure(figsize=(12, 4))
plt.plot(t, ug, color="green", linewidth=0.3, label="ground displacement")
plt.plot(t, ub, color="blue", linewidth=0.3, label="building displacement")
plt.plot([0, t[-1]], [0, 0], color="black", linewidth=0.2, linestyle=":")
plt.xlim(0, t[-1])
plt.legend()
plt.xlabel("time [s]")
plt.ylabel("displacement [cm]")
plt.show()

Task 2.5

After interpretation of the previous results, you are requested to perform similar calculations for other parameter combinations. Vary \(k\), \(c\) and \(m\) (factor 10 larger (and smaller), by editing the cell at the top of this notebook, compared to their original value) and see how they influence the results. Think about whether your observations make physical sense. If ‘oscillations’ occur, try explaining them.

Take notes on paper – you will need them later for the report.

Solution 2.5

Answers to these questions can be obtained by varying parameters \(k\), \(c\) and \(m\) in the range [1, 10, 100] by editing the cell at the top of this notebook, and plotting the results:

• Higher stiffness makes the building following the ground motion more closely. It also causes higher vibration/oscillation frequencies.

• Higher damping also makes the building follow the ground motion more closely. It also damps out vibrations/oscillations more.

• Higher mass makes the building following the ground motion more slowly. It also causes lower vibration/oscillation frequencies.

• (For the cases where \(k\) = \(c\) = \(m\) the solution is (approximately) the same, irrespective of the actual values).

• Oscillations that might be observed with some parameter combinations could, in general, have a physical or numerical origin:

  • Physical: If the earthquake signal approximates the eigenfrequency of the system (based on the specific combination of parameters), oscillations could be a result of resonance.

  • Numerical: If the time step is too large, oscillations could be a result of numerical instability (when using a Forward Difference scheme). This will be discussed in more detail in Week 1.3.

u = []
values = [1, 10, 100]

idx = 0
for k in values:
    for c in values:
        for m in values:
            u.append(solution(k, c, m, num))
            idx += 1

fig, axes = plt.subplots(9, 3, figsize=(15, 30))
axes = axes.flatten()

for idx in range(27):
    k_val = values[idx // 9]
    c_val = values[(idx // 3) % 3]
    m_val = values[idx % 3]

    axes[idx].plot(t, u[idx], linewidth=0.5)
    axes[idx].plot(t, ug, color="green", linewidth=0.3, label="ground displacement")
    axes[idx].set_xlabel("time [s]")
    axes[idx].set_ylabel("displacement [cm]")
    axes[idx].set_title(f"k={k_val}, c={c_val}, m={m_val}")
    axes[idx].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

By Ronald Brinkgreve and Anna Störiko, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.