University of Florida/Eml4507/s13 Team4Report1

Main Code
Below is the main page of MATLAB for reproducing all examples.

EDU>> A = [1 2 3;4 5 6;7 8 9]

A =

1    2     3     4     5     6     7     8     9

EDU>> A = [ 1 2 3 4 5 6 7 8 9 ]

A =

1    2     3     4     5     6     7     8     9

EDU>> A = [1 2;3 4] + i*[5 6;7 8]

A =

Column 1

1.0000 + 5.0000i 3.0000 + 7.0000i

Column 2

2.0000 + 6.0000i 4.0000 + 8.0000i

EDU>> A = [1+5i 2+6i;3+7i 4+8i]

A =

Column 1

1.0000 + 5.0000i 3.0000 + 7.0000i

Column 2

2.0000 + 6.0000i 4.0000 + 8.0000i

EDU>> A = [1 2 3;4 5 6;7 8 9]

A =

1    2     3     4     5     6     7     8     9

EDU>> B = [A, zeros(3,2); zeros(2,3), eye(2)]

B =

1    2     3     0     0     4     5     6     0     0     7     8     9     0     0     0     0     0     1     0     0     0     0     0     1

EDU>> x = [];

EDU>> n = 3;

EDU>> for i = 1:n

x = [x,i^2] end

x =

1

x =

1    4

x =

1    4     9

EDU>> x = [];

EDU>> for i = n:-1:1

x = [x,i^2] end

x =

9

x =

9    4

x =

9    4     1

EDU>> m = 4

m =

4

EDU>> for i = 1:m

for j = 1:n

H(i, j) = 1/(i+j-1);

end

end

EDU>> H

H =

1.0000   0.5000    0.3333    0.5000    0.3333    0.2500    0.3333    0.2500    0.2000    0.2500    0.2000    0.1667

EDU>> s = 0;

EDU>> for c = A

s = s + sum(c);

end

EDU>> s

s =

45

EDU>> n = 0;

EDU>> a = 12;

EDU>> while 2^n < a

n = n + 1;

end

EDU>> n

n =

4

EDU>> if n < 0

parity = 0;

elseif rem(n,2) == 0

parity = 2;

else parity = 1;

end

EDU>> parity

parity =

2

EDU>> A

A =

1    2     3     4     5     6     7     8     9

EDU>> y = eig(A)

y =

16.1168  -1.1168   -0.0000

EDU>> [U,D] = eig(A)

U =

-0.2320  -0.7858    0.4082   -0.5253   -0.0868   -0.8165   -0.8187    0.6123    0.4082

D =

16.1168        0         0         0   -1.1168         0         0         0   -0.0000

EDU>> m = 2; n = 3; x = 0:.01:2*pi; y = sin(m*x); z = cos(n*x); plot(x,y,x,z)



EDU>> C = 0.2:0.2:1.2

C =

Columns 1 through 3

0.2000   0.4000    0.6000

Columns 4 through 6

0.8000   1.0000    1.2000

EDU>> C = 5:-1:1

C =

5    4     3     2     1

EDU>> x = [0.0:0.1:2.0]';

EDU>> y = sin(x);

EDU>> [x y]

ans =

0        0    0.1000    0.0998    0.2000    0.1987    0.3000    0.2955    0.4000    0.3894    0.5000    0.4794    0.6000    0.5646    0.7000    0.6442    0.8000    0.7174    0.9000    0.7833    1.0000    0.8415    1.1000    0.8912    1.2000    0.9320    1.3000    0.9636    1.4000    0.9854    1.5000    0.9975    1.6000    0.9996    1.7000    0.9917    1.8000    0.9738    1.9000    0.9463    2.0000    0.9093

EDU>> A(1:3,2)

ans =

2    5     8

EDU>> A(:,3)

ans =

3    6     9

EDU>> A(1:3,:)

ans =

1    2     3     4     5     6     7     8     9

EDU>> A(:,[2 3])

ans =

2    3     5     6     8     9

EDU>> z = randint(4,5)

z =

1    9     3     7     0     2     9     9     3     1     4     4     3     2     9     0     4     1     4     9

EDU>> xm = stat(z)

xm =

Columns 1 through 3

1.7500   6.5000    4.0000

Columns 4 through 5

4.0000   4.7500

EDU>> M = magic(6);

EDU>> E = zeros(6,50);

EDU>> for j = 1:50

E(:,j) = eig(M^i);

end

EDU>> E

E =

1.0e+08 *

Columns 1 through 3

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 4 through 6

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 7 through 9

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 10 through 12

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 13 through 15

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 16 through 18

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 19 through 21

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 22 through 24

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 25 through 27

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 28 through 30

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 31 through 33

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 34 through 36

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 37 through 39

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 40 through 42

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 43 through 45

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 46 through 48

1.5181   1.5181    1.5181    0.0000    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001    0.0001    0.0001

Columns 49 through 50

1.5181   1.5181    0.0000    0.0000    0.0053    0.0053    0.0053    0.0053    0.0001    0.0001    0.0001    0.0001

EDU>> s = 'This is a test'

s =

This is a test

EDU>> disp('this message is hereby displayed')

this message is hereby displayed

EDU>> error('Sorry, the matrix must be symmetric')

??? Sorry, the matrix must be

symmetric EDU>> inter = input('Enter the number of iterations:  ')

Enter the number of iterations:  32

inter =

[]

ans =

32

EDU>> x = -4:0.01:4; y = sin(x); plot(x,y)



EDU>> x = -1.5:0.01:1.5; y = exp(-x.^2); plot(x,y)



EDU>> t = 0:0.001:2*pi; x=cos(3*t); y=sin(2*t); plot(x,y)

EDU>> title ('Best Least Squares Fit')

EDU>> axis([-1,1,-1,1])

EDU>> x=0:0.01:2*pi;y1=sin(x) ;y2=sin(2*x);y3=sin(4*x);plot(x,y1,x,y2,x,y3)



EDU>> x=0:0.01:2*pi; Y=[sin(x)', sin(2*x)', sin(4*x)']; plot(x,Y)

EDU>> x=0:0.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x);

EDU>> plot(x,y1,'--',x,y2,':',x,y3,'+')



EDU>> t=0.01:0.01:20*pi; x=cos(t); y=sin(t); z=t.^3; plot3(x,y,z)



EDU>> mesh(eye(10))



EDU>> surf(eye(10))



EDU>> xx = -2:0.2:2;

EDU>> yy = xx;

EDU>> [x,y] = meshgrid(xx,yy);

EDU>> z = exp(-2.^2-y.^2);

EDU>> mesh(z)



EDU>> F = floor(10*rand(6)); F = triu(tril(F,1),-1);

EDU>> S = sparse(F)

S =

(1,1)       5   (2,2)        1   (3,2)        6   (2,3)        2   (3,3)        7   (4,3)        1   (3,4)        7   (5,4)        9   (4,5)        3   (5,5)        5   (6,5)        5   (5,6)        8   (6,6)        5

EDU>> F = full(S)

F =

5    0     0     0     0     0     0     1     2     0     0     0     0     6     7     7     0     0     0     0     1     0     3     0     0     0     0     9     5     8     0     0     0     0     5     5

EDU>> m = 6; n = 6; e = ones(n,1); d = -2*e;

EDU>> T = spdiags([e,d,e],[-1,0,1],m,n)

T =

(1,1)      -2   (2,1)        1   (1,2)        1   (2,2)       -2   (3,2)        1   (2,3)        1   (3,3)       -2   (4,3)        1   (3,4)        1   (4,4)       -2   (5,4)        1   (4,5)        1   (5,5)       -2   (6,5)        1   (5,6)        1   (6,6)       -2

EDU>> i = [1 2 3 4 4 4]; j = [1 2 3 1 2 3]; s = [5 6 7 8 9 10];

EDU>> S = sparse(i,j,s,4,3),full(S)

S =

(1,1)       5   (4,1)        8   (2,2)        6   (4,2)        9   (3,3)        7   (4,3)       10

ans =

5    0     0     0     6     0     0     0     7     8     9    10

EDU>> sparse(i,j,s,m,n)

ans =

(1,1)       5   (4,1)        8   (2,2)        6   (4,2)        9   (3,3)        7   (4,3)       10

EDU>> n = 6; e = floor(10*rand(n-1,1)); E = sparse(2:n,1:n-1,e,n,n)

E =

(2,1)       3   (3,2)        9   (4,3)        8   (5,4)        5   (6,5)        6

EDU>> n=20;e=ones(n,1);d=-2*e; T=spdiags([e,d,e],[-1,0,1],n,n); A=full(T);

EDU>> b=ones(n,1);s=sparse(b);tic,T\s;sparsetime=toc, tic,A\b;fulltime=toc

sparsetime =

0.0011

fulltime =

0.0394

Bisect.m
Bisect returns a zero of the function. fun is a string containing the name of a real-valued MATLAB function of a single real variable. x is a starting guess. Optional 3rd input is tol sets the accuracy of the result.

function [b, steps] = bisect(fun, x, tol) if nargin < 3, tol = eps; end trace = (nargout == 2); if x ~= 0, dx = x/20; else, dx = 1/20; end a = x - dx; fa = feval(fun,a); b = x + dx; fb = feval(fun,b); while (fa > 0) == (fb > 0) dx = 2.0*dx; a = x - dx; fa = feval(fun,a); if (fa > 0) ~= (fb > 0), break, end b = x + dx; fb = feval(fun,b); end if trace, steps = [a fa; b fb]; end while abs(b - a) > 2.0*tol*max(abs(b) ,1.0) c = a + 0.5*(b-a);fc = feval(fun,c); if trace, steps = [steps; [c fc]]; end if (fb > 0) == (fc > 0) b = c; fb = fc; else a = c; fa = fc; end end

Randint.m
Randint returns an m by n matrix with random entries between a and b. If left empty, a=0 and b=9.

function a = randint(m,n,a,b) if nargin < 3, a = 0; b = 9; end a = floor((b-a+1)*rand(m,n)) + a;

Stat.m
For a vector x, the stat function returns the mean and standard deviation of x.

function [mean, stdev] = stat(x) [m n] = size(x); if m == 1 m = n;   end mean = sum(x)/m; stdev = sqrt(sum(x.^2)/m - mean.^2);

Torus.m
Torus generates a plot of a torus with central radius a and lateral radius r. n controls the number of facets on the surface.

function [x,y,z] = torus(r,n,a) if nargin < 3, a = 1; end if nargin < 2, n = 30; end if nargin < 1, r = 0.5 end theta = pi*(0:2:2*n)/n; phi = 2*pi*(0:2:n)'/n; xx = (a + r*cos(phi))*cos(theta); yy = (a + r*cos(phi))*sin(theta); zz = r*sin(phi)*ones(size(theto)); if nargout == 0 surf(xx,yy,zz) ar = (a + r)/sqrt(2); axis([-ar,ar,-ar,ar,-ar,ar]) else x = xx; y = yy; z = zz; end

R.1.2a
$$

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

Problem R1.3
Description:	For the model 2-D truss with 2 inclined elements, provide the numerical values for the element stiff matrices and find the coordinates of the global nodes 1 and 3, with the global node 2 at the origin of the global coordinate system.



Img1.3.1: Obtained from the Eml4500.f08.1.djvu p.13



Img1.3.2: Obtained from the Eml4500.f08.1.djvu p.17



Img1.3.3: Obtained from the Eml4500.f08.1.djvu p.16

$$

\begin{array}{lcr} \Theta_{1} = 30^{0} \\

\Theta_{2} = \frac{\pi}{4}

\end{array} $$

>> l1 = cos(30*pi/180);

>> m1 = sin(30*pi/180);

>> k1 = zeros(4,4)

>> k1(1,1) = l1^2;

>> k1(2,2) = m1^2;

>> k1(3,3) = l1^2;

>> k1(4,4) = m1^2;

>> k1(1,2) = l1*m1;

>> k1(1,3) = -(l1^2);

>> k1(1,4) = -(l1*m1);

>> k1(2,3) = -(l1*m1);

>> k1(2,4) = -(m1^2);

>> k1(3,4) = l1*m1;

>> for i = 1:4

for j = 1:4

if k1(i,j)==0

k1(i,j)=k1(j,i)

end

end

end

>> k1

k1 =

0.7500   0.4330   -0.7500   -0.4330    0.4330    0.2500   -0.4330   -0.2500   -0.7500   -0.4330    0.7500    0.4330   -0.4330   -0.2500    0.4330    0.2500

$$

\begin{array}{lcr} k^{(1)}= \frac{E^{(1)}*A^{(1)}}{L^{(1)}}=\frac{(3)*(1)}{(4)}=0.75 \end{array} $$

>> k1 = k1*0.75

k1 =

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

>> l2 = cos(pi/4);

>> m2 = sin(pi/4);

>> k2(1,1) = l2^2;

>> k2(2,2) = m2^2;

>> k2(3,3) = l2^2;

>> k2(4,4) = m2^2;

>> k2(1,2) = l2*m2;

>> k2(1,3) = -(l2^2);

>> k2(1,4) = -(l2*m2);

>> k2(2,3) = -(l2*m2);

>> k2(2,4) = -(m2^2);

>> k2(3,4) = l2*m2;

>> for i = 1:4

for j = 1:4

if k2(i,j)==0

k2(i,j)=k2(j,i)

end

end

end

k2 =

0.5000   0.5000   -0.5000   -0.5000    0.5000    0.5000   -0.5000   -0.5000   -0.5000   -0.5000    0.5000    0.5000   -0.5000   -0.5000    0.5000    0.5000

$$

\begin{array}{lcr} k^{(2)}= \frac{E^{(2)}*A^{(2)}}{L^{(2)}}=\frac{(5)*(2)}{(2)}=5 \end{array} $$

>> k2 = k2*5

k2 =

2.5000   2.5000   -2.5000   -2.5000    2.5000    2.5000   -2.5000   -2.5000   -2.5000   -2.5000    2.5000    2.5000   -2.5000   -2.5000    2.5000    2.5000

The coordinated of global nodes one and three can be represented as seen below with the global node two being the origin.

$$

\begin{array}{lcr} node_{2} = (0 \hat{\underline{i}}, 0 \hat{\underline{j}}) \\

node_{3} = 4*(cos(45^{^{0}}) \hat{\underline{i}}, -sin(45^{^{0}}) \hat{\underline{j}}) = (2.8284 \hat{\underline{i}}, -2.8284 \hat{\underline{j}}) \\

node_{1} = 2*(-cos(30^{^{0}}) \hat{\underline{i}}, -sin(30^{^{0}}) \hat{\underline{j}}) = (-1.7320 \hat{\underline{i}}, -1 \hat{\underline{j}})

\end{array} $$

Problem R1.4
“Description:	Consider a plane truss consisting of three bars with the properties E = 200GPa, A1 = 6.0 x 10-4 m2, A2 = 3.0 x 10-4 m2, A3 = 10.0 x 10-4 m2, and loaded by a single force P = 80 kN. The corresponding finite element model consists of three elements and eight degrees of freedom.” (calfem p.234)



Img1.4: Obtained from the calfem manual (calfem.pdf) p.234

The first step is to define a topology matrix: Edof

>> Edof = [1,1,2,5,6;2,5,6,7,8;3,3,4,5,6]

Edof =

1    1     2     5     6     2     5     6     7     8     3     3     4     5     6

K is the stiffness matrix and f is the force matrix

>> K = zeros(8);

>> f = zeros(8,1);

>> f(6)=-80e3

f =

0          0           0           0           0      -80000           0           0

E is the modulus of elasticity. A1, A2, and A3 are the corresponding areas of each truss. Combined, we make the property matrices for each truss in the system, defined as ep1, ep2, and ep3.

>> E=2.0e11;

>> A1=6.0e-4;

>> A2=3.0e-4;

>> A3=10.0e-4;

>> ep1=[E A1]

ep1 =

1.0e+011 *

2.0000   0.0000

>> ep2=[E A2];

>> ep3=[E A3];

The coordinates of each connection of the truss system are defined below.

>> ex1=[0,1.6];

>> ex2=[1.6,1.6];

>> ex3=[0,1.6];

>> ey1=[0,0];

>> ey2=[0,1.2];

>> ey3=[1.2,0];

Using the calfem34 designed function bar2e, we are able to construct the element stiffness matrices as seen below.

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

Ke1 =

1.0e+007 *

7.5000        0   -7.5000         0         0         0         0         0   -7.5000         0    7.5000         0         0         0         0         0

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

Ke2 =

1.0e+007 *

0        0         0         0         0    5.0000         0   -5.0000         0         0         0         0         0   -5.0000         0    5.0000

>> Ke3=bar2e(ex3,ey3,ep3)

Ke3 =

64000000  -48000000   -64000000   48000000   -48000000    36000000    48000000  -36000000   -64000000    48000000    64000000  -48000000    48000000   -36000000   -48000000   36000000

The next step is to create the global stiffness matrix by using the topology and assem function.

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

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

>> K=assem(Edof(3,:),K,Ke3);

>> K

K =

1.0e+008 *

0.7500        0         0         0   -0.7500         0         0         0         0         0         0         0         0         0         0         0         0         0    0.6400   -0.4800   -0.6400    0.4800         0         0         0         0   -0.4800    0.3600    0.4800   -0.3600         0         0   -0.7500         0   -0.6400    0.4800    1.3900   -0.4800         0         0         0         0    0.4800   -0.3600   -0.4800    0.8600         0   -0.5000         0         0         0         0         0         0         0         0         0         0         0         0         0   -0.5000         0    0.5000

Using the physical displacements for this problem, the solveq function is used to solve the system of equations for the displacement vector a and the support force vector r.

>> bc= [1,0;2,0;3,0;4,0;7,0;8,0];

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

a =

0        0         0         0   -0.0004   -0.0012         0         0

r =

1.0e+004 *

2.9845        0   -2.9845    2.2383    0.0000    0.0000         0    5.7617

The normal forces for section 1, 2, and 3 are calculated using the bar2s function that came with the calfem34 toolbox.

>> ed1=extract(Edof(1,:),a);

>> N1=bar2s(ex1,ey1,ep1,ed1)

N1 =

-2.9845e+004

>> ed2=extract(Edof(2,:),a);

>> N2=bar2s(ex2,ey2,ep2,ed2)

N2 =

5.7617e+004

>> ed3=extract(Edof(3,:),a);

>> N3=bar2s(ex3,ey3,ep3,ed3)

N3 =

3.7306e+004

Problem R1.5














Problem R1.6


BC is a zero force member D= 2 in σ=F/A=F/(πr^2 )

Summing the moments about joint A and solving for the reaction force at joint F $$ \sum M_a= -F_{fy} (36)+400(9)-1200(24)=0 $$ $$ F_{fy} (36)=32400$$ $$ F_{fy}=900 lb$$

Now we can solve for the reaction forces at joint A $$ \sum F_y= F_{ay}-1200+900=0$$ $$ F_{ay}=300 lb $$ $$ \sum F_x= F_{ax}-400=0$$ $$ F_{ax}=400 lb$$

Making a FBD of joint A,F,D,E and summing forces in the x and y direction we can solve for the axial forces in the two force members. Using this force and the known cross sectional area we can solve for the stiffness.

Joint A FBD: $$ \sum F_x = -400 + F_{ac} - F_{ab} \frac{12}{15} = 0$$ $$ \sum F_y = 300 - F_{ab} \frac{9}{15} = 0$$ $$ F_{ab} =300 \frac{9}{15} =500 lb$$ $$ F_{ac} = 500 \frac{12}{15} +400 = 800$$ $$ \sigma_{ab} = \frac{500lb}{\pi(1in)^2}=159.1549 \frac{lb}{in^2}= 159.2psi$$ $$ \sigma_{ac} =\frac{800}{\pi(1in)^2} =254.6479psi$$

Joint F FBD: $$ \Sigma F_x = -F_{ef} +F_{df} \frac{12}{15} = 0$$ $$ \Sigma F_y = 900 - F_{df}\frac{9}{15} = 0$$ $$ F_{df} =900\frac{15}{9} =1500lb$$ $$ F_{ef}=1500\frac{12}{15} =1200lb$$ $$ \sigma_{df} =\frac{1500lb}{\pi(1in)^2}=477.5psi$$ $$ \sigma_{ef}=\frac{1200}{\pi(1in)^2}=382psi$$

Joint D FBD: $$ \Sigma F_x = 400-F_{bd} -1500\frac{12}{15} = 0$$ $$ \Sigma F_y = -F_{de}+1500\frac{9}{15} = 0$$ $$ F_{bd} =1500\frac{12}{15}-400 =800lb$$ $$ F_{de} =1500\frac{9}{15} =900lb$$ $$ \sigma_{bd}=\frac{800lb}{\pi(1in)^2}=254.6psi $$ $$ \sigma_{de} =\frac{900lb}{\pi(1in)^2}=286.7psi$$

Joint E FBD: $$ \Sigma F_x = 1200 -F_{ce} -F_{be}\frac{12}{15} = 0$$ $$ \Sigma F_y =900 -1200 +F_{be}\frac{9}{15}=0$$ $$ F_{be} =(-900+1200)\frac{15}{9} =500lb$$ $$ F_{ce} =-500\frac{12}{15} +1200 =800lb$$ $$ \sigma_{be} = \frac{500lb}{\pi(1in)^2} = 159.2psi$$ $$ \sigma_{ce} = \frac{800}{\pi(1in)^2}=254.6psi$$

Array showing all finalized answers $$ \begin{bmatrix} Member & \sigma (psi) & \\ AB &159.2 &Compression \\AC &254.6 &Tension \\BD &254.6 &Tension \\BE &159.2 &Tension \\BC &0 &n/a \\ CE &254.6  &Tension \\DF &477.5 &Compression \\DE &286.7 &Tension \\EF &382.0 &Tension \end{bmatrix}$$

The following MATLAB code automatically solves for the same values as above. This code verifies our hand calculations were correct.

Coord=[0 0;12 0;24 0;36 0;12 9;24 9]; Dof=[1 2;3 4;5 6;7 8;9 10;11 12]; 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]; [Ex,Ey]=coordxtr(Edof,Coord,Dof,2); K=zeros(12); F=zeros(12,1); F(6)=-1200;F(11)=400; ep=[10000 pi]; for i=1:9 Ke=bar2e(Ex(i,:),Ey(i,:),ep); K=assem(Edof(i,:),K,Ke); end bc=[1 0;2 0;8 0]; Q=solveq(K,F,bc); Ed=extract(Edof,Q); for i=1:9 N(i)=bar2s(Ex(i,:),Ey(i,:),ep,Ed(i,:)); end eldraw2(Ex,Ey); plotpar=[1 3 0];scale=1; eldisp2(Ex,Ey,Ed,plotpar,scale); r=1; A=pi*r^2; stressAB=N(4)/A stressAC=N(1)/A stressBD=N(7)/A stressBE=N(6)/A stressBC=0 stressCE=N(2)/A stressDF=N(9)/A stressDE=N(8)/A stressEF=N(3)/A The negative and positive terms denote compression and tension respectively. stressAB = -159.1549 stressAC = 254.6479 stressBD = 254.6479 stressBE = 159.1549 stressBC =    0 stressCE = 254.6479 stressDF = -477.4648 stressDE = 286.4789 stressEF = 381.9719

--EML4507.s13.team4ever.Wulterkens (discuss • contribs) 22:13, 30 January 2013 (UTC)