Configure the plotting machinery
In [6]:
%pylab inline
%config InlineBackend.figure_format = 'svg'
import json
s = json.load( open("mplrc.json") )
matplotlib.rcParams.update(s)
matplotlib.rcParams['figure.figsize'] = 9,4
black="#404060" # plots containing "real black" elements look artificial
Configure the rendering of this notebook with CSS
In [7]:
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
Out[7]:
In [36]:
subplot(1,2,1);plot((0,1,3),(0,40,0));xticks((0,1,3)); yticks((0,40))
xlabel('t/s');ylabel('p(t)/kN');grid();
subplot(1,2,2);plot((0,1,3,8),(0,40,0,0));xticks((0,1,3,8)); yticks((0,40))
xlabel('t/s');ylabel('p(t)/kN');grid();
figsize(9.0,0.1);
We are to replicate the solution obtained for a triangular load using the Duhamel Integral.
Our load is 3s long, but we need a stretch of zeros to damp out the response and simulate rest initial conditions. We add zeros to the end of the load up to a total duration of 8s. The period of our loading is hence, rather arbitrarily, $T=8\,{}$s.
In our exercise, we are free to choose the number of samples per second, so we chose 512 sps.
How many samples are there? $N = 8\times512=4096$, note that $N$ is a power of 2.
In [9]:
T = 8
sps = 512
print "The Nyquist frequency is ", sps/2.0,"Hz"
The array t contains the times at which our signal was sampled, the load p is computed starting with an empty list and appending a new value for each element in the vector t
In [10]:
t = arange(0,T,1./sps)
p = []
for time in t:
if time<1.0:
p.append(time*40000)
elif time<3.0:
p.append(60000-20000*time)
else:
p.append(0.0)
Am I sure that the list p contains the values of the loading?
Let's try to plot p vs t...
In [11]:
matplotlib.rcParams['figure.figsize'] = 9,4
plot(t, p, black) ; xlabel("t/s") ; ylabel("p(t)/N")
ylim((-5000,45000))
print len(t), len(p)
It looks OK...
Now, the fast Fourier transform of the sequence p is computed, and given a name, P.
It is customary to denote Fourier pairs by the same letter, the small letter for the time domain representation and the capital letter for the frequency domain representation.
In [12]:
P = fft.fft(p)
iP = fft.ifft(P)
I have computed also the inverse FFT of the FFT of the loading, naming it iP, it is a sequence of complex numbers and here we plot the real and the imaginary part of each component versus time.
In [13]:
plot(t,real(iP),black,t,1*imag(iP),'y') ; xlabel("t/s") ; ylabel("p(t)/N") ;
It seems OK...
Next, we use a convenience function to compute a sequence of frequencies (in Hertz!) associated with the components of P, the FFT of p. The parameters are the number of points and the sampling interval..
Note that the sequence of frequencies has a discontinuity when the Nyquist frequency is reached, i.e., the next frequency is the most negative one.
In [14]:
f = fft.fftfreq(T*sps, 1./sps)
f, f[2046:2051]
Out[14]:
The x axis is streching over the interval $-f_\text{Ny}$, $+f_\text{Ny}$
In [15]:
plot(f,real(P), black, f,imag(P),'b')
xlim(-256,256) ; xticks(range(-256,257,64))
xlabel("f/Hz") ; ylabel("P(f)") ;
The plot above is not much clear, because the frequency components are significantly different from zero only in a narrow range of frequencies around the origin.
In the next 3 plots we zoom near the origin of the frequency axis to have a bit more of detail. There are 3 plots, first the absolute value of P vs f, then the real part and finally the imaginary part of P, versus f.
In [16]:
plot(f,abs(P), black)
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("abs(P(f))")
xlim(-4,4) ; ylim(-0.2E7,3.3E7);
Not afwully nice, this last plot...the baseline and the missing line at the left of the zero are artifacts, due do the particular sequence with which the positive and negative frequencies are arranged in the DFT output.
To obviate these problems we can use the function fftshift, that reorders (shifts)
the elements in an array such that the sequence goes from the most negative
frequency to the most positive.
In [17]:
plot(fftshift(f),fftshift(abs(P)), black)
axhline(0,color=black, linewidth=0.25)
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("abs(P(f))")
xlim(-4,4) ; ylim(-0.2E7,3.3E7);
and now the other two plots I promised,
In [18]:
plot(fftshift(f),fftshift(real(P)), black)
axhline(0,color=black, linewidth=0.25)
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("real(P(f))")
xlim(-4,4) ; ylim(-3.3E7,3.3E7);
In [19]:
plot(fftshift(f),fftshift(imag(P)), black)
axhline(0,color=black, linewidth=0.25)
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("imag(P(f))")
xlim(-4,4) ; ylim(-3.3E7,3.3E7);
Until now, we did without the SDOF, now it's time to describe it and derive its response function.
All the parameters are the same as in the excel example, we compute k because we need it to normalize the response.
In [29]:
z = 0.1; fn = 1/0.6 ; m =6E5 ; wn = fn*2*pi ; k = m*wn**2
def H(f):
b = f/fn
return 1./((1-b*b)+1j*(2*z*b))
As usual, we plot the response function, or rather the absolute value of, against a short span of the frequency axis, centered about the origin, to show the details of the response function itself.
In [30]:
plot(fftshift(f),fftshift(abs(H(f)))) ; xlabel("f/Hz") ; xlim(-8,8) ; ylabel("H(f)") ;
The FFT of the response is computed multiplying, term by term, P by the response or transfer function, then we compute the IFFT of X to obtain x, the time domain representation of the response domain.
In [31]:
X = [ P_*H(f_) for P_, f_ in zip(P,f)]
x = ifft(X)/k
Note that the response function is periodic with period $T=8\,{}$s.
In the end, we remain with the task of plotting the response function, that is the real part of x. Just to be certain we plot also the imaginary part of x, so we can be sure that it is negligible with respect to the real part
In [32]:
plot(t,1000*real(x))
xlabel("t/s"); ylabel(r"$\Re(x)/$mm");
axhline(0,color=black, linewidth=0.25)
show()
In [33]:
plot(t,1E18*imag(x),linewidth=0.33)
axhline(0,color=black, linewidth=0.25)
xlabel("t/s"); ylabel(r"$\Im(x)/$am");
show()
The importance of the zero trail to adjust for initial rest condition cannot be underestimated. The length required depends, of course, on how much damped our system is, the lesser the damping, the longer the time required to damp out the response.
Lets try to see what happens if we go from $\zeta=0.10$ to $\zeta=0.01$:
In [34]:
z = 0.01
X = [ P_*H(f_) for P_, f_ in zip(P,f)]
x = ifft(X)/kfunction
plot(t,1000*real(x))
xlabel("t/s"); ylabel(r"$\Re(x)/$mm");
axhline(0,color=black, linewidth=0.25)
show()
In [ ]: