User:Eml4507.s13.team4.vcb

Table of Assignments R7

Problem Assignments R5
Problem #	Solved&Typed by	Reviewed by
1a	Vernon Babich, Chad Colocar	All
1b	Vernon Babich, Tyler Wulterkens	All
2	David Bonner, Chad Colocar	All

Problem R7.1a: Verify the dim of the ${\underline {k}}_{6x6}^{(e)}{\underline {d}}_{6x1}^{(e)}$ matrices (fead.f08.mtgs.[37-41] pg. 2)

On our honor, we did this assignment on our own.

Given: The desired dimensions for k and d with index of 6

${\underline {k}}_{6x6}^{(e)}{\underline {d}}_{6x1}^{(e)}$

Find:

1.Use an index of 6 to verify the dimensions of k and d

Obtain:

${\underline {k}}_{6x6}^{(e)}{\underline {d}}_{6x1}^{(e)}={\underline {f}}_{6x1}^{(e)}$

With supporting: ${\underline {k}}_{6x6}^{(e)}={\underline {\tilde {T}}}_{6x6}^{(e)T}{\underline {\tilde {k}}}_{6x6}^{(e)}{\underline {\tilde {T}}}_{6x6}^{(e)}$

From:

${\underline {\tilde {k}}}_{6x6}^{(e)}{\underline {\tilde {d}}}_{6x1}^{(e)}={\underline {\tilde {f}}}_{6x1}^{(e)}$

Solution: Verify the dimensions of k and d

The constructed 6x6 matrix shown below represents ${\underline {\tilde {T}}}$

${\begin{Bmatrix}{\tilde {d}}_{1}\\{\tilde {d}}_{2}\\{\tilde {d}}_{3}\\{\tilde {d}}_{4}\\{\tilde {d}}_{5}\\{\tilde {d}}_{6}\end{Bmatrix}}={\begin{Bmatrix}R_{11}&R_{12}&0&0&0&0\\R_{21}&R_{22}&0&0&0&0\\0&0&1&0&0&0\\0&0&0&R_{11}&R_{12}&0\\0&0&0&R_{21}&R_{22}&0\\0&0&0&0&0&1\end{Bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\end{Bmatrix}}$

The following matrix represents the element stiffness matrix in local coordinates.

${\underline {\tilde {k}}}_{6x6}^{(e)}={\begin{bmatrix}{\frac {EA}{L}}&0&0&{\frac {-EA}{L}}&0&0\\0&{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}&0&{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}&0&{\frac {-6EI}{L^{2}}}&{\frac {2EI}{L}}\\{\frac {-EA}{L}}&0&0&{\frac {EA}{L}}&0&0\\0&{\frac {-12EI}{L^{3}}}&{\frac {-6EI}{L^{2}}}&0&{\frac {12EI}{L^{3}}}&{\frac {-6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}&0&{\frac {-6EI}{L^{2}}}&{\frac {4EI}{L}}\end{bmatrix}}$

Using the above two matricies and the transpose of T, you can construct the following equation:

${\underline {k}}_{6x6}^{(e)}={\underline {\tilde {T}}}_{6x6}^{(e)T}{\underline {\tilde {k}}}_{6x6}^{(e)}{\underline {\tilde {T}}}_{6x6}^{(e)}$

Multiplying ${\underline {\tilde {k}}}_{6x6}^{(e)}$ by ${\underline {\tilde {d}}}_{6x1}^{(e)}$ will yield the required ${\underline {\tilde {f}}}_{6x1}^{(e)}$

Problem R7.1b: Solve 2 element frame system (fead.f08.mtgs.[37-41] pg. 3)

On our honor, we did this assignment on our own.

Given: Information on the two element truss system

Assume a square cross section. Also use same data from fea.f08.mtgs.p5-4.

Figure 1: Two member system (fig.JPG)

Element length: (1) =4

Element length: (2) =2

Young's modulus: (1) =3

Young's modulus: (2) =5

Cross section area: (1) =1

Cross section area: (2) =2

Inclination angle: (1) = 30 deg

Inclination angle: (2) = -45 deg

Solution:

1. Plot Undeformed Shape

Undeformed Shape.

Figure 1: Two member system (2bnodef.JPG)

View of element with single force member view.

Figure 1: Two member system (2b.JPG)

2. Plot Deformed Shape 2-bar Truss

Plotting the deformed shape 2-bar truss for the given problem from p5-2 of the fead08 lecture notes. The deformation can be seen the truss using the forces given from p5-2.

Figure 1: Two member system (2bdef2.JPG)

3. Plot Deformed Shape 2-bar Frame

Plotting the deformed shape 2-bar frame for the given problem from p5-2 of the fead08 lecture notes. The constraints are preserved during the analysis in this mode to keep the frame.

Figure 1: Two member system (2bdef.JPG)

ref

for i=1:2:5
    x_1((i+1)/2)=P(1,((i+1)/2))+V(i,5);
end

for i=2:2:6
    y_1((i/2))=P(2,(i/2))+V(i,1);
end

%Required deformation plotting
hold on
for i=1:2
    l1=N(1,i);
    l2=N(2,i);
    x1=[x_1(l1),x_1(l2)];
    y1=[y_1(l1),y_1(l2)];
    axis([-1 5 -1 5])
    plot(x1,y1,'b')
    title(['Constrained Shape'])
    hold on
end

%Required deformation plotting
for i=1:2:5
    x_o((i+1)/2)=P(1,((i+1)/2));
end

for i=2:2:6
    y_o((i/2))=P(2,(i/2));
end
 for i=1:2
    l1=N(1,i);
    l2=N(2,i);
    x1=[x_o(l1),x_o(l2)];
    y1=[y_o(l1),y_o(l2)];
    axis([-1 5 -1 5])
    plot(x1,y1,'-.or')
    hold on
 end
end

Problem R7.2

On my honor, I have neither given nor recieved unauthorized aid in doing this assignment.

Description

We are to resolve problem 5.7. We are to solve for motion of the truss using modal superposition using the three lowest eigenvalues.

Solution

The following code was used for problem 5.7 to obtain the K and M matrices and to obtain the eigenvalues and the eigenvectors:

function R7P2a

p = 2;
E = 5;
A = 0.5;
L = 1;
Ld = 1*sqrt(2);

%degree of freedom for each element
dof = zeros(10,5);
dof(1,:) = [1 1 2 3 4];
dof(2,:) = [2 1 2 5 6];
dof(3,:) = [3 3 4 5 6];
dof(4,:) = [4 5 6 9 10];
dof(5,:) = [5 5 6 7 8];
dof(6,:) = [6 3 4 9 10];
dof(7,:) = [7 3 4 7 8];
dof(8,:) = [8 7 8 9 10];
dof(9,:) = [9 9 10 11 12];
dof(10,:) = [10 7 8 11 12];

%position of each node
pos = zeros(2,6);
pos(:,1) = [0;0];
pos(:,2) = [L;0];
pos(:,3) = [L;L];
pos(:,4) = [2*L;0];
pos(:,5) = [2*L;L];
pos(:,6) = [3*L;0];

%Connections
conn = zeros(2,10);
conn(:,1) = [1;2];
conn(:,2) = [1;3];
conn(:,3) = [2;3];
conn(:,4) = [3;5];
conn(:,5) = [3;4];
conn(:,6) = [2;5];
conn(:,7) = [2;4];
conn(:,8) = [4;5];
conn(:,9) = [5;6];
conn(:,10) = [4;6];

%seperates into x and y coordinates
x = zeros(10,2);
y = zeros(10,2);
for i = 1:10
    x(i,:) = [pos(1,conn(1,i)),pos(1,conn(2,i))];
    y(i,:) = [pos(2,conn(1,i)),pos(2,conn(2,i))];
end

%Set up K and M matrix
K = zeros(12);
M = zeros(12);
for i = 1:10
    xm = x(i,2)-x(i,1);
    ym = y(i,2)-y(i,1);
    L = sqrt(xm^2+ym^2);
    l = xm/L;
    m = ym/L;
    k = 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*p;
    m = [m/2 0 0 0;0 m/2 0 0;0 0 m/2 0;0 0 0 m/2];
    Dof = dof(i,2:5);
    K(Dof,Dof) = K(Dof,Dof)+k;
    M(Dof,Dof) = M(Dof,Dof)+m;
end
bc = [1;2;12];
K
M
F = [0;0;0;0;0;0;0;-5;0;0;0;0];

[L,X] = eigen(K,M,bc)

From this, we obtain M and K matrices which are:

K =

 Columns 1 through 9

   3.3839    0.8839   -2.5000         0   -0.8839   -0.8839         0         0         0
   0.8839    0.8839         0         0   -0.8839   -0.8839         0         0         0
  -2.5000         0    5.8839    0.8839         0         0   -2.5000         0   -0.8839
        0         0    0.8839    3.3839         0   -2.5000         0         0   -0.8839
  -0.8839   -0.8839         0         0    4.2678         0   -0.8839    0.8839   -2.5000
  -0.8839   -0.8839         0   -2.5000         0    4.2678    0.8839   -0.8839         0
        0         0   -2.5000         0   -0.8839    0.8839    5.8839   -0.8839         0
        0         0         0         0    0.8839   -0.8839   -0.8839    3.3839         0
        0         0   -0.8839   -0.8839   -2.5000         0         0         0    4.2678
        0         0   -0.8839   -0.8839         0         0         0   -2.5000         0
        0         0         0         0         0         0   -2.5000         0   -0.8839
        0         0         0         0         0         0         0         0    0.8839

 Columns 10 through 12

        0         0         0
        0         0         0
  -0.8839         0         0
  -0.8839         0         0
        0         0         0
        0         0         0
        0   -2.5000         0
  -2.5000         0         0
        0   -0.8839    0.8839
   4.2678    0.8839   -0.8839
   0.8839    3.3839   -0.8839
  -0.8839   -0.8839    0.8839

M =

 Columns 1 through 9

   1.2071         0         0         0         0         0         0         0         0
        0    1.2071         0         0         0         0         0         0         0
        0         0    2.2071         0         0         0         0         0         0
        0         0         0    2.2071         0         0         0         0         0
        0         0         0         0    2.4142         0         0         0         0
        0         0         0         0         0    2.4142         0         0         0
        0         0         0         0         0         0    2.2071         0         0
        0         0         0         0         0         0         0    2.2071         0
        0         0         0         0         0         0         0         0    2.4142
        0         0         0         0         0         0         0         0         0
        0         0         0         0         0         0         0         0         0
        0         0         0         0         0         0         0         0         0

 Columns 10 through 12

        0         0         0
        0         0         0
        0         0         0
        0         0         0
        0         0         0
        0         0         0
        0         0         0
        0         0         0
        0         0         0
   2.4142         0         0
        0    1.2071         0
        0         0    1.2071

This will also output the lowest eigenpairs:

$\gamma _{1},\phi _{1}$

$\gamma _{2},\phi _{2}$

$\gamma _{3},\phi _{3}$

Then, with these, we will find what the modal equations are according to the equation:

$z^{''}+\gamma _{i}z=\phi _{i}^{T}F(t)$

This will yield 3 unique differential equations. We solve for the complete solution to these differential equation using the boundary conditions from:

$z_{i}(0)={\bar {\phi }}_{i}^{T}Md(0)$ and

$z_{i}^{'}(0)={\bar {\phi }}_{i}^{T}Md^{'}(0)$

With these solutions for z, we can then find and plot the actual displacements from modal superposition using:

$d(t)=\Sigma _{j=1}^{3}z_{j}\phi _{j}$

Problem 6.4 Deformation in 3D Space

On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Problem Statement

There is a 3D truss system with applied forces. We need to optimize each beam to meet a factor of safety of 1.5.

Matlab Solution

Set up degrees of freedom and coordinates

%Constants
E = 30000000;   %Youngs Modulus
L1 = 3;         %Short lengths
L2 = 8;         %Long lengths
v = 0.3;        %poissons ratio
YS = 37000;     %Yield Strength
p = 0.00073;    %density
FS = 1.5;        %Factor of Safety
%coordinates
pos = zeros(3,10);
pos(:,1) = [-L1/2;0;8];
pos(:,2) = [L1/2;0;8];
pos(:,3) = [-L1/2;L1/2;4];
pos(:,4) = [L1/2;L1/2;4];
pos(:,5) = [L1/2;-L1/2;4];
pos(:,6) = [-L1/2;-L1/2;4];
pos(:,7) = [-L2/2;L2/2;0];
pos(:,8) = [L2/2;L2/2;0];
pos(:,9) = [L2/2;-L2/2;0];
pos(:,10) = [-L2/2;-L2/2;0];
%Connections
conn = zeros(2,25);
conn(:,1) = [1;2];
conn(:,2) = [1;4];
conn(:,3) = [2;3];
conn(:,4) = [1;5];
conn(:,5) = [2;6];
conn(:,6) = [2;4];
conn(:,7) = [2;5];
conn(:,8) = [1;3];
conn(:,9) = [1;6];
conn(:,10) = [3;6];
conn(:,11) = [4;5];
conn(:,12) = [3;4];
conn(:,13) = [5;6];
conn(:,14) = [3;10];
conn(:,15) = [6;7];
conn(:,16) = [4;9];
conn(:,17) = [5;8];
conn(:,18) = [4;7];
conn(:,19) = [3;8];
conn(:,20) = [5;10];
conn(:,21) = [6;9];
conn(:,22) = [6;10];
conn(:,23) = [3;7];
conn(:,24) = [4;8];
conn(:,25) = [5;9];
%Degree of Freedom
dof = zeros(25,7);
D = zeros(10,3);
    %assign nodes degree of freedoms
for i = 1:10
    j = (i-1)*2+i;
    D(i,:) = [j,j+1,j+2];
end
    %Assign members their DOF's
for i = 1:25
    dof(i,:) = [i,D(conn(1,i)',:),D(conn(2,i)',:)];
end
%seperates into x and y coordinates
x = zeros(25,2);
y = zeros(25,2);
z = zeros(25,2);
for i = 1:25
    x(i,:) = [pos(1,conn(1,i)),pos(1,conn(2,i))];
    y(i,:) = [pos(2,conn(1,i)),pos(2,conn(2,i))];
    z(i,:) = [pos(3,conn(1,i)),pos(3,conn(2,i))];
end

Set up K matrix to be used for calculations

%Set up K and M matrix
K = zeros(30);
M = zeros(30);

A = zeros(25,1);
for i = 1:25
    A(i) = pi()*(2/2)^2;   %Area of all members
end
for i = 1:25
    xm = x(i,2)-x(i,1);
    ym = y(i,2)-y(i,1);
    zm = z(i,2)-z(i,1);
    L = sqrt(xm^2+ym^2+zm^2);
    l = xm/L;
    m = ym/L;
    n = zm/L;
    k = E*A(i)/L*[l^2 l*m l*n -l^2 -l*m -l*n;l*m m^2 m*n -l*m -m^2 -m*n;l*n m*n n^2 -l*n -m*n -n^2;-l^2 -l*m -l*n l^2 l*m l*n;-l*m -m^2 -m*n l*m m^2 m*n;-l*n -m*n -n^2 l*n m*n n^2];
    m = L*A(i)*p;
    m = [m/2 0 0 0 0 0;0 m/2 0 0 0 0;0 0 m/2 0 0 0;0 0 0 m/2 0 0;0 0 0 0 m/2 0;0 0 0 0 0 m/2];
    Dof = dof(i,2:7);
    K(Dof,Dof) = K(Dof,Dof)+k;
    M(Dof,Dof) = M(Dof,Dof)+m;
end

Reduce matrices and make initial conditions

%Makes an alternate K and M matrix to calncel rows and columns for boundary
%conditions.
Kr = K;
Mr = M;
Kr((19:30),:) = [];
Kr(:,(19:30)) = [];
Mr((19:30),:) = [];
Mr(:,(19:30)) = [];

%Applied Forces
F = zeros(18,1);
F(2) = 60000;
F(5) = 60000;

Determine the deformation of each DOF

d = Kr\F;
d = [d;zeros(12,1)];

Set up new coordinates and determine stress and strain

posn = zeros(3,10);
%new position vectors
for i = 1:10
    j = i*3;
    posn(1,i) = pos(1,i) + d(j-2);
    posn(2,i) = pos(2,i) + d(j-1);
    posn(3,i) = pos(3,i) + d(j);
end

xn = zeros(25,2);
yn = zeros(25,2);
zn = zeros(25,2);
%Put in new coordinates
for i = 1:25
    xn(i,:) = [posn(1,conn(1,i)),posn(1,conn(2,i))];
    yn(i,:) = [posn(2,conn(1,i)),posn(2,conn(2,i))];
    zn(i,:) = [posn(3,conn(1,i)),posn(3,conn(2,i))];
end

%Finding stress and FOS
e = zeros(25,1);
o = zeros(25,1);
for i = 1:25
    xm = x(i,2)-x(i,1);
    ym = y(i,2)-y(i,1);
    zm = z(i,2)-z(i,1);
    L = sqrt(xm^2+ym^2+zm^2);
    xl = xn(i,2)-xn(i,1);
    yl = yn(i,2)-yn(i,1);
    zl = zn(i,2)-zn(i,1);
    Ll = sqrt(xl^2+yl^2+zl^2);
    e(i) = (Ll-L)/L;                %Strain
    o(i) = e(i)*E;                  %Stress
end

Determine factor of safety and reduce area

FOS = zeros(25,1);
on = 1.5*YS;
for i = 1:25
    A(i) = A(i)*abs(o(i))/on;
    FOS(i) = YS/abs(o(i));
    if A(i) <= pi()*(0.1/2)^2
       A(i) = pi()*(0.1/2)^2;
    elseif A(i) >= pi()*(2.5/2)^2
        A(i) = pi()*(2.5/2)^2;
    end

end

Conclusion

The final area of each member was able to be reduced to its minimum of 0.1 inches and stay above a factor of safety of 1.5. Final stresses are in psi:

1. 13.94
2. -8487.15
3. -8487.15
4. 8611.63
5. 8611.63
6. -12871.54
7. 13045.99
8. -12871.54
9. 13045.99
10. 63.42
11. 63.42
12. 788.95
13. -778.48
14. -5171.02
15. 5174.49
16. -5171.02
17. 5174.49
18. -7787.74
19. -7787.74
20. 7792.39
21. 7792.39
22. 14680.46
23. -14679.76
24. -14679.76
25. 14680.46

CALFEM Verification

The Matlab calculations were confirmed using the CALFEM toolbox in Matlab.

The following presents the constants used and options for different designs

E = 30000000;   %modulus
p = 0.00073;    %density
L1 = 3;         %Short lengths
L2 = 8;         %Long lengths
v = 0.3;        %poissons ratio
YS = 37000;     %Yield Strength
FS = 1.5;        %Factor of Safety
A = pi()*(2/2)^2;
A = 0.1 %used if reducing area to 0.1 in
ep=[E A p*A];

This Edof matrix uses column 1 as the as the element number. The next three columns are the DOFs for the first node and the next three columns are for the second node since this is a three dimensional problem.

  Edof =[1     1     2     3     4     5     6
     2     1     2     3    10    11    12
     3     4     5     6     7     8     9
     4     1     2     3    13    14    15
     5     4     5     6    16    17    18
     6     4     5     6    10    11    12
     7     4     5     6    13    14    15
     8     1     2     3     7     8     9
     9     1     2     3    16    17    18
    10     7     8     9    16    17    18
    11    10    11    12    13    14    15
    12     7     8     9    10    11    12
    13    13    14    15    16    17    18
    14     7     8     9    28    29    30
    15    16    17    18    19    20    21
    16    10    11    12    25    26    27
    17    13    14    15    22    23    24
    18    10    11    12    19    20    21
    19     7     8     9    22    23    24
    20    13    14    15    28    29    30
    21    16    17    18    25    26    27
    22    16    17    18    28    29    30
    23     7     8     9    19    20    21
    24    10    11    12    22    23    24
    25    13    14    15    25    26    27];

Coord contains the x, y, and z coordinates for each node.

Coord =[-1.5000         0    8.0000
    1.5000         0    8.0000
   -1.5000    1.5000    4.0000
    1.5000    1.5000    4.0000
    1.5000   -1.5000    4.0000
   -1.5000   -1.5000    4.0000
   -4.0000    4.0000         0
    4.0000    4.0000         0
    4.0000   -4.0000         0
   -4.0000   -4.0000         0];

Dof =[1     2     3
     4     5     6
     7     8     9
    10    11    12
    13    14    15
    16    17    18
    19    20    21
    22    23    24
    25    26    27
    28    29    30];

The bar3e is used to make the stiffness matrix for a 3D system. The bar3s is used to find the member forces in a 3D system. One can see that the 3 denotes the dimension of the system being observed. The stresses can be found by dividing the forces by the areas.

Ex, Ey, Ez]=coordxtr(Edof, Coord, Dof, 2);

K=zeros(30);F=zeros(30,1);F(1)=60,000; F(2)=60,000;

for i=1:25
    Ke=bar3e(Ex(i,:),Ey(i,:),Ez(i,:),ep);
    K=assem(Edof(i,:),K,Ke);
end

bc=[7 0;8 0;9 0;10 0;];

Q=solveq(K,F,bc);

Ed=extract(Edof,Q);

for i=1:25
    N(i)=bar3s(Ex(i,:),Ey(i,:),Ez(i,:),ep,Ed(i,:));
end

S=N/A

Problem R5.7: Consider the free vibration of truss system p.53-22b (fead.f13.sec53b.[21])

On our honor, we did this assignment on our own.

Given: Truss system under free vibration with applied force and zero initial velocity

Consider the truss system. The initial force is applied at node four and the force equals five.

Use standard node and element naming conventions.

Figure 1: Truss system (fead.f13.sec53b p.53-22b). ([1])

Solution: Solve for the truss motion

K =

 Columns 1 through 5

   3.3839    0.8839   -2.5000         0   -0.8839
   0.8839    0.8839         0         0   -0.8839
  -2.5000         0    5.8839    0.8839         0
        0         0    0.8839    3.3839         0
  -0.8839   -0.8839         0         0    4.2678
  -0.8839   -0.8839         0   -2.5000         0
        0         0   -2.5000         0   -0.8839
        0         0         0         0    0.8839
        0         0   -0.8839   -0.8839   -2.5000
        0         0   -0.8839   -0.8839         0
        0         0         0         0         0
        0         0         0         0         0

 Columns 6 through 10

  -0.8839         0         0         0         0
  -0.8839         0         0         0         0
        0   -2.5000         0   -0.8839   -0.8839
  -2.5000         0         0   -0.8839   -0.8839
        0   -0.8839    0.8839   -2.5000         0
   4.2678    0.8839   -0.8839         0         0
   0.8839    5.8839   -0.8839         0         0
  -0.8839   -0.8839    3.3839         0   -2.5000
        0         0         0    4.2678         0
        0         0   -2.5000         0    4.2678
        0   -2.5000         0   -0.8839    0.8839
        0         0         0    0.8839   -0.8839

 Columns 11 through 12

        0         0
        0         0
        0         0
        0         0
        0         0
        0         0
  -2.5000         0
        0         0
  -0.8839    0.8839
   0.8839   -0.8839
   3.3839   -0.8839
  -0.8839    0.8839

p = 2;
E = 5;
A = 0.5;
L = 1;
Ld = 1*sqrt(2);

%degree of freedom for each element
dof = zeros(10,5);
dof(1,:) = [1 1 2 3 4];
dof(2,:) = [2 1 2 5 6];
dof(3,:) = [3 3 4 5 6];
dof(4,:) = [4 5 6 9 10];
dof(5,:) = [5 5 6 7 8];
dof(6,:) = [6 3 4 9 10];
dof(7,:) = [7 3 4 7 8];
dof(8,:) = [8 7 8 9 10];
dof(9,:) = [9 9 10 11 12];
dof(10,:) = [10 7 8 11 12];

%position of each node
pos = zeros(2,6);
pos(:,1) = [0;0];
pos(:,2) = [L;0];
pos(:,3) = [L;L];
pos(:,4) = [2*L;0];
pos(:,5) = [2*L;L];
pos(:,6) = [3*L;0];

%Connections
conn = zeros(2,10);
conn(:,1) = [1;2];
conn(:,2) = [1;3];
conn(:,3) = [2;3];
conn(:,4) = [3;5];
conn(:,5) = [3;4];
conn(:,6) = [2;5];
conn(:,7) = [2;4];
conn(:,8) = [4;5];
conn(:,9) = [5;6];
conn(:,10) = [4;6];

%seperates into x and y coordinates
x = zeros(10,2);
y = zeros(10,2);
for i = 1:10
    x(i,:) = [pos(1,conn(1,i)),pos(1,conn(2,i))];
    y(i,:) = [pos(2,conn(1,i)),pos(2,conn(2,i))];
end

%Set up K and M matrix
K = zeros(12);
M = zeros(12);
for i = 1:10
    xm = x(i,2)-x(i,1);
    ym = y(i,2)-y(i,1);
    L = sqrt(xm^2+ym^2);
    l = xm/L;
    m = ym/L;
    k = 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*p;
    m = [m/2 0 0 0;0 m/2 0 0;0 0 m/2 0;0 0 0 m/2];
    Dof = dof(i,2:5);
    K(Dof,Dof) = K(Dof,Dof)+k;
    M(Dof,Dof) = M(Dof,Dof)+m;
end

Table of Assignments R5

Problem Assignments R5
Problem #	Solved&Typed by	Reviewed by
1	David Bonner, Vernon Babich	All
2	Chad Colocar, David Bonner	All
3	Bryan Tobin, Tyler Wulterkens	All
4	Tyler Wulterkens, Vernon Babich	All
5	Vernon Babich, Tyler Wulterkens	All
6	David Bonner, Vernon Babich	All
7	Tyler Wulterkens, Chad Colocar, Vernon Babich	All

Problem R5.5a: Eigen vector plot for zero evals of the 2 bar truss system p.21-2 (fead.f08.mtgs.[21])

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

Given: Two member zero evals

Consider a 2 member system. Plot the eigen vectors corresponding to zero evals and interpret the results.

Use standard node and element naming conventions.

Figure 1: Two member system (fig.JPG)

Solution: Plot the eigen values and interpret

This is the general stiffness matrix for a two bar system as given.

${\vec {K}}={\begin{bmatrix}&1&2&3&4&5&6\\1&(l^{(1)})^{2}*K^{(1)}&l^{(1)}m^{(1)}*K^{(1)}&-(l^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&0&0\\2&l^{(1)}m^{(1)}*K^{(1)}&(m^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&-(m^{(1)})^{2}*K^{(1)}&0&0\\3&-(l^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&(l^{(1)})^{2}*K^{(1)}+(l^{(2)})^{2}*K^{(2)}&l^{(1)}m^{(1)}*K^{(1)}+l^{(2)}m^{(2)}*K^{(2)}&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}\\4&-l^{(1)}m^{(1)}*K^{(1)}&-(m^{(1)})^{2}*K^{(1)}&l^{(1)}m^{(1)}*K^{(1)}+l^{(2)}m^{(2)}*K^{(2)}&(m^{(1)})^{2}*K^{(1)}+(m^{(2)})^{2}*K^{(2)}&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}\\5&0&0&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}&(l^{(2)})^{2}*K^{(2)}&l^{(2)}m^{(2)}*K^{(2)}\\6&0&0&-l^{(2)}m^{(2)}*K^{(2)}&-(m^{(2)})^{2}*K^{(2)}&l^{(2)}m^{(2)}*K^{(2)}&(m^{(2)})^{2}*K^{(2)}\\\end{bmatrix}}$

The global stiffness matrix in numerical form match to this specific problem is as follows.

${\vec {K}}={\begin{bmatrix}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.0625&-2.1752&-2.5000&2.5000\\-0.3248&-0.1875&-2.1752&2.6875&2.5000&-2.5000\\0&0&-2.5000&2.5000&2.5000&-2.5000\\0&0&2.5000&-2.5000&-2.5000&2.5000\\\end{bmatrix}}$

The eigen vector matrix is as follows.

${\vec {V}}={\begin{bmatrix}-0.1118&0.5043&-0.0000&-0.5931&0.6174&-0.0139\\-0.0814&-0.8634&0.0000&-0.3476&0.3565&-0.0080\\-0.4628&0.0089&0.0000&-0.4803&-0.5409&0.5123\\0.5266&-0.0053&-0.0000&-0.5429&-0.4330&-0.4904\\-0.4947&0.0071&-0.7071&0.0313&-0.0765&-0.4984\\0.4947&-0.0071&-0.7071&-0.0313&0.0765&0.4984\\\end{bmatrix}}$

${\vec {D}}={\begin{bmatrix}0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&0&0&1.4705&0\\0&0&0&0&0&10.0295\\\end{bmatrix}}$

The eigen vectors correspond to the columns of the eigen vector matrix shown above. The eigen value mode shapes are a linear combination of the pure mode shapes (pure rigid body and pure mechanism).

Figure 1: The Vectors plotted from the eigen vector matrix ([2])

function R5_5e2

K= [9/16 (3*sqrt(3))/16 -9/16 -(3*sqrt(3))/16 0 0;
(3*sqrt(3))/16 3/16 -(3*sqrt(3))/16 -3/16 0 0;
-9/16 -(3*sqrt(3))/16 49/16 ((3*sqrt(3)-40)/16) -5/2 5/2;
-(3*sqrt(3))/16 -3/16 ((3*sqrt(3)-40)/16) 43/16 5/2 -5/2;
0 0 -5/2 5/2 5/2 -5/2;
0 0 5/2 -5/2 -5/2 5/2]

Eval=eig(K)

[V,D]=eig(K)

V
D=abs(D)

%Each column of V is a mode

 for i=1:4
   figure
   plot(V(:,i));
   title(['Mode ',num2str(i)])
 end

end

As a thought experiment, we plotted the position plus the deformation derived from the eigen vector matrix. This will supply 4 shapes but they are not constrained to the boundary conditions. Applying the boundary conditions to the K matrix reduces the number of non zero eigen vector outputs to two. The results are shown below. The figure plots came from the eigen vectors in the fifth and sixth columns of the constrained eigen vector matrix. Using the other columns produced a figure that was not constrained to the fixed supports.

Figure 1: Mode shape constrained to the boundary conditions. The dotted line represents the original shape. ([3])

Figure 1: Mode shape constrained to the boundary conditions. The dotted line represents the original shape. ([4])

function R5_5e

K= [9/16 (3*sqrt(3))/16 -9/16 -(3*sqrt(3))/16 0 0;
(3*sqrt(3))/16 3/16 -(3*sqrt(3))/16 -3/16 0 0;
-9/16 -(3*sqrt(3))/16 49/16 ((3*sqrt(3)-40)/16) -5/2 5/2;
-(3*sqrt(3))/16 -3/16 ((3*sqrt(3)-40)/16) 43/16 5/2 -5/2;
0 0 -5/2 5/2 5/2 -5/2;
0 0 5/2 -5/2 -5/2 5/2]

%K_red=K([3 4],[3 4]);

K1=[0 0 0 0 0 0;
    0 0 0 0 0 0;
    0 0 K(3,3) K(3,4) 0 0;
    0 0 K(4,3) K(4,4) 0 0;
    0 0 0 0 0 0;
    0 0 0 0 0 0]

[V D] = eig(K1)

N=[1 2;
    2 3];

P=[0 cos(pi/6)*4 cos(pi/6)*4+cos(pi/4)*2;
   0 sin(pi/6)*4 sin(pi/6)*4-sin(pi/4)*2];

for i=1:2:5
    x_1((i+1)/2)=P(1,((i+1)/2))+V(i,1);
end

for i=2:2:6
    y_1((i/2))=P(2,(i/2))+V(i,1);
end

%Plot deformed system with constraints
hold on
for i=1:2
    l1=N(1,i);
    l2=N(2,i);
    x1=[x_1(l1),x_1(l2)];
    y1=[y_1(l1),y_1(l2)];
    axis([-1 5 -1 5])
    plot(x1,y1,'b')
    title(['Constrained Shape'])
    hold on
end

%Plot undeformed system
for i=1:2:5
    x_o((i+1)/2)=P(1,((i+1)/2));
end

for i=2:2:6
    y_o((i/2))=P(2,(i/2));
end
 for i=1:2
    l1=N(1,i);
    l2=N(2,i);
    x1=[x_o(l1),x_o(l2)];
    y1=[y_o(l1),y_o(l2)];
    axis([-1 5 -1 5])
    plot(x1,y1,'-.or')
    hold on
 end
end

Problem R5.5b: Eigen vector plot for zero evals of the square 3 bar truss system p.21-3 (figure a) (fead.f08.mtgs.[21])

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

Given: Square 3 bar truss system

Consider a 3 member system. Plot the eigen vectors corresponding to zero evals and interpret the results.

a=b=1 E=2 A=3

Figure 1: Two member system (b.JPG)

Find:

1.Plot the eigen values

Use standard node and element naming conventions.

Solution: Plot the eigen values

${\vec {V}}={\begin{bmatrix}-0.7071&0&0&-0.7071\\0&1.0000&0&0\\-0.7071&0&0&0.7071\\0&0&1.0000&0\\\end{bmatrix}}$

${\vec {D}}={\begin{bmatrix}0&0&0&0\\0&6&0&0\\0&0&6&0\\0&0&0&12\\\end{bmatrix}}$

Figure 1: The only shape corresponding to zero evals. ([5])

function R5_5ec
E=2;
A=3;
L=1;
k=(E*A)/L;
k1=k*[0 0 0 0;
    0 1 0 -1;
    0 0 0 0;
    0 -1 0 1]
k3=k*[0 0 0 0;
    0 1 0 -1;
    0 0 0 0;
    0 -1 0 1]
k2=k*[1 0 -1 0;
    0 0 0 0;
    -1 0 1 0;
    0 0 0 0]

K=zeros(8,8)

K(1:4,1:4)=k1
K(5:8,5:8)=k3

Kt=zeros(8,8)
Kt(3:6,3:6)=k2

Kg=K+Kt

%boundary conditions applied
K1=[0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0;
    0 0 Kg(3,3) Kg(3,4) Kg(3,5) Kg(3,6) 0 0;
    0 0 Kg(4,3) Kg(4,4) Kg(4,5) Kg(4,6) 0 0;
    0 0 Kg(5,3) Kg(5,4) Kg(5,5) Kg(5,6) 0 0;
    0 0 Kg(6,3) Kg(6,4) Kg(6,5) Kg(6,6) 0 0;
    0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0]

Kred=[Kg(3,3) Kg(3,4) Kg(3,5) Kg(3,6);
     Kg(4,3) Kg(4,4) Kg(4,5) Kg(4,6);
     Kg(5,3) Kg(5,4) Kg(5,5) Kg(5,6);
     Kg(6,3) Kg(6,4) Kg(6,5) Kg(6,6)]

[V D] = eig(Kred)

 for i=1:6
   figure
   plot(V(:,i));
   title(['Mode ',num2str(i)])
 end

end

Table of Assignments

Problem Assignments
Problem #	Solved&Typed by	Reviewed by
1	David Bonner, Bryan Tobin	All
2	Bryan Tobin, David Bonner	All
3	Vernon Babich, Tyler Wulterkens	All
4	Tyler Wulterkens, Chad Colocar, Vernon Babich, David Bonner	All

Problem R4.3: Post processing of member forces and reactions p.11-3 (fead.f08.mtgs.[10-18])

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

Given: Two member system with initial force P

Consider a 2 member system with an applied force. The displacements are assumed to be known. This allows for the post processing for the reactions both globally and for the members.

Figure 1: Two member system (fig.JPG)

Solution: Post process for member forces and reactions

Both methods first require the displacements which are assumed as given since we are post processing for the member forces.

This is the general stiffness matrix for a two bar system as given.

${\vec {K}}={\begin{bmatrix}&1&2&3&4&5&6\\1&(l^{(1)})^{2}*K^{(1)}&l^{(1)}m^{(1)}*K^{(1)}&-(l^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&0&0\\2&l^{(1)}m^{(1)}*K^{(1)}&(m^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&-(m^{(1)})^{2}*K^{(1)}&0&0\\3&-(l^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}&(l^{(1)})^{2}*K^{(1)}+(l^{(2)})^{2}*K^{(2)}&l^{(1)}m^{(1)}*K^{(1)}+l^{(2)}m^{(2)}*K^{(2)}&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}\\4&-l^{(1)}m^{(1)}*K^{(1)}&-(m^{(1)})^{2}*K^{(1)}&l^{(1)}m^{(1)}*K^{(1)}+l^{(2)}m^{(2)}*K^{(2)}&(m^{(1)})^{2}*K^{(1)}+(m^{(2)})^{2}*K^{(2)}&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}\\5&0&0&-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}&(l^{(2)})^{2}*K^{(2)}&l^{(2)}m^{(2)}*K^{(2)}\\6&0&0&-l^{(2)}m^{(2)}*K^{(2)}&-(m^{(2)})^{2}*K^{(2)}&l^{(2)}m^{(2)}*K^{(2)}&(m^{(2)})^{2}*K^{(2)}\\\end{bmatrix}}$

The first two columns and the last two columns are struck out due to fixed positions of nodes 1 and 3. Since we are post processing for the reactions forces and member forces, we know the displacements at node 2. These displacements are used along with the global stiffness matrix to solve for the reaction forces. The reaction forces will be parallel to the members since the truss elements are two force members.

${\begin{Bmatrix}F_{1x}\\F_{1y}\\F_{3x}\\F_{3y}\end{Bmatrix}}={\begin{Bmatrix}-(l^{(1)})^{2}*K^{(1)}&-l^{(1)}m^{(1)}*K^{(1)}\\-(l^{(1)})m^{(1)}*K^{(1)}&-(m^{(1)})^{2}*K^{(1)}\\-(l^{(2)})^{2}*K^{(2)}&-l^{(2)}m^{(2)}*K^{(2)}\\-(l^{(2)})m^{(2)}*K^{(2)}&-(m^{(2)})^{2}*K^{(2)}\end{Bmatrix}}{\begin{Bmatrix}u2\\v2\end{Bmatrix}}$

Method 1 (square root of sum of squares)

Once the x and y components of the global forces on node 1 have been found, the square root of the sum of squares of the component can be used to find the member force magnitude. This is the member force since we are observing a two force member.

$\mathbf {Element1force} ={\sqrt {F_{1x}^{2}+F_{1y}^{2}}}$

$\mathbf {Element2force} ={\sqrt {F_{3x}^{2}+F_{3y}^{2}}}$

Method 2 (transformation matrix)

The local items are denoted by the overhead bar for each term. The global items do not have this accent.

This first step describes how you obtain the local displacements from the global system.

${\begin{bmatrix}\mathbf {T} \end{bmatrix}}{\begin{Bmatrix}\mathbf {q} \end{Bmatrix}}={\begin{Bmatrix}\mathbf {\overline {q}} \end{Bmatrix}}$

${\begin{bmatrix}{\overline {u_{1}}}\\{\overline {v_{1}}}\\{\overline {u_{2}}}\\{\overline {v_{2}}}\end{bmatrix}}={\begin{bmatrix}l&m&0&0\\-l&m&0&0\\0&0&l&m\\0&0&-l&m\end{bmatrix}}{\begin{bmatrix}u_{1}\\v_{1}\\u_{2}\\v_{2}\end{bmatrix}}$

${\begin{bmatrix}\mathbf {T} \end{bmatrix}}{\begin{Bmatrix}\mathbf {f} \end{Bmatrix}}={\begin{Bmatrix}\mathbf {\overline {f}} \end{Bmatrix}}$

${\begin{bmatrix}f_{1{\overline {x}}}\\f_{1{\overline {y}}}\\f_{2{\overline {x}}}\\f_{2{\overline {y}}}\end{bmatrix}}={\frac {EA}{L}}{\begin{bmatrix}1&0&-1&0\\0&0&0&0\\-1&0&1&0\\0&0&0&0\end{bmatrix}}{\begin{bmatrix}u_{1}\\v_{1}\\u_{2}\\v_{2}\end{bmatrix}}$

Discussion of Methods

As you could see the method 1 is most simple for easy systems. However, getting into larger systems, the second method will be more useful for automating with an application such as Matlab. The systematic approach of method two allows for large sets of data to be systematically crunched based on the global system matrix. It can be seen that the computational efficiency could vary based on the application. I would choose method two for most approaches since it is more generally applied to large systems.

Problem Pb-53.5: Draw the FBDs and derive the dif. eq. of motion from p.53-13 (2 DOFs spring-mass-damper system)

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

Given a Linearly Elastic Spring-Mass_Damper System With 2 DOFs

Consider a system of three linearly elastic springs and dampers with two masses.

Figure 1: Spring system obtained from EML4507 sec.53.b, Pb-53.5 (Dr. Vu-Quoc) ([6])

Find:

2.Derive the differential equation of motion and the coefficient matrices

$\mathbf {M} {\ddot {\mathbf {d} }}+\mathbf {C} {\dot {\mathbf {d} }}+\mathbf {K} {\mathbf {d} }=\mathbf {F} (t)$

$\mathbf {d} ={\begin{Bmatrix}d_{1}\\d_{2}\end{Bmatrix}}_{2\times 1}$

$\mathbf {F} ={\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}_{2\times 1}$

$\mathbf {M} ={\begin{bmatrix}m_{1}&0\\0&m_{2}\end{bmatrix}}_{2\times 2}$

$\mathbf {C} ={\begin{bmatrix}(c_{1}+c_{2})&-c_{2}\\-c_{2}&(c_{2}+c_{3})\end{bmatrix}}_{2\times 2}$

$\mathbf {K} ={\begin{bmatrix}(k_{1}+k_{2})&-k_{2}\\-k_{2}&(k_{2}+k_{3})\end{bmatrix}}_{2\times 2}$

Solution: FBDs of all components

Construct the first FBD (mass 1)

This FBD is of the left mass (mass 1).

Figure 1: Original Image: FBD of mass 1. ([7])

Construct the second FBD (mass 2)

This FBD is of the right mass (mass 2).

Figure 1: Original Image: FBD of mass 2. ([8])

Construct the third FBD (left wall support)

This FBD is of the left wall support (item 3).

Figure 1: Original Image: FBD of left wall support. ([9])

Construct the fourth FBD (right wall support)

This FBD is of the right wall support (item 4).

Figure 1: Original Image: FBD of right wall support. ([10])

Solution: Derive the differential equation of motion and the coefficient matrices

Approach

The position of the masses in the system are described by d1 and d2. The walls are fixed. This means that the values for d3 and d4, for the left and right walls respectively, are zero.

Force Derivations

The springs are assumed to be at their unstretched length (equilibrium) in the figure. A positive d1 causes tension in spring 1 and compression in spring 2. This observation is made from holding m2 stationary. Viewing the last spring from equilibrium, if the mass 2 moves to the right and d2 becomes positive, the third spring will compress.

The differential equations of motion will show that the system is coupled.

$F_{s1}=k_{1}d_{1}$

$F_{d1}=c_{1}{\dot {d}}_{1}$

$F_{s2}=k_{2}(d_{2}-d_{1})$

$F_{d2}=c_{2}({\dot {d}}_{2}-{\dot {d}}_{1})$

$F_{s3}=k_{3}d_{2}$

$F_{d3}=c_{3}{\dot {d}}_{2}$

Summing the Forces

$m_{1}{\ddot {d}}_{1}+(c_{1}+c_{2}){\dot {d}}_{1}-c_{2}{\dot {d}}_{2}+(k_{1}+k_{2})d_{1}-k_{2}d_{2}=F_{1}$

$m_{2}{\ddot {d}}_{2}-c_{2}{\dot {d}}_{1}+(c_{2}+c_{3}){\dot {d}}_{2}-k_{2}d_{1}+(k_{2}+k_{3})d_{2}=F_{2}$

Writing the Two Differential Equations in Matrix Form

Grouping the dependent equations and writing the coefficients in matrix form yields the results shown below.

$\mathbf {M} {\ddot {\mathbf {d} }}+\mathbf {C} {\dot {\mathbf {d} }}+\mathbf {K} {\mathbf {d} }=\mathbf {F} (t)$

$\mathbf {d} ={\begin{Bmatrix}d_{1}\\d_{2}\end{Bmatrix}}_{2\times 1}$

$\mathbf {F} ={\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}_{2\times 1}$

$\mathbf {M} ={\begin{bmatrix}m_{1}&0\\0&m_{2}\end{bmatrix}}_{2\times 2}$

$\mathbf {C} ={\begin{bmatrix}(c_{1}+c_{2})&-c_{2}\\-c_{2}&(c_{2}+c_{3})\end{bmatrix}}_{2\times 2}$

$\mathbf {K} ={\begin{bmatrix}(k_{1}+k_{2})&-k_{2}\\-k_{2}&(k_{2}+k_{3})\end{bmatrix}}_{2\times 2}$

Contribution

Problem Assignments
Problem #	Solved&Typed by	Reviewed by
1	Chad Colocar, Bryan Tobin	All
2	Vernon Babich, Tyler Wulterkens	All
3	Bryan Tobin, Chad Colocar	All
4	David Bonner, Vernon Babich	All
5	Was Dean Pickett, David Bonner	All
6	Was Dean Pickett, Tyler Wulterkens	All
7	Tyler Wulterkens, Chad Colocar	All
8	David Bonner, Vernon Babich	All

Problem 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

Figure 1: Spring System from CALFEM Manual 3.4 p9.3-2 ([11])

Use the following diagram for the node locations

Figure 1: Spring System node numbering from CALFEM Manual 3.4 p9.3-2 ([12])

Solution: Results of CALFEM

Construct the topology matrix

This matrix contains the node numbers and degrees of freedom.

${\begin{bmatrix}1&1&2\\2&2&3\\3&2&3\end{bmatrix}}$

Construct the global stiffness matrix and load vector f

$\mathbf {K} ={\begin{bmatrix}0&0&0\\0&0&0\\0&0&0\end{bmatrix}}$

The load vector is in position 2 with F=100

$\mathbf {f} ={\begin{bmatrix}0\\100\\0\end{bmatrix}}$

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

$\mathbf {Ke1} ={\begin{bmatrix}1500&-1500\\-1500&1500\end{bmatrix}}$

$\mathbf {Ke2} ={\begin{bmatrix}3000&-3000\\-3000&3000\end{bmatrix}}$

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:

$\mathbf {K} ={\begin{bmatrix}3000&-3000&0\\-3000&3000&0\\0&0&0\end{bmatrix}}$

Assembling the first element at the second row:

$\mathbf {K} ={\begin{bmatrix}3000&-3000&0\\-3000&4500&-1500\\0&-1500&1500\end{bmatrix}}$

Assembling the second element at the third row:

$\mathbf {K} ={\begin{bmatrix}3000&-3000&0\\-3000&7500&-4500\\0&-4500&4500\end{bmatrix}}$

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.

$\mathbf {bc} ={\begin{bmatrix}1&0\\3&0\end{bmatrix}}$

$\mathbf {a} ={\begin{bmatrix}0\\0.0133\\0\end{bmatrix}}$

$\mathbf {a} ={\begin{bmatrix}-40\\0\\-60\end{bmatrix}}$

Evaluate the element forces from the element displacements

$\mathbf {ed1} ={\begin{bmatrix}0&0.0133\end{bmatrix}}$

$\mathbf {ed2} ={\begin{bmatrix}0.0133&0\end{bmatrix}}$

$\mathbf {ed3} ={\begin{bmatrix}0.0133&0\end{bmatrix}}$

Evaluate the spring forces from function spring1s

$\mathbf {es1} ={\begin{bmatrix}40\end{bmatrix}}$

$\mathbf {es2} ={\begin{bmatrix}-20\end{bmatrix}}$

$\mathbf {es3} ={\begin{bmatrix}-40\end{bmatrix}}$

Solution: Results of verification

This section verifies the CALFEM code by manual calculations

Nodal equilibrium equations

$f_{1}^{(1)}=F_{1}=R_{1}$

$f_{2}^{(1)}+f_{2}^{(2)}+f_{2}^{(3)}=F=100$

$f_{3}^{(2)}+f_{3}^{(3)}=F_{2}=R_{2}$

Setting up the equation to solve

${\begin{bmatrix}K_{s}\end{bmatrix}}{\begin{Bmatrix}Q_{s}\end{Bmatrix}}={\begin{Bmatrix}F_{s}\end{Bmatrix}}$

$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}}$

Implementing the boundary conditions

$1000{\begin{bmatrix}7.5\end{bmatrix}}{\begin{Bmatrix}u_{2}\end{Bmatrix}}={\begin{Bmatrix}100\end{Bmatrix}}$

Inverting the global stiffness matrix to find displacements

$u_{2}=0.01333$

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.

$P=k(u_{j}-u_{i})$

$P^{(1)}=3000(0.01333-0)=40$

$P^{(2)}=1500(0-0.01333)=-20$

$P^{(3)}=3000(0-0.01333)=-40$

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.

R.1.2

R.1.2a

Problem diagram for the mass spring damper system. Obtained from FEAD:sec.53a, p.53-2b, Pb-53.1.

Free body diagram observing rectilinear motion for the mass spring damper system. Original image.

${\begin{array}{lcr}{\mbox{Inertial Frame fixed at t=0 shown on diagram}}\\{\underline {E}}_{y}:to\ the\ right\\{\underline {E}}_{x}:upward\\{\underline {E}}_{z}:{\underline {E}}_{y}\times {\underline {E}}_{x}\\\\{\mbox{Kinematics}}\\{\mbox{Observing rectilinear motion}}\\_{}^{F}{\textrm {\underline {r}}}=y{\underline {E}}_{y}\\_{}^{F}{\textrm {\underline {v}}}={\dot {y}}{\underline {E}}_{y}={\mbox{v}}{\underline {E}}_{y}\\\\{\mbox{Kinetics}}\\{\underline {F}}_{s}=-k(l-l_{o}){\underline {u}}_{s}\\l=l_{o}+y\\{\underline {F}}_{s}=-k(y){\underline {E}}_{y}\\{\underline {F}}_{d}=-c{\frac {d{\underline {r}}}{d{\underline {t}}}}=-c{\mbox{v}}{\underline {E}}_{y}\\{\underline {f}}(t)=f(t){\underline {E}}_{y}\\{\underline {F}}=-k(y){\underline {E}}_{y}-c{\mbox{v}}{\underline {E}}_{y}+f(t){\underline {E}}_{y}=m{\ddot {y}}{\underline {E}}_{y}\\{\mbox{The components in the y-direction yield a differential equation}}\\m{\ddot {y}}+c{\dot {y}}+ky={\underline {f}}(t)\end{array}}$\\

R.1.2b

Problem diagram for the mass spring damper system. Obtained from IEA S12:p.1-4

Free body diagram observing rectilinear motion for the mass spring damper system. Original image.

${\begin{array}{lcr}{\mbox{Kinematics}}\\y=y_{k}+y_{c}(1)\\\\{\mbox{Kinetics}}\\m{y}''+f_{1}=f(t)(2)\\f_{1}=f_{k}=f_{c}(3)\\\\{\mbox{Constitutive Equations}}\\f_{k}=ky_{k}(4)\\f_{c}=c{y_{c}}'(5)\\(1):y''={y_{k}}''+{y_{c}}''(6)\\(3):f_{k}=f_{c}\Rightarrow k{y_{k}}=c{y_{c}}'(7)\\\Rightarrow {y_{c}}'={\frac {k}{c}}{y_{k}}(8)\\{\mbox{Using}}(1)and(8)\\{y}''={y_{k}}''+{({y_{c}}')}'={y_{k}}''+{({\frac {k}{c}}{y_{k}})}'={y_{k}}''+{\frac {k}{c}}{{y_{k}}'}(9)\\{\mbox{Using}}(9)and(2)-(3)\\m({y_{k}}''+{\frac {k}{c}}{{y_{k}}'})+k{y_{k}}=f(t)(10)\end{array}}$\\

R.1.5

Solution plotted for roots equal to r = -0.5 and r = -1.5

Solution plotted for both roots equal to r = -0.5

Solution plotted for roots equal to r = -0.5 + 2i and r = -0.5 - 2i

Plot of all solutions