In [1]:
from sympy import *
%load_ext sympyprinting
%config InlineBackend.figure_format = 'svg'
from IPython.core.display import display, Image, Math, Latex
rcParams['savefig.dpi'] = 120
ndof = 3
In [2]:
Image(filename='figures/structure.png')
Out[2]:
First we need a number of symbols, to be used in the following.
Note that the symbols function is rather smart when dealing with pedices and greek letters...
In [3]:
x, t, k, m, EJ, L, wn, w0, w, delta, Lambda = symbols(
'x, t, k, m, EJ, L, omega_n, omega_0, omega, delta, Lambda')
w2, w20, w2n = map(lambda xxx: xxx*xxx, (w,w0,wn))
psi = Matrix(symbols('psi_1 psi_2 psi_3'))
A, B, C = symbols('A B C')
beta = symbols('beta_n')
influence = Matrix((1, 0, 1))
Next, using some of the symbols previously defined we assing a name to an algebraic expression, so defining the ground displacement and later, by double differentiation with respect to t, the ground acceleration
In [4]:
ug = delta*sin(w*t) ; ag = diff(ug,t,2)
display(Eq(symbols('u_g'),ug), Eq(symbols('a_g'),ag))
The flexibility matrix for the 4 degrees of freedom system is
In [5]:
F0 = Matrix((( 24, 16, 9, -32),
( 16, 16, 5, -33),
( 9, 5, 4, -10),
( -32, -33, -10, 72)))/6
The flexibility matrix for the 3 DOF system can be computed using the "slice" notation, the stiffness matrix is computed inverting the flexibility and the mass matrix is simply the unit matrix. All structural matrices are normalized.
While we are at it, we derive the matrix that gives the reaction $f_4$ as a linear combination of the nodal loads.
In [6]:
F = Matrix(F0[:-1,:-1] - F0[:-1,-1:]*F0[-1:,-1:].inv()*F0[-1:,:-1])
R4_from_fs = F0[-1:,-1:].inv()*F0[-1:,:-1]
K = F.inv() ; M = eye(3)
display(Math(r'\mathbf{F} = \frac1{432k}' + latex(F*432)))
display(Math(r'\mathbf{K} = \frac{k}{209}' + latex(K*209)))
display(Math(r'\mathbf{M} = m' + latex(M)))
Now we can write the equation of free vibrations, denoting by Free the matrix of coefficients of the
equation, we can write the scalar equations
In [7]:
Free = K - Lambda * M
dummy = [display(Math(latex(line)+'=0')) for line in Free*psi]
To have non trivial solution, we need that the determinant of Free is zero
In [8]:
char_eq = Free.det().expand()
display(Eq(char_eq))
The eigenvalues are normalized with respect to $\omega_0^2 = k/m$, that is $\omega_i^2 = \Lambda_i × \omega_0^2$.
The eigenvalues are the solutions of the characteristic equation,
In [9]:
Lambdas = [(r).evalf().subs(I,0) for r in solve(char_eq)]
Lambdas.sort()
Omegas = map(sqrt, Lambdas)
Omegas = [om.evalf(n=8) for om in Omegas]
display(Math(r'\begin{align*}' +
',\\\ '.join([
r'\frac{\omega^2_%d}{k/m} &='%(i+1,)+str(Lambdas[i]) for i in range(ndof)]) +
r'.\end{align*}'))
Now we solve the equation of free vibration once for each eigenvalue, to obtain the corresponding eigenvectors.
In [10]:
Psi = eye(3)
for i, r in enumerate(Lambdas):
equ = [Eq(r,) for r in (Free.subs(Lambda, r)*psi)[1:]]
print "The linearly indipendent equations of free vibrations for i=%d"%(i+1,)
display(equ)
d = solve(equ,psi[1:]) ; d[psi[0]] = 1
unnorm = psi.subs(d)
factor2 = (unnorm.transpose()*M*unnorm)[0]
print "The non-normalized eigenvector and the modal mass for i=%d"%(i+1,)
display((unnorm, factor2))
Psi[:,i] = (unnorm/sqrt(factor2))
print "---------------------------------------------------------\n"
print "The matrix of normalized eigenvectors:"
display(Psi)
In [11]:
display(Latex(r'''The modal loads$$$$
$$\mathbf{p}^* = -\Psi^T M E a_g ='''+latex(Psi.transpose()*M*influence)+"\\;"+latex(-ag)+'$$$$$$'))
In [12]:
q123 = []
for i in range(3):
wi = symbols('omega_%d'%(i+1,))
betai = symbols('beta_%d'%(i+1,))
qi = symbols('q_%d'%(i+1,))
display(Latex("The particular integral for mode %d, \
$\\xi(t)=C\\times \\sin(\\omega t)$ and the corresponding equation of equilibrium" % (i+1,)))
xi = C*sin(w*t)
p = -(Psi[:,i].transpose()*influence*ag)[0]
display(Eq(xi*symbols('omega_%d'%(i+1,))**2+xi.diff(t,2),p))
xi = xi.subs(C, solve(xi.diff(t,2)+wn*wn*xi-p,C)[0])
# homogeneous solution, displacements and velocities
h_d = A*sin(wn*t) + B*cos(wn*t)
display(Latex(r'$$q_{%d}(t) = '%(i+1,) + latex((h_d+xi).subs(wn,wi)) + "$$"))
h_v = h_d.diff(t)
d = solve(map(lambda x: x.subs(t,0),(h_d + xi, h_v + xi.diff(t))), A,B)
display(Latex(r'''Imposing the rest initial conditions we have
\begin{align*} A & = %s, & (%s & = %s)\\ B & = %s. \end{align*}''' % (
latex(d[A].subs(w,wn*betai).simplify()), latex(betai), latex(w/wi),
latex(d[B]))))
d[w] = 2.75 ; d[w2n] = wn*wn ; d[wn] = Omegas[i]
q = (h_d+xi).subs(d).simplify()
d[t] = w0*t ; d[delta] = 1
q123.append(q.subs(d))
display(Latex(
'''Substituting the numerical values, expressing the frequencies
as multiples of $\\omega_0$ and'''))
display(Latex('''normalizing with respect to
$\\delta$ we have'''+
r'$$\frac{%s(t)}{\delta} = %s.$$'%(latex(qi),latex(q.subs(d).evalf(n=8)))))
a = display("----+----|"*7+'--') if i<ndof-1 else 1
What is the definition of $x_3$?
In [13]:
display(Math(r'x_3(t) = %s + %s + %s' % tuple(r'\psi_{3%d}\,q_{%d}(t)'%(i,i) for i in range(1,ndof+1))))
We
In [14]:
display(Psi[2,:]*Matrix(symbols('q_1 q_2 q_3')))
x3 = (Psi[2,:]*Matrix(q123))[0].simplify()
display(Math(r"\frac{x_3(t)}{\delta} = %s"%(latex(x3.evalf(n=6)),)))
The x3 expression is not suitable for numerical evaluation, as in plotting,
so we have to take an extra step to transform the symbolic expression in a numerical function.
In [15]:
x3t = lambdify(t, x3.subs(w0,1), modules='numpy')
It's time to plot our results
In [16]:
time = linspace(0.0,10.0,201)
plot(time,x3t(time))
xlabel(r'$\omega_0t$', fontsize='x-large')
ylabel(r'$\frac{x_3}{\delta}$', fontsize='xx-large', rotation=0)
axis(xmin=0, xmax=10, ymin=-4, ymax=+4)
grid()
To compute the reaction $f_4$
In [17]:
x = Psi*Matrix(q123)
fs = K*x
R4 = R4_from_fs*fs ; R4 = R4[0]
display(Math('f_4 = k\delta\;(%s)'%(latex(R4.evalf(n=5)),)))
As before, the sybolic expression must be transformed in a numerical function
In [18]:
R4t = lambdify(t, R4.subs(w0,1), modules='numpy')
In this case the resonant behaviour of the system is particularly apparent.
In [19]:
plot(time,R4t(time))
xlabel(r'$\omega_0t$', fontsize='x-large')
ylabel(r'$f_4/(k\delta)$', fontsize='x-large')
grid()
In [20]:
k = mat(((912,-684,-1482),(-684,2064,864),(-1482,864,2928)))/209.
m = mat(eye(3)*1.0)
r = mat('1;0;1')
h = 0.002 ; dur = 10.0
w = 2.75 ; w2 = w*w
def load(t):
return -r*float(-w2*sin(w*t))
A = m*3 ; V = m*6/h
K = k + m*6/h/h
KI = K.I
MI = m # m being the unit matrix
t = 0
x = mat('0;0;0')
v = mat('0;0;0')
p = load(t)
cont = []
while t<dur+h/2:
a = p-k.dot(x)
cont.append([t, x[2,0]])
t = t+h
dp = load(t)-p
dp_star = dp+A*a+V*v
dx = KI*dp_star
dv = (dx/h-v)*3-a*h/2
x, v, p = x+dx, v+dv, p+dp
cont = array(cont)
Our last task will be to plot, with a continuous line, the numerically computed response and then superimpose on the plot dots obtained by sampling the analytical response function,
In [21]:
time = linspace(0,10,101)
plot(cont[:,0], cont[:,1], '-k', label='Numerical Response',
linewidth=0.5)
plot(time, x3t(time), 'ok', label='Analytical Response',
alpha=0.5, markersize=2)
plot((0,10),(0,0),'-', color='lightgray', linewidth=1)
axis(xmin=0, xmax=10, ymin=-4, ymax=4)
title('Comparison of Numerical and Analytical Solutions')
xlabel(r'$\omega_0 t$', fontsize='x-large')
ylabel(r'$\frac{x_3}{\delta}$', fontsize='xx-large', rotation=0)
legend(fancybox=True, shadow=True, numpoints=5)
setp(gca().get_legend().get_texts(), fontsize='small')
setp(gca().get_legend().get_frame(), color='#FFF040')
setp(gca().get_legend().get_frame(), edgecolor='black')
grid()
The agreement between the two solutions is indeed very good.
In [21]: