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

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…

Driveby contribution to Python Cryptography

While at PyConAU 2016 I attended the Monday sprints and spent some time looking at a proposed feature I hoped would soon be part of cryptography. As most readers of this blog will know, cryptography is a very respected project within the Python ecosystem and it was an interesting experience to see how such a prominent open source project handles contributions and reviews.

The feature in question is the Diffie-Hellman Key Exchange algorithm used in many cryptography applications. Diffie-Helman Key Exchange is a way of generating a shared secret between two parties where the secret can't be determined by an eavesdropper observing the communication. DHE is extremely common - it is one of the primary methods used to provide "perfect forward secrecy" every time you initiate a TLS connection to an HTTPS website. Mathematically it is extremely elegant and the inventors were the recipients of the 2015 Turing award.

I wanted to write about this particular contribution because man…

Python, Virtualenv and Docker

Unsurprisingly I use some very popular Scientific Python packages like Numpy, Scipy and Scikit Learn. These packages don't get on that well with virtualenv and pip as they take a lot of external dependencies to build. These dependencies can be optional libraries like libblas and libatlas which if present will make Numpy run faster, or required dependencies like a fortran compiler.

Back in the good old days you wouldn't pin all your dependency versions down and you'd end up with a precarious mix of apt-get installed and pip installed packages. Working with other developers, especially on different operating system update schedules could be a pain. It was time to update your project when it breaks because of a dependency upgraded by the operating system.

Does virtualenv fully solve this? No, not when you have hard requirements on the binaries that must be installed at a system level.

Docker being at a lower level gives you much more control without adding too much extra comp…