View source for University of Florida/Eml4507/Fead-s13-team4-R2

==Problem R2.1 Model 2-D Truss With 2 Inclined Elements and Verify Assembly==
===Given 2-D Truss Model and Constants===
Consider the 2D truss system

[[image:R2.1.jpg|frame|center|Figure 1: 2D Truss system analyzed in this problem. ([http://en.wikiversity.org/wiki/File:R2.1.jpg])]]


Use the following constants to very calculations

<math>l^{(1)}=4</math>

<math>l^{(2)}=2</math>

<math>E^{(1)}=3</math>

<math>E^{(2)}=5</math>

<math>A^{(1)}=1</math>

<math>A^{(2)}=2</math>

<math>\O ^{(1)}=30^{\circ}</math>

<math>\O ^{(2)}=135^{\circ}</math>

<math>P=7</math>





===Find:===
====1) Explain the Assembly Process of the Above 2D Truss System====
====2) Verify With CALFEM====
===Solution: Assembly Process for 2d Truss System===
====Construct local stiffness matrices====
The stiffness matrices will be constructed for individual elements 1 and 2. The rows and columns will be numbered to represent the force and degree of freedom. 
For element 1

<math>
\vec{K^{(1)}}=\begin{bmatrix}
&1&2&3&4\\
1&(l^{(1)} )^{2} & l^{(1)}m^{(1)}  & -(l^{(1)} )^{2} & -l^{(1)}m^{(1)}\\ 
2&l^{(1)}m^{(1)} & (m^{(1)} )^{2} & -l^{(1)}m^{(1)} & -(m^{(1)} )^{2}\\ 
3&-(l^{(1)} )^{2}& -l^{(1)}m^{(1)} & (l^{(1)} )^{2} & l^{(1)}m^{(1)}\\ 
4&-l^{(1)}m^{(1)} & -(m^{(1)} )^{2} & l^{(1)}m^{(1)} & (m^{(1)} )^{2}\\ 
\end{bmatrix}*K^{(1)}
</math>

For element 2

<math>
\vec{K^{(2)}}=\begin{bmatrix}
&3&4&5&6\\
3&(l^{(2)} )^{2} & l^{(2)}m^{(2)}  & -(l^{(2)} )^{2} & -l^{(2)}m^{(2)}\\ 
4&l^{(2)}m^{(2)} & (m^{(2)} )^{2} & -l^{(2)}m^{(2)} & -(m^{(2)} )^{2}\\ 
5&-(l^{(2)} )^{2}& -l^{(2)}m^{(2)} & (l^{(2)} )^{2} & l^{(2)}m^{(2)}\\ 
6&-l^{(2)}m^{(2)} & -(m^{(2)} )^{2} & l^{(2)}m^{(2)} & (m^{(2)} )^{2}\\ 
\end{bmatrix}*K^{(2)}
</math>

The variable K can be determined by

<math>K^{(e)}=\frac{EA}{L}</math>

If the variables are plugged in, we get the following matrices:

For element 1

<math>
\vec{K^{(1)}}=\begin{bmatrix}
&1&2&3&4\\
1&0.5625 & 0.3248  & -0.5625 & -0.3248\\ 
2&0.3248 & 0.1875 & -0.3248 & -0.1875\\ 
3&-0.5625& -0.3248 & 0.5625& 0.3248\\ 
4&-0.3248 & -0.1875 & 0.3248 & 0.1875\\ 
\end{bmatrix}
</math>

For element 2

<math>
\vec{K^{(2)}}=\begin{bmatrix}
&3&4&5&6\\
3&2.5&-2.5&-2.5&2.5\\ 
4&-2.5&2.5&2.5&-2.5\\ 
5&-2.5&2.5&2.5&-2.5\\ 
6&2.5&-2.5&-2.5&2.5\\ 
\end{bmatrix}
</math>


It is clear that if we wish to combine the matrices, it will place coordinate 3,3 of element 2 on the location of 3,3 of element 1 due to them referencing the same point in the same direction. When combining the matrices, addition will be used since the matrices represent directional forces and forces always combine in the same direction by addition. 

====Combining local stiffness matrices into a global matrix====
The local matrices from element 1 and 2 will be combined as described above to give

<math>
\vec{K}=\begin{bmatrix}
&1&2&3&4&5&6\\
1&0.5625 & 0.3248  & -0.5625 & -0.3248&0&0\\ 
2&0.3248 & 0.1875 & -0.3248 & -0.1875&0&0\\ 
3&-0.5625& -0.3248 & 3.0625& -2.1752&-2.5&2.5\\ 
4&-0.3248 & -0.1875 & -2.1752 & 2.6875&2.5&-2.5\\ 
5&0&0&-2.5&2.5&2.5&-2.5\\
6&0&0&2.5&-2.5&-2.5&2.5\\
\end{bmatrix}
</math>


===Solution: Verify With CALFEM===


====MatLab Code====
EDU>> Edof = [1  1 2 3 4

2  3 4 5 6];

EDU>> K = zeros(6);

EDU>> f = zeros(6,1);

EDU>> f(4) = -7;

EDU>> E1 = 3; E2 = 5;

EDU>> L1 = 4;L2 = 2;

EDU>> A1 = 1;A2 = 2;

EDU>> ep1 = [E1 A1]

ep1 =

     3     1


EDU>> ep2 = [E2 A2]

ep2 =

     5     2


EDU>> ex1 = [-3.464 0]; ex2 = [0 1.414];

EDU>> ey1 = [-2 0]; ey2 = [0 -1.414];

EDU>> Ke1 = bar2e(ex1,ey1,ep1)

Ke1 =

    0.5625    0.3248   -0.5625   -0.3248
    0.3248    0.1875   -0.3248   -0.1875
   -0.5625   -0.3248    0.5625    0.3248
   -0.3248   -0.1875    0.3248    0.1875


EDU>> Ke2 = bar2e(ex2,ey2,ep2)

Ke2 =

    2.5004   -2.5004   -2.5004    2.5004
   -2.5004    2.5004    2.5004   -2.5004
   -2.5004    2.5004    2.5004   -2.5004
    2.5004   -2.5004   -2.5004    2.5004


EDU>> K = assem(Edof(1,:),K,Ke1);

EDU>> K = assem(Edof(2,:),K,Ke2)

K =

    0.5625    0.3248   -0.5625   -0.3248         0         0
    0.3248    0.1875   -0.3248   -0.1875         0         0
   -0.5625   -0.3248    3.0629   -2.1756   -2.5004    2.5004
   -0.3248   -0.1875   -2.1756    2.6879    2.5004   -2.5004
         0         0   -2.5004    2.5004    2.5004   -2.5004
         0         0    2.5004    -2.5004  -2.5004    2.5004


EDU>> bc = [1 0;2 0;5 0;6 0];

EDU>> [a,r]=solveq(K,f,bc)

a =

         0
         0
   -4.3519
   -6.1268
         0
         0


r =

    4.4378
    2.5622
         0
   -0.0000
   -4.4378
    4.4378

===Conclusion===
By comparing the value from the hand calculation and the value from the CALFEM program, we can see that the two matrices are identical. This verifies the assembly process of using addition when combining local matrices by addition.


== Problem R2.2 ==

===Disccription:===
We are asked to consider the L2-ODE-CC system referenced in (1) p.53-2 to have the following complex conjugate roots:

<math> \lambda_{1,2} = -0.5 \pm 2i </math>

We are asked to find the homogeneous solution in standard form and to find the damping ratio, <math> \zeta </math>

We are then to consider the initial conditions

<math> y(0) = 1, y'(0) = 0 </math>

We are to plot this solution and then use that plot to find the log decrement and compare it to the log decrement that would be obtained using (6) p.53-5b.

Finally, we're to find the quality factor, Q, and the loss factor, <math>  \eta </math>




===Solution:===

From R1.5, we have already obtained that for the above conditions,
 
<math> y(t) = e^{-0.5 t} cos(2t) + (\frac{1}{4}) e^{-0.5 t} sin(2t) </math>

with coefficients to the characteristic equation of 

<math> a = 1</math> and <math> b = \frac{17}{4}</math>

From (1) and (2) p.53-5, we have that

<math> \omega^{2} = b </math>

<math> 2 \omega \zeta = a </math>

Knowing both a and b we get

<math> \omega = \sqrt{ \frac{17}{4}} = 2.0615 </math>

<math> \zeta = \frac{1}{2 \omega} = \frac{1}{2*2.0615} = 0.2425 </math>

We then have from equation (1) p.53-5d

<math> \delta = \frac{1}{3}[\delta^{(1)} + \delta^{(2)} + \delta^{(3)}]  </math>

[[File:Pr 5 third graph.jpg]]

Using the above graph of the homogeneous solution, accurate estimates were able to be made as to the various values of y at the peaks

Using equation (3) p.53-5c we are able to find

<math> \delta^{(1)} = log(1/0.2075) = 0.68235 </math><br>
<math> \delta^{(2)} = log(0.2075/0.0432) = 0.68216 </math><br>
<math> \delta^{(3)} = log(0.0432/0.009) = 0.68124 </math>

===Conclusion===

We then arrive at

<math> \delta = \frac{1}{3}[0.68235 + 0.68216 + 0.68124] = 0.68192  </math>

For comparison, the equation (6) p.53-5b gives

<math> \delta = \frac{2 \pi \zeta}{\sqrt{1 - \zeta^{2}}} = \frac{2 \pi 0.2425}{\sqrt{1 - 0.2425^{2}}} = 1.5705 </math>

We now find the quality factor, Q.  From (2) p.53-5d, we have

<math> Q = \frac{\delta}{\pi} = \frac{\sqrt{1 - \zeta^{2}}}{2 \zeta} = \frac{\sqrt{1 - 0.2425^{2}}}{2 *0.2425} = 2 </math>

From (5) p.53-5d, we know that the loss factor is just the inverse of the quality factor so

<math> \eta = \frac{1}{Q} = \frac{1}{2} = 0.5 </math>

==Problem R2.3==
=== Problem Statement ===
KS 2008 p.103 sec.2.7 pb.21
write a matlab program to assemble the global stiffness matrix, 
compute the unknown disp dofs, the reaction forces, the member forces; 
compare the results with exact hand calculation 
Statics, Mechanics of Materials).

run CALFEM to verify the results, but now do this problem 
in a different way (as in commercial FE codes). 

first, establish 2 arrays: 
(1) the global node coordinate array “coord”, and 
(2) the elem connectivity array “conn”.  next, write a 
matlab code to loop over the number of element to call 
the function “bar2e” of CALFEM to compute the element 
stiffness matrices.  with this method, you let your 
matlab code construct the arrays ex1, ey1, ex2, ey2, etc. 
for you, using the info from the arrays “coord” and “conn”.



=== Solution ===
=====MATLAB Code=====
The following is a MATLAB code that assembles the global stiffness matrix, computes the unknown displacement dofs, reaction forces, and member forces
    function FEA_R2p3
    E = 10000; % Young's modulus
    A = pi*2^2/4; % Area
    % Lengths of each element
    L1 = 12;
    L2 = 12;
    L3 = 12;
    L4 = sqrt((12^2)+(9^2));
    L5 = 9;
    L6 = sqrt((12^2)+(9^2));
    L7 = 12;
    L8 = 9;
    L9 = sqrt((12^2)+(9^2));
    % topology matrix
    % [element# node1x node1y node2x node2y]
    Edof = [1 1  2  3  4;2  3  4 5 6;3  5  6 7 8;
            4 1  2  9 10;5  9 10 3 4;6  9 10 5 6;
            7 9 10 11 12;8 11 12 5 6;9 11 12 7 8];
    % element stiffness matrix in local coordinates for elements 1 through 9
    esmloc1 = (E*A/L1)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
    esmloc2 = (E*A/L2)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
    esmloc3 = (E*A/L3)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
    esmloc4 = (E*A/L4)*[(L1/L4)^2 (L1/L4)*(L5/L4) -(L1/L4)^2 -(L1/L4)*(L5/L4);
                        (L1/L4)*(L5/L4) (L5/L4)^2 -(L1/L4)*(L5/L4) -(L5/L4)^2;
                        -(L1/L4)^2 -(L1/L4)*(L5/L4) (L1/L4)^2 (L1/L4)*(L5/L4);
                        -(L1/L4)*(L5/L4) -(L5/L4)^2 (L1/L4)*(L5/L4) (L5/L4)^2];
                        % [cos^2 cossin -cos^2 -cossin;cossin sin^2 -cossin -sin^2;
                        %  -cos^2 -cossin cos^2 cossin;-cossin -sin^2 cossin sin^2]
    esmloc5 = (E*A/L5)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1];
    esmloc6 = (E*A/L6)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4);
                        -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2;
                        -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4);
                         (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2];
                        % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2;
                        %  -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2]
    esmloc7 = (E*A/L7)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
    esmloc8 = (E*A/L8)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1];
    esmloc9 = (E*A/L9)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4);
                        -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2;
                        -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4);
                         (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2];
                        % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2;
                        %  -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2]
    % element nodal coordinates
    % ex1 and ey1 are the x and y coordinates to the nodes of element 1
    ex1 = [0 12];  ey1 = [0 0];
    ex2 = [12 24]; ey2 = [0 0];
    ex3 = [24 36]; ey3 = [0 0];
    ex4 = [0 12];  ey4 = [0 9];
    ex5 = [12 12]; ey5 = [9 0];
    ex6 = [12 24]; ey6 = [9 0];
    ex7 = [12 24]; ey7 = [9 9];
    ex8 = [24 24]; ey8 = [9 0];
    ex9 = [24 36]; ey9 = [9 0];
    K = zeros(12);
    F = zeros(12,1); 
    F(6)=-1200; F(11)=400;
    ep = [E A];
    % assembly of the global stiffness matrix
    K =  assem(Edof(1,:),K,esmloc1);
    K =  assem(Edof(2,:),K,esmloc2);
    K =  assem(Edof(3,:),K,esmloc3);
    K =  assem(Edof(4,:),K,esmloc4);
    K =  assem(Edof(5,:),K,esmloc5);
    K =  assem(Edof(6,:),K,esmloc6);
    K =  assem(Edof(7,:),K,esmloc7);
    K =  assem(Edof(8,:),K,esmloc8);
    K =  assem(Edof(9,:),K,esmloc9);
    % boundary conditions
    bc = [1 0;2 0;8 0];
    % reaction force vector
    Q = solveq(K,F,bc);
    % finding element displacements
    ed1 = extract(Edof(1,:),Q); ed2 = extract(Edof(2,:),Q);
    ed3 = extract(Edof(3,:),Q); ed4 = extract(Edof(4,:),Q);
    ed5 = extract(Edof(5,:),Q); ed6 = extract(Edof(6,:),Q);
    ed7 = extract(Edof(7,:),Q); ed8 = extract(Edof(8,:),Q);
    ed9 = extract(Edof(9,:),Q);
    disp = [ed1;ed2;ed3;ed4;ed5;ed6;ed7;ed8;ed9];
    % reaction forces at each element; member forces (stress)
    N1 = bar2s(ex1,ey1,ep,ed1); MF1 = N1/A;
    N2 = bar2s(ex2,ey2,ep,ed2); MF2 = N2/A;
    N3 = bar2s(ex3,ey3,ep,ed3); MF3 = N3/A;
    N4 = bar2s(ex4,ey4,ep,ed4); MF4 = N4/A;
    N5 = bar2s(ex5,ey5,ep,ed5); MF5 = N5/A;
    N6 = bar2s(ex6,ey6,ep,ed6); MF6 = N6/A;
    N7 = bar2s(ex7,ey7,ep,ed7); MF7 = N7/A;
    N8 = bar2s(ex8,ey8,ep,ed8); MF8 = N8/A;
    N9 = bar2s(ex9,ey9,ep,ed9); MF9 = N9/A;
    RF = [N1;N2;N3;N4;N5;N6;N7;N8;N9];
    MF = [MF1;MF2;MF3;MF4;MF5;MF6;MF7;MF8;MF9];

With this, the output is as follows:
    GlobalStiffnessMatrix = K 
    dispDofs = disp
    reactionForces = RF
    memberForces = MF

 GlobalStiffnessMatrix = 1.0e+003 *
    3.9584    1.0053   -2.6180         0         0         0         0         0   -1.3404   -1.0053         0         0
    1.0053    0.7540         0         0         0         0         0         0   -1.0053   -0.7540         0         0
   -2.6180         0    5.2360         0   -2.6180         0         0         0         0         0         0         0
         0         0         0    3.4907         0         0         0         0         0   -3.4907         0         0
         0         0   -2.6180         0    6.5764   -1.0053   -2.6180         0   -1.3404    1.0053         0         0
         0         0         0         0   -1.0053    4.2446         0         0    1.0053   -0.7540         0   -3.4907
         0         0         0         0   -2.6180         0    3.9584   -1.0053         0         0   -1.3404    1.0053
         0         0         0         0         0         0   -1.0053    0.7540         0         0    1.0053   -0.7540
   -1.3404   -1.0053         0         0   -1.3404    1.0053         0         0    5.2988         0   -2.6180         0
   -1.0053   -0.7540         0   -3.4907    1.0053   -0.7540         0         0         0    4.9986         0         0
         0         0         0         0         0         0   -1.3404    1.0053   -2.6180         0    3.9584   -1.0053
         0         0         0         0         0   -3.4907    1.0053   -0.7540         0         0   -1.0053    4.2446

 dispDofs =
         0         0    0.3056   -1.4992
    0.3056   -1.4992    0.6112   -2.1836
    0.6112   -2.1836    1.0695         0
         0         0    0.8260   -1.4992
    0.8260   -1.4992    0.3056   -1.4992
    0.8260   -1.4992    0.6112   -2.1836
    0.8260   -1.4992    0.5204   -1.9258
    0.5204   -1.9258    0.6112   -2.1836
    0.5204   -1.9258    1.0695         0

 reactionForces = 1.0e+003 *
    0.8000
    0.8000
    1.2000
   -0.5000
         0
    0.5000
   -0.8000
    0.9000
   -1.5000

 memberForces =
  254.6479
  254.6479
  381.9719
 -159.1549
         0
  159.1549
  254.6479
  286.4789
 -477.4648

It can be seen [http://en.wikiversity.org/wiki/Team4Report1#Problem_R1.6 here] that the problem done by hand corresponds to the force output of the MATLAB program.<br><br><br>

=====CALFEM MATLAB =====
    % CALFEM
    % global coordinate matrix
    Coord = [0 0;12 0;24 0;36 0;12 9;24 9];
    % element connectivity array
    conn = Edof;
    % global nodal dof matrix
    Dof = [1 2;3 4;5 6;7 8;9 10;11 12];
    % contructing element arrays
    [ex,ey]=coordxtr(conn,Coord,Dof,2);
    K=zeros(12);
    EA=[E A]; % Young's modulus, Area
    % compute element stiffness matrix Ke and global stiffness matrix
    for i=1:9
       Ke = bar2e(ex(i,:),ey(i,:),EA); % element stiffness matrix
       K = assem(conn(i,:),K,Ke);      % global stiffness matrix
    end
    % element displacement vector
    ed=extract(conn,Q);
    % reaction forces
    for i=1:9
       N(i)=bar2s(ex(i,:),ey(i,:),ep,ed(i,:));
    end
    % member forces
    MF1C = N(1)/A;  % stress in AC
    MF2C = N(2)/A;  % stress in CE
    MF3C = N(3)/A;  % stress in EF
    MF4C = N(4)/A;  % stress in AB
    MF5C = 0;  % stress in BC
    MF6C = N(6)/A;  % stress in BE
    MF7C = N(7)/A;  % stress in BD
    MF8C = N(8)/A;  % stress in DE
    MF9C = N(9)/A;  % stress in DF
    MFC = [MF1C;MF2C;MF3C;MF4C;MF5C;MF6C;MF7C;MF8C;MF9C];

The output to the CALFEM code is as follows:
    CGlobalStiffnessMatrix = K 
    CdispDofs = ed
    CreactionForces = N'
    CmemberForces = MFC

 CGlobalStiffnessMatrix = 1.0e+003 *
    3.9584    1.0053   -2.6180         0         0         0         0         0   -1.3404   -1.0053         0         0
    1.0053    0.7540         0         0         0         0         0         0   -1.0053   -0.7540         0         0
   -2.6180         0    5.2360         0   -2.6180         0         0         0         0         0         0         0
         0         0         0    3.4907         0         0         0         0         0   -3.4907         0         0
         0         0   -2.6180         0    6.5764   -1.0053   -2.6180         0   -1.3404    1.0053         0         0
         0         0         0         0   -1.0053    4.2446         0         0    1.0053   -0.7540         0   -3.4907
         0         0         0         0   -2.6180         0    3.9584   -1.0053         0         0   -1.3404    1.0053
         0         0         0         0         0         0   -1.0053    0.7540         0         0    1.0053   -0.7540
   -1.3404   -1.0053         0         0   -1.3404    1.0053         0         0    5.2988         0   -2.6180         0
   -1.0053   -0.7540         0   -3.4907    1.0053   -0.7540         0         0         0    4.9986         0         0
         0         0         0         0         0         0   -1.3404    1.0053   -2.6180         0    3.9584   -1.0053
         0         0         0         0         0   -3.4907    1.0053   -0.7540         0         0   -1.0053    4.2446

 CdispDofs =
         0         0    0.3056   -1.4992
    0.3056   -1.4992    0.6112   -2.1836
    0.6112   -2.1836    1.0695         0
         0         0    0.8260   -1.4992
    0.8260   -1.4992    0.3056   -1.4992
    0.8260   -1.4992    0.6112   -2.1836
    0.8260   -1.4992    0.5204   -1.9258
    0.5204   -1.9258    0.6112   -2.1836
    0.5204   -1.9258    1.0695         0

 CreactionForces = 1.0e+003 *
    0.8000
    0.8000
    1.2000
   -0.5000
         0
    0.5000
   -0.8000
    0.9000
   -1.5000

 CmemberForces =
  254.6479
  254.6479
  381.9719
 -159.1549
         0
  159.1549
  254.6479
  286.4789
 -477.4648

The results of all methods verifies their accuracy.

== Problem R2.4 ==

===Description:===
We are to find an expression for the excitation referenced in (1) p.53-6 and then find the actual solution and the amplification factor, A for the following conditions:

<math> f_{0}/m = r_{0} = 1 </math><br>
<math> \bar{w} = 0.9 w </math>

We are then to find the complete solution.


===Solution:===

<math> y'' + y' +(\frac{17}{4})y = \frac{f_{0}}{m} cos(\bar{w} t) </math>

Let us guess the solution is:

<math> y_{p}(t) = A sin(\bar{w} t) + B cos(\bar{w} t) </math>

then

<math> y_{p}(t) = A sin(\bar{w} t) + B cos(\bar{w} t) </math><br>
<math> y'_{p}(t) = A \bar{w} cos(\bar{w} t) - B \bar{w} sin(\bar{w} t) </math><br>
<math> y''_{p}(t) = - A \bar{w}^{2} cos(\bar{w} t) - B \bar{w}^{2} sin(\bar{w} t) </math>

We then write

<math> [ - A \bar{w}^{2} - B \bar{w} + (\frac{17}{4}) A ] sin(\bar{w} t) + [ - B \bar{w}^{2} + A \bar{w} + (\frac{17}{4}) B] cos(\bar{w} t) = \frac{f_{0}}{m} cos(\bar{w} t) </math>

It follows that

<math> [ - A \bar{w}^{2} - B \bar{w} + (\frac{17}{4}) A ] = 0 </math> and that <math> [ - B \bar{w}^{2} + A \bar{w} + (\frac{17}{4}) B] = \frac{f_{0}}{m} </math>

We solve for A and B to be

<math> A = \frac{f_{0}/m}{[\bar{w} + (\frac{17}{4})^{2}(\frac{1}{\bar{w}}) - \frac{34 \bar{w}}{4} + \bar{w}^{3} ]} </math><br>
<math> B = \frac{ \frac{f_{0}}{m} [ \frac{17}{4} - \bar{w}^{2}] }{[[\bar{w} + (\frac{17}{4})^{2}(\frac{1}{\bar{w}}) - \frac{34 \bar{w}}{4} + \bar{w}^{3} ] (\bar{w})]} </math>

Once given the conditions

<math> f_{0}/m = r_{0} = 1 </math><br>
<math> \bar{w} = 0.9 w </math>

we can find

<math> \bar{w} = 0.9 w  =  \bar{w} = 0.9 * 2.0615 = 1.855 </math> 

Plugging these values in yields

<math> A = 0.4528 </math><br>
<math> B = 0.1975 </math>

The particular solution is then

<math> y_{p}(t) = 0.4528 sin(1.855 t) + 0.1975 cos(1.855 t) </math>

Equations (3) p. 53-8 and (2) p. 53-8 give us

<math> \rho = \frac{\bar{w}}{w} = \frac{1.855}{2.0615} = 0.9</math><br>
<math> A = \frac{1}{[(1 - \rho^{2})^{2} + (2 \zeta \rho)^{2}]^{\frac{1}{2}}} =  \frac{1}{[(1 - 0.9^{2})^{2} + (2 *0.2425 *0.9)^{2}]^{\frac{1}{2}}} = 2.628 </math>

The following is the graph of the homogeneous solution:

[[File:Pr 5 third graph.jpg]]

The following is the graph of the particular solution:

[[File:Particular sol.jpg]]

The following is the graph of the total solution, <math> y(t) = y_{h}(t) + y_{p}(t) </math>:

[[File:Total sol.jpg]]

A graph of all three plots is produced below:

[[File:All graphs.jpg]]

==Problem R2.5 CALFEM 3.4 Manual exs1 on a Spring System and Verification==
===Given a Linearly Elastic Spring System===
Consider a system of 3 linearly elastic springs 
[[image:R.2.5 Fig 1.JPG|frame|center|Figure 1: Spring System from CALFEM Manual 3.4 p9.3-2 ([http://en.wikiversity.org/wiki/File:R.2.5_Fig_1.JPG])]]

Use the following diagram for the node locations
[[image:R.2.5 fig 2.JPG|frame|center|Figure 1: Spring System node numbering from CALFEM Manual 3.4 p9.3-2 ([http://en.wikiversity.org/wiki/File:R.2.5_fig_2.JPG])]]

===Find: ===
====1.Find the results when running exs1 from CALFEM====
====2.Verify the results of CALFEM to manual calculations====
===Solution: Results of CALFEM===
====Construct the topology matrix====
This matrix contains the node numbers and degrees of freedom.

<math>
\begin{bmatrix}
1 & 1 & 2 \\
2 & 2 & 3 \\
3 & 2 & 3
\end{bmatrix}
</math>

====Construct the global stiffness matrix and load vector f====
<math>
\mathbf{K}=
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
</math>

The load vector is in position 2 with F=100

<math>
\mathbf{f}=
\begin{bmatrix}
0 \\
100 \\
0 
\end{bmatrix}

</math>

====Generating the element stiffness matrices====
The element stiffness matrices are generated using the CALFEM
function spring1e.  The matrices are generated using k and 2k where k=1500

<math>
\mathbf{Ke1}=
\begin{bmatrix}
1500 & -1500 \\
-1500 & 1500
\end{bmatrix}
</math>

<math>
\mathbf{Ke2}=
\begin{bmatrix}
3000 & -3000 \\
-3000& 3000
\end{bmatrix}
</math>


====Assembling the element stiffness matrices into the global stiffness matrix====
The assembly is done using the assem fucntion from CALFEM.

Assembling the second element at the first row:

<math>
\mathbf{K}=
\begin{bmatrix}
3000 & -3000 & 0 \\
-3000 & 3000 & 0 \\
0 & 0 & 0
\end{bmatrix}
</math>

Assembling the first element at the second row:

<math>
\mathbf{K}=
\begin{bmatrix}
3000 & -3000 & 0 \\
-3000 & 4500 & -1500 \\
0 & -1500 & 1500
\end{bmatrix}
</math>

Assembling the second element at the third row:

<math>
\mathbf{K}=
\begin{bmatrix}
3000 & -3000 & 0 \\
-3000 & 7500 & -4500 \\
0 & -4500 & 4500
\end{bmatrix}
</math>



====Solve the system of equations using boundary conditions====
The CALFEM fucntion solveq is used to solve the system of equations subject to certain boundary conditions.

<math>
\mathbf{bc}=
\begin{bmatrix}
1 & 0 \\
3 & 0
\end{bmatrix}
</math>

<math>
\mathbf{a}=
\begin{bmatrix}
0 \\
0.0133 \\
0
\end{bmatrix}
</math>

<math>
\mathbf{a}=
\begin{bmatrix}
-40 \\
0 \\
-60
\end{bmatrix}
</math>


====Evaluate the element forces from the element displacements====

<math>
\mathbf{ed1}=
\begin{bmatrix}
0 & 0.0133
\end{bmatrix}
</math>

<math>
\mathbf{ed2}=
\begin{bmatrix}
0.0133 & 0
\end{bmatrix}
</math>

<math>
\mathbf{ed3}=
\begin{bmatrix}
0.0133 & 0
\end{bmatrix}
</math>

====Evaluate the spring forces from function spring1s====

<math>
\mathbf{es1}=
\begin{bmatrix}
40
\end{bmatrix}
</math>

<math>
\mathbf{es2}=
\begin{bmatrix}
-20
\end{bmatrix}
</math>

<math>
\mathbf{es3}=
\begin{bmatrix}
-40
\end{bmatrix}
</math>

===Solution: Results of verification===
This section verifies the CALFEM code by manual calculations

====Nodal equilibrium equations====

<math>f_1^{(1)}=F_1 = R_1</math>

<math>f_2^{(1)}+ f_2^{(2)}+ f_2^{(3)}=F = 100 </math> 

<math>f_3^{(2)}+f_3^{(3)}=F_2=R_2</math>

====Setting up the equation to solve====

<math>
\begin{bmatrix}K_s\end{bmatrix}\begin{Bmatrix}Q_s\end{Bmatrix}=\begin{Bmatrix}F_s \end{Bmatrix}
</math>


<math> 
1000\begin{bmatrix}3 &-3  &-2  &0 \\  -3&7.5  &-4.5 \\  0&-4.5  &4.5 \end{bmatrix}\begin{Bmatrix}u_1\\ u_2\\ u_3\end{Bmatrix}=\begin{Bmatrix}R_1\\ 100\\ R_2\end{Bmatrix}</math>

====Implementing the boundary conditions====

<math>
1000\begin{bmatrix}7.5 \end{bmatrix}\begin{Bmatrix} u_2 \end{Bmatrix}=\begin{Bmatrix}100 \end{Bmatrix}
</math>

====Inverting the global stiffness matrix to find displacements====

<math> u_2=0.01333 </math>

Displacements 1 and 3 are zero due to boundary conditions.

====Calculate forces in springs from equation====

The forces in the springs are found by the following equations. Negative forces display compressions and positive tensions.
 
<math>P=k(u_j-u_i)</math>

<math> P^{(1)}=3000(0.01333-0)=40</math>

<math>P^{(2)}=1500(0-0.01333)=-20</math>

<math>P^{(3)}=3000(0-0.01333)=-40</math>

====Discuss if the solution is verified====

The solution by CALFEM is verified since both the displacements and the forces are the same by both methods.

==Problem R2.6==
===Problem Statement===
Find the displacements of the nodes as well as the forces in the springs(tension/compression). Also determine the reaction forces at the walls.
[[File:R.2.6 fig 1.png]]
 

===Nodal Equilibrium Equations===
  <math>f_1^{(1)}+ f_1^{(4)}=F_1 = R_1</math>
  <math>f_2^{(1)}+ f_2^{(2)}+ f_2^{(3)}=F_2 = 0 </math> 
  <math>f_3^{(3)}+ f_3^{(4)}+ f_3^{(5)}=F_3 = 1000</math>
  <math>f_4^{(2)}+ f_4^{(5)}+ f_4^{(6)}=F_4 = 0</math>
  <math>f_5^{(6)}=F_5=R_5</math>

 <math>\begin{bmatrix}K_s\end{bmatrix}\begin{Bmatrix}Q_s\end{Bmatrix}=\begin{Bmatrix}F_s \end{Bmatrix}</math>
 

 <math> 100\begin{bmatrix}7 &-5  &-2  &0  &0 \\  -5&15  &-6  &-4  &0 \\  -2&-6  &12  &-4  &0 \\  0&-4&-4  &11  &-3 \\ 0&0  &0  &-3  &3 \end{bmatrix}\begin{Bmatrix}u_1\\ u_2\\ u_3\\ u_4\\ u_5\end{Bmatrix}=\begin{Bmatrix}R_1\\ 0\\ 1000\\ 0\\ R_5\end{Bmatrix}</math>
 <math>100\begin{bmatrix}15  &-6  &-4\\ -6  &12  &-4\\-4&-4  &11\\ \end{bmatrix}\begin{Bmatrix} u_2\\ u_3\\ u_4\end{Bmatrix}=\begin{Bmatrix}0\\ 1000\\ 0\end{Bmatrix}</math>


'''Using Matlab'''<br>
'''Solutions'''<br>
'''Displacements'''
[[File:R.2.6 fig 2.png]]
 
These are the displacements at <math>u_2</math>,   <math>u_3</math>, and    <math>u_4</math> respectively. The negative sign shows that the nodes are moving in the direction of the applied force which is against the convention of positive displacement left to right. 

Now it is possible to find the forces in the spring. Negative forces denote compression and positive forces denote tension.<br>
===Forces in Springs===
 <math>P=k(u_j-u_i)</math>
 <math> P^{(1)}=500(-0.8542-0)=-427.1N</math>
 <math>P^{(2)}=400(-0.8750+0.8542)=-8.32N</math>
 <math>P^{(3)}=600(-1.5521+0.8542)=-419N</math>
 <math>P^{(4)}=200(-0.8750-0)=-310N</math>
 <math>P^{(5)}=400(-0.8750+1.5521)=270.8N</math>
 <math>P^{(6)}=300(0+0.8750)=262.5N</math>
Now,
===Reaction Forces===
 <math> f_1^{(1)}+ f_1^{(4)}=F_1 = R_1</math>

 <math>R_1=-427N-310N=-737N</math>

 <math> f_5^{(6)}=F_5=R_5</math>

 <math>R_5=262.5N</math>

===Verification with Calfem===

'''a''' values are the nodal displacements 1 through 5.<br>
'''r''' values are the reaction forces at the walls.<br>
'''es1-es6''' are the forces in the springs. Negative values denote compression and positive values denote tension.
These answers match those previously calculated.


<syntaxhighlight lang="matlab">
===Code===
Edof=[1 1 2;2 2 4;3 2 3;4 1 3;5 3 4;6 4 5];
ep1=500;
ep2=400;
ep3=600;
ep4=200;
ep5=400;
ep6=300;
Ke1=spring1e(ep1);
Ke2=spring1e(ep2);
Ke3=spring1e(ep3);
Ke4=spring1e(ep4);
Ke5=spring1e(ep5);
Ke6=spring1e(ep6);
K=zeros(5,5);
K=assem(Edof(1,:),K,Ke1);
K=assem(Edof(2,:),K,Ke2);
K=assem(Edof(3,:),K,Ke3);
K=assem(Edof(4,:),K,Ke4);
K=assem(Edof(5,:),K,Ke5);
K=assem(Edof(6,:),K,Ke6);
bc= [1 0; 5 0];
f=zeros(5,1);
f(3)=-1000;
[a,r]=solveq(K,f,bc)
ed1=extract(Edof(1,:),a);
ed2=extract(Edof(2,:),a);
ed3=extract(Edof(3,:),a);
ed4=extract(Edof(4,:),a);
ed5=extract(Edof(5,:),a);
ed6=extract(Edof(6,:),a);
es1=spring1s(ep1,ed1)
es2=spring1s(ep2,ed2)
es3=spring1s(ep3,ed3)
es4=spring1s(ep4,ed4)
es5=spring1s(ep5,ed5)
es6=spring1s(ep6,ed6)

'''Solution'''
a =

         0
   -0.8542
   -1.5521
   -0.8750
         0


r =

  737.5000
   -0.0000
    0.0000
    0.0000
  262.5000


es1 =

 -427.0833


es2 =

   -8.3333


es3 =

 -418.7500


es4 =

 -310.4167


es5 =

  270.8333


es6 =

  262.5000
</syntaxhighlight>






--[[User:EML4507.s13.team4ever.Wulterkens|EML4507.s13.team4ever.Wulterkens]] ([[User talk:EML4507.s13.team4ever.Wulterkens|discuss]] • [[Special:Contributions/EML4507.s13.team4ever.Wulterkens|contribs]]) 07:21, 6 February 2013 (UTC)

{{CourseCat}}