Application Center - Maplesoft

App Preview:

3-D oscillator with two masses, coupled by elastic springs

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


Image 

3-D Oscillator with Two Masses, Coupled by Elastic Springs 

? 2001 Harald Kammerer, GERB Schwingungsisolierungen, Germany
www.gerb.com 

 

Plot_2d 

Introduction 

This worksheet shows the calculation of the motion of a system with 12 degrees of freedom excited by ground motion.  

The following Maple techniques are highlighted: 

 

  • Developing a general model
 

 

 

  • Reading in specific data from a file
 

 

 

  • Calculation of system matrices
 

 

 

  • Creating a 3D animation of the system
 

 

Problem Description 

Image 

Consider a system of two masses. Both masses are described by rectangular prisms. The support for the two masses are elastic.  The ground motion is known and is stored in the file "motion.dat".
Note: The file "motion.dat" includes four columns. First the discrete time steps, second the values of the ground acceleration in x-direction, third in y-direction and fourth in z-direction. 

Contents 

 

  • 1) In the first section of the worksheet the mechanical system is described. All data points are in m, kg, N and sec.
 

 

 

  • 2) In the second section, the system matrices are calculated. The stiffness matrix is calculated according [1].
 

 

 

  • 3) In section three, the eigenvalue problem is solved and the modal matrix is calculated.
 

 

 

  • 4) In section four, the time history of the system movement is calculated numerically. The NEWMARK-algorithm is used to solve the incremental equations of motion [2].
 

 

 

  • 5) In the last section, the solutions are displayed. First, the eigenvalues and eigenmodes are listed. Then the time history of the motion in every direction relative to the ground are plotted. Finally, the system motion is illustrated in the form of a small animation (which is duplicated at the top of this worksheet).
 

 

The lower mass and the lower spring in general are described by the index "l", the upper ones by the index "o".  

The coordinates 1, 2 and 3 are the translation of the lower mass in x-, y- and z-direction,  4, 5 and 6 the corresponding rotations. The coordinates 7, 8 and 9 are the translations of the upper mass, 10, 11 and 12 the rotations. 

Initialization 

CodeEditor-Button#initialization 

(Klick for restart and load packages) 

Description of the Mechanical Model 

Dimensions of the Bodies 

The two bodies are cuboids and have the following dimensions and masses (lower body; upper body): 

x-direction: ll; lo [m] 

y-direction: bl; bo [m] 

z-direction: hl; ho [m] 

mass: ml; mo [kg] 

lower body 

length: 

> `assign`(ll, 10.0); -1
 

width: 

> `assign`(bl, 8.0); -1
 

height: 

> `assign`(hl, 2.0); -1
 

weight: 

> `assign`(ml, 0.1e7); -1
 

The coordinates of two diagonal positioned edges are later needed for the animation of the bodies. Here they are defined relative to the c.o.g. of the body. 

End points of 

the diagonal: 

> `assign`(ENDl1, [VectorCalculus:-`-`(VectorCalculus:-`*`(ll, `/`(1, 2))), VectorCalculus:-`-`(VectorCalculus:-`*`(bl, `/`(1, 2))), VectorCalculus:-`-`(VectorCalculus:-`*`(hl, `/`(1, 2)))]); -1
 

> `assign`(ENDl2, [VectorCalculus:-`*`(ll, `/`(1, 2)), VectorCalculus:-`*`(bl, `/`(1, 2)), VectorCalculus:-`*`(hl, `/`(1, 2))]); -1
 

upper body 

length: 

> `assign`(lo, 8.0); -1
 

width: 

> `assign`(bo, 7.0); -1
 

height: 

> `assign`(ho, 2.0); -1
 

weight: 

> `assign`(mo, 0.1e7); -1
 

Again we need define the end points of one diagonal for later animations 

End points of 

the diagonal: 

> `assign`(ENDo1, [VectorCalculus:-`-`(VectorCalculus:-`*`(lo, `/`(1, 2))), VectorCalculus:-`-`(VectorCalculus:-`*`(bo, `/`(1, 2))), VectorCalculus:-`-`(VectorCalculus:-`*`(ho, `/`(1, 2)))]); -1
 

> `assign`(ENDo2, [VectorCalculus:-`*`(lo, `/`(1, 2)), VectorCalculus:-`*`(bo, `/`(1, 2)), VectorCalculus:-`*`(ho, `/`(1, 2))]); -1
 

distance between the centers of gravity of both bodies 

this is necessary to be known for the animation later 

> `assign`(DSlo, [0, 0, VectorCalculus:-`+`(hl, ho)]); -1
 

General Spring Characteristics 

The general characteristics of the springs are their stiffness in horizontal and vertical direction. We assume in this calculation that the stiffness in all horizontal directions is the same. Further we define one set of values for the horizontal stiffness and the vertical stiffness for all the lower springs and one set of values for all the upper springs. This general data are modified in the following subsection for each individual spring. 


coh/clh: stiffness of the upper/lower spring  in horizontal direction [N/m] 

cov/clv: stiffness of the upper/lower spring  in vertical direction [N/m] 

 

The rotational stiffness of the springs are neglected.
 

lower springs 

