Skip to main content

Analysis of a 3 degree of freedom building in an earthquake with scipy

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 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 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.

Popular posts from this blog

Matplotlib in Django

The official django tutorial is very good, it stops short of displaying
data with matplotlib - which could be very handy for dsp or automated
testing. This is an extension to the tutorial. So first you must do the
official tutorial!
Complete the tutorial (as of writing this up to part 4).

Adding an image to a view

To start with we will take a static image from the hard drive and
display it on the polls index page.
Usually if it really is a static image this would be managed by the
webserver eg apache. For introduction purposes we will get django to
serve the static image. To do this we first need to change the

Change the template
At the moment poll_list.html probably looks something like this:

<h1>Django test app - Polls</h1> {% if object_list %} <ul> {% for object in object_list %} <li><a href="/polls/{{}}">{{ object.question }}</a></li> {% endfor %} </ul> {% else %} <p>No polls are available.</p> …

My setup for downloading & streaming movies and tv

I recently signed up for Netflix and am retiring my headless home media pc. This blog will have to serve as its obituary. The box spent about half of its life running FreeNAS, and half running Archlinux. I’ll briefly talk about my experience with FreeNAS, the migration, and then I’ll get to the robust setup I ended up with.

The machine itself cost around $1000 in 2014. Powered by an AMD A4-7300 3.8GHz cpu with 8GB of memory. A SilverStone DS380 case is both functional, quiet and looks great. The hard drives have been updated over the last two years until it had a full compliment of 6 WD Green 4TiB drives - all spinning bits of metal though.

Initially I had the BSD based FreeNAS operating system installed. I had a single hard drive in its own ZFS pool for TV and Movies, and a second ZFS pool comprised of 5 hard drives for documents and photos.

FreeNAS is straight forward to use and setup, provided you only want to do things supported out of the box or by plugins. Each plugin is install…

Python and Gmail with IMAP

Today I had to automatically access my Gmail inbox from Python. I needed the ability to get an unread email count, the subjects of those unread emails and then download them. I found a library on sourceforge, but it actually opened the normal gmail webpage and site scraped the info. I wanted something much faster, luckily gmail can now be accessed with both pop and imap.

After a tiny amount of research I decided imap was the better albiet slightly more difficult protocol. Enabling imap in gmail is straight forward, it was under labs.

The address for gmail's imap server is:

Python has a library module called imaplib, we will make heavy use of that to access our emails. I'm going to assume that we have already defined two globals - username and password. To connect and login to the gmail server and select the inbox we can do: