Computational Seismology
Finite Differences - Seismometer Equation

This notebook is part of the supplementary material to Computational Seismology: A Practical Introduction, Oxford University Press, 2016.

Authors:¶

This exercise covers the following aspects:

• Solving seismometer equation with finite difference method
• Getting familiar with seismometer response function

Basic Equations¶

Please refer to the Exercise 4.19 from the book.

The seismometer equation is given by $$\ddot{x} + 2\epsilon \dot{x} + \omega^2_0 x = - \ddot{u}$$

Where,

$\epsilon$ is the damping parameter,

$\omega_0$ is the eigenfrequency,

$\ddot{u}(t)$ is the ground motion by which the seismometer is excited, and

$x(t)$ is the motion of the seismometer mass.

We replace the time derivative by centered finite-differentiation $$\dot{x} \ \approx \ \frac{x (t + \mathrm{d}t) - x ( t- \mathrm{d}t)} {2\mathrm{d}t}$$

$$\ddot{x} \ \approx \ \frac{x( t+ \mathrm{d}t)-2x(t) + x( t- \mathrm{d}t)} {\mathrm{d}t^2}$$

Now, solving for $x(t+\mathrm{d}t)$ the extrapolation scheme is

$$x(t+\mathrm{d}t) = \frac{ - \ddot{u}(t) \mathrm{d}t^2 + (2-\omega^2_0 \mathrm{d}t^2) x(t) + (\epsilon \mathrm{d}t-1) x(t-\mathrm{d}t)} {(1+\epsilon \mathrm{d}t)}$$

Exercise¶

Part 1

While running the following cells frequency of forcing (the frequency of ground motion) and the damping parameter will be asked to enter. First try using undamped seismometer (i.e. h = 0) for some forcing frequency (eg. 0.1 Hz, 1 Hz, 2Hz, 3Hz, 4Hz, 5Hz, etc.) and interpret the results.

Part 2

Now try frequency of forcing of your choice (eg. 1 HZ) and try to search for suitable damping parameter (h).

Message: Once you become familiar with all the codes below you can go to the Cell tab on the toolbar and click Run All.

In [ ]:
#Configuration step (Please run it before the simulation code!)

import numpy as np
import matplotlib
# Show Plot in The Notebook
matplotlib.use("nbagg")
import matplotlib.pyplot as plt
matplotlib.rcParams['figure.facecolor'] = 'w'                 # remove grey background

In [ ]:
#Initialization of parameters

f0 = 1.                                                       # eigenfrequency of seismometer (hertz)
w = 2. * np.pi * f0                                           # in radians per second
dt = .01                                                      # time increment for numerical scheme
isnap = 2                                                     # snapshot frequency for visualization

# Central frequency of forcing or ground motion (will be asked)
fu0 = 1.0
# Uncomment for interactivity.
# fu0 = float(input('Give frequency of forcing (e.g. f=1 Hz) : '))

# Damping parameter of seismometer (will be asked)
h = 0.5
# Uncomment for interactivity.
# h = float(input('Give damping parameter (e.g. h=0.5) : '))

# Initialize ground motion

# Initialize parameters for ground motion
p = 1. / fu0                                                  # period
nts = int(2. * p / dt)                                        # time steps
uii = np.zeros(nts)                                           # ground motion
t0 = p / dt                                                   # time (S)
a = 4. / p                                                    # half-width (so called sigma)

# we use derivative of a Gaussian as our ground motion
for it in range(nts):
t = (it - t0) * dt
uii[it] = -2 * a * t * np.exp(-(a * t) ** 2)

nt = int(round(5. * 1. / fu0 / dt))                           # total number of time steps
src = np.zeros(nt)                                            # initialize an source array of zeros
src[0:len(uii)] = uii                                         # make derivative of gaussian as source

# End initialization of ground motion

# Initial conditions
eps = h * w                                                   # damping factor
x = np.zeros(nt)
xnow = 0
xold = 0
x_vector = np.zeros(nt)

In [ ]:
# Extrapolation scheme and the plots

# Initialization of plots

# lines1 will plot the seismometer response and lines2 for source function

lines1 = plt.plot(np.dot((np.arange(1, nt+1)), dt), x_vector[0:nt] /
np.max(np.abs(x[0:nt])), color = "red", lw = 1.5, label="Seismometer response")
lines2 = plt.plot(np.zeros(nt), color = 'blue',lw = 1.5,  label="Ground Acceleration")

plt.title("At rest")
plt.axis([0, nt*dt, -1, 1])
plt.xlabel("Time (s)")
plt.ylabel("Displacement")
plt.legend(loc="upper right")
plt.ion()
plt.show()

# Begin extrapolation and update the plot

# Extrapolation scheme (Centered finite-difference)

for i in range(nt):
if i == 0:
xold = xnow

xnew = (-src[i] * dt ** 2 + (2 - w ** 2 * dt ** 2) * xnow + (eps * dt - 1) * xold) / (1 + eps * dt)

xold = xnow                                               # for next loop present will be past
xnow = xnew                                               # for next step future will be present
x[i] = xnow
x_vector[i] = x[i]

# Updating the plots
if not i % isnap:
for l in lines1:
l.remove()
del l
for k in lines2:
k.remove()
del k

lines1  = plt.plot(np.dot((np.arange (1, i+1)), dt), x_vector[0:i] /
np.max(np.abs (x[0:nt])),color = "red",lw = 1.5, label="Seismometer response")
lines2 = plt.plot(np.dot((np.arange (1, i+1)), dt), src[0:i] /
np.max(src[0:nt]), color = 'blue',lw = 1.5, label="Ground Acceleration")
plt.title("F0 = 1Hz, SRC = %.2f Hz, h = %.2f " % (fu0, h))
plt.gcf().canvas.draw()

plt.ioff()
plt.show()