> `assign`(clh, VectorCalculus:-`*`(5., `^`(10, 6))); -1`assign`(clv, VectorCalculus:-`*`(5., `^`(10, 6))); -1
 

upper springs 

> `assign`(coh, VectorCalculus:-`*`(1., `^`(10, 6))); -1`assign`(cov, VectorCalculus:-`*`(1., `^`(10, 6))); -1
 

Position of the Bodies Supporting Springs and their Individual Characteristics 

All positions are given in metres (m). They are given with respect to the c.o.g. of the bodies. The stiffness of the individual springs can be manipulated in the following, see e.g. 1. spring of lower mass.
 

lower body 


: position of the lower spring i in direction j with respect to the lower bodies c.o.g.
: stiffness of the lower spring i in direction j 

where i = 1..4 is the number of the spring and j = x, y, z is the direction in which the stiffness is defined.
 

Spring 1 

The coordinates for the position of the springs are defined with respect to the c.o.g. of the lower mass. 

Note: The stiffness are varied from the general stiffness defined previously. 

position in x-direction (horizontal) 

> `assign`(a[1, x], VectorCalculus:-`-`(VectorCalculus:-`*`(ll, `/`(1, 2))))
 

position in y-direction (horizontal) 

> `assign`(a[1, y], VectorCalculus:-`-`(VectorCalculus:-`*`(bl, `/`(1, 2))))
 

position in z-direction (vertical) 

> `assign`(a[1, z], VectorCalculus:-`-`(VectorCalculus:-`*`(hl, `/`(1, 2))))
 

spring stiffness in x-direction (horizontal) 

> `assign`(cl[1, x], VectorCalculus:-`*`(.9, clh))
 

spring stiffness in y-direction (horizontal) 

> `assign`(cl[1, y], VectorCalculus:-`*`(.85, clh))
 

spring stiffness in z-direction (vertical) 

> `assign`(cl[1, z], VectorCalculus:-`*`(.8, clv))
 

Spring 2 

Notice that we assume the stiffness (in the z direction) of this spring is significantly different from that of the others.  This could be caused by damage to the spring. 

position in x-direction (horizontal) 

> `assign`(a[2, x], VectorCalculus:-`*`(ll, `/`(1, 2)))
 

position in y-direction (horizontal) 

> `assign`(a[2, y], VectorCalculus:-`-`(VectorCalculus:-`*`(bl, `/`(1, 2))))
 

position in z-direction (vertical) 

> `assign`(a[2, z], VectorCalculus:-`-`(VectorCalculus:-`*`(hl, `/`(1, 2))))
 

spring stiffness in x-direction (horizontal) 

> `assign`(cl[2, x], clh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(cl[2, y], clh)
 

spring stiffness in z-direction (vertical) 

> `assign`(cl[2, z], VectorCalculus:-`*`(.1, clv))
 

 

Spring 3 

position in x-direction (horizontal) 

> `assign`(a[3, x], VectorCalculus:-`*`(ll, `/`(1, 2)))
 

position in y-direction (horizontal) 

> `assign`(a[3, y], VectorCalculus:-`*`(bl, `/`(1, 2)))
 

position in z-direction (vertical) 

> `assign`(a[3, z], VectorCalculus:-`-`(VectorCalculus:-`*`(hl, `/`(1, 2))))
 

spring stiffness in x-direction (horizontal) 

> `assign`(cl[3, x], clh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(cl[3, y], clh)
 

spring stiffness in z-direction (vertical) 

> `assign`(cl[3, z], clv)
 

Spring 4 

position in x-direction (horizontal) 

> `assign`(a[4, x], VectorCalculus:-`-`(VectorCalculus:-`*`(ll, `/`(1, 2))))
 

position in y-direction (horizontal) 

> `assign`(a[4, y], VectorCalculus:-`*`(bl, `/`(1, 2)))
 

position in z-direction (vertical) 

> `assign`(a[4, z], VectorCalculus:-`-`(VectorCalculus:-`*`(hl, `/`(1, 2))))
 

spring stiffness in x-direction (horizontal) 

> `assign`(cl[4, x], clh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(cl[4, y], clh)
 

spring stiffness in z-direction (vertical) 

> `assign`(cl[4, z], clv)
 

Overall Parameters 

Here, we will display the two sets of parameters in matrix form 

> a = Matrix([seq([seq(a[j, i], `in`(i, [x, y, z]))], j = 1 .. 4)])
 

> cl = Matrix([seq([seq(cl[j, i], `in`(i, [x, y, z]))], j = 1 .. 4)])
 

upper mass 

 

: position of the upper spring i in direction j
: position of the upper spring i in direction j
: vertical position of the upper springs lower end (which is connected with the lower mass) with respect to the c.o.g. of the lower mass
Again i = 1..4 is the number of the spring and j = x, y, z is the direction in which the stiffness is defined.
 

Spring 1 

Notice that the upper springs here need additional information on the position of the lower mass support. 

position in x-direction (horizontal) 

> `assign`(b[1, x], VectorCalculus:-`-`(VectorCalculus:-`*`(lo, `/`(1, 2))))
 

position in y-direction (horizontal) 

> `assign`(b[1, y], VectorCalculus:-`-`(VectorCalculus:-`*`(bo, `/`(1, 2))))
 

position in z-direction (vertical) 

> `assign`(b[1, z], VectorCalculus:-`-`(VectorCalculus:-`*`(ho, `/`(1, 2))))
 

position in z-direction against the lower mass 

> `assign`(bq[1, z], VectorCalculus:-`*`(hl, `/`(1, 2)))
 

spring stiffness in x-direction (horizontal) 

> `assign`(co[1, x], coh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(co[1, y], coh)
 

spring stiffness in z-direction (vertical) 

> `assign`(co[1, z], cov)
 

Spring 2 

position in x-direction (horizontal) 

> `assign`(b[2, x], VectorCalculus:-`*`(lo, `/`(1, 2)))
 

position in y-direction (horizontal) 

> `assign`(b[2, y], VectorCalculus:-`-`(VectorCalculus:-`*`(bo, `/`(1, 2))))
 

position in z-direction (vertical) 

> `assign`(b[2, z], VectorCalculus:-`-`(VectorCalculus:-`*`(ho, `/`(1, 2))))
 

position in z-direction against the lower mass 

> `assign`(bq[2, z], VectorCalculus:-`*`(hl, `/`(1, 2)))
 

spring stiffness in x-direction (horizontal) 

> `assign`(co[2, x], coh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(co[2, y], coh)
 

spring stiffness in z-direction (vertical) 

> `assign`(co[2, z], cov)
 

Spring 3 

position in x-direction (horizontal) 

> `assign`(b[3, x], VectorCalculus:-`*`(lo, `/`(1, 2)))
 

position in y-direction (horizontal) 

> `assign`(b[3, y], VectorCalculus:-`*`(bo, `/`(1, 2)))
 

position in z-direction (vertical) 

> `assign`(b[3, z], VectorCalculus:-`-`(VectorCalculus:-`*`(ho, `/`(1, 2))))
 

position in z-direction against the lower mass 

> `assign`(bq[3, z], VectorCalculus:-`*`(hl, `/`(1, 2)))
 

spring stiffness in x-direction (horizontal) 

> `assign`(co[3, x], coh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(co[3, y], coh)
 

spring stiffness in z-direction (vertical) 

> `assign`(co[3, z], cov)
 

Spring 4 

position in x-direction (horizontal) 

> `assign`(b[4, x], VectorCalculus:-`-`(VectorCalculus:-`*`(lo, `/`(1, 2))))
 

position in y-direction (horizontal) 

> `assign`(b[4, y], VectorCalculus:-`*`(bo, `/`(1, 2)))
 

position in z-direction (vertical) 

> `assign`(b[4, z], VectorCalculus:-`-`(VectorCalculus:-`*`(ho, `/`(1, 2))))
 

position in z-direction against the lower mass 

> `assign`(bq[4, z], VectorCalculus:-`*`(hl, `/`(1, 2)))
 

spring stiffness in x-direction (horizontal) 

> `assign`(co[4, x], coh)
 

spring stiffness in y-direction (horizontal) 

> `assign`(co[4, y], coh)
 

spring stiffness in z-direction (vertical) 

> `assign`(co[4, z], cov)
 

Overall Parameters 

Here, we will display the two sets of parameters in matrix form 

> b = Matrix([seq([seq(b[j, i], `in`(i, [x, y, z]))], j = 1 .. 4)])
 

> co = Matrix([seq([seq(co[j, i], `in`(i, [x, y, z]))], j = 1 .. 4)])
 

Calculate the System Matrices 

System Stiffness 

The system stiffness is calculated in some steps one after the other. 

The basic idea is to get the springs stiffness from the relation of the potential energy of the deformed springs with the stiffness matrix K and the vector of the deformation q. This relation is valid for the total system as well as for every single spring with every main direction (x, y and z) for itselves. The potential energy of the total system is the sum of the potential energy of all single springs with respect to every main directions. 

 

And this is the method to get the stiffness matrix. We consider the potential energy Uij = cijuj2 for spring i in direction j = x, y, z. First we do this in general manner for the lower and the upper springs. 

 

Next we extract the entries of the global stiffness matrix from the this expression for the potential energy for each direction. This can easy be done. Just find the coefficient in the relation for U of all global coordinates xi xj with x = q, p (lower mass and upper mass) and i, j = 1..6, see the set xlist. 

 

So the first step to get the stiffness matrix is to calculate the deformation of the springs and the square of all this. 

 

deformation of the lower springs 

Here, we describe the deformation of the lower springs in general manner in terms of the position and orientation of the lower mass. 

 

q1, q2,...,q6 are the degrees of freedom (q1, q2, q3 being the linear motion in x, y, z direction respectively and q4, q5 and q6 being the rotation about the x, y, z axis with respect to the c.o.g.) of the lower body. 

 

ak,x, ak,y, ak,z are the general descriptions of the lower springs positions with respect to the c.o.g. of the lower mass. They are substituted later by the real spring positions of every spring one after the other when the stiffness matrix is generated. 

 

u1, u2 and u3 are the spring deformations in the x, y and z direction respectively. Note that the springs rotational stiffness is neglected. 

 

> `assign`(u1, VectorCalculus:-`+`(VectorCalculus:-`+`(q1, VectorCalculus:-`*`(a[k, z], q5)), VectorCalculus:-`-`(VectorCalculus:-`*`(a[k, y], q6)))); -1`assign`(u2, VectorCalculus:-`+`(VectorCalculus:-`+`(q2, VectorCalculus:-`-`(VectorCalculus:-`*`(a[k, z], q4))), VectorCalculus:-`*`(a[k, x], q6))); -1`assign`(u3, VectorCalculus:-`+`(VectorCalculus:-`+`(q3, VectorCalculus:-`*`(a[k, y], q4)), VectorCalculus:-`-`(VectorCalculus:-`*`(a[k, x], q5)))); -1
 

> `assign`(uq1, expand(simplify(`*`(`^`(u1, 2))))); 1`assign`(uq2, expand(simplify(`*`(`^`(u2, 2))))); 1`assign`(uq3, expand(simplify(`*`(`^`(u3, 2)))))
 

deformation of the upper springs 

For the upper springs, the deformation is described in terms of the position and orientation of the upper and the lower body. 

 

p1, p2,...,p6 are the degrees of freedom of the upper body (similar to the q defined for the lower body). 

 

bk,x, bk,y, bk,z are the general descriptions of the upper end of the upper springs positions with respect to the c.o.g. of the upper mass, bqk,z is the general description of the lower end of the upper springs positions with respect to the c.o.g. of the lower mass. They are substituted later by the real spring positions of every spring one after the other when the stiffness matrix is generated. 

 

o1, o2 and o3 are the spring deformations in the x, y, and z direction respectively. 

 

> `assign`(o1, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`-`(q1), VectorCalculus:-`-`(VectorCalculus:-`*`(bq[k, z], q5))), Vecto...`assign`(o2, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`-`(q2), VectorCalculus:-`*`(bq[k, z], q4)), VectorCalculus:-`-`(Vector...`assign`(o3, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`-`(q3), VectorCalculus:-`-`(VectorCalculus:-`*`(b[k, y], q4))), Vector...
 

> `assign`(oq1, expand(simplify(`*`(`^`(o1, 2))))); 1`assign`(oq2, expand(simplify(`*`(`^`(o2, 2))))); 1`assign`(oq3, expand(simplify(`*`(`^`(o3, 2)))))
 

general form of stiffness matrix for individual spring directions 


Using the above deformation equations, we can extract the general form of the spring stiffness components for each direction from the potential energy according to [1]. 

cl1, cl2, cl3: 

12 x 12 stiffness matrices for the lower springs in x, y, and z direction respectively. 

co1, co2, co3: 

12 x 12 stiffness matrices for the upper springs in x, y, and z direction respectively. 

> `assign`(xlist, [seq(q || i, i = 1 .. 6), seq(p || i, i = 1 .. 6)]); 1`assign`(x2list, [seq(`*`(`^`(xlist[i], 2)), i = 1 .. nops(xlist))])
 

Now we calculate the potential energy of the springs in every direction one after the other. The coefficients of the expression cspring number,direction qi qj (same for pi and pj) yield the entry i,j of the stiffness matrix of each direction of the lower and the upper springs. We do this in two steps. First get the diagonal matrix with the coefficients if j=i using x2list, next we add the entries for i≠j

> for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
for l to 3 do `assign`(lxyz, [x, y, z]); `assign`(cl || l, evalm(VectorCalculus:-`+`(DiagonalMatrix([seq(coeff(VectorCalculus:-`*`(cl[k, lxyz[l]], uq || l), x2list[i]), i = 1 .. 12)]), VectorCalculus:...
 

stiffness matrix for every spring in every direction 

With the above general form of the stiffness matrices in the x, y, and z direction, we can generate the stiffness matrix associated with each spring by substituting the parameters corresponding to each of the springs. 

 

The following index notion are used for the resulting table of stiffness matrices: 

 

first index (): 1, 2, 3 and 4 for lower springs; 5, 6, 7 and 8 for upper springs, 

second index (): for direction x, y and z 

 

> for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
for s_index to 4 do `assign`(n, 1); for d_index in [x, y, z] do `assign`(k[s_index, d_index], Matrix(eval(cl || n, {seq(a[k, i] = a[s_index, i], `in`(i, {x, y, z})), seq(b[k, i] = b[s_index, i], `in`(...
 

overall stiffness matrix 

 

Finally, to get the overall system stiffness matrix, we add up all of the individual stiffness matrices. 

 

> `assign`(K, Matrix(12, 12, 0)); -1for i to 8 do `assign`(K, VectorCalculus:-`+`(K, VectorCalculus:-`+`(VectorCalculus:-`+`(k[i, x], k[i, y]), k[i, z]))) end do; -1
 

> K
 

(right mouse click on the place holder and choose " Browse " to see the entries of the matrix) 

System Mass 

matrix of one body 

 

To obtain the system mass matrix, we first define the mass matrix for one body. 

 

> `assign`(Mr, DiagonalMatrix([m, m, m, VectorCalculus:-`*`(VectorCalculus:-`*`(m, VectorCalculus:-`+`(`*`(`^`(b, 2)), `*`(`^`(h, 2)))), `/`(1, 12)), VectorCalculus:-`*`(VectorCalculus:-`*`(m, VectorCal...
 

general mass matrix 

 

Then we combine them to get the overall system mass matrix. 

 

> `assign`(M1, evalm(subs({b = bl, h = hl, l = ll, m = ml}, evalm(Mr)))); -1
 

> `assign`(M2, evalm(subs({b = bo, h = ho, l = lo, m = mo}, evalm(Mr)))); -1
 

> `assign`(M, DiagonalMatrix([M1, M2]))
 

(right mouse click on the place holder and choose " Browse " to see the entries of the matrix) 

System Matrix 

 

Finally, the system matrix is the product of the inverse of the mass matrix and the system stiffness matrix. 

 

> `assign`(SM, Typesetting:-delayDotProduct(MatrixInverse(M), K)); 1
 

(right mouse click on the place holder and choose " Browse " to see the entries of the matrix) 

Eigenvalue Problem Solution 

 

With the system matrix, we can now look at the eigenvalue problem. 

 

Eigenvalues and Eigenmodes 

 

First, we compute the eigenvalues and eigenvectors of the system matrix using the Eigenvectors procedure from the LinearAlgebra package. 

 

> `assign`(lambda, vecs, Eigenvectors(SM))
 

Here, we split the eigenvector matrix into individual vectors so that we can sort them based on the associated eigenvalues. 

 

> for i to 12 do `assign`(v || i, convert(SubMatrix(vecs, [1 .. 12], [i]), Vector)) end do; -1
 

 

Sorting with respect to eigenvalues: 

 

> `assign`(lstsrt, sort([seq(Re(lambda[i]), i = 1 .. 12)]))
 

In this calculation we consider no damping, so we know that the result have to be real and we suppress the imaginary part 0. I 

> for i to 12 do for j to 12 do if lstsrt[i] = lambda[j] then `assign`(num[i], j); if `and`(`>`(i, 1), `<>`(num[i], num[VectorCalculus:-`+`(i, VectorCalculus:-`-`(1))])) then `assign`(j, 12) end if end ...
for i to 12 do for j to 12 do if lstsrt[i] = lambda[j] then `assign`(num[i], j); if `and`(`>`(i, 1), `<>`(num[i], num[VectorCalculus:-`+`(i, VectorCalculus:-`-`(1))])) then `assign`(j, 12) end if end ...
for i to 12 do for j to 12 do if lstsrt[i] = lambda[j] then `assign`(num[i], j); if `and`(`>`(i, 1), `<>`(num[i], num[VectorCalculus:-`+`(i, VectorCalculus:-`-`(1))])) then `assign`(j, 12) end if end ...
for i to 12 do for j to 12 do if lstsrt[i] = lambda[j] then `assign`(num[i], j); if `and`(`>`(i, 1), `<>`(num[i], num[VectorCalculus:-`+`(i, VectorCalculus:-`-`(1))])) then `assign`(j, 12) end if end ...
for i to 12 do for j to 12 do if lstsrt[i] = lambda[j] then `assign`(num[i], j); if `and`(`>`(i, 1), `<>`(num[i], num[VectorCalculus:-`+`(i, VectorCalculus:-`-`(1))])) then `assign`(j, 12) end if end ...
 

> for i to 12 do `assign`(j, num[i]); `assign`(ev || i, v || j) end do; -1
 

> `assign`(ew, Vector([seq(lstsrt[i], i = 1 .. 12)]))
 

(right mouse click on the place holder and choose " Browse " to see the entries of the vector). 

 

Now numi is the sequence of the position number of the eigenvalues in the vector λ in increasing order, ew is the vector of the the eigenvalues with increasing order and the ev||i are the corresponding eigenvectors. It is in general not neccessary to sort the eigenvalues and eigenvectors but it is usual to do so. 

> seq(num[i], i = 1 .. 12)
 

Eigenfrequencies 

 

Now we extract the eigenfrequencies from the eigenvalues. 

 

> `assign`(omega, [seq(sqrt(ew[i]), i = 1 .. 12)]); 1; `assign`(f, [seq(evalf(VectorCalculus:-`*`(omega[i], `/`(1, `*`(VectorCalculus:-`*`(2, Pi))))), i = 1 .. 12)]); 1
 

>
 

Modal Matrix 

 

We obtain the modal matrix by normalizing the eigenvectors. 

 

> for i to 12 do `assign`(phi || i, Normalize(ev || i)) end do; -1
 

> `assign`(phi, Matrix([seq(phi || i, i = 1 .. 12)]))
 

> `assign`(`φt`, Transpose(phi)); -1
 

(right mouse click on the place holder and choose " Browse " to see the entries of the matrix) 

Modal Mass and Modal Stiffness 

The modal mass matrix and the modal stiffness matrix are obtained by the following calculation. The results have to be diagonal matrices, so we set the nondiagonal elements to zero to remove any numeric trash. 

> `assign`(Ms, Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(`φt`, M), phi)); -1`assign`(Ms, DiagonalMatrix([seq(Ms[i, i], i = 1 .. 12)]))
 

> `assign`(Ks, Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(`φt`, K), phi)); -1`assign`(Ks, DiagonalMatrix([seq(Ks[i, i], i = 1 .. 12)]))
 

(right mouse click on the place holder and choose " Browse " to see the entries of the matrices) 

Numeric Solution of the Equation of Motion 

 

Here, we will obtain the numeric solution to the equation of motion using the NEWMARK algorithm [2] 

 

Time-History of the Ground Motion 

Reading in data from file 

 

First, we will read in the ground motion from the data file motion.dat. We make some rearrangements to this data.  

 

Recall that the data is listed in 4 columns: 

1. time, 2. x-acceleration, 3. y-acceleration, 4. z-acceleration 

 

> `assign`(loaddat, Array(readdata(
 

(right mouse click on the place holder and choose " Browse " to see the entries) 

 

Extract the time step from the data 

> `assign`(dt, VectorCalculus:-`+`(loaddat[4][1], VectorCalculus:-`-`(loaddat[3][1])))
 

 

Getting the number of samples 

> `assign`(ndt, max(seq(i, i = ArrayDims(loaddat)[1])))
 

 

Setting the total time 

> `assign`(maxT, VectorCalculus:-`*`(VectorCalculus:-`+`(ndt, VectorCalculus:-`-`(1)), dt))
 

 

Setting up the discrete time variables for use later. 

> for i from 0 to ndt do `assign`(t || i, VectorCalculus:-`*`(i, dt)) end do; -1
 

 

Assigning the direction data into separate variables. 

> `assign`(a0x, [seq(loaddat[i][2], i = 1 .. ndt)]); -1`assign`(a0y, [seq(loaddat[i][3], i = 1 .. ndt)]); -1`assign`(a0z, [seq(loaddat[i][4], i = 1 .. ndt)]); -1
 

 

For the animation later, the ground acceleration is numerically integrated twice to get the ground displacement. For the numeric integration,  it is assumed that the acceleration between two discrete time steps is linear. Additionally, we assume that the velocities and the displacement at the beginning of the time period is zero. And because we only need the motion in discrete time steps it is exact enough to do the integration by the trapezoidal rule. 

> `assign`(v0x, Vector[column](ndt, 0)); -1`assign`(v0y, Vector[column](ndt, 0)); -1`assign`(v0z, Vector[column](ndt, 0)); -1
 

> `assign`(x0x, Vector[column](ndt, 0)); -1`assign`(x0y, Vector[column](ndt, 0)); -1`assign`(x0z, Vector[column](ndt, 0)); -1
 

> for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
for j from 2 to ndt do for i in [x, y, z] do `assign`(v0 || i[j], VectorCalculus:-`+`(v0 || i[VectorCalculus:-`+`(j, VectorCalculus:-`-`(1))], VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`...
 

Plotting the input data 

Setting up the line colour. 

> `assign`(clax, red); -1`assign`(clay, blue); -1`assign`(claz, green); -1
 

Ground Acceleration 

x - direction 

> `assign`(SSax, [seq([t || j, a0x[j]], j = 1 .. ndt)]); -1`assign`(Pax, curve(SSax, color = clax)); -1
 

y - direction 

> `assign`(SSay, [seq([t || j, a0y[j]], j = 1 .. ndt)]); -1`assign`(Pay, curve(SSay, color = clay)); -1
 

z-direction 

> `assign`(SSaz, [seq([t || j, a0z[j]], j = 1 .. ndt)]); -1`assign`(Paz, curve(SSaz, color = claz)); -1
 

Ground Displacement 

x - direction 

> `assign`(SSdx, [seq([t || j, x0x[j]], j = 1 .. ndt)]); -1`assign`(Pdx, curve(SSdx, color = clax)); -1
 

y - direction 

> `assign`(SSdy, [seq([t || j, x0y[j]], j = 1 .. ndt)]); -1`assign`(Pdy, curve(SSdy, color = clay)); -1
 

z - direction 

> `assign`(SSdz, [seq([t || j, x0z[j]], j = 1 .. ndt)]); -1`assign`(Pdz, curve(SSdz, color = claz)); -1
 

Plots 

Accelerations 

> display({Pax, Pay, Paz}, title =
display({Pax, Pay, Paz}, title =
 

Displacement 

> display({Pdx, Pdy, Pdz}, title =
display({Pdx, Pdy, Pdz}, title =
 

Computing the excitation forces from the ground motion 

 

The equation of motion is 

 

with the vector of the relativ motion . The vector contains the influece of the ground acceleration on the degrees of freedom of the structure. This vector is now prepared for all discrete time steps. 

The  

> for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
for i to ndt do `assign`(a0vec[i], Vector(12, 0)); `assign`(a0vec[i][1], a0x[i]); `assign`(a0vec[i][7], a0x[i]); `assign`(a0vec[i][2], a0y[i]); `assign`(a0vec[i][8], a0y[i]); `assign`(a0vec[i][3], a0z...
 

See that there are no ground rotations in this example, this is normally assumed in earthquake calculations. Take the term to the right side of the equation and substitute allows to use the standard formulation, here we use the NEWMARK-allgorithm, for solve the equation of motion for excitation by forces: 

 

So we calculate now the vector of the substitute excitation force for all discrete time steps. 

> for i to ndt do `assign`(F[i], VectorCalculus:-`-`(Typesetting:-delayDotProduct(M, a0vec[i]))) end do; -1
 

 

After all we transform the substitute force into modal coordinates 

 

> for i to ndt do `assign`(Fs[i], Typesetting:-delayDotProduct(`φt`, F[i])) end do; -1
 

NEWMARK-Algorithm 

 

Now, apply the NEWMARK algorithm. This algorithm is assumed to be known and is not described in detail in this document. 

For the vector of the motion in modal coordinates we use the variable Q, the corresponding velocity Qp and the corresponding acceleration Qpp

 

First, setup the NEWMARK-constants 

> `assign`(alpha, .25); 1`assign`(delta, .5)
 

 

Then, we compute the modified mass matrix 

> `assign`(Mn, VectorCalculus:-`+`(Ms, VectorCalculus:-`*`(Typesetting:-delayDotProduct(Ks, alpha, true), `*`(`^`(dt, 2)))))
 

Next initialize the start vectors. 

> `assign`(RS[1], Vector(12, 0)); -1; `assign`(Q[1], Vector(12, 0)); -1; `assign`(Qp[1], Vector(12, 0)); -1; `assign`(Qpp[1], Vector(12, 0)); -1
`assign`(RS[1], Vector(12, 0)); -1; `assign`(Q[1], Vector(12, 0)); -1; `assign`(Qp[1], Vector(12, 0)); -1; `assign`(Qpp[1], Vector(12, 0)); -1
 

Then compute the output one step at a time. 

> for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
for i from 2 to ndt do `assign`(RS[i], Vector([seq(VectorCalculus:-`+`(Fs[i][j], VectorCalculus:-`-`(VectorCalculus:-`*`(Ks[j, j], VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Q[VectorC...
 

Finally, we can transform the solution from the modal coordinates Q back into real coordinates X. Remember: This coordinates are relative to the ground motion. 

> for i to ndt do `assign`(X[i], map(Re, Typesetting:-delayDotProduct(phi, Q[i]))) end do; -1
 

Viewing the Solution 

Graphical View of the Time History of the Motion 

Translations 

x - direction (relative to ground motion) 

Setup the data points. 

> `assign`(SSXl, [seq([t || j, X[j][1]], j = 1 .. ndt)]); -1`assign`(SSXo, [seq([t || j, X[j][7]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(PXl, curve(SSXl, color = red)); -1`assign`(PXo, curve(SSXo, color = blue)); -1
 

> display({PXl, PXo}, title =
 

y - direction (relative to ground motion) 

Setup the data points. 

> `assign`(SSYl, [seq([t || j, X[j][2]], j = 1 .. ndt)]); -1; `assign`(SSYo, [seq([t || j, X[j][8]], j = 1 .. ndt)]); -1
`assign`(SSYl, [seq([t || j, X[j][2]], j = 1 .. ndt)]); -1; `assign`(SSYo, [seq([t || j, X[j][8]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(PYl, curve(SSYl, color = red)); -1`assign`(PYo, curve(SSYo, color = blue)); -1
 

> display({PYl, PYo}, title =
 

z -direction (relative to ground motion) 

Setup the data points. 

> `assign`(SSZl, [seq([t || j, X[j][3]], j = 1 .. ndt)]); -1`assign`(SSZo, [seq([t || j, X[j][9]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(PZl, curve(SSZl, color = red)); -1`assign`(PZo, curve(SSZo, color = blue)); -1
 

> display({PZl, PZo}, title =
 

Rotations 

Rotation about x 

Setup data points. 

> `assign`(SRXl, [seq([t || j, X[j][4]], j = 1 .. ndt)]); -1`assign`(SRXo, [seq([t || j, X[j][10]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(RXl, curve(SRXl, color = red)); -1`assign`(RXo, curve(SRXo, color = blue)); -1
 

> display({RXl, RXo}, title =
 

Rotation about y 

Setup data points. 

> `assign`(SRYl, [seq([t || j, X[j][5]], j = 1 .. ndt)]); -1`assign`(SRYo, [seq([t || j, X[j][11]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(RYl, curve(SRYl, color = red)); -1`assign`(RYo, curve(SRYo, color = blue)); -1
 

> display({RYl, RYo}, title =
 

Rotation about z 

Setup data points. 

> `assign`(SRZl, [seq([t || j, X[j][6]], j = 1 .. ndt)]); -1`assign`(SRZo, [seq([t || j, X[j][12]], j = 1 .. ndt)]); -1
 

Setup the plot. 

> `assign`(RZl, curve(SRZl, color = red)); -1`assign`(RZo, curve(SRZo, color = blue)); -1
 

> display({RZl, RZo}, title =
 

Animation 

 

To reduce execution time, only some time steps are used for the animation. If you have a fast machine, you can increase this number. 

> `assign`(nmx, 100); 1
 

 

Setting the number of plots used for the animation. 

> if `<`(nmx, ndt) then `assign`(nplot, nmx) else `assign`(nplot, ndt) end if
 

 

Animation stepsize 

> `assign`(ds, floor(VectorCalculus:-`*`(ndt, `/`(1, `*`(nplot)))))
 

 

Maximal motion 

> `assign`(maxm, max(seq(seq(X[j][i], i = 1 .. 12), j = 1 .. ndt)))
 

 

Calculate the normalization factor of the motion for the animation (using the Maximal motion magnitude). 

> `assign`(norml, VectorCalculus:-`*`(DSlo[3], `/`(1, `*`(VectorCalculus:-`*`(maxm, 10)))))
 

Setting up the basic elements 

Nominal ground plate location 

> `assign`(GROUND, translate(display(cuboid([VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll, lo), `/`(1, 2)), 1)), VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll,...
`assign`(GROUND, translate(display(cuboid([VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll, lo), `/`(1, 2)), 1)), VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll,...
`assign`(GROUND, translate(display(cuboid([VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll, lo), `/`(1, 2)), 1)), VectorCalculus:-`-`(VectorCalculus:-`+`(VectorCalculus:-`*`(max(ll,...
 

Nominal Diagonal Points for each prism at each animation time 

Here we need the absolute motion. So we add the ground displacement and the relative motion. 

> `assign`(iplot, 0)
 

> for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(xg[i], x0x[j]); `assign`(yg[i], x0y[j]); `assign`(zg[i], x0z[j]); `assign`(rotl[i], VectorCalculus:-`...
 

Creating the plot objects. 

> for j to nplot do `assign`(CUBlc[j], cuboid(ENDl1t[j], ENDl2t[j], color = red)); `assign`(CUBoc[j], cuboid(ENDo1t[j], ENDo2t[j], color = blue)) end do; -1
for j to nplot do `assign`(CUBlc[j], cuboid(ENDl1t[j], ENDl2t[j], color = red)); `assign`(CUBoc[j], cuboid(ENDo1t[j], ENDo2t[j], color = blue)) end do; -1
for j to nplot do `assign`(CUBlc[j], cuboid(ENDl1t[j], ENDl2t[j], color = red)); `assign`(CUBoc[j], cuboid(ENDo1t[j], ENDo2t[j], color = blue)) end do; -1
for j to nplot do `assign`(CUBlc[j], cuboid(ENDl1t[j], ENDl2t[j], color = red)); `assign`(CUBoc[j], cuboid(ENDo1t[j], ENDo2t[j], color = blue)) end do; -1
 

Translate and Rotation the basic elements 

Ground and the Cuboids 

Create the plots of the translated cuboids and the ground. 

The plot of the upper cuboid must be translated by the constant distance between both bodies in the initial position, because all motions are calculated with respect of the bodies' center of gravity. Finally, the cuboids are rotated around the rotation vectors. 

> `assign`(G[1], GROUND)
 

> `assign`(CUBl[1], CUBlc[1]); -1`assign`(CUBo[1], translate(CUBoc[1], DSlo[1], DSlo[2], DSlo[3])); -1
 

> for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
for j from 2 to nplot do `assign`(G[j], translate(GROUND, VectorCalculus:-`*`(xg[j], norml), VectorCalculus:-`*`(yg[j], norml), VectorCalculus:-`*`(zg[j], norml))); `assign`(CUBl[j], rotate(CUBlc[j], ...
 

Spring Connection Positions 

Next, we calculate the spring connections between the cuboids. We use the spring deformation relationships derived above in the calculation. 

 

> `assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSgx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(x0x[j], norml), a[n, x]...
 

 

> `assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
 

 

> `assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FTlx[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(eva...
 

 

> `assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
`assign`(iplot, 0); -1; for j by ds to ndt do `assign`(iplot, VectorCalculus:-`+`(iplot, 1)); `assign`(i, iplot); `assign`(FSox[i], [seq(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(Vec...
 

Connecting the Bodies 

For the animation the springs are plotted by lines. 

For the lower end of the lower springs the end points are the connection between the lower springs and the ground, for the upper end of the lower springs the end points are the connection between the lower springs and the lower cuboid. 

> for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
for j to nplot do `assign`(SL[j], line([FSgx[j][1], FSgy[j][1], FSgz[j][1]], [FSlx[j][1], FSly[j][1], FSlz[j][1]], color = green, thickness = 3), line([FSgx[j][2], FSgy[j][2], FSgz[j][2]], [FSlx[j][2]...
 

For the lower end of the upper springs the end points are the connection between the upper springs and the lower cuboid, for the upper end of the upper springs the end points are the connection between the upper springs and the upper cuboid. 

> for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
for j to nplot do `assign`(SO[j], line([FSox[j][1], FSoy[j][1], FSoz[j][1]], [FTlx[j][1], FTly[j][1], FTlz[j][1]], color = green, thickness = 3), line([FSox[j][2], FSoy[j][2], FSoz[j][2]], [FTlx[j][2]...
 

Displaying the Animation 

> `assign`(ANIM, seq(display({G[j], SL[j], SO[j], CUBl[j], CUBo[j]}), j = 1 .. nplot)); -1
 

> display(ANIM, insequence = true, style = patch, scaling = constrained, axes = none, orientation = [50, 80])
display(ANIM, insequence = true, style = patch, scaling = constrained, axes = none, orientation = [50, 80])
 

>  
 

References  

[1] E. Luz, "Schwingungsprobleme im Bauwesen", (Expert Verlag, 1992) 

[2] K.-J. Bathe; E. L. Wilson, "Numerical Methods In Finite Elemente Analysis", (Prentice-Hall, 1976) 


Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.

 

Image