• Introduction:
• Section I: The animation of the first mode
• Section 2: The animation of the second mode
  • Section 3: The animation of the third mode
  • Section 4: The animation of a subway car overcoming a bump
• Conclusion:

Modelling and Animation of Subway Cars

by J. Angeles, McGill University, Dept. of Mech. Eng. and Centre of Intelligent Machines, Canada
2000 J. Angeles, angeles@cim.mcgill.ca, www.cim.mcgill.ca/~rmsl/Angeles_html,

NOTE: This worksheet demonstrates the use of Maple for the visualization of the modal analysis of a three-degree-of-freedom system, in the area of vertical vibration of mass-transport systems. It illustrates how to compute: (a) the natural frequencies of the system; (b) the modal vectors of the same; (c) the time response of the system to a pulse-type of input; and (d) the visualization of the above items.

Introduction:

Modelling and animation of the three natural modes of the vibrations of the subway cars of Montreal, Paris, Mexico City, etc., running on pneumatic tires. The system is fully described in Example 3.5.3 of

Angeles, J., 2000, 305-315A Dynamics of Vibrations Lecture Notes, Department of Mechanical Engineering, McGill University, Montreal.

The numerical values of the system parameters are taken from manufacturer's information, as available in

Angeles, J. and Espinosa, I., 1981, ``Suspension-system synthesis for mass transport vehicles with prescribed dynamic behavior,'' ASME Paper 81-DET-44, Proc. 1981 ASME Design Engineering Technical Conference, Sept. 20--23, Hartford.

> restart; with(plots): with(plottools): with(linalg):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

In the sequel, we will need the auto procedure to order the eigenvalues from smallest to largest, since Maple does not do this.

Procedure auto computes the eigenvalues and eigenvectors of a 3X3 matrix A, that is assumed symmetric, and returns the eigenvalues alpha_1 <= alpha_2 <= alpha_3 as well as their corresponding eigenvectors ordered in array E:=[e1 e2 e3]. Arguments i, j, k are to be assigned any integer values when calling the procedure. They are set by auto to one of the three values 1, 2, 3, with 1 associated to the smallest and 3 to the largest eigenvalue of A.

> auto:=proc(A)
global alpha,Ev;
local lambda,EE,e1,e2,e3,ee,i,j,k,t;
EE:=eigenvectors(A):
lambda:=<EE[1][1],EE[2][1],EE[3][1]>; print(lambda);
ee[1]:=EE[1][3]: ee[2]:=EE[2][3]: ee[3]:=EE[3][3]:
print(ee[1],ee[2],ee[3]); #lists of eigenvectors of 1st, 2nd, and 3rd eigenvalues
if lambda[1]<=lambda[2] then i:=1: j:=2 else i:=2: j:=1 fi;
print(i,j); #lambda_i is smaller of lambda_1 and lambda_2
if lambda[3]>lambda[j] then k:=3 else k:= j: j:=3:
print(j,k); fi; #lambda_k is larger of lambda_j and lambda_3
if lambda[i]>lambda[j] then t:=j: j:=i: i:=t fi;
print(i,j,k); #if lambda_i is larger than lambda_j, then i and j are exchanged
alpha:=<lambda[i],lambda[j],lambda[k]>; print(alpha);
e1:=ee[i][1]: e2:=ee[j][1]: e3:=ee[k][1]:
print(e1,e2,e3); #[1] indicates first of the list of eigenvectors of ith, jth, or kth eigenvalue. We assume here that each eigenvalue has geometric multiplicity of 1.
Ev:=concat(e1,e2,e3); #entered blockwise
end:

The model is a three-dof undamped system with mass and stiffness matrices given below:

> m_1:=1971: m_2:=3256: m_3:=15780: #masses in kg
M:=matrix([[m_2,0,0],[0,m_1,0],[0,0,m_3]]):
k_1:=4900000: k_2:=3430000: k_3:=837000: k_4:=1783000: #stiffness values in N/m
k_11:=4*(2*k_1+k_2+k_4): k_12:=-4*(2*k_1+k_2):
k_22:=2*(4*k_1+2*k_2+k_3): k_23:=-2*k_3: k_33:=-k_23:
K:=matrix([[k_11,k_12,0],[k_12,k_22,k_23],[0,k_23,k_33]]):
N:= matrix([[evalf(sqrt(m_2)),0,0],[0,evalf(sqrt(m_1)),0],[0,0,evalf(sqrt(m_3))]]):
N_inv:=inverse(N):

> Omega2:=evalm(N_inv&*K&*N_inv): #frequency matrix-squared
auto(Omega2):
omega:=<sqrt(alpha[1]),sqrt(alpha[2]),sqrt(alpha[3])>:
E:=evalm(Ev): ET:=transpose(E):

Therefore, the velocities at which these frequencies are excited because of wheel unbalances are, for a wheel diameter of 0.960 m:

> d:=0.960: v[1]:=omega[1]*d/2: v[2]:=omega[2]*d/2:
v[3]:=omega[3]*d/2: #velocities in m/s
vv[1]:=3.6*v[1]: vv[2]:=3.6*v[2]:
vv[3]:=3.6*v[3]: #foregoing velocities in km/h
Phi:=evalm(N_inv&*E):

Test orthonormality of modal vectors

> P:=evalm(transpose(Phi)&*M&*Phi):

Apparently, modal vectors are orthonormal w.r.t. mass matrix. Normalize the modal vectors by dividing them by their Euclidean norm:

> phi_1:=vector([Phi[1,1],Phi[2,1],Phi[3,1]]):
phi_2:=vector([Phi[1,2],Phi[2,2],Phi[3,2]]):
phi_3:=vector([Phi[1,3],Phi[2,3],Phi[3,3]]):
norm1:=norm(phi_1,2): norm2:=norm(phi_2,2):
norm3:=norm(phi_3, 2):
psi_1:=evalm(phi_1/norm1): psi_2:=evalm(phi_2/norm2): psi_3:=evalm(phi_3/norm3): #normalized modal vectors
Delta:=evalm(ET&*Omega2&*E): #Omega^2 in diagonal form

Roundoff nondiagonal entries to zero and write square root of resulting diagonal matrix

> Del:=matrix([[sqrt(Delta[1,1]),0,0],[0,sqrt(Delta[2,2]), 0],[0,0,sqrt(Delta[3,3])]]): #frequency matrix in diagonal form
Omega:=evalm(E&*Del&*ET): #frequency matrix back to original coordinates
U:=inverse(Omega):
h:=evalf((1/4)*2*Pi/omega[3]): #time increment defined as 1/4 of shortest natural period, in s
Oh:=evalm(Omega*h): #preparations for simulation start here
Delh:=evalm(Del*h):
DC:=matrix([[cos(Delh[1,1]),0,0],[0,cos(Delh[2,2]),0],
[0,0,cos(Delh[3,3])]]): #cos(Omega*h) in diagonal form
DS:=matrix([[sin(Delh[1,1]),0,0],[0,sin(Delh[2,2]),0],
[0,0,sin(Delh[3,3])]]): #sin(Omega*h) in diagonal form
C:=evalm(E&*DC&*ET): #cos(Omega*h) in original coordinates
S:=evalm(E&*DS&*ET): #sin(Omega*h) in original coordinates
CS:=concat(C,S): SC:=concat(-S,C):
F:=stackmatrix(CS,SC): #F matrix for simulation, entered blockwise

B:=0.5: #B, given in m, is used to dimension the boxes representing the masses, for animation purposes.

Section I: The animation of the first mode

> v[1]:=omega[1]*d/2: v0:=evalm(v[1]*psi_1): #initial velocity, of magnitude v[1], in m/s, and in direction of first mode
ww:=evalm(N&*v0):
w:=linsolve(Omega,ww):
z:=vector([0,0,0,w[1],w[2],w[3]]): #initial state vector
n:=250: golden:=evalf((1+sqrt(5))/2): #"golden" is golden-section ratio, which we use to proportion the mass boxes, whose areas are proportional to the masses of the elements that they represent

xronos:= 0: v:=v[1]: # v is travelling velocity
one:=B/2: two:=2*B: three:=two+B/golden:
four:= evalf(sqrt(m_1/m_2)*B/2): five:=4*B:
six:=evalf(five+sqrt(m_1/m_2)*B/golden):
seven:=evalf(sqrt(m_3/m_2)*B/2): eight:=6*B:
nine:=evalf(eight+sqrt(m_3/m_2)*B/golden):
ten:=evalf(sqrt(m_3/m_2)*B/2): eleven:=B/4:
#parameters needed to dimension the boxes

Simulation loop follows. This is based on a zero-input discrete-time model, as per (Angeles, 2000):

z [k+1] = F z [k], z [0] = initial state

> for i from 1 to n do ti:=xronos; vi:=v*ti;
z:=evalm(F&*z); displ:=vector([z[1],z[2],z[3]]);
x:= evalm(N_inv&*displ); #displacement of actual model
o:=[vi,0];
p:=point(o,symbol=circle,color=blue);
m2:=rectangle([vi-one,three+x[1]],[vi+one,two+x[1]],color=blue):
m1:=rectangle([vi-four,six+x[2]],[vi+four,five+x[2]],color=red):
m3:=rectangle([vi-seven,nine+x[3]],[vi+seven,eight+x[3]],color=green);
di:=disk([vi,0],eleven,color=red);
frms||i:=display({p,m1,m2,m3,di},color=red);
xronos:=ti+h; next;
od:
display(seq(frms||i,i=1..n),insequence=true,scaling=constrained);

Section 2: The animation of the second mode

> v[1]:=omega[1]*d/2: v0:=evalm(v[1]*psi_2): #initial velocity, of magnitude v[1], in m/s, and in direction of second mode
ww:=evalm(N&*v0):
w:=linsolve(Omega,ww):
z:=vector([0,0,0,w[1],w[2],w[3]]): #initial state vector
n:=250: golden:=evalf((1+sqrt(5))/2): #"golden" is golden-section ratio, which we use to proportion the mass boxes

xronos:= 0: v:=v[1]: #v is travelling velocity
one:=B/2: two:=2*B: three:=two+B/golden:
four:=evalf(sqrt(m_1/m_2)*B/2): five:=4*B:
six:=evalf(five+sqrt(m_1/m_2)*B/golden):
seven:=evalf(sqrt(m_3/m_2)*B/2): eight:=6*B:
nine:=evalf(eight+sqrt(m_3/m_2)*B/golden):
ten:=evalf(sqrt(m_3/m_2)*B/2): eleven:=B/4: #parameters needed to dimension the boxes, whose areas are proportional to the masses

Simulation loop follows. This is based on the zero-input discrete-time model:

z [k+1] = F z [k], z [0] = initial state

> for i from 1 to n do ti:=xronos; vi:= v*ti;
z:= evalm(F&*z); displ:=vector([z[1], z[2], z[3]]);
x:= evalm(N_inv&*displ);#displacement of actual model
o:=[vi, 0];
p:=point(o, symbol=circle, color=blue);
m2:=rectangle([vi-one,three+x[1]],[vi+one,two+x[1]],color=blue):
m1:=rectangle([vi-four,six+x[2]],[vi+four,five+x[2]],color=red):
m3:=rectangle([vi-seven,nine+x[3]],[vi+seven,eight+x[3]],color=green);
di:= disk([vi,0],eleven,color=red);
frms||i:=display({p,m1,m2,m3,di},color=red);
xronos:= ti+h; next;
od:
display(seq(frms||i,i=1..n),insequence=true,scaling=constrained);

Section 3: The animation of the third mode

> v[1]:=omega[1]*d/2: v0:=evalm(2*v[1]*psi_3): #initial velocity, of magnitude 2v[1], in m/s, and in direction of third mode
ww:=evalm(N&*v0):
w:=linsolve(Omega,ww):
z:=vector([0,0,0,w[1],w[2],w[3]]): #initial state vector
n:=250:

xronos:= 0: v:=v[1]: #v is travelling velocity
one:= B/2: two:=2*B: three:=two+B/golden:
four:=evalf(sqrt(m_1/m_2)*B/2): five:=4*B:
six:=evalf(five+sqrt(m_1/m_2)*B/golden):
seven:=evalf(sqrt(m_3/m_2)*B/2): eight:=6*B:
nine:=evalf(eight+sqrt(m_3/m_2)*B/golden):
ten:=evalf(sqrt(m_3/m_2)*B/2): eleven:=B/4: #parameters needed to dimension the boxes, whose areas are proportional to the masses

Simulation loop follows. This is based on zero-input discrete-time model:

z [k+1] = F z[ k], z [0] = initial state

> for i from 1 to n do ti:=xronos; vi:= v*ti;
z:=evalm(F&*z); displ:=vector([z[1],z[2],z[3]]);
x:=evalm(N_inv&*displ);#displacement of actual model
o:=[vi,0];
p:=point(o,symbol=circle,color=blue);
m2:=rectangle([vi-one,three+x[1]],[vi+one,two+x[1]],color=blue):
m1:=rectangle([vi-four,six+x[2]],[vi+four,five+x[2]],color=red):
m3:=rectangle([vi-seven,nine+x[3]],[vi+seven,eight+x[3]],color=green);
di:= disk([vi,0],eleven,color=red);
frms||i:=display({p,m1,m2,m3,di},color=red);
xronos:= ti+h; next;
od:
display(seq(frms||i,i=1..n),insequence=true,scaling=constrained);

Section 4: The animation of a subway car overcoming a bump

Animation is based on simulation of a subway car overcoming a sinusoidal bump with height equal to the radius of its wheel. Discrete-time model takes the form derived using the zero-order hold to sample the input force.
z [k+1] = F z [k] + G u [k], z [0] = initial state
where z [k] is the state vector at instant t_k, F & G are constant matrices and u [k] = ( N W )^(-1) f [k], with f [k] = f (t_k) being the excitation: f (t_k) = [-4k_4 0 0]^T. Note that W is the frequency matrix and k_4 is the stiffness of the (rubber) tires, its numerical value being given below.

Define input:
Bump is plotted here

> B:=0.25: lambda:=1: v[2]:=omega[2]*d/2: v:=v[2]: T:=evalf(lambda/v): #B in m; lambda in m; v in m/s is the travelling speed. Change parameters as needed
omega:= evalf(2*Pi*v/lambda):

bump :=proc(xronos)global B,omega,T,height:
height:=B*(sin(omega*xronos)*Heaviside(xronos)+sin(omega*(xronos- T/2))*Heaviside(xronos-T/2)):
end:
plot(bump(t),t =-T/2..T,thickness=2,axes=BOXED);

> h:= T/20: #time increment defined as 1/10 of time it takes to overcome the bump, in s
Oh:=evalm(Omega*h): Delh:=evalm(Del*h):

DC:=matrix([[cos(Delh[1,1]),0,0],[0,cos(Delh[2,2]),0],[0,0,cos(Delh[3,3])]]): #cos(Omega*h) in diagonal form
DS:=matrix([[sin(Delh[1,1]),0,0],[0,sin(Delh[2,2]),0], [0,0,sin(Delh[3,3])]]): #sin(Omega*h) in diagonal form
C:=evalm(E&*DC&*ET): #cos(Omega*h) in original coordinates
S:=evalm(E&*DS&*ET): #sin(Omega*h) in original coordinates
CS:=concat(C,S): SC:=concat(-S,C):
F:=stackmatrix(CS,SC): #F matrix for simulation, entered blockwise
ID:=Matrix(3,3,shape=identity):
GU:=matadd(ID,-C): #Upper block of G matrix for simulation
G:=stackmatrix(GU,S): #G matrix for simulation, entered blockwise
k:=vector([-4*k_4,0,0]): #auxiliary vector for simulation
psi:=linsolve(K,k):
phi:=evalm(N&*psi): #Defining input vector for simulation. Vector phi must be multiplied by the bump function
z:=vector([0,0,0,0,0,0]): xronos:=-T/2:
n:=120:
one:=2*B: two:=8*B: three:=two+4*B/golden:
four:=evalf(sqrt(m_1/m_2)*2*B): five:=16*B:
six:=evalf(five+sqrt(m_1/m_2)*4*B/golden):
seven:=evalf(sqrt(m_3/m_2)*2*B): eight:=24*B:
nine:=evalf(eight+sqrt(m_3/m_2)*4*B/golden):
ten:=evalf(sqrt(m_3/m_2)*2*B): #parameters needed to dimension the boxes, whose areas are proportional to the masses

for i from 1 to n do ti:=xronos; vi:=v*ti;
if ti=0 or ti=T/2 then b:=0 else b:=bump(ti) fi;
force:=evalm(phi*b); #input to math model in normal form
z:=evalm(F&*z+G&*force);
displ:=vector([z[1],z[2],z[3]]);
x:=evalm(N_inv&*displ); #displacement of actual model
o:=[vi,0];
p:=point(o,symbol=circle,color=blue);
m2:=rectangle([vi-one,three+x[1]],[vi+one,two+x[1]], color=blue):
m1:=rectangle([vi-four,six+x[2]],[vi+four,five+x[2]], color=red):
m3:=rectangle([vi-seven,nine+x[3]],[vi+seven,eight+x[3]], color=green);
di:=disk([vi,b],one,color=red);
frms||i:=display({p,m1,m2,m3,di},color=red);
xronos:=ti+h; next;
od:
display(seq(frms||i,i=1..n),insequence=true,scaling=constrained);

>

Conclusion:

Maple was used to: a) perform the modal analysis of the three-degree-of-freedom (three-dof) model of a real-life physical system under vertical vibration; b) illustrate the procedure to compute the frequency matrix of a multi-dof mechanical system; c) visualize the natural modes of vibration of the system; and d) simulate and visualize the dynamic response of the system when perturbed by a bump found on its way.

Disclaimer: While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.