This is my building, I want to know how much each floor moves in an earthquake... Also it would be nice to know the acceleration experienced by each floor.

First off we need some data. I have an earthquake data file containing ground acceleration data from the Kobe earthquake of 1995. It is a matlab file, so the data will have to be extracted into a numpy format. Scipy has an

Firstly I'll make a small helper function that enforces complex symmetry, useful for after the frequency domain analysis:

First off we need some data. I have an earthquake data file containing ground acceleration data from the Kobe earthquake of 1995. It is a matlab file, so the data will have to be extracted into a numpy format. Scipy has an

*io*module which contains a*matlab*submodule. Since we want to visualize this data somehow, the*pylab*package will also be used.Firstly I'll make a small helper function that enforces complex symmetry, useful for after the frequency domain analysis:

from scipy.io import matlab as mio from pylab import plot, show, ylabel, xlabel, title, figure, legend, annotate import numpy as np

def enforce_complex_sym(array, nyquist): '''Enforce complex symmetry''' array[:,nyquist+1:] = np.conj(array[:,nyquist-1:0:-1])Next I'll write two plotting functions, that will annotate the maximum of some time indexed data:

def analyse_time_data(t, a, plot_title='', variable='', plot_label='', yunits='m/s^{2}'): def plot_responses(t, arrayN, var, unit, N=3):

So with all the setup done, we can get the data from the matlab file:

# Load earthquake data from matlab file kobe_data = mio.loadmat('Kobe.mat') a, t, dt = [kobe_data[index][0] for index in ['f', 't', 'dt']]

So this opens the matlab file "Kobe.mat" which contains matrices 'f', 't', and 'dt'. Once the file is open we can load the data into numpy arrays a, t and dt. At this point I used my analyse_time_data function to plot the raw ground acceleration.

Now we have the data, and it looks very earthquake like, lets create our system:

# data on each story: m = 10000 k = 1600000 c = 13000 # Construct System matricies: M = np.matrix(np.diag(3*[m])) C = np.matrix( [[2*c, -c, 0], [ -c, 2*c, -c], [ 0, -c, c]]) K = np.matrix( [[2*k, -k, 0], [-k, 2*k, -k], [0, -k, k]])

M, C and K are all (3 x 3) matrices, F is the force on each floor, N is the number of sample points:

F = ( -M*np.matrix(np.ones(3)).transpose() ) * a N = F.shape[1] nyquist = (N/2)

Now we could numerically integrate this, or we could jump into the frequency domain, I like my fft's so lets do that. We then iterate over the frequency domain data applying a transfer function.

# Transform into frequency domain Fs = np.fft.fft(F, axis=1) w = np.array(2 * np.pi * np.fft.fftfreq(N, dt)) # Calculate and apply the transfer function upto nyquist frequency Vs, dVs, ddVs = [np.zeros([3,N],np.complex) for i in range(3)] for i in xrange(nyquist): Gs = ((K - (w[i])**2 *M) + 1j*w[i]*C).I Vs[:,i] = (Gs * np.c_[Fs[:,i]]).reshape(3) dVs[:,i] = Vs[:,i] * 1j * w[i] ddVs[:,i] = Vs[:,i] * w[i]**2 [enforce_complex_sym(array, nyquist) for array in [Vs, dVs, ddVs]]

The reason we only go up to the nyquist frequency is because the data is all contained in the first half of the complex signal, we simply enforce complex symmetry then convert back to the time domain:

# Convert back to the time domain vt, dvt, ddvt = [np.fft.ifft(array, axis=1).real for array in [Vs, dVs, ddVs]] # Plot the response plot_responses(t, vt, 'Displacement', 'm') plot_responses(t, dvt, 'Velocity', 'm/s') plot_responses(t, ddvt, 'Acceleration', 'm/s^{2}')

And now for some plots:

It can be seen that the max displacement for floor 1 is about half that of the 3rd floor. The maximum the top story ends up moving is 0.33m.

Acceleration shows a similar pattern, the top story suffers the worst, reaching over 1G laterally. Possibly in a future post I will do the same analysis the numerical integration way.