# University of Florida/Eml4507/s13.team4ever.Wulterkens.R1
## Problem R1.1
### 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 =
EDU>> A = [ 1 2 3 4 5 6 7 8 9 ]
A =
EDU>> A = [1 2;3 4] + i*[5 6;7 8]
A =
EDU>> A = [1+5i 2+6i;3+7i 4+8i]
A =
EDU>> A = [1 2 3;4 5 6;7 8 9]
A =
EDU>> B = [A, zeros(3,2); zeros(2,3), eye(2)]
B =
EDU>> x =;
EDU>> n = 3;
EDU>> for i = 1:n
x = [x,i^2] end
x =
x =
x =
EDU>> x =;
EDU>> for i = n:-1:1
x = [x,i^2] end
x =
x =
x =
EDU>> m = 4
m =
EDU>> for i = 1:m
for j = 1:n
H(i, j) = 1/(i+j-1);
end
end
EDU>> H
H =
EDU>> s = 0;
EDU>> for c = A
s = s + sum(c);
end
EDU>> s
s =
EDU>> n = 0;
EDU>> a = 12;
EDU>> while 2^n < a
n = n + 1;
end
EDU>> n
n =
EDU>> if n < 0
parity = 0;
elseif rem(n,2) == 0
parity = 2;
else parity = 1;
end
EDU>> parity
parity =
EDU>> A
A =
EDU>> y = eig(A)
y =
EDU>> [U,D] = eig(A)
U =
D =
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 =
EDU>> C = 5:-1:1
C =
EDU>> x = [0.0:0.1:2.0]';
EDU>> y = sin(x);
EDU>> [x y]
ans =
EDU>> A(1:3,2)
ans =
EDU>> A(:,3)
ans =
EDU>> A(1:3,:)
ans =
EDU>> A(:,[2 3])
ans =
EDU>> z = randint(4,5)
z =
EDU>> xm = stat(z)
xm =
EDU>> M = magic(6);
EDU>> E = zeros(6,50);
EDU>> for j = 1:50
E(:,j) = eig(M^i);
end
EDU>> E
E =
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 =
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 =
EDU>> F = full(S)
F =
EDU>> m = 6; n = 6; e = ones(n,1); d = -2*e;
EDU>> T = spdiags([e,d,e],[-1,0,1],m,n)
T =
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 =
ans =
EDU>> sparse(i,j,s,m,n)
ans =
EDU>> n = 6; e = floor(10*rand(n-1,1)); E = sparse(2:n,1:n-1,e,n,n)
E =
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 =
fulltime =

### Functions
#### 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.

#### Randint.m

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

#### Stat.m

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

#### 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.

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

## 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 =
${\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 =
>> 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 =
${\begin{array}{lcr}k^{(2)}={\frac {E^{(2)}*A^{(2)}}{L^{(2)}}}={\frac {(5)*(2)}{(2)}}=5\end{array}}$
>> k2 = k2*5
k2 =
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 =
K is the stiffness matrix and f is the force matrix
>> K = zeros(8);
>> f = zeros(8,1);
>> f(6)=-80e3
f =
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 =
>> 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 =
>> Ke2=bar2e(ex2,ey2,ep2)
Ke2 =
>> Ke3=bar2e(ex3,ey3,ep3)
Ke3 =
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 =
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 =
r =
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 =
>> ed2=extract(Edof(2,:),a);
>> N2=bar2s(ex2,ey2,ep2,ed2)
N2 =
>> ed3=extract(Edof(3,:),a);
>> N3=bar2s(ex3,ey3,ep3,ed3)
N3 =

## Problem R1.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

## 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
Now we can solve for the reaction forces at joint A
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:
Joint F FBD:
Joint D FBD:
Joint E FBD:
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.
The negative and positive terms denote compression and tension respectively.