In [1]:
import numpy as np
import matplotlib.pyplot as plt
plt.style.use(['fivethirtyeight', './00_mplrc'])
%matplotlib inline
from IPython.display import HTML
HTML(open("00_custom.css", "r").read())
Out[1]:
Our data is as follows
In [2]:
E = 30E9
R = 1.20
r = 0.90
m = 300E3
m1 = m
m2 = m+m
L = 45.0
We have to compute the moment of inertia of the annular beam and the flexibility coefficents ($x_1$ is at $L/2$, $x_2$ at $L$).
In [3]:
J = np.pi*(R**4-r**4)/4
EJ = E*J
f11 = (L/2)**3/(3*EJ)
f22 = (L/1)**3/(3*EJ)
In [4]:
f12 = (5*L**3) / (48*EJ)
We can construct the flexibility matrix and, by inversion, the stiffness matrix of our system, and eventually the mass matrix too.
In [5]:
F = np.array(((f11,f12),(f12,f22)))
F
Out[5]:
In [6]:
K = np.linalg.inv(F)
K
Out[6]:
In [7]:
M = np.array(((m1,0),(0,m2)))
M
Out[7]:
Next, we have to compute the eigenvalues, our system has 2 DOF so we can write explicitely the canonical equation,
(k11 - w m1) (k22 - w m2) - k12^2 = 0
w^2 m1 m2 - (k11 m2 + m1 k22) w + (k11 k22 - k12 k21) = 0
We write the coefficients of the algebraic equation and the corresponding solution.
In [8]:
k11, k12, k21, k22 = K.flatten()
a = m1*m2 ; b = - k11*m2 - k22*m1 ; c = k11*k22-k12**2
print a, b, c
print 1.0, b/a, c/a
d = np.sqrt(b*b-4*a*c)
L1, L2 = (-b-d)/2/a, (-b+d)/2/a
print L1, L2
print 2*np.pi/np.sqrt(L1), 2*np.pi/np.sqrt(L2)
Having the eigenvalues, we can compute the eigenvectors... first we compute the second component from the first equation, then we normalize with respect to the mass matrix.
In [9]:
evc = np.zeros((2,2))
evl = np.array((L1,L2))
period = 2*np.pi/np.sqrt(evl)
for i, Lambda in enumerate(evl):
#A = K - Lambda*M
#print np.linalg.det(A)
psi1 = (-1)**i
psi2 = -psi1*(k11-Lambda*m1)/k12
print psi1, psi2
f = np.sqrt(psi1**2*m1+psi2**2*m2)
evc[:,i] = psi1/f, psi2/f
print evl
print period
print evc
We define the force modulus, the time step according to the text of the exercise and the varios constants used by the linear acceleration algorithm
In [10]:
def force(t):
return 80000.0*np.sin(t) if t <= 2*np.pi else 0.0
load = np.array((1.,0.))
h = period[1]/6.0
ks = K + 6*M/h/h
fs = np.linalg.inv(ks)
ma = 3*M
mv = 6*M/h
mi = np.linalg.inv(M)
T, X, V, A = [], [], [], []
The initial values of our system
In [11]:
t0 = 0.0
x0 = np.zeros((2,))
v0 = np.zeros((2,))
p0 = load*force(0)
The numerical integration, sensu stricto,
In [12]:
while t0 < 4*np.pi:
a0 = np.dot(mi, p0 - np.dot(K,x0))
for l, i in ((T, t0), (X, x0), (V, v0), (A, a0)): l.append(i)
t1 = t0 + h
p1 = load*force(t1)
dp = p1 - p0 + np.dot(ma, a0) + np.dot(mv, v0)
dx = np.dot(fs, dp)
dv = dx*3/h - 3*v0 - a0*h/2
x1 = x0 + dx
v1 = v0 + dv
t0, x0, v0, p0 = t1, x1, v1, p1
T, X, V, A = map(np.array, (T, X, V, A))
In [13]:
%matplotlib inline
plt.figure(figsize=(10,4))
plt.plot(T/np.pi,X[:,1])
plt.xlabel('$t/\\pi s$')
plt.ylabel('$x_2/m$')
plt.title('System response, larger mass deflections')
plt.grid(1)
plt.savefig('numerical.pdf')
To check our results, we use the analytical solution for the forced response phase, the expressions should be obvious except maybe for the slicing tricks used to build the modal response array...
In [14]:
pstar = evc.T.dot(load)*80000.0
C = pstar/(evl-1.0)
A = -C/np.sqrt(evl)
t = np.linspace(0, np.pi*2, 101)
q = A[:,None]*np.sin(np.sqrt(evl)[:,None]*t) + C[:,None]*np.sin(t)
x = evc.dot(q)
plt.plot(t/np.pi,x[1]);