In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from scipy import *
import matplotlib.pylab as pl
from matplotlib import rcParams
import json ; s = json.load(open("mplrc.json")) ; del json
rcParams.update(s)
rcParams['figure.figsize'] = 10,4
In [2]:
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
Out[2]:
The system definition
In [3]:
m = 1200. ; mass = m
k = 800000.
T = 0.2500
Compute the other characteristics of the system, in particular the damping ratio and the damping coefficient.
In [4]:
wd = 2*pi/T
wd2 = wd*wd
wn2 = k/m
wn = sqrt(wn2)
z = sqrt(1-wd2/wn2)
is_it_T = 2*pi/(wn*sqrt(1.0-z*z))
print "z = %6.3f%%"%(z*100,)
print "T =", is_it_T
damp = 2*z*sqrt(k*m)
The characteristics of the loading
In [5]:
n = 6.25
wf = n*pi/T
po = 80000.0
We compute the particular integral, $\xi = S \sin \omega_ft + C \cos \omega_f t$, determining the costants $S$ and $C$
In [6]:
M = matrix(( (wn**2-wf**2, -2*z*wn*wf),
( +2*z*wn*wf, wn**2-wf**2) ))
print M
sc = matrix((wn*wn*po/k, 0.0)).T
print sc
SC = M.I*sc
S = SC[0,0]
C = SC[1,0]
We define a set of functions that return the values of the particular integral and the values of its time derivatives
In [7]:
def r0(t): return S*sin(wf*t)+C*cos(wf*t)
def r1(t): return wf*(S*cos(wf*t)-C*sin(wf*t))
def r2(t): return -wf**2*r0(t)
Just to be sure, we plot the difference between the load and the value that we obtain substituting the particular integral in the left member of the eq. of dynamic equilibrium.
In [8]:
def error(t): return mass*r2(t)+damp*r1(t)+k*r0(t)-p(t)
t = linspace(0,T/n,101)
def p(t): return po*sin(wf*t)
pl.plot(t,error(t));pl.xlim(0,.04);pl.show()
pl.plot(t,p(t));pl.xlim(0,.04);
Having $\xi(t)$ and $\dot\xi(t)$ we can determine the constants of integration in the homogeneous integral imposing the respect of the initial rest conditions.
In [9]:
# exp(-z wn t) * (A sin(wd t) + B cos(wd t))
# for t=0: 1*(A*0+B*1) = B
# wd exp() * (A cos(wd t) - B sin(wd t) - z wn exp() * (A sin(wd t) + B cos(wd t))
# for t=0: wd*1*(A*1-B*0) -z*wn*1*(A*0+B*1) = wd A - z wn B
B = - r0(0)
A = (z*wn*B-r1(0))/wd
def xforced(t,A=A,B=B): return exp(-z*wn*t)*(A*sin(wd*t)+B*cos(wd*t))+r0(t)
def vforced(t,A=A,B=B): return -z*wn*exp(-z*wn*t)*(A*sin(wd*t)+B*cos(wd*t)) + wd*exp(-z*wn*t)*(-B*sin(wd*t)+A*cos(wd*t)) + r1(t)
Let's plot the forced response in terms of displacements and velocities
In [10]:
pl.plot(t, xforced(t), label = r'$x(t)$')
pl.plot(t, vforced(t)/100, label = r'$v(t)/100$')
pl.legend(loc=0)
pl.xlim(0,.04);
Let's see the values of the response at the beginning and at the end of the external loading
In [11]:
print xforced(0), xforced(0.04)
print vforced(0), vforced(0.04)
For an undamped system (our system is heavily damped) with maximum after half-sine loading, the max load is given by $$x_0 = \Delta_s \frac{2\beta}{\beta^2-1}\cos\left(\frac{\pi}{2\beta}\right)$$
In [12]:
beta = T/2/0.04
print beta
print 80000/k*2*beta/(beta**2-1)*cos(pi/2/beta)
Use a function to compute the constants of integration of the h. response for $t>t_0$
In [13]:
def coeff(T,x,v):
M = matrix((
( sin(wd*T), cos(wd*T) ),
( -z*wn*sin(wd*T) + wd*cos(wd*T), -z*wn*cos(wd*T) - wd*sin(wd*T))
))*exp(-z*wn*T)
xvT = matrix(((x(T),),
(v(T),)
))
AB = M.I*xvT
return AB[0,0], AB[1,0]
Use the function above to determine the constants of integration and define
the free response and a (vector) function that returns either the forced or the free response according to the value of t
In [14]:
Af, Bf = coeff(T/n,xforced,vforced)
def xfree(t):
return exp(-z*wn*t)*(Af*sin(wd*t)+Bf*cos(wd*t))
def x12(t):
return where(t<T/n, xforced(t), xfree(t))
Now we plot the response in terms of displacements and forces
In [15]:
time = linspace(0.0, T, 4001)
pl.plot(time, x12(time));
pl.show()
pl.plot(time, k*x12(time));
In [16]:
# 1
h = 0.001
h2 = h*h
# 2
stop = 0.500+h/2
# 3
k_ = k+2*damp/h+4*mass/h2
a0_coeff = 2*mass
v0_coeff = 4*mass/h+damp+damp
# 4
def p(t): return where(t<0.04, 80000*sin(pi*t/0.04), 0)
# 5
t0 =0 ; x0 = 0 ; v0 = 0 ; p0 = 0 ; a0 = 0
# 6
t= [] ; X = [] ; V = [] ; A = []
The actual loop that iterates the system state according to the constant acceleration algorithm
In [17]:
while t0<stop:
t.append(t0) ; X.append(x0) ; V.append(v0) ; A.append(a0)
t1 = t0+h
p1 = p(t1)
dp = p1 - p0 + a0_coeff*a0 + v0_coeff*v0
dx = dp/k_
dv = 2*(dx/h-v0)
t0 = t1
p0 = p1
x0 = x0+dx
v0 = v0+dv
a0 = (p0 - k*x0 - damp*v0)/mass
Lets plot our results near the maximum, in terms of spring force, against the exact response we previously found.
In [18]:
pl.plot(t,[k*d for d in X], label='Numerical')
time = linspace(0,0.5, 5001)
pl.plot(time, k*x12(time),label='Analytical')
pl.xlim(0,0.5)
pl.legend(loc=0)
if 1:
pl.ylim(35000,37500)
pl.xlim(0.06,0.09)
Repeat using the linear acceleration algorithm
In [19]:
h = 0.001
h2 = h*h
stop = 0.500+h/2
k_ = k + 3*damp/h + 6*mass/h2
a0_coeff = 3*mass + damp*h/2
v0_coeff = 6*mass/h + 3*damp
dx2dv = 3.0/h ; v2dv = 3.0 ; a2dv = h/2.0
def p(t): return where(t<0.04, 80000*sin(pi*t/0.04), 0)
t0 =0 ; x0 = 0 ; v0 = 0 ; p0 = 0 ; a0 = 0
t= [] ; X = [] ; V = [] ; A = []
while t0<stop:
t.append(t0) ; X.append(x0) ; V.append(v0) ; A.append(a0)
t1 = t0+h
p1 = p(t1)
dp = p1 - p0 + a0_coeff*a0 + v0_coeff*v0
dx = dp/k_
dv = dx*dx2dv - v0*v2dv - a0*a2dv
t0 = t1
p0 = p1
x0 = x0+dx
v0 = v0+dv
a0 = (p0 - k*x0 - damp*v0)/mass
and plot the exact force response and numerical force response
In [20]:
pl.plot(t,[k*d for d in X], label='Numerical')
time = linspace(0,0.5, 5001)
pl.plot(time, k*x12(time),label='Analytical')
pl.xlim(0,0.5)
pl.legend(loc=0)
if 1:
pl.ylim(35000,37500)
pl.xlim(0.06,0.09)
In [21]:
pl.plot(t,X)
pl.xlabel('Time/s')
pl.ylabel('x(t)/m')
pl.hlines(0,0,0.5,alpha=0.2)
pl.xlim(0,0.5)
pl.title('Linear acceleration algorithm')
pl.show()
time = linspace(0,0.5, 5001)
pl.plot(time, x12(time))
pl.xlabel('Time/s')
pl.ylabel('x(t)/m')
pl.hlines(0,0,0.5,alpha=0.2)
pl.xlim(0,0.5)
pl.title('Exact integration');
In [21]: