# University of Florida/Eml4507/Team7 Report2

## Problem 1

#### Find

Verify the equilibrium of node 2; Explain the assembly process; Using CALFEM verify the results of the 2-D truss system with 2 inclined members [1]

#### Given

${\displaystyle P=7;f_{3}^{(1)}=4.4378f_{4}^{(1)}=2.5622;f_{1}^{(2)}=-4.4378;f_{2}^{(2)}2=4.4578}$

#### Solution

##### Equilibrium of Node 2
 {\displaystyle {\begin{aligned}&\Sigma F_{x}=0\\\end{aligned}}} (1.1)
 {\displaystyle {\begin{aligned}&\left(-f_{3}^{(1)}\right)+0+\left(-f_{1}^{(2)}\right)=0\\\end{aligned}}} (1.2)

{\displaystyle {\begin{aligned}&-4.4378+0+4.4378=0\\\end{aligned}}}

 {\displaystyle {\begin{aligned}&\Sigma F_{y}=0\\\end{aligned}}} (1.3)
 {\displaystyle {\begin{aligned}&\left(-f_{4}^{(1)}\right)+P+\left(-f_{2}^{(2)}\right)=0\\\end{aligned}}} (1.4)

{\displaystyle {\begin{aligned}&-2.5622+7-4.4578=0\\\end{aligned}}}

##### Explanation for the Assembly Process

To assemble the global stiffness matrix for the structure; sum the force-displacement matrices in global coordinates:

 ${\displaystyle {\begin{Bmatrix}F_{1}\\F_{2}\\F_{3}\\F_{4}\\F_{5}\\F_{6}\\\end{Bmatrix}}={\begin{bmatrix}f_{1}^{(1)}\\f_{2}^{(1)}\\f_{3}^{(1)}\\f_{4}^{(1)}\\f_{5}^{(1)}\\f_{6}^{(1)}\\\end{bmatrix}}+{\begin{bmatrix}f_{1}^{(2)}\\f_{2}^{(2)}\\f_{3}^{(2)}\\f_{4}^{(2)}\\f_{5}^{(2)}\\f_{6}^{(2)}\\\end{bmatrix}}}$ (1.5)

Substituting for the force displacement matrix for each element:

 ${\displaystyle {\begin{Bmatrix}F_{1}\\F_{2}\\F_{3}\\F_{4}\\F_{5}\\F_{6}\\\end{Bmatrix}}={\begin{bmatrix}k_{11}^{(1)}&k_{12}^{(1)}&k_{13}^{(1)}&k_{14}^{(1)}&k_{15}^{(1)}&k_{16}^{(1)}\\k_{21}^{(1)}&k_{22}^{(1)}&k_{23}^{(1)}&k_{24}^{(1)}&k_{25}^{(1)}&k_{26}^{(1)}\\k_{31}^{(1)}&k_{32}^{(1)}&k_{33}^{(1)}&k_{34}^{(1)}&k_{35}^{(1)}&k_{36}^{(1)}\\k_{41}^{(1)}&k_{42}^{(1)}&k_{43}^{(1)}&k_{44}^{(1)}&k_{45}^{(1)}&k_{46}^{(1)}\\k_{51}^{(1)}&k_{52}^{(1)}&k_{53}^{(1)}&k_{54}^{(1)}&k_{55}^{(1)}&k_{56}^{(1)}\\k_{61}^{(1)}&k_{62}^{(1)}&k_{63}^{(1)}&k_{64}^{(1)}&k_{65}^{(1)}&k_{66}^{(1)}\\\end{bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\\\end{Bmatrix}}+{\begin{bmatrix}k_{11}^{(2)}&k_{12}^{(2)}&k_{13}^{(2)}&k_{14}^{(2)}&k_{15}^{(2)}&k_{16}^{(2)}\\k_{21}^{(2)}&k_{22}^{(2)}&k_{23}^{(2)}&k_{24}^{(2)}&k_{25}^{(2)}&k_{26}^{(2)}\\k_{31}^{(2)}&k_{32}^{(2)}&k_{33}^{(2)}&k_{34}^{(2)}&k_{35}^{(2)}&k_{36}^{(2)}\\k_{41}^{(2)}&k_{42}^{(2)}&k_{43}^{(2)}&k_{44}^{(2)}&k_{45}^{(2)}&k_{46}^{(2)}\\k_{51}^{(2)}&k_{52}^{(2)}&k_{53}^{(2)}&k_{54}^{(2)}&k_{55}^{(2)}&k_{56}^{(2)}\\k_{61}^{(2)}&k_{62}^{(2)}&k_{63}^{(2)}&k_{64}^{(2)}&k_{65}^{(2)}&k_{66}^{(2)}\\\end{bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\\\end{Bmatrix}}}$ (1.6)

Substituting zero for where there is no stiffness for each element and changing to local coordinates:

 ${\displaystyle {\begin{Bmatrix}F_{1}\\F_{2}\\F_{3}\\F_{4}\\F_{5}\\F_{6}\\\end{Bmatrix}}={\begin{bmatrix}k_{11}^{(1)}&k_{12}^{(1)}&k_{13}^{(1)}&k_{14}^{(1)}&0&0\\k_{21}^{(1)}&k_{22}^{(1)}&k_{23}^{(1)}&k_{24}^{(1)}&0&0\\k_{31}^{(1)}&k_{32}^{(1)}&k_{33}^{(1)}&k_{34}^{(1)}&0&0\\k_{41}^{(1)}&k_{42}^{(1)}&k_{43}^{(1)}&k_{44}^{(1)}&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\\end{bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\\\end{Bmatrix}}+{\begin{bmatrix}0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&k_{11}^{(2)}&k_{12}^{(2)}&k_{13}^{(2)}&k_{14}^{(2)}\\0&0&k_{21}^{(2)}&k_{22}^{(2)}&k_{23}^{(2)}&k_{24}^{(2)}\\0&0&k_{31}^{(2)}&k_{32}^{(2)}&k_{33}^{(2)}&k_{34}^{(2)}\\0&0&k_{41}^{(2)}&k_{42}^{(2)}&k_{44}^{(2)}&k_{44}^{(2)}\\\end{bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\\\end{Bmatrix}}}$ (1.7)

Summing the matrices:

 ${\displaystyle {\begin{Bmatrix}F_{1}\\F_{2}\\F_{3}\\F_{4}\\F_{5}\\F_{6}\\\end{Bmatrix}}={\begin{bmatrix}k_{11}^{(1)}&k_{12}^{(1)}&k_{13}^{(1)}&k_{14}^{(1)}&0&0\\k_{21}^{(1)}&k_{22}^{(1)}&k_{23}^{(1)}&k_{24}^{(1)}&0&0\\k_{31}^{(1)}&k_{32}^{(1)}&\left(k_{33}^{(1)}+k_{11}^{(2)}\right)&\left(k_{34}^{(1)}+k_{12}^{(2)}\right)&k_{13}^{(2)}&k_{14}^{(2)}\\k_{41}^{(1)}&k_{42}^{(1)}&\left(k_{43}^{(1)}+k_{21}^{(2)}\right)&\left(k_{44}^{(1)}+k_{22}^{(2)}\right)&k_{23}^{(2)}&k_{24}^{(2)}\\0&0&k_{31}^{(2)}&k_{32}^{(2)}&k_{33}^{(2)}&k_{34}^{(2)}\\0&0&k_{41}^{(2)}&k_{42}^{(2)}&k_{43}^{(2)}&k_{44}^{(2)}\\\end{bmatrix}}{\begin{Bmatrix}d_{1}\\d_{2}\\d_{3}\\d_{4}\\d_{5}\\d_{6}\\\end{Bmatrix}}}$ (1.8)
##### CALFEM for 2-D Truss System

MATLAB script file:


Edof=[1 1 2 3 4;
2 3 4 5 6;];
bc=[1 0;2 0;5 0;6 0;]
f=zeros(6,1);
f(4)=7;
K=zeros(6);
x1=-2*sqrt(3); y1=-2;
x2=0; y2=0;
x3=sqrt(2); y3=-sqrt(2);
A1=1; A2=2;
E1=3; E2=5;
ex1=[x1 x2]; ey1=[y1 y2]; ep1=[E1 A1];
ex2=[x2 x3]; ey2=[y2 y3]; ep2=[E2 A2];

Ke1=bar2e(ex1,ey1,ep1)
Ke2=bar2e(ex2,ey2,ep2)

K=assem(Edof(1,:),K,Ke1);
K=assem(Edof(2,:),K,Ke2)

a=solveq(K,f,bc);

ed1=extract(Edof(1,:),a)
ed2=extract(Edof(2,:),a)

f1=Ke1*ed1'
f2=Ke2*ed2'


MATLAB output:


>> r2p1

bc =

1     0
2     0
5     0
6     0

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

Ke2 =

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

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

ed1 =

0         0    4.3520    6.1271

ed2 =

4.3520    6.1271         0         0

f1 =

-4.4378
-2.5622
4.4378
2.5622

f2 =

-4.4378
4.4378
4.4378
-4.4378


## Problem 2

#### Find

Homogenous L2-ODE-CC, damping ratio zeta, plot of the solution, log decrement delta, quality factor Q and loss factor eta.

#### Given

Consider the spring-mass damper system on p.53-2.

${\displaystyle \lambda _{1}=-0.5+i2,\lambda _{2}=-0.5-i2}$


Initial Conditions:

${\displaystyle y(0)=1,y'(0)=0}$


#### Solution

Case: Underdamped

${\displaystyle \Delta =a^{2}-4b<0}$
${\displaystyle \alpha =-0.5,\beta =2}$
${\displaystyle (\lambda -(-0.5+2i))(\lambda -(-0.5-2i))=0}$
${\displaystyle \lambda ^{2}+1\lambda +2.25=0}$

${\displaystyle a=1,b=2.25}$

 ${\displaystyle \omega ^{2}=b}$ (2.1)
 ${\displaystyle 2\omega \zeta =a}$ (2.2)

Solving equation (2.1) and (2.2)
${\displaystyle \omega =1.5}$
${\displaystyle \zeta =0.333}$

Therefore L2-ODE-CC: ${\displaystyle y_{h}''+y_{h}'+2.25y_{h}=0}$

${\displaystyle y=C_{1}e^{(\alpha x)}cos(\beta x)+C_{2}e^{(\alpha x)}sin(\beta x)}$

 ${\displaystyle y=C_{1}e^{(-0.5x)}cos(2x)+C_{2}e^{(-0.5x)}sin(2x)}$ (2.3)
 ${\displaystyle y'=-0.5C_{1}e^{(-0.5x)}cos(2x)-2C_{1}e^{(-0.5x)}sin(2x)-0.5C_{2}e^{(-0.5x)}sin(2x)+2C_{2}e^{(-0.5x)}cos(2x)}$ (2.4)

Using the two initial conditions and solving equation (2.3) and (2.4) the coefficients are as follows: ${\displaystyle C_{1}=1,C_{2}=.25}$

Therefore,
${\displaystyle y=e^{(-0.5x)}cos(2x)+0.25e^{(-0.5x)}sin(2x)}$

The next step in analyzing the model is comparing and evaluating Log Decrements. This can be done numerically and by a formula. The log decrement in this case was done both ways and then compared to verify the value.
Numerical method: ${\displaystyle \delta ^{(n)}=log[(y_{h}^{(n)})/(y_{h}^{(n+1)})]}$

The log decrement was calculated by finding the average value: ${\displaystyle \delta =[1/3](\delta ^{(1)}+\delta ^{(2)}+\delta ^{(3)})}$

Verifying the log decrement value by using the formula: ${\displaystyle \delta =(2\pi \zeta )/{\sqrt {1-\zeta ^{2}}}}$

Comparing the two values: ${\displaystyle \delta =2.21}$

Finding the quality factor : ${\displaystyle Q=\delta /\pi }$
${\displaystyle Q=0.703}$

Finding the loss factor: ${\displaystyle \eta =1/Q}$
${\displaystyle \eta =1.42}$

## Problem 3

#### Find

Determine the member force and axial stress in each member of the truss shown in the figure using one of the FE analysis programs in the Appendix. Assume that Young’s modulus is 104 psi and all cross-sections are circular with a diameter of 2 in. Compare the results with the exact solutions that are obtained from the free-body diagram.[2]

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

##### By Hand Calculations
 ${\displaystyle \sum F_{x}=0=400lb+T_{1x}}$ (3.1)

${\displaystyle T_{1x}=-400lb}$

 ${\displaystyle \sum M_{1}=0=-1200lb(24')-400lb(9')+T_{4y}(36')}$ (3.2)

${\displaystyle T_{4y}=900lb}$

 ${\displaystyle \sum F_{y}=0=-1200lb+900lb+T_{1y}}$ (3.3)

${\displaystyle T_{1y}=300lb}$

 ${\displaystyle \sum M_{6}=0=-300lb(12')-400lb(9')+F_{8}(9')}$ (3.4)

${\displaystyle F_{8}=800lb}$

 ${\displaystyle \sum M_{3}=0=-300lb(24')+F_{2}(9')}$ (3.5)

${\displaystyle F_{8}=800lb}$

 ${\displaystyle \sum F_{y}=0=300lb-{\frac {3}{5}}F_{5}}$ (3.6)

${\displaystyle F_{5}=500lb}$

##### By MATLAB

MATLAB script file:


%X= [element number; EA/L; l; m];
X = [1:9; 174.44 218.06 174.44 290.74 174.44 290.74 218.06 218.06 218.06; 0.80 -1.00 -0.80 0.00 -0.80 0.00 1.00 1.00 1.00;
0.60 0.00 0.60 1.00 0.60 1.00 0.00 0.00 0.00]
i=1;
for i=1:9
str=fprintf('For element %d\n',i);
k=X(2,i).*[X(3,i)^2 X(3,i)*X(4,i) -(X(3,i)^2) -X(3,i)*X(4,i); X(3,i)*X(4,i) X(4,i)^2 -X(3,i)*X(4,i) -(X(4,i)^2);
-(X(3,i)^2)-X(3,i)*X(4,i) X(3,i)^2 X(3,i)*X(4,i); -X(3,i)*X(4,i) -(X(4,i)^2) X(3,i)*X(4,i) X(4,i)^2]
i=i+1;
end


MATLAB output:

>> report2

X =

   1.0000    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000    8.0000    9.0000
174.4400  218.0600  174.4400  290.7400  174.4400  290.7400  218.0600  218.0600  218.0600
0.8000   -1.0000   -0.8000         0   -0.8000         0    1.0000    1.0000    1.0000
0.6000         0    0.6000    1.0000    0.6000    1.0000         0         0         0


For element 1

k =

 111.6416   83.7312 -111.6416  -83.7312
83.7312   62.7984  -83.7312  -62.7984
-111.6416  -83.7312  111.6416   83.7312
-83.7312  -62.7984   83.7312   62.7984


For element 2

k =

 218.0600         0 -218.0600         0
0         0         0         0
-218.0600         0  218.0600         0
0         0         0         0


For element 3

k =

 111.6416  -83.7312 -111.6416   83.7312
-83.7312   62.7984   83.7312  -62.7984
-111.6416   83.7312  111.6416  -83.7312
83.7312  -62.7984  -83.7312   62.7984


For element 4

k =

        0         0         0         0
0  290.7400         0 -290.7400
0         0         0         0
0 -290.7400         0  290.7400


For element 5

k =

 111.6416  -83.7312 -111.6416   83.7312
-83.7312   62.7984   83.7312  -62.7984
-111.6416   83.7312  111.6416  -83.7312
83.7312  -62.7984  -83.7312   62.7984


For element 6

k =

        0         0         0         0
0  290.7400         0 -290.7400
0         0         0         0
0 -290.7400         0  290.7400


For element 7

k =

 218.0600         0 -218.0600         0
0         0         0         0
-218.0600         0  218.0600         0
0         0         0         0


For element 8

k =

 218.0600         0 -218.0600         0
0         0         0         0
-218.0600         0  218.0600         0
0         0         0         0


For element 9

k =

 218.0600         0 -218.0600         0
0         0         0         0
-218.0600         0  218.0600         0
0         0         0         0


## Problem 4

### Part 1

#### Find

For each case, derive the corresponding standard L2-ODE-CC, find the overall homogeneous solution, and plot the solution in a separate figure; then plot all 3 solutions in the same figure for comparison.

#### Given

Two Roots: ${\displaystyle \lambda _{1}=-0.5,\lambda _{2}=-1.5}$

Initial conditions y(0)=1, y'(0)=0

#### Solution

Since \lambda_1 and \lambda_2 are unequal and negative, the system is overdamped. Currently we have:

${\displaystyle (\lambda +0.5)(\lambda +1.5)=0}$

Standard form

${\displaystyle (\lambda ^{2}+2*\lambda +0.75)=0}$

The standard homogeneous L2-ODE-CC, or the equation of motion, is therefore:

${\displaystyle y_{h}''+2y_{h}'+0.75y_{h}=0}$

The general solution is given by:

${\displaystyle y_{h}(t)=c_{1}e^{(-0.5t)}+c_{2}e^{(-1.5t)}}$

Using the given initial conditions, we solve for the constants:

${\displaystyle y^{'}(0)=0=-0.5c_{1}-1.5c_{2}}$

${\displaystyle c_{1}={\frac {3}{4}},c_{2}=-{\frac {1}{4}}}$

Plugging these in to find the overall homogeneous solution:

${\displaystyle y_{p}(t)={\frac {3}{4}}e^{(-0.5t)}-{\frac {1}{4}}e^{(-1.5t)}}$

The graph of the solution is:

### Part 3

#### Find

Derive the corresponding standard L2-ODE-CC, find the overall homogeneous solution, and plot the solution in a separate figure; then plot all 3 solutions in the same figure for comparison.

#### Given

Two roots: ${\displaystyle \lambda _{1}=-0.5+i2,\lambda _{2}=-1.5-i2}$

Initial conditions:

${\displaystyle y(0)=1,y'(0)=v_{0}=0}$

#### Solution

From the characteristic equation:

${\displaystyle (\lambda -\lambda _{1})(\lambda -\lambda _{2})=0}$

${\displaystyle (\lambda +0.5-i2)(\lambda +0.5+i2)=0}$

Multiplying out to the standard form:

${\displaystyle \lambda ^{2}+\lambda +4.25=0}$

${\displaystyle \lambda ^{2}+a\lambda +b=0}$

Where:

${\displaystyle a=c/m=1andb=k/m=4.25}$

The standard homogeneous L2-ODE-CC, or the equation of motion, is therefore:

${\displaystyle y_{h}''+y_{h}'+4.25y_{h}=0}$

The general solution is given by:

${\displaystyle y_{h}(t)=Ae^{-at/2}cos(\omega _{d}-\phi )}$

The variables for this equation are found to be:

${\displaystyle \omega _{n}={\sqrt {b}}=2.06[itex]rad/s[itex]\zeta =a/(2\omega _{n})=0.243}$

${\displaystyle \omega _{d}=\omega _{n}{\sqrt {1-\zeta ^{2}}}=2.00}$ rad/s

${\displaystyle A={\sqrt {((v_{0}+\zeta \omega _{n}y_{0})^{2}+(y_{0}\omega _{d})^{2})/(\omega _{d})^{2}}}=1.03}$

${\displaystyle \phi =tan^{-1}(y_{0}\omega _{d})/(v_{0}+\zeta \omega _{n}y_{0}))=75.9}$

Therefore:

${\displaystyle y_{h}(t)=1.03e^{-0.5t}cos(2t-75.9)}$

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

## Problem 5

#### Find

The general procedure in linear finite element calculations.

Generate element matrices, assemble element matrices into the Global System of equations. Solve the global system of equations and evaluate element forces.

#### Solution

Using Matlab, we first define the topology matrix by Edof containing element numbers and global degrees of freedom

>> Edof=[1 1 2; 2 2 3; 3 2 3]
1     1     2
2     2     3
3     2     3


Now the global stiffness matrix K ${\displaystyle (3x3)}$ of zeros.

>> K=zeros(3,3)
0     0     0
0     0     0
0     0     0


The load vector ${\displaystyle f(3x1)}$ with load ${\displaystyle F=100}$ and load it in position 2.

>> f=zeros(3,1); f(2)=100
0
100
0


Stiffness matrices are generated by using the function spring1e
The function is given by:

function [Ke]=spring1e(ep);
k = ep;
Ke = [ k -k;
-k  k];


Now we can set the element property as ep for the springs with their respective spring constants ${\displaystyle k}$ and ${\displaystyle 2k}$ Where ${\displaystyle k=1500}$

Commands in matlab are as follows:

>> k=1500; ep1=k; ep2=2*k;
>> Ke1=spring1e(ep1)
Ke1 =
1500       -1500
-1500        1500
>> Ke2=spring1e(ep2)
Ke2 =
3000       -3000
-3000        3000


With the assem function we can assemble stiffness matrices into the global stiffness matrix K. Which leads too:

assem function:

function [K,f]=assem(edof,K,Ke,f,fe)
[nie,n]=size(edof);
t=edof(:,2:n);
for i = 1:nie
K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke;
if nargin==5
f(t(i,:))=f(t(i,:))+fe;
end
end


Stiffness matrix K:

>> K=assem(Edof(1,:),K,Ke2)
K =
3000       -3000           0
-3000        3000           0
0           0           0
>> K=assem(Edof(2,:),K,Ke1)
K =
3000       -3000           0
-3000        4500       -1500
0       -1500        1500
>> K=assem(Edof(3,:),K,Ke2)
K =
3000       -3000           0
-3000        7500       -4500
0       -4500        4500


The Global system of equations is solved with boundary condtions given in bc as well as with function solveq.

Function solveq is as follows:

function [d,Q]=solveq(K,f,bc)
if nargin==2 ;
d=K\f ;
elseif nargin==3;
[nd,nd]=size(K);
fdof=[1:nd]';
d=zeros(size(fdof));
Q=zeros(size(fdof));
pdof=bc(:,1);
dp=bc(:,2);
fdof(pdof)=[];
s=K(fdof,fdof)\(f(fdof)-K(fdof,pdof)*dp);
d(pdof)=dp;
d(fdof)=s;
end
Q=K*d-f;



With Matlab commands it results in:

>> bc=[1 0; 3 0];
>> [a,r]=solveq(K,f,bc)
a =
0
0.0133
0
r =
-40.0000
0
-60.0000


Element forces are now computed from the element displacements. These are obtained from global displacements a using the function extract.

Function extract:

function [ed]=extract(edof,a)
[nie,n]=size(edof);
t=edof(:,2:n);
for i = 1:nie
ed(i,1:(n-1))=a(t(i,:))';
end


Now the corresponding Matlab commands are as follows:

>> ed1=extract(Edof(1,:),a)
ed1 =
0    0.0133
>> ed2=extract(Edof(2,:),a)
ed2 =
0.0133         0
>> ed3=extract(Edof(3,:),a)
ed3 =
0.0133         0


The last step is to evaluate the spring forces by using the function spring1s.

We have function spring1s as written below:

function [es]=spring1s(ep,ed)
k = ep;
u=ed;
es = k*(u(2)-u(1));


We can now finally evaluate the spring forces. Using Matlab the forces are as follows:

>> es1=spring1s(ep2,ed1)
es1 =
40
>> es2=spring1s(ep1,ed2)
es2 =
-20
>> es3=spring1s(ep2,ed3)
es3 =
-40


• Note: Instructions and commands were all followed from CALFEM 3.4 Manual.

## Problem 6

#### Find

Find the displacements of the three bodies and the forces ( tensile/compressive) in the springs. What are the reactions at the walls? Assume only horizontal translation.

#### Solution

The nodes and elements are labeled in the problem. Create a connectivity table.

Construct Elemental equilibrium equations for all 6 elements.

Here is a sample equilibrium equation for element 3:

 ${\displaystyle {\binom {f_{2}^{(3)}}{f_{3}^{(3)}}}=k_{3}{\begin{bmatrix}1&-1\\-1&1\\\end{bmatrix}}{\begin{bmatrix}u_{2}\\u_{3}\\\end{bmatrix}}}$ (6.1)

The notation for  : The superscript parenthesis signifies the element number, while the subscript signifies the associated node number.

With element forces calculated, a node equilibrium can be established. Start with a free-body diagram of each node, and then create the five nodal equations:

Make sure all forces are positive and modeled in tension.

Combine your elemental equilibrium equations with your nodal equilibrium equations for the global stiffness correlation:

 ${\displaystyle [F_{s}]=[Q_{s}][K_{s}]}$ (6.2)
 ${\displaystyle {\begin{bmatrix}F_{1}=R_{1}\\F_{2}=0\\F_{3}=1000\\F_{4}=0\\F_{5}=R_{5}\end{bmatrix}}={\begin{bmatrix}k_{1}+k_{4}&-k_{1}&-k_{4}&0&0\\-k_{1}&k_{1}+k_{2}+k_{3}&-k_{3}&-k_{2}&0\\-k_{4}&-k_{3}&k_{3}+k_{4}+k_{5}&-k_{5}&0\\0&-k_{4}&-k_{4}&k_{2}+k_{5}+k_{6}&0\\0&0&0&-k_{6}&k_{6}\end{bmatrix}}{\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\\u_{5}\\\end{bmatrix}}}$ (6.3)

##### Solving via MATLAB
%Spring constants - in N/mm
k1 = 500;
k2 = 400;
k3 = 600;
k4 = 200;
k5 = 400;
k6 = 300;

%Declaration of Global Displacement Variables
u1 = 0; %due to being at reactions
u2 = 1; %placeholder
u3 = 1; %placeholder
u4 = 1; %placeholder
u5 = 0; %due to being at reactions

%Combining the Nodal equilibrium and the Element Equilibrium, we can find
%the displacements

kglobal =[k1+k4   -k1     -k4     0     0;
-k1  k1+k2+k3  -k3    -k2    0;
-k4    -k3   k3+k4+k5 -k5    0;
0     -k2     -k5 k2+k5+k6 -k6;
0      0       0    -k6     k6;]

%Rows 1 and 5 can be excluded for now, because the displacements at node 1
%and 5 are zero. This makes the corresponding rows and columns irrelevant

kglobalsimplified = [k1+k2+k3  -k3    -k2;
-k3     k3+k4+k5 -k5;
-k2       -k5  k2+k5+k6;]

%kglobalsimplified*displacements = External Forces

Fsimplified = [0; %no external force
1000; %given
0;]; %no external force

Displacementsimplified = kglobalsimplified\Fsimplified %u2, u3, u4

u2 = Displacementsimplified(1)
u3 = Displacementsimplified(2)
u4 = Displacementsimplified(3)

klocal = [1 -1;
-1 1;]; %Local Equilibrium equations are multiplied by this

%Element Equilibrium Equations
k1ef = k1*klocal*[u1;u2]; %k1ef means K1 Element Forces
k2ef = k2*klocal*[u2;u4];
k3ef = k3*klocal*[u2;u3];
k4ef = k4*klocal*[u1;u3];
k5ef = k5*klocal*[u3;u4];
k6ef = k6*klocal*[u4;u5];

%The above matrices are 2x1s
%We need to access the specific values within each matrix
%So we'lll reassign the matrices to new variables

f11 = k1ef(1); %k11 = k(node 1)(element 1)
f21 = k1ef(2);

f22 = k2ef(1);
f42 = k2ef(2);

f23 = k3ef(1);
f33 = k3ef(2);

f14 = k4ef(1);
f34 = k4ef(2);

f35 = k5ef(1);
f45 = k5ef(2);

f46 = k6ef(1);
f56 = k6ef(2);

%With our element force variables, we can create the nodal equilibirum eqns
F3 = 1000;

R1 = f11 + f14 %Reaction Force
F2 = f21 + f23 + f22; %unknown force
F3 = f33 + f34 + f35; %Given force
F4 = f42 + f45 + f46; %unknown force
R5 = f56 %Reaction Force


Resulting in the answers (in units mm for the displacements and N for the Reaction):

>> run r26.m

u2 =

0.8542

u3 =

1.5521

u4 =

0.8750

R1 =

-737.5000

R5 =

-262.5000

##### Solving via CALFEM
ep1 = 500; %spring constants
ep2 = 400;
ep3 = 600;
ep4 = 200;
ep5 = 400;
ep6 = 300;

%create element spring stiffness matrix
k1 = spring1e(ep1);
k2 = spring1e(ep2);
k3 = spring1e(ep3);
k4 = spring1e(ep4);
k5 = spring1e(ep5);
k6 = spring1e(ep6);

%define degrees of freedom for each element. All elements have 2 dof.
Edof = [1 1 2; % el1 is element 1
2 2 4;
3 2 3;
4 1 3;
5 3 4;
6 4 5;];
fe = zeros(5,1); %declare f, and then populate f with extractions
fe(1) = spring1s(ep1,Edof(1,:));
fe(2) = spring1s(ep2,Edof(2,:));
fe(3) = spring1s(ep3,Edof(3,:));
fe(4) = spring1s(ep4,Edof(4,:));
fe(5) = spring1s(ep5,Edof(5,:));

%Assembling structural stiffness matrix
K = zeros(5);
[K,f] = assem(Edof(1,:),K,k1,f,fe);
[K,f] = assem(Edof(2,:),K,k2,f,fe);%adding element 2 to the structural matrix
[K,f] = assem(Edof(3,:),K,k3,f,fe);%adding element 3 to the structural matrix
[K,f] = assem(Edof(4,:),K,k4,f,fe);
[K,f] = assem(Edof(5,:),K,k5,f,fe);
[K,f] = assem(Edof(6,:),K,k6,f,fe);

a = solveq(K,f); %solves K*a = f

Ed = extract(Edof,a) %extracts element displacements from global
%solution vector a

Ed(3)
Ed(4)
Ed(5)

f(1)
f(5)


Resulting in the answers (in units mm for the displacements and N for the Reaction):

>> run calfem26.m

Ed(3) =

0.8542

Ed(4) =

1.5521

Ed(5) =

0.8750

f(1) =

-737.5000

f(5) =

-262.5000


## References

2. KS 2008 p.103 sec.2.7 pb.21
3. KS 2008 p.103 sec.2.7 pb.21
5. KS 2008 p.98 sec.2.7 pb.21