University of Florida/Eml4507/Team 7 Report 5

From testwiki
Jump to navigation Jump to search

Problem 1

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

Solve the general eigenvalue problem for the spring-mass-damper system given on p.53-13 of 'Fead.s13.sec.53b', compare results obtained in Pb-53-6 and verify mass orthogonality of eigenvectors.

Given

M=(3002)

K=(30202035)

Solution

Solving the general eigenvalue problem

Kx=λMx
λMxKx=0
(λMK)x=0

λMK=(3λ3020202λ35)

det(λMK)=0

Calculating the determinant
(3λ30)(2λ35)400=0
6λ2165λ+1050=0

Therefore, the Eigenvalues are: λ1=17.5,λ2=10

The next step is finding the Eigenvectors
λ1=17.5
(λ1MK)X1=(22.520200)(x11x21)=(00)

Letting x21=1 and calculating x12

 (x11x21)=(0.8891) 


λ1=10
(λ1MK)X1=(0202015)(x12x22)=(00)

Letting x12=1 and calculating x22

 (x12x22)=(11.333) 


|π—πŸπ—πŸ|=[x11x12x21x22]=[0.889111.333]


These results compare favorably with those from problem 5.6 using CALFEM.

Mass Orthogonality: xiTMxj=0
The eigenvectors are therefore proven to be mass-orthogonal.

Problem 2

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

Do all work with your own matlab FE code, then verify the results with CALFEM; display both results clearly for easy comparison.

Solve the generalized eigenvalue problem of the above truss. Display the results for the lowest 3 eigenpairs. Plot and animate the lowest 3 mode shapes. Verify the mass orthogonality of these 3 eigenvectors.

Now consider the same truss, but with 2 missing braces:

Solve the generalized eigenvalue for this truss, and plot the mode shapes with zero frequency.

Given

L23=L45=1,L12=L24=L46=1,A=1/2,E=5,ρ=2

Solution

The lumped mass matrix, MΒ― , is shown below:

>> M=[1 0 0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 1];

K=zeros(10,10)

Because they are identical, we can use the CALFEM function β€œassem” for elements 2 and 6:

>> edof=[2 1 2 5 6; 6 3 4 9 10]; >> l=cos(45); m=sin(45); >> Ke=1.768*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   0.4879    0.7903   -0.4879   -0.7903
   0.7903    1.2801   -0.7903   -1.2801
  -0.4879   -0.7903    0.4879    0.7903
  -0.7903   -1.2801    0.7903    1.2801

>> [K]=assem(edof,K,Ke); For elements 1,4,7, and 10: >> edof=[1 1 2 3 4; 4 5 6 9 10; 7 3 4 7 8; 10 7 8 11 12]; >> l=1; >> m=0; >> Ke=2.5*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   2.5000         0   -2.5000         0
        0         0         0         0
  -2.5000         0    2.5000         0
        0         0         0         0

>> [K]=assem(edof,K,Ke);

For the elements 3 and 8:

>> edof=[3 3 4 5 6; 8 7 8 9 10]; l=0; m=1; Ke=2.5*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

        0         0         0         0
        0    2.5000         0   -2.5000
        0         0         0         0
        0   -2.5000         0    2.5000

>> [K]=assem(edof,K,Ke);

For elements 5 and 9:

>> edof=[5 5 6 7 8; 9 9 10 11 12]; >> l=cos(135); >> m=sin(135); >> Ke=1.768*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   1.7542   -0.1556   -1.7542    0.1556
  -0.1556    0.0138    0.1556   -0.0138
  -1.7542    0.1556    1.7542   -0.1556
   0.1556   -0.0138   -0.1556    0.0138

>> [K]=assem(edof,K,Ke)

K =

 Columns 1 through 4
   2.9879    0.7903   -2.5000         0
   0.7903    1.2801         0         0
  -2.5000         0    5.4879    0.7903
        0         0    0.7903    3.7801
  -0.4879   -0.7903         0         0
  -0.7903   -1.2801         0   -2.5000
        0         0   -2.5000         0
        0         0         0         0
        0         0   -0.4879   -0.7903
        0         0   -0.7903   -1.2801
        0         0         0         0
        0         0         0         0
 Columns 5 through 8
  -0.4879   -0.7903         0         0
  -0.7903   -1.2801         0         0
        0         0   -2.5000         0
        0   -2.5000         0         0
   4.7421    0.6347   -1.7542    0.1556
   0.6347    3.7939    0.1556   -0.0138
  -1.7542    0.1556    6.7542   -0.1556
   0.1556   -0.0138   -0.1556    2.5138
  -2.5000         0         0         0
        0         0         0   -2.5000
        0         0   -2.5000         0
        0         0         0         0
 Columns 9 through 12
        0         0         0         0
        0         0         0         0
  -0.4879   -0.7903         0         0
  -0.7903   -1.2801         0         0
  -2.5000         0         0         0
        0         0         0         0
        0         0   -2.5000         0
        0   -2.5000         0         0
   4.7421    0.6347   -1.7542    0.1556
   0.6347    3.7939    0.1556   -0.0138
  -1.7542    0.1556    4.2542   -0.1556
   0.1556   -0.0138   -0.1556    0.0138

In order to solve the generalized eigenvalue problem we used the function β€œeigen”, which returns the eigenvalues, λ, and the eigenvectors, β€˜x’:

Kx=Ξ»Mx >> [L,X]=eigen(K,M) L =

   0.0000
   0.0000
   0.0043
   0.3370
   1.4874
   2.1472
   3.9323
   5.0164
   6.2876
   6.9486
   7.2731
  10.7101


X =

 Columns 1 through 4
  -0.4065   -0.0382    0.0125   -0.2771
   0.0382   -0.4065   -0.2166    0.6103
  -0.4065   -0.0382    0.0122   -0.1890
   0.0382   -0.4065   -0.1923    0.0564
  -0.4065   -0.0382   -0.0217    0.1776
   0.0382   -0.4065   -0.1948    0.1689
  -0.4065   -0.0382    0.0102   -0.0012
   0.0382   -0.4065   -0.1524   -0.4836
  -0.4065   -0.0382   -0.0412    0.1877
   0.0382   -0.4065   -0.1539   -0.4109
  -0.4065   -0.0382    0.0280    0.1020
   0.0382   -0.4065    0.9099    0.0588
 Columns 5 through 8
  -0.5124   -0.3701   -0.1840    0.2470
   0.2249   -0.5457    0.0289   -0.0011
  -0.1249   -0.3416    0.1335   -0.2501
  -0.3404    0.4774   -0.0870   -0.2533
   0.2127    0.0780   -0.5547   -0.1457
  -0.4841    0.0930    0.1690    0.2458
   0.1231   -0.0039    0.3146   -0.2511
   0.4246   -0.0295   -0.0415    0.5429
   0.1230    0.3502   -0.3333    0.2405
   0.1826    0.0003   -0.0315   -0.5382
   0.1785    0.2873    0.6240    0.1593
  -0.0076    0.0046   -0.0379    0.0040
 Columns 9 through 12
   0.3808   -0.2946   -0.0054    0.1736
   0.2043   -0.0207    0.1439    0.0439
  -0.2599    0.4888    0.2687   -0.4541
   0.0756   -0.2008    0.5664   -0.1042
  -0.2724   -0.3927   -0.2785   -0.3392
  -0.3957    0.1521   -0.5052   -0.0067
  -0.4206   -0.2130    0.1168    0.6482
  -0.2138   -0.0829    0.2081   -0.0524
   0.2109    0.5935    0.0403    0.3022
   0.3342    0.1351   -0.4180    0.1104
   0.3612   -0.1820   -0.1419   -0.3307
  -0.0045    0.0171    0.0047    0.0091

For the first eigenpair, the eigenvalue and corresponding eigenvector are as follows:

L= 0 X= {-0.406457927237019 0.0381918846465145 -0.406457927237019 0.0381918846465150 -0.406457927237019 0.0381918846465148 -0.406457927237019 0.0381918846465148 -0.406457927237019 0.0381918846465148 -0.406457927237020 0.0381918846465028}

These eigenvectors were used to calculate the displacement of each node on the truss and plotted using the MATLAB code below: % truss original position % start nodes (x/y) / end nodes (x/y)

sn=[     0     0
         0     0
         1     0
         1     0
         1     0
         1     1
         1     1
         2     0
         2     0
         2     1];
en=[     1     0
         1     1
         1     1
         2     0
         2     1
         2     0
         2     1
         2     1
         3     0
         3     0];

% the engine

    x=[sn(:,1),en(:,1)].'; % <- note: transpose is important!
    y=[sn(:,2),en(:,2)].';
    hold on
    plot(x,y,'b','linewidth',4);
    axis([-0.5 3.5 -0.5 1.5])
    
    % start nodes (x then y) and end nodes (x then y) in 10x2 matrices
sn=[     0 0 %1
         0 0 %1
         0.89354 0.0381918846465150 %2
         0.89354 0.0381918846465150%2
         0.89354 0.0381918846465150%2
         0.89354 1.038192 %3
         0.89354 1.038192 %3
         1.89354 0.038192%4
         1.89354 0.038192 %4
         1.89354 1.038192]; %5
en=[     0.89354 0.0381918846465150 %2
         0.89354 1.038192  %3
         0.89354 1.038192  %3
         1.89354 0.038192%4
         1.89354 1.038192%5
         1.89354 0.038192%4
         1.89354 1.038192%5
         1.89354 1.038192%5
         2.89354 0 %6
         2.89354 0]; %6

% the engine

    x1=[sn(:,1),en(:,1)].'; % <- note: transpose is important!
    y1=[sn(:,2),en(:,2)].';
    plot(x1,y1,'r','linewidth',3);
    axis([-0.5 3.5 -0.5 1.5])
    xlabel('x')
    ylabel('y')
    title('Mode Shape 1')

File:Modeshape1.jpg


For the second eigenpair:

L=0 X={-0.0381918846465119 -0.406457927237036 -0.0381918846465119 -0.406457927237034 -0.0381918846465145 -0.406457927237035 -0.0381918846465120 -0.406457927237031 -0.0381918846465160 -0.406457927237031 -0.0381918846465106 -0.406457927236949}

File:Modeshape2 2.jpg

For the third eigenpair:

L=0.0043 X={0.0124740073031415 -0.216589482612068 0.0122217899111943 -0.192300781266623 -0.0216971729033128 -0.194763453030673 0.0102232821955044 -0.152379327406532 -0.0411885014874805 -0.153869406652153 0.0279665949809443 0.909902450968240}

File:Modeshape3 1.jpg

Analyzing the second truss, we get the lumped mass matrix, MΒ― , shown below: >> M=[1 0 0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 1];

For element 2 we use the CALFEM function β€œassem”: >> edof=[2 1 2 5 6]; >> l=cos(45); >>m=sin(45); >> Ke=1.768*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   0.4879    0.7903   -0.4879   -0.7903
   0.7903    1.2801   -0.7903   -1.2801
  -0.4879   -0.7903    0.4879    0.7903
  -0.7903   -1.2801    0.7903    1.2801

>> [K]=assem(edof,K,Ke);

For elements 1,4,7, and 6:

>> edof=[1 1 2 3 4; 4 5 6 9 10; 7 3 4 7 8; 6 7 8 11 12]; >> l=1; >> m=0; >> Ke=2.5*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   2.5000         0   -2.5000         0
        0         0         0         0
  -2.5000         0    2.5000         0
        0         0         0         0

>> [K]=assem(edof,K,Ke);

For the elements 3 and 8:

>> edof=[3 3 4 5 6; 8 7 8 9 10]; l=0; m=1; Ke=2.5*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

        0         0         0         0
        0    2.5000         0   -2.5000
        0         0         0         0
        0   -2.5000         0    2.5000

>> [K]=assem(edof,K,Ke); For element 5:

>> edof=[5 5 6 7 8]; >> l=cos(135); >> m=sin(135); >> Ke=1.768*[l^2 l*m -l^2 -l*m; l*m m^2 -l*m -m^2; -l^2 -l*m l^2 l*m; -l*m -m^2 l*m m^2]

Ke =

   1.7542   -0.1556   -1.7542    0.1556
  -0.1556    0.0138    0.1556   -0.0138
  -1.7542    0.1556    1.7542   -0.1556
   0.1556   -0.0138   -0.1556    0.0138

>> [K]=assem(edof,K,Ke)

>> K

K =

 Columns 1 through 4
   2.9879    0.7903   -2.5000         0
   0.7903    1.2801         0         0
  -2.5000         0    5.0000         0
        0         0         0    2.5000
  -0.4879   -0.7903         0         0
  -0.7903   -1.2801         0   -2.5000
        0         0   -2.5000         0
        0         0         0         0
        0         0         0         0
        0         0         0         0
        0         0         0         0
        0         0         0         0
 Columns 5 through 8
  -0.4879   -0.7903         0         0
  -0.7903   -1.2801         0         0
        0         0   -2.5000         0
        0   -2.5000         0         0
   2.9879    0.7903         0         0
   0.7903    3.7801         0         0
        0         0    5.0000         0
        0         0         0    2.5000
  -2.5000         0         0         0
        0         0         0   -2.5000
        0         0   -2.5000         0
        0         0         0         0
 Columns 9 through 12
        0         0         0         0
        0         0         0         0
        0         0         0         0
        0         0         0         0
  -2.5000         0         0         0
        0         0         0         0
        0         0   -2.5000         0
        0   -2.5000         0         0
   2.5000         0         0         0
        0    2.5000         0         0
        0         0    2.5000         0
        0         0         0         0

In order to solve the generalized eigenvalue problem we used the function β€œeigen”, which returns the eigenvalues, λ, and the eigenvectors, β€˜x’:

Kx=Ξ»Mx >> [L,X]=eigen(K,M)

L =

        0
   0.0000
   0.0000
   0.0000
   0.0000
   1.1650
   2.1580
   5.0000
   5.0000
   5.0000
   6.5940
   8.6191


X =

 Columns 1 through 4
        0   -0.0507    0.0194    0.4831
        0   -0.3353    0.1932   -0.2386
        0   -0.0507    0.0194    0.4831
        0   -0.4522    0.3724    0.0226
        0    0.1385   -0.2708    0.0601
        0   -0.4522    0.3724    0.0226
        0   -0.0507    0.0194    0.4831
        0    0.4637    0.5182    0.0235
        0    0.1385   -0.2708    0.0601
        0    0.4637    0.5182    0.0235
        0   -0.0507    0.0194    0.4831
   1.0000         0         0         0
 Columns 5 through 8
   0.0167   -0.4448    0.4308   -0.0000
   0.4831    0.4811    0.5118    0.0000
   0.0167   -0.0991    0.3317   -0.0000
   0.1279   -0.3136   -0.4503    0.3241
   0.5920   -0.1034   -0.0380    0.5250
   0.1279   -0.1675   -0.0616   -0.3241
   0.0167    0.2928   -0.0537   -0.0000
   0.1262   -0.0000    0.0000    0.3454
   0.5920   -0.1936   -0.2780   -0.5250
   0.1262    0.0000    0.0000   -0.3454
   0.0167    0.5482   -0.3928   -0.0000
        0         0         0         0
 Columns 9 through 12
  -0.0612    0.4647    0.2458   -0.3159
  -0.0000    0.0000    0.2514   -0.0574
   0.0612   -0.4647    0.0067    0.6510
   0.2072   -0.1840    0.3943   -0.0396
   0.2745    0.1666   -0.3986    0.0599
  -0.2072    0.1840   -0.6456    0.0970
   0.0612   -0.4647   -0.2501   -0.6265
  -0.6117   -0.0806    0.0000   -0.0000
  -0.2745   -0.1666    0.2434   -0.0245
   0.6117    0.0806   -0.0000    0.0000
  -0.0612    0.4647    0.1528    0.2560
        0         0         0         0

For the first eigenpair, the eigenvalue and corresponding eigenvector are as follows:

L= 0 X= {0 0 0 0 0 0 0 0 0 0 0 1}

The truss was plotted in MATLAB using the following code:

% truss original position % start nodes (x/y) / end nodes (x/y)

sn=[     0     0
         0     0
         1     0
         1     0
         1     1
         2     0
         2     0
         2     1];
en=[     1     0
         1     1
         1     1
         2     0
         2     1
         2     1
         3     0
         3     0];

% the engine

    x=[sn(:,1),en(:,1)].'; % <- note: transpose is important!
    y=[sn(:,2),en(:,2)].';
    hold on
    plot(x,y,'b','linewidth',4);
    axis([-0.5 3.5 -0.5 1.5])
    
    % start nodes (x then y) and end nodes (x then y) in 10x2 matrices
sn=[     0 0 %1
         0 0 %1
         1 0 %2
         1 0%2
         1 1%6-3
         2 0%4
         2 0%4
         2 1]; %5
en=[     1 0 %2
         1 1 %3
         1 1 %3
         2 0%4
         2 1%5
         2 1%5
         3 0%6
         3 0];%6

% the engine

    x1=[sn(:,1),en(:,1)].'; % <- note: transpose is important!
    y1=[sn(:,2),en(:,2)].';
    plot(x1,y1,'r','linewidth',3);
    axis([-0.5 3.5 -0.5 1.5])
    xlabel('x')
    ylabel('y')
    title('Mode Shape 1')

File:Modeshape1 2.png

Problem 3

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

Find, plot and compare the eigenvectors.
Compare the mode shapes using two different assumptions.
Animate each format in the given sine wave.

Given

|π—πŸπ—πŸ|=[x11x12x21x22]

(2.1)

Solve for the eigenvectors under two differing assumptions:

x11=x12=1
x21=x22=1

(2.2)

Plot mode shapes according to y=𝐗𝐒sin(wit)

where wit is the circular frequency, wit=γi,

according to the eigenvalues γ1=4+5 and γ2=45.

Solution

Firstly, the eigenvectors must be found.
The x21=x22=1 assumption is in effect.
This value was found with γ1=4+5

[kγ1I]𝐗=[0]

(2.3)

[1+5221+5][x11]=[00]

(2.4)

2x1+1+5=0
x1=(1+5)/2

(2.5)


and this value was found with γ2=45

[kγ2I]𝐗=[0]

(2.6)

[152215][x11]=[00]

(2.7)

2x1+15=0
x1=(15)/2

(2.8)

With these values solved for, construct the eigenvector:

[π—πŸπ—πŸ]=[(1+5)/2(15)/211]


(2.9)

Now assuming x11=x12=1 the eigenvectors are found using the same eigenvalues γ1,2=4±5

[kγ1I]𝐗=[0]

(2.10)

[1+5221+5][1x2]=[00]

(2.11)

1+52x2=0
x2=(1+5)/2

(2.12)

and this value was found with γ2=45

[kγ2I]𝐗=[0]

(2.13)

[152215][1x2]=[00]

(2.14)

152x2=0
x2=(15)/2

(2.15)

With these values solved for, construct the eigenvector:

[π—πŸπ—πŸ]=[11(1+5)/2(15)/2]

(2.16)

For both assumptions, the mode shape are plotted, not according to scale.

[Row 2 = 1 Assumption] - Mode 1


[Row 2 = 1 Assumption] - Mode 2


[Row 1 = 1 Assumption] - Mode 1


[Row 1 = 1 Assumption] - Mode 2

Comparing the assumptions:
The constant 1 switches to the other point, upon switching the Row = 1 assumption.
Going from mode 1 to mode 2 within an assumption leaves the solved variable with a -2.236 reduction in magnitude.
Going from one assumption to another, within the same mode, reduces the magnitude of the solved variable by 1.


Modes animated according to y=𝐗𝐒sin(wit) where wit is the circular frequency, wit=γi, according to the eigenvalues γ1=4+5 and γ2=45:

M1 =[

  1.618 -.618;
    1       1;
];
%quiver(0,M1(1,1));
%hold on;
%quiver(0,M1(2,1));

M2 =[

   1      1;
 0.618  -1.618;
];
%quiver(0,0,M2(1,1),M2(2,1));
%hold on;
%quiver(0,0,M2(1,2),M2(2,2));


%% Sine wave:
X=[1 1;
   (-1+sqrt(5))/2 (-1-sqrt(5))/2];
la=[4+sqrt(5);4+sqrt(5)];
ex=[0 3];
y=[0;0];

for kk=1:2
    eigenvalx=la(kk);
    eigenvecx=X(:,kk);
    w=sqrt(eigenvalx);
    T=2*pi/w;
    dt=T/100;
    t=0:dt:T;
    figure(kk)
    
    if kk==1;
        filename = 'mode1_3.gif';
    elseif kk==2;
        filename = 'mode2_3.gif';
    end
    
    for jj = 1:100
        tt=t(jj);
        d=eigenvecx*sin(w*tt);
        d=d';
        xd=(ex+d)';
        s1=['o' , 'k'];
        axis manual
        plot(xd,y,s1)
        set(gca, 'ylim', [-1.5 1.5], 'xlim', [-3 5]);
        drawnow
        frame = getframe(1);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        if jj == 1;
            imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
        else
            imwrite(imind,cm,filename,'gif','WriteMode','append');
        end
    end
end
File:R5p3 mode1 top.gif
Mode 1 Animation


File:R5p3 mode2 top.gif
Mode 2 Animation

Problem 4

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

The Lowest 3 eigenpairs (ωj,ϕj)forj=1,2,3 and plot the 3 lowest mode shapes in a gif image.

Given

File:25mt.PNG
Figure 6.1: 25 member truss with L=0.3 m

E=100GPa,A=1.0cm2,L=0.3m,ρ=5000kgm3

Solution

 
%GIVENS
L=0.3;
E=100e9;
A=1e-4;
rho=5000;
ep=[E A];

%ELEMENT DOF
Edof=zeros(25,5);
for i=1:6
    j=2*(i-1);
    Edof(i,:)=[i,1+j,2+j,3+j,4+j];
end
for i=7:12
    j=2*(i+1);
    Edof(i,:)=[i,j-1,j,j+1,j+2];
end
for i=13:19
    j=2*(i-13);
    Edof(i,:)=[i,1+j,2+j,15+j,16+j];
end
for i=20:25
    j=2*(i-20);
    Edof(i,:)=[i,1+j,2+j,17+j,18+j];
end

%COORDINATES
ex=zeros(14,2);
ey=zeros(14,2);
for i=1:25
    if i<7
        ex(i,:)=[L*(i-1),L*i];
        ey(i,:)=[0,0];
    elseif i>6&&i<13
        ex(i,:)=[L*(i-7),L*(i-6)];
        ey(i,:)=[L,L];
    elseif i>12&&i<20
        ex(i,:)=[L*(i-13),L*(i-13)];
        ey(i,:)=[0,L];
    else
        ex(i,:)=[L*(i-20),L*(i-19)];
        ey(i,:)=[0,L];
    end
end

%STIFFNESS AND MASS MATRIX
K=zeros(28);
M=zeros(28);
for i=1:25
    E=ep(1);A=ep(2);
    xt=ex(i,2)-ex(i,1);
    yt=ey(i,2)-ey(i,1);
    L=sqrt(xt^2+yt^2);
    l=xt/L; m=yt/L;
    ke=E*A/L*[l^2 l*m -l^2 -l*m;
              l*m m^2 -l*m -m^2;
              -l^2 -l*m l^2 l*m;
              -l*m -m^2 l*m m^2;];
    m=L*A*rho;
    me=[m/2 0 0 0;
        0 m/2 0 0;
        0 0 m/2 0;
        0 0 0 m/2];
    edoft=Edof(i,2:5);
    K(edoft,edoft)=K(edoft,edoft)+ke;
    M(edoft,edoft)=M(edoft,edoft)+me;
end

%FIND THE EIGEN VALUES AND VECTORS
nd=size(K,1);
fdof=[1:nd]';
b=[1 15 16]';
fdof(b(:))=[];
[X,D]=eig(K(fdof,fdof),M(fdof,fdof));
nfdof=size(X,1);
eigenval=diag(D);
eigenvec=zeros(nd,nfdof);
eigenvec(fdof,:)=X;

%PLOT THE GIF
for kk=1:3
    eigenvalx=eigenval(kk);
    eigenvecx=eigenvec(:,kk);
    w=sqrt(eigenvalx);
    T=2*pi/w;
    dt=T/100;
    t=0:dt:T;
    figure(kk)
    
    if kk==1;
        filename = 'mode1.gif'
    elseif kk==2;
        filename = 'mode2.gif'
    elseif kk==3;
        filename = 'mode3.gif'
    end
    
    for jj = 1:100
        tt=t(jj);
        d=eigenvecx*sin(w*tt);
        d=d';
        rr=Edof(:,2:5);
        ed=zeros(25,4);
        for i=1:25;
            ed(i,1:4)=d(rr(i,:));
        end
        xd=(ex+ed(:,[1 3]))';
        yd=(ey+ed(:,[2 4]))';
        s1=['-' , 'k'];
        axis manual
        plot(xd,yd,s1)
        set(gca, 'ylim', [-1.5 1.5], 'xlim', [0 3]);
        drawnow
        frame = getframe(1);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        if jj == 1;
            imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
        else
            imwrite(imind,cm,filename,'gif','WriteMode','append');
        end
    end
end

The resulting egienpairs are:

ω1=1.6890×105ϕ1=[00.00820.02790.07570.04790.19820.06050.35810.06700.53840.06900.72440.06900.9052000.03640.06770.06450.19060.08470.35140.09740.53320.10380.72130.10580.9045]ω2=2.6676×106ϕ2=[00.07450.02270.41760.09730.65120.18840.62520.26030.31370.29410.18250.29760.6937000.07320.35410.07140.61490.01530.62280.06030.33680.12010.15530.14550.6853]ω3=7.0219×106ϕ3=[00.03520.21780.04160.39070.10230.51910.12830.61150.07510.67220.02530.69410.0971000.12040.00890.26560.07800.42210.11950.56630.08090.67100.01540.71780.0940]

The mode shapes for the previous eigenpairs are"

File:Mode1 team7.gif
Figure 4.2: First eigen shape
File:Mode2 team7.gif
Figure 4.3: Second eigen shape
File:Mode3 team7.gif
Figure 4.4: Third eigen shape

Problem 5

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

Do the HW problem on p.21-1 (2-bar truss, unconstrained, no supports, use 6x6 stiffness matrix ) and the HW problem on p.21-3 (4-bar truss, with supports, use reduced stiffness matrix ). For problem on p. 21-1, plot the eigenvectors corresponding to zero evals of 2 bar truss system and interpret results. For problem on p. 21-3, plot eigenvectors corresponding to zero evals for case (a).

Given

The HW problems on p.21-1 and p.21-3.

Solution

For problem 21-1:



K x u = 0 x u = 0

Σα x u = W

K x W = K x (Σα x u) = Σα x (K x u) = Σα x 0 = 0



For problem 21-3: K x v = Ξ» x v



where a = b = 1m, E=2, and A=3

Problem 6

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

File:MDOF.JPG

Using the method of Modal Superposition, find the governing equations of motion for the system shown above.

Given

m1=3,m2=2
c1=1/2,c2=1/4,c3=1/3
k1=10,k2=20,k3=15
F1(t)=0,F2(t)=0
d1(0)=1,d2(0)=2
d1(0)=0,d2(0)=0

Solution

First, with the information given we can substitute them into the equation with the form:

𝐌dΒ¨+𝐂dΛ™+𝐊d=𝐅(t)

Where M, C, and K are the mass, damper and stiffness matrices.

𝐌=[m100m2]

𝐂=[(c1+c2)c2c2(c2+c3)]

𝐊=[(k1+k2)k2k2(k2+k3)]

Substituting the known values, we now have the equation as matrix form as:

[3002]dΒ¨+[3/41/41/47/12]dΛ™+[30202035]d=𝐅(t)

Eq. 1

We can now guess a solution that will solve Equation 1. In this case the following will be used:

d(t)=uβ†’ejωt

d(t)Λ™=jωuβ†’ejωt

d(t)Β¨=j2ω2uβ†’ejωt

Where u→ represents the eigen vectors and j=(1) and u→0

The equation can now be written as:

𝐌(j2ω2uβ†’ejωt)+𝐂jωuβ†’ejωt+𝐊uβ†’ejωt=𝐅(t)

Eq. 2

The following matrix can be used in order to represent the modes of vibration of the system by its eigenvalues and eigenvectors:

[M]1[K]=[3002]1[30202035]=[106.6671017.5]

Where the inverse matrix of M is determined by:

[M]1=adj(M)det(M)


To find the solutions to the system we use the guesses that we made (after plugging them back into the matrix equation, Equation 2) and obtain the following worked out equation:

(ω2[M]1[K])d0

When substituting we obtain:

(ω2[106.6671017.5])d0

Eq. 3

And d0 are the eigen vectors or mode shapes.

To find ω we require the determinant for Equation 3:

det(ω2[106.6671017.5])=det[ω2106.66710ω217.5]

Eq. 4

Taking the determinant, we now have the equation:

ω427.5ω2+175

Eq. 5

By using the quadratic formula we can find the respective ω values.

ω2=b±b24ac2a

Eq. 6

And now,ω12=17.5 and ω22=10
Using these values we can now calculate the mode shapes by first substituting them into the matrix in Equation 4 for each ω value.
For ω12=17.5 we substitute and obtain:

[296.256.66710288.75]{d1d2}={00}

Eq. 7

solving for d1:

296.25d1+6.667d2=0

and

10d1+288.75d2=0

Form these equations we obtain:

d1={.022528.875}

Now we substitute ω22=10 into the matrix in equation 4 and repeat the above steps to obtain d2 we have:

[906.6671082.5]{d1d2}={00}

Corresponding equations are now:

90d1+6.667d2=0

and

10d1+82.5d2=0

The modal vector is now :

d2={13.5.1212}

The modal matrix can now be written as:

ϕ=[.022513.528.875.1212]

Now the method of modal superposition can be implemented by the following equation:

{d1d2}[.022513.528.875.1212]={z1z2}

Where z1 and z2 are the modal coordinates.
For mass matrix:

[ϕ]T[M][ϕ]=[.022528.87513.5.1212][3002][.022513.528.875.1212]=[1667.57.97.9546.8]

For damping matrix:

[ϕ]T[C][ϕ]=[.022528.87513.5.1212][3/41/41/47/12][.022513.528.875.1212]=[486.038895.1895.18135.878]

For Stiffness matrix:

[ϕ]T[K][ϕ]=[.022528.87513.5.1212][30202035][.022513.528.875.1212]=[29156766576655403]

With the results the equation of motion for the system can be represented by modal coordinates in two 1-degree of freedom equations, with modal superposition we obtain:

1667.5z1Β¨+486.0388z1Λ™+29156z1=0

Eq. 8

and

546.8z2Β¨+135.878z2Λ™+5403z2=0

Eq. 8

*Reference: http://mecsys.mecc.polimi.it/Courses/msd_lc/ta/Lesson3.pdf
 The above link was used as a guide to help me work through this problem. The steps are similar, however, some of the steps were       
explained in more detail than the reference. Derivations and variables were explained as well.

Problem 7

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Find

Given

Solution

References


Contributing Team Members

Problem Assignments
Problem # Solved & Typed by Reviewed by
1 Zeyn Kermani All
2 Spencer Herran All
3 Joshua Plicque All
4 Matthew Gidel All
5 Kristin Howe All
6 Kevin Li All
7 Brandon Wright All

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Template:CourseCat