# University of Florida/Eml4507/Team7 Report3

## Problem 1

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


### Find

For the 3-bar truss example in K.2008.2 on pg. 43, solve this problem, but with non-zero prescribed disp dofs, using your own FE matlab code(not CALFEM). Using the same "connectivity array" given in the same problem statement, provide all of the results, including all of the unknown global disp dofs, force components(reactions), and member forces. Plot the deformed shape superposed on the undeformed shape. Verify the results with CALFEM.

### Given

d3= 2 cm, d4= -1 cm,d5= -3 cm, d6= 5 cm

### Solution

MATLAB script file:


d3=2; d4=-1; d5=-3; d6=5;
theta=45;
%X= [element number; EA/L; l; m];
X = [1:3; 174.44 218.06 174.44; 0.80 -1.00 -0.80;
0.60 0.00 0.60]
i=1;
for i=1:3
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:

>> report3

X =

   1.0000    2.0000    3.0000
174.4400  218.0600  174.4400
0.8000   -1.0000   -0.8000


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


F1=

 45.32


F2=

 53.68


F3=

 41.84


F=

 211.45


## Problem 2

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

### Find

1. Draw the FBDs of all components in the spring-mass-damper system in Figure 1 below, with the known displacement degree of freedoms being at the left and right supports. (${\displaystyle d_{3}=0}$ and ${\displaystyle d_{4}=0}$)

2. Derive equations 1, 2, 3, and 4 from the figure above.

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

${\displaystyle \mathbf {d} ={\begin{Bmatrix}d_{1}\\d_{2}\end{Bmatrix}}}$
${\displaystyle \mathbf {F} ={\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}}$

 ${\displaystyle \mathbf {M} ={\begin{bmatrix}m_{1}&0\\0&m_{2}\end{bmatrix}}}$ (2.2)
 ${\displaystyle \mathbf {C} ={\begin{bmatrix}(c_{1}+c_{2})&-c_{2}\\-c_{2}&(c_{2}+c_{3})\end{bmatrix}}}$ (2.3)
 ${\displaystyle \mathbf {K} ={\begin{bmatrix}(k_{1}+k_{2})&-k_{2}\\-k_{2}&(k_{2}+k_{3})\end{bmatrix}}}$ (2.4)

### Solution

FBD of Mass 1:

FBD of Mass 2:

#### Equation Derivations

In order to derive Equation 3.2.1, we can first apply Newton's second law of motion in order to obtain the equations of motion from the two masses.
From the Free body diagram of Mass 1(Figure 2), we can use Newton's second law in order to obtain:

 ${\displaystyle 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}}$ (2.5)

Applying Newton's second law again to Mass 2(Figure 3), we also obtain:

 ${\displaystyle 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}}$ (2.6)

From equations 2.5 and 2.6, the motion of the masses will influence each other. Where, the motion of mass 1 will influence mass 2 and vice versa. Equations 2.5 and 2.6 can be rewritten in matrix form. Therefore, we can now write them as:

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

Which is the same as equation 2.1. Where ${\displaystyle \mathbf {M} }$ is the mass matrix, ${\displaystyle \mathbf {C} }$ is the damping matrix, and ${\displaystyle \mathbf {K} }$ is the stiffness matrix and ${\displaystyle d(t)}$ and ${\displaystyle F(t)}$ are the displacement and force vectors,respectively.Using equations 2.5 and 2.6 the mass, stiffness, and damping matrices can be determined.
The vectors for displacement and force can be represented by:
${\displaystyle \mathbf {d} ={\begin{Bmatrix}d_{1}\\d_{2}\end{Bmatrix}}}$
and

 ${\displaystyle \mathbf {F} ={\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}}$ (2.8)

As for the mass matrix, we see that in Equation 3.2.5 there is only mass 1 present and in equation 3.2.6 only mass 2 is pressent, therefore we have:

 ${\displaystyle \mathbf {M} ={\begin{bmatrix}m_{1}&0\\0&m_{2}\end{bmatrix}}}$ (2.9)

For the damping matrix, it is determined by the coefficients in front of${\displaystyle {\dot {d}}_{1}}$ and ${\displaystyle {\dot {d}}_{2}}$. So we now have:

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

Lastly we determine the stiffness matrix by the coefficients in front of ${\displaystyle d_{1}}$ and ${\displaystyle d_{2}}$. Which is:

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

Essentially, from Equations 2.5 and 2.6 the mass, stiffness, and damping matrices were determined. Equations 2.8 through 2.11 are the same as Equations 2.1 through 2.4, respectively. Without the equation of motions for both masses, neither the simplified matrix form of the equation nor the matrices can be derived.

## Problem 3

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

### Find

FE formulation leads to a system of coupled ODEs in time of the form:

${\displaystyle Mx''+Cx'+Kx=F(t)}$

Integrate the governing system of L2-ODEs-CC in order to generate the time histories for the displacement degrees of freedom. Once the displacement matrix has been found, plot the time histories for the displacement degrees of freedom.

### Given

${\displaystyle m_{1}=3,m_{2}=2}$
${\displaystyle c_{1}=1/2,c_{2}=1/4,c_{3}=1/3}$
${\displaystyle k_{1}=10,k_{2}=20,k_{3}=15}$
${\displaystyle F_{1}(t)=0,F_{2}(t)=0}$
${\displaystyle d_{1}(0)=-1,d_{2}(0)=2}$
${\displaystyle d_{1}'(0)=0,d_{2}'(0)=0}$

### Solution

The CALFEM toolbox is used to solve the system of equations and integrate the equations of motion. First the CALFEM command “eigen” is used to solve for the eigenvalues, ${\displaystyle \lambda }$, and eigenvectors, ${\displaystyle x}$, of the generalized eigenvalue problem:

${\displaystyle Kx=\lambda Mx}$

After entering the matrices into MATLAB, the “eigen” function is run. The MATLAB script file for “eigen” is:

function [L,X]=eigen(K,M,b)
% [L]=eigen(K,M)
% [L]=eigen(K,M,b)
% [L,X]=eigen(K,M)
% [L,X]=eigen(K,M,b)
%-------------------------------------------------------------
% PURPOSE
%  Solve the generalized eigenvalue problem
%  [K-LM]X = 0, considering boundary conditions.
%
% INPUT:
%    K : global stiffness matrix, dim(K)= nd x nd
%    M : global mass matrix, dim(M)= nd x nd
%    b : boundary condition matrix
%        dim(b)= nb x 1
% OUTPUT:
%    L : eigenvalues stored in a vector with length (nd-nb)
%    X : eigenvectors dim(X)= nd x nfdof, nfdof : number of dof's
%-------------------------------------------------------------

% Copyright (c)  Division of Structural Mechanics and
%                Department of Solid Mechanics.
%                Lund Institute of Technology
%-------------------------------------------------------------
[nd,nd]=size(K);
fdof=[1:nd]';
%
if nargin==3
pdof=b(:);
fdof(pdof)=[];
if nargout==2
[X1,D]=eig(K(fdof,fdof),M(fdof,fdof));
[nfdof,nfdof]=size(X1);
for j=1:nfdof;
mnorm=sqrt(X1(:,j)'*M(fdof,fdof)*X1(:,j));
X1(:,j)=X1(:,j)/mnorm;
end
d=diag(D);
[L,i]=sort(d);
X2=X1(:,i);
X=zeros(nd,nfdof);
X(fdof,:)=X2;
else
d=eig(K(fdof,fdof),M(fdof,fdof));
L=sort(d);
end
else
if nargout==2
[X1,D]=eig(K,M);
for j=1:nd;
mnorm=sqrt(X1(:,j)'*M*X1(:,j));
X1(:,j)=X1(:,j)/mnorm;
end
d=diag(D);
[L,i]=sort(d);
X=X1(:,i);
else
d=eig(K,M);
L=sort(d);
end
end
%--------------------------end--------------------------------


The function takes in the K and M matrices and outputs the eigenvalues, L, and eigenvectors, X, shown below:

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

L =

8.2884
39.2116

X =

-0.5458   -0.8379
-0.5925    0.3859


The two eigenvalues are listed in increasing order, the smaller of which is used to find the fundamental circular frequency of the system, ${\displaystyle \omega _{1}}$ , using the formula:

${\displaystyle \lambda _{1}=(\omega _{1})^{2}}$
${\displaystyle \omega _{1}=2.88rad/s}$

The fundamental period of the system, T1, is found through the relation:

${\displaystyle T_{1}=2\pi /\omega _{1}=2.18s}$

The fundamental period is then used to find the time step size, dt:

${\displaystyle dt=T_{1}/10=0.218s}$

The CALFEM command “step2” is then utilized to numerically integrate the equations of motion for time, t, in the following interval:

${\displaystyle 0\leq t\leq 10.91s}$

The function “step2” is an algorithm for dynamic solution of second-order finite elemental equations considering boundary conditions. Its script is shown here:

function [Dsnap,D,V,A]=step2(K,C,M,d0,v0,ip,f,pdisp)
%Dsnap=step2x(K,C,M,d0,v0,ip,f,pdisp)
%[Dsnap,D,V,A]=step2x(K,C,M,d0,v0,ip,f,pdisp)
%-------------------------------------------------------------
% PURPOSE
%  Algorithm for dynamic solution of second-order
%  FE equations considering boundary conditions.
%
% INPUT:
%    K : stiffness matrix, dim(K)= nd x nd
%    C : damping matrix, dim(C)= nd x nd
%    M : mass matrix, dim(M)= nd x nd
%    d0 : initial vector d(0), dim(f)= nd x 1
%    v0 : initial vector v(0), dim(f)= nd x 1
%    ip : [dt tottime alfa delta [nsnap nhist t(i) ...  dof(i) ... ]] :
%         integration and output parameters
%    f : load vector, dim(f)= n x nstep+1,
%        If dim(f)= n x 1 it is supposed that the values
%        are kept constant during time integration.
%    pdisp : boundary condition matrix, dim(pdisp)= nbc x nstep+2,
%            where nbc = number of boundary conditions
%            (constant or time dependent)
%            The first column contains the degrees-of-freedom with  prescribed values
%            and the subsequent columns contain there the time history.
%            If dim(pdisp)= nbc x 2 it is supposed that the values
%            are kept constant during time integration.
%
% OUTPUT:
%    Dsnap : displacement snapshots at 'nsnap' timesteps, specified in ip.
%          the time are also specified in ip. dim(Dsnap)=nd x nsnap
%    D :   solution matrix containing time history displacement
%          d at the selected dof's. dim(D)=nhist x nstep+1
%    V :   solution matrix containing time history of the first time derivative
%          of d at the selected dof's. dim(V)=nhist x nstep+1
%    A :   solution matrix containing time history of the second time derivative
%          of d at the selected dof's. dim(A)=nhist x nstep+1
%-------------------------------------------------------------

% Copyright (c)  Division of Structural Mechanics and
%                Department of Solid Mechanics.
%                Lund Institute of Technology
%-------------------------------------------------------------
% Calfem_com='step2x(K,C,M,d0,v0,ip,f,pdisp)';  Calfem_arg=[7,8];
%   errtest;   if Calfem_fail==1 return; end;
%-------------------------------------------------------------
[nd,nd]=size(K);
[ndc,ndc]=size(C); if (ndc==0); C=zeros(nd,nd);  end
dt=ip(1);  tottime=ip(2);   alfa=ip(3);   delta=ip(4);
b1 = dt*dt*0.5*(1-2*alfa);  b2 = (1-delta)*dt;
b3 = delta*dt;              b4 = alfa*dt*dt;

nstep=1;
[nr nc]=size(f);    if (nc>1);     nstep = nc-1; end
if nargin==7;                      bound=0;      end
if nargin==8
[nr nc]=size(pdisp); if (nc>2);  nstep = nc-2; end
bound=1;             if (nc==0); bound=0;      end
end
ns=tottime/dt;   if (ns < nstep | nstep==1); nstep=ns;     end

[nr nc]=size(f);
tf = zeros(nd,nstep+1);
if (nc==1);     tf(:,:)=f(:,1)*ones(1,nstep+1);  end
if (nc>1);      tf = f;        end

[nr ncip]=size(ip);
if (ncip >= 5);
nsnap=ip(5);             nhist=ip(6);
lists=ip(7:6+nsnap);     listh=ip(7+nsnap:ncip);
if (nhist > 0);
[nr nc]=size(listh);
D=zeros(nc,nstep+1);   V=zeros(nc,nstep+1);   A=zeros(nc,nstep+1);
end
if (nsnap > 0);          Dsnap=zeros(nd,nsnap); end
end
a0=M\(tf(:,1)-C*v0-K*d0);
if (nhist > 0);
D(:,1) = d0(listh);   V(:,1) = v0(listh);   A(:,1) = a0(listh);
end

tempd=zeros(nd,1);    tempv=zeros(nd,1);    tempa=zeros(nd,1);
fdof=[1:nd]';
if (bound==1);
[nr nc]=size(pdisp);
if (nc==2);
pd=pdisp(:,2)*ones(1,nstep+1);      pv=zeros(nr,nstep+1);
end
if (nc>2);
pd=pdisp(:,2:nstep+2);      pv(:,1)=(pd(:,2)-pd(:,1))/dt;
%size(pd), size(pdisp),size(pv),
pv(:,2:nstep+1)=(pd(:,2:nstep+1)-pd(:,1:nstep))/dt;
end
pdof=pdisp(:,1); fdof(pdof)=[];
Keff = M(fdof,fdof)+b3*C(fdof,fdof)+b4*K(fdof,fdof);
end
if (bound==0);  Keff = M+b3*C+b4*K;         end
[L,U]=lu(Keff);

dnew=d0(fdof);    vnew=v0(fdof);    anew=a0(fdof);
isnap=1;
for j = 1:nstep;
time=dt*j;
dold=dnew;      vold=vnew;      aold=anew;
dpred=dold+dt*vold+b1*aold;     vpred=vold+b2*aold;
if (bound==0); reff=tf(:,j+1)-C*vpred-K*dpred;  end
if (bound==1);
pdeff=C(fdof,pdof)*pv(:,j+1)+K(fdof,pdof)*pd(:,j+1);
reff=tf(fdof,j+1)-C(fdof,fdof)*vpred-K(fdof,fdof)*dpred-pdeff;
end
y=L\reff;         anew=U\y;  dnew=dpred+b4*anew;  vnew=vpred+b3*anew;
if (nhist > 0 | nsnap > 0);
if (bound==1);  tempd(pdof)=pd(:,j+1);  tempv(pdof)=pv(:,j+1); end
tempd(fdof)=dnew;         tempv(fdof)=vnew;         tempa(fdof)=anew;
if (nhist > 0);
D(:,j+1) = tempd(listh);  V(:,j+1) = tempv(listh);  A(:,j+1) = tempa(listh);
end
if (nsnap > 0  & isnap <= nsnap );
if (time >= lists(isnap)); Dsnap(:,isnap) = tempd; isnap=isnap+1; end
end
end
end
%--------------------------end--------------------------------


After plugging in all the corresponding values and running the "step2" function, a plot was made for each mass. The plots of the displacement degrees of freedom, ${\displaystyle d_{1}(t)}$ and ${\displaystyle d_{2}(t)}$, over the time interval are shown below:

${\displaystyle d_{1}(t)}$

${\displaystyle d_{2}(t)}$

## Problem 4

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

### Find

Using MATLAB find the displacement at the center and reactions ${\displaystyle R_{L}}$ and ${\displaystyle R_{R}}$

### Given

${\displaystyle E=100GPa,F=10,000N,}$ area of cross sections of the three portions shown are, starting from the left,${\displaystyle A_{1}=10^{-4}m^{2},A_{2}=2*10^{-4}m^{2},}$ and${\displaystyle A_{3}=10^{-4}m^{2}}$

### Solution

Since forces can only be applied at nodes, divide the middle bar into two different elements. There are 4 elements and 5 nodes.

MATLAB script file:

Variables are created for the constants f changes for each case:

E=100e9;

F=10000;
L1 = 0.3; %in meters
A1 = 10e-4; %in meters^2
L2 = 0.2;
A2 = 2*10e-4;
L3 = 0.3;
A3 = 10e-4;

%Element 1
k1=(E*A1/L1)*[1 -1;
-1  1];
%Element 2
k2=(E*A2/L2)*[1 -1;
-1  1];
%Element 3
k3=(E*A2/L2)*[1 -1;
-1  1];
%Element 4
k4=(E*A3/L3)*[1 -1;
-1  1];

Fsimp=[0 F 0]';
Ksimp=[k1(2,2)+k2(1,1) k2(1,2) 0;
k2(2,1) k2(2,2)+k3(1,1) k3(1,2);
0 k3(2,1) k3(2,2)+k4(1,1)];

innerdisplacements=Ksimp\Fsimp;
Center_displacement = innerdisplacements(2)%outer displacements are zero

R_L = k1(1,2)* innerdisplacements(1)
R_R = k4(2,1)* innerdisplacements(3)%u1 and u5 are zero so no reason to include

%Undisplaced Elements
Bar1x1=[0 .1 .2 .3];
Bar1y1=[.4 .4 .4 .4];
Bar1y2=[.3 .3 .3 .3];
Bar2x1=[.3 .35 .4 .45 .5];
Bar2y1=[.6 .6 .6 .6 .6];
Bar2y2=[.2 .2 .2 .2 .2];
Bar3x1=[.5 .55 .6 .65 .7];
Bar3y1=[.6 .6 .6 .6 .6];
Bar3y2=[.2 .2 .2 .2 .2];
Bar4x1=[.7 .8 .9 1];
Bar4y1=[.4 .4 .4 .4];
Bar4y2=[.3 .3 .3 .3];

%Displaced
dBar1x1=[0 .1 .2 .300015];
dBar2x1=[.300015 .35 .4 .45 .500035];
dBar3x1=[.500035 .55 .6 .65 .70005];
dBar4x1=[.70005 .8 .9 1.00005];

plot(Bar1x1,Bar1y1,'b',Bar2x1,Bar2y1,'b',Bar3x1,Bar3y1,'b',Bar4x1,Bar4y1,'b',...
Bar1x1,Bar1y2,'b',Bar2x1,Bar2y2,'b',Bar3x1,Bar3y2,'b',Bar4x1,Bar4y2,'b')
line([.3 ; .3],[.2 ; .6])
line([.7 ; .7],[.2 ; .6])
line([.300015 ; .300015],[.2 ; .6])
line([.70005 ; .70005],[.2 ; .6])

plot(dBar1x1,Bar1y1,'r--',dBar2x1,Bar2y1,'r--',dBar3x1,Bar3y1,'r--',dBar4x1,Bar4y1,'r--')


MATLAB output:

Center_displacement =

  2.0000e-05


R_L =

 -5.0000e+03


R_R =

 -5.0000e+03


The displacements were so insignificant that the displaced lines didn't appear to render:

## Problem 5

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

### Find

Determine displacements at Node 1 and axial forces in Element 1 and 2. Use Von Mises' theory to determine if the elements will yield or not. Use Euler buckling load (Pcr=pi^2*E*I/L^2) to determine if the elements under compressive loads will buckle.

### Given

Outer Dimension= 0.012 meters

Inner Dimension= 0.009 meters

E=70GPa

Force at Node 1= 1000N

### Solution

function FEA_R3p5
E=70000000000;
ao=0.012;
ai=0.009;
sig=70000000;
F=1000;
A=ao^2-ai^2;
phi1=26.6;
phi2=-90;
L1=sqrt(5);
L2=2;
I=(ao^4-ai^4)/12;
%Element 1
k1=(E*A/L1)*[cos(phi1)^2 cos(phi1)*sin(phi1) -cos(phi1)^2 -cos(phi1)*sin(phi1);
cos(phi1)*sin(phi1) sin(phi1)^2 -cos(phi1)*sin(phi1) -sin(phi1)^2;
-cos(phi1)^2 -cos(phi1)*sin(phi1) cos(phi1)^2 cos(phi1)*sin(phi1);
-cos(phi1)*sin(phi1) -sin(phi1)^2 cos(phi1)*sin(phi1) sin(phi1)^2];
Pcr1=pi^2*E*I/L1^2
%Element 2
k2=(E*A/L2)*[cos(phi2)^2 cos(phi2)*sin(phi2) -cos(phi2)^2 -cos(phi2)*sin(phi2);
cos(phi2)*sin(phi2) sin(phi2)^2 -cos(phi2)*sin(phi2) -sin(phi2)^2;
-cos(phi2)^2 -cos(phi2)*sin(phi2) cos(phi2)^2 cos(phi2)*sin(phi2);
-cos(phi2)*sin(phi2) -sin(phi2)^2 cos(phi2)*sin(phi2) sin(phi2)^2];
Pcr2=pi^2*E*I/L2^2
U=[0 F]';
R2=[k1(3,3)+k2(1,1) k1(3,4)+k2(1,2) ;
k1(4,3)+k2(2,1) k1(4,4)+k2(2,2) ];

Displacements_at_Node1_y_x_in_Meters=R2\U

f=[k1(1:2,3:4) ;k2(3:4,1:2)];
Forces_at_Node2_and_Node3_in_Newtons=f*Displacements_at_Node1_y_x_in_Meters

%Undisplaced Elements
Element1x1=[0 .5 1 1.5 2];
Element1y1=-.5*Element1x1+1;
Element2x1=[0 .5 1 1.5 2];
Element2y1=[0 0 0 0 0];
%Displaced Elements
Element1x2=[0 .5 1 1.5 2.0009];
Element1y2=[1 .75 .5 .25 -.002];
Element2x2=[0 .5 1 1.5 2.0009];
Element2y2=[0 0 0 -.0005 -.002];

% plot(Element1x1,Element1y1,'g',Element2x1,Element2y1,'r')
% axis([0,2.5,-.1,1.5])

plot(Element1x1,Element1y1,'g',Element2x1,Element2y1,'r',Element1x2,Element1y2, 'b--',Element2x2,Element2y2,'k--')
axis([1.85,2.01,-.01,0.1]);


EDU>> FEA_R3p5

Pcr1 =

 163.2186


Pcr2 =

 204.0232


Displacements_at_Node1_y_x_in_Meters =

  -0.0020
0.0009


Forces_at_Node2_and_Node3_in_Newtons =

 1.0e+003 *

  -0.1311
-1.2615
0.1311
0.2615


## Problem 6

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

### Find

Using MATLAB find the deflections and element forces and plot the displacement for the following load cases:

Load Case A: F13 = F27 = 10,000 N

Load Case A: F14 = F28 = 10,000 N

Load Case C: F13 = 10,000 N and F27 = -10,000 N

### Given

${\displaystyle E=100GPa,A=1.0cm^{2},L=0.3m}$

### Solution

#### MATLAB for 2-D Truss System

MATLAB script file:

Variables are created for the constants f changes for each case:


%GIVENS
L=0.3;
E=100e9;
A=1e-4;
ep=[E A];
f=zeros(28,1);
f(14)=10000;
f(28)=10000;


Build the connectivity array:


%ELEMENT DOF
Edof=zeros(25,5);
for i=1:6
j=2*(i-1);
Edof(i,:)=[i,1+j,2+j,3+j,4+j];
end
for i=7:12
j=2*(i+1);
Edof(i,:)=[i,j-1,j,j+1,j+2];
end
for i=13:19
j=2*(i-13);
Edof(i,:)=[i,1+j,2+j,15+j,16+j];
end
for i=20:25
j=2*(i-20);
Edof(i,:)=[i,1+j,2+j,17+j,18+j];
end


ex and ey are matrices of the coordinates of the nodes:


%COORDINATES
ex=zeros(14,2);
ey=zeros(14,2);
for i=1:25
if i<7
ex(i,:)=[L*(i-1),L*i];
ey(i,:)=[0,0];
elseif i>6&&i<13
ex(i,:)=[L*(i-7),L*(i-6)];
ey(i,:)=[L,L];
elseif i>12&&i<20
ex(i,:)=[L*(i-13),L*(i-13)];
ey(i,:)=[0,L];
else
ex(i,:)=[L*(i-20),L*(i-19)];
ey(i,:)=[0,L];
end
end


Solve for the displacements, plot, and find the forces


%STIFFNESS MATRIX
K=zeros(28);
for i=1:25
E=ep(1);A=ep(2);
xt=ex(i,2)-ex(i,1);
yt=ey(i,2)-ey(i,1);
L=sqrt(xt^2+yt^2);
l=xt/L; m=yt/L;
ke=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;];
edoft=Edof(i,2:5);
K(edoft,edoft)=K(edoft,edoft)+ke;
end


Build the global stiffness matrix:



%SOLVE
bc=[1 0; 15 0; 16 0];
elm=[1:28]';
elm(bc(:,1))=[];
d=zeros(1);
d(elm)=K(elm,elm)\f(elm);

%PLOT
x=ex' ;
y=ey';
xc=x ;yc=y;
s1=['-' , 'k'];
s2='ko';
hold on
axis equal
plot(xc,yc,s1)
plot(x,y,s2)
hold off
ed=zeros(25,4);
n=5;
rr=Edof(:,2:n);
for i=1:25;
ed(i,1:4)=d(rr(i,:));
end
x=(ex+ed(:,[1 3]))';
y=(ey+ed(:,[2 4]))';
xc=x ;yc=y;
s1=['-' , 'm'];
s2='k.';
hold on
axis equal
plot(xc,yc,s1)
plot(x,y,s2)
hold off

%ELEMENT FORCES
for i=1:25
E=ep(1);A=ep(2);
xt=ex(i,2)-ex(i,1);
yt=ey(i,2)-ey(i,1);
L=sqrt(xt^2+yt^2);
l=xt/L; m=yt/L;
T=[l m 0 0;
0 0 l m];
N(i)=E*A/L*[-1 1]*T*ed(i,:)';
end

table=zeros(25,6);
for i=1:25
table(i,1)=i;
table(i,2:5)=ed(i,:);
table(i,6)=N(i);
end
format shortG
table


MATLAB output:

           1            0  -1.3553e-18       0.0003      -0.0003        10000
2       0.0003      -0.0003       0.0006      -0.0006        10000
3       0.0006      -0.0006       0.0009      -0.0009        10000
4       0.0009      -0.0009       0.0012      -0.0012        10000
5       0.0012      -0.0012       0.0015      -0.0015        10000
6       0.0015      -0.0015       0.0018      -0.0018        10000
7            0            0       0.0003      -0.0003        10000
8       0.0003      -0.0003       0.0006      -0.0006        10000
9       0.0006      -0.0006       0.0009      -0.0009        10000
10       0.0009      -0.0009       0.0012      -0.0012        10000
11       0.0012      -0.0012       0.0015      -0.0015        10000
12       0.0015      -0.0015       0.0018      -0.0018        10000
13            0  -1.3553e-18            0            0   4.5178e-11
14       0.0003      -0.0003       0.0003      -0.0003   4.7294e-11
15       0.0006      -0.0006       0.0006      -0.0006   4.0018e-11
16       0.0009      -0.0009       0.0009      -0.0009   2.9104e-11
17       0.0012      -0.0012       0.0012      -0.0012   1.4552e-11
18       0.0015      -0.0015       0.0015      -0.0015            0
19       0.0018      -0.0018       0.0018      -0.0018   -7.276e-12
20            0  -1.3553e-18       0.0003      -0.0003  -6.2755e-11
21       0.0003      -0.0003       0.0006      -0.0006  -6.8212e-11
22       0.0006      -0.0006       0.0009      -0.0009  -4.9113e-11
23       0.0009      -0.0009       0.0012      -0.0012  -2.0009e-11
24       0.0012      -0.0012       0.0015      -0.0015  -1.0914e-11
25       0.0015      -0.0015       0.0018      -0.0018   1.0914e-11


           1            0       0.0006        0.003    0.0064971        1e+05
2        0.003    0.0064971       0.0054     0.018394        80000
3       0.0054     0.018394       0.0072     0.035091        60000
4       0.0072     0.035091       0.0084     0.055388        40000
5       0.0084     0.055388        0.009     0.078085        20000
6        0.009     0.078085        0.009      0.10168            0
7            0            0      -0.0036    0.0058971     -1.2e+05
8      -0.0036    0.0058971      -0.0066     0.017794       -1e+05
9      -0.0066     0.017794       -0.009     0.034491       -80000
10       -0.009     0.034491      -0.0108     0.054788       -60000
11      -0.0108     0.054788       -0.012     0.077485       -40000
12       -0.012     0.077485      -0.0126      0.10138       -20000
13            0       0.0006            0            0       -20000
14        0.003    0.0064971      -0.0036    0.0058971       -20000
15       0.0054     0.018394      -0.0066     0.017794       -20000
16       0.0072     0.035091       -0.009     0.034491       -20000
17       0.0084     0.055388      -0.0108     0.054788       -20000
18        0.009     0.078085       -0.012     0.077485       -20000
19        0.009      0.10168      -0.0126      0.10138       -10000
20            0       0.0006      -0.0036    0.0058971        28284
21        0.003    0.0064971      -0.0066     0.017794        28284
22       0.0054     0.018394       -0.009     0.034491        28284
23       0.0072     0.035091      -0.0108     0.054788        28284
24       0.0084     0.055388       -0.012     0.077485        28284
25        0.009     0.078085      -0.0126      0.10138        28284


           1            0   4.3966e-18       0.0003       0.0003        10000
2       0.0003       0.0003       0.0006       0.0012        10000
3       0.0006       0.0012       0.0009       0.0027        10000
4       0.0009       0.0027       0.0012       0.0048        10000
5       0.0012       0.0048       0.0015       0.0075        10000
6       0.0015       0.0075       0.0018       0.0108        10000
7            0            0      -0.0003       0.0003       -10000
8      -0.0003       0.0003      -0.0006       0.0012       -10000
9      -0.0006       0.0012      -0.0009       0.0027       -10000
10      -0.0009       0.0027      -0.0012       0.0048       -10000
11      -0.0012       0.0048      -0.0015       0.0075       -10000
12      -0.0015       0.0075      -0.0018       0.0108       -10000
13            0   4.3966e-18            0            0  -1.4655e-10
14       0.0003       0.0003      -0.0003       0.0003  -1.4916e-10
15       0.0006       0.0012      -0.0006       0.0012  -1.3824e-10
16       0.0009       0.0027      -0.0009       0.0027  -1.0186e-10
17       0.0012       0.0048      -0.0012       0.0048  -2.9104e-11
18       0.0015       0.0075      -0.0015       0.0075   5.8208e-11
19       0.0018       0.0108      -0.0018       0.0108            0
20            0   4.3966e-18      -0.0003       0.0003   2.0736e-10
21       0.0003       0.0003      -0.0006       0.0012   2.1464e-10
22       0.0006       0.0012      -0.0009       0.0027   1.7826e-10
23       0.0009       0.0027      -0.0012       0.0048    7.276e-11
24       0.0012       0.0048      -0.0015       0.0075   9.4587e-11
25       0.0015       0.0075      -0.0018       0.0108  -8.7311e-11


#### CALFEM for 2-D Truss System

MATLAB script file:


L=0.3;
E=100e9;
A=1e-4;
ep=[E A];
f=zeros(28,1);
f(13)=10000;
f(27)=10000;

Edof=zeros(25,5);
for i=1:6
j=2*(i-1);
Edof(i,:)=[i,1+j,2+j,3+j,4+j];
end
for i=7:12
j=2*(i+1);
Edof(i,:)=[i,j-1,j,j+1,j+2];
end
for i=13:19
j=2*(i-13);
Edof(i,:)=[i,1+j,2+j,15+j,16+j];
end
for i=20:25
j=2*(i-20);
Edof(i,:)=[i,1+j,2+j,17+j,18+j];
end

ex=zeros(14,2);
ey=zeros(14,2);
for i=1:25
if i<7
ex(i,:)=[L*(i-1),L*i];
ey(i,:)=[0,0];
elseif i>6&&i<13
ex(i,:)=[L*(i-7),L*(i-6)];
ey(i,:)=[L,L];
elseif i>12&&i<20
ex(i,:)=[L*(i-13),L*(i-13)];
ey(i,:)=[0,L];
else
ex(i,:)=[L*(i-20),L*(i-19)];
ey(i,:)=[0,L];
end
end

K=zeros(28);
for i=1:25
Ke=bar2e(ex(i,:),ey(i,:),ep);
K=assem(Edof(i,:),K,Ke);
end

bc=[1 0; 15 0; 16 0];
a=solveq(K,f,bc);
Ed=extract(Edof,a);
eldraw2(ex,ey)
plotpar=[1 3 0];scale=10;
eldisp2(ex,ey,Ed,plotpar,scale);


MATLAB output: For Load Case A:

           1            0  -1.3553e-18       0.0003      -0.0003        10000
2       0.0003      -0.0003       0.0006      -0.0006        10000
3       0.0006      -0.0006       0.0009      -0.0009        10000
4       0.0009      -0.0009       0.0012      -0.0012        10000
5       0.0012      -0.0012       0.0015      -0.0015        10000
6       0.0015      -0.0015       0.0018      -0.0018        10000
7            0            0       0.0003      -0.0003        10000
8       0.0003      -0.0003       0.0006      -0.0006        10000
9       0.0006      -0.0006       0.0009      -0.0009        10000
10       0.0009      -0.0009       0.0012      -0.0012        10000
11       0.0012      -0.0012       0.0015      -0.0015        10000
12       0.0015      -0.0015       0.0018      -0.0018        10000
13            0  -1.3553e-18            0            0   4.5178e-11
14       0.0003      -0.0003       0.0003      -0.0003   4.7294e-11
15       0.0006      -0.0006       0.0006      -0.0006   4.0018e-11
16       0.0009      -0.0009       0.0009      -0.0009   2.9104e-11
17       0.0012      -0.0012       0.0012      -0.0012   1.4552e-11
18       0.0015      -0.0015       0.0015      -0.0015            0
19       0.0018      -0.0018       0.0018      -0.0018   -7.276e-12
20            0  -1.3553e-18       0.0003      -0.0003  -6.2755e-11
21       0.0003      -0.0003       0.0006      -0.0006  -6.8212e-11
22       0.0006      -0.0006       0.0009      -0.0009  -4.9113e-11
23       0.0009      -0.0009       0.0012      -0.0012  -2.0009e-11
24       0.0012      -0.0012       0.0015      -0.0015  -1.0914e-11
25       0.0015      -0.0015       0.0018      -0.0018   1.0914e-11


           1            0       0.0006        0.003    0.0064971        1e+05
2        0.003    0.0064971       0.0054     0.018394        80000
3       0.0054     0.018394       0.0072     0.035091        60000
4       0.0072     0.035091       0.0084     0.055388        40000
5       0.0084     0.055388        0.009     0.078085        20000
6        0.009     0.078085        0.009      0.10168            0
7            0            0      -0.0036    0.0058971     -1.2e+05
8      -0.0036    0.0058971      -0.0066     0.017794       -1e+05
9      -0.0066     0.017794       -0.009     0.034491       -80000
10       -0.009     0.034491      -0.0108     0.054788       -60000
11      -0.0108     0.054788       -0.012     0.077485       -40000
12       -0.012     0.077485      -0.0126      0.10138       -20000
13            0       0.0006            0            0       -20000
14        0.003    0.0064971      -0.0036    0.0058971       -20000
15       0.0054     0.018394      -0.0066     0.017794       -20000
16       0.0072     0.035091       -0.009     0.034491       -20000
17       0.0084     0.055388      -0.0108     0.054788       -20000
18        0.009     0.078085       -0.012     0.077485       -20000
19        0.009      0.10168      -0.0126      0.10138       -10000
20            0       0.0006      -0.0036    0.0058971        28284
21        0.003    0.0064971      -0.0066     0.017794        28284
22       0.0054     0.018394       -0.009     0.034491        28284
23       0.0072     0.035091      -0.0108     0.054788        28284
24       0.0084     0.055388       -0.012     0.077485        28284
25        0.009     0.078085      -0.0126      0.10138        28284


           1            0   4.3966e-18       0.0003       0.0003        10000
2       0.0003       0.0003       0.0006       0.0012        10000
3       0.0006       0.0012       0.0009       0.0027        10000
4       0.0009       0.0027       0.0012       0.0048        10000
5       0.0012       0.0048       0.0015       0.0075        10000
6       0.0015       0.0075       0.0018       0.0108        10000
7            0            0      -0.0003       0.0003       -10000
8      -0.0003       0.0003      -0.0006       0.0012       -10000
9      -0.0006       0.0012      -0.0009       0.0027       -10000
10      -0.0009       0.0027      -0.0012       0.0048       -10000
11      -0.0012       0.0048      -0.0015       0.0075       -10000
12      -0.0015       0.0075      -0.0018       0.0108       -10000
13            0   4.3966e-18            0            0  -1.4655e-10
14       0.0003       0.0003      -0.0003       0.0003  -1.4916e-10
15       0.0006       0.0012      -0.0006       0.0012  -1.3824e-10
16       0.0009       0.0027      -0.0009       0.0027  -1.0186e-10
17       0.0012       0.0048      -0.0012       0.0048  -2.9104e-11
18       0.0015       0.0075      -0.0015       0.0075   5.8208e-11
19       0.0018       0.0108      -0.0018       0.0108            0
20            0   4.3966e-18      -0.0003       0.0003   2.0736e-10
21       0.0003       0.0003      -0.0006       0.0012   2.1464e-10
22       0.0006       0.0012      -0.0009       0.0027   1.7826e-10
23       0.0009       0.0027      -0.0012       0.0048    7.276e-11
24       0.0012       0.0048      -0.0015       0.0075   9.4587e-11
25       0.0015       0.0075      -0.0018       0.0108  -8.7311e-11


## Problem 7

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

### Find

 Verify the element stiffness transformation from local coordinates to global coordinates by: ${\displaystyle \mathbf {k^{(e)}} =\mathbf {T} ^{(e)T}\,\mathbf {\hat {k}} ^{(e)}\mathbf {T} ^{(e)}}$ (7.1)
 Verify the element stiffness transformation from local coordinates to global coordinates by: ${\displaystyle \mathbf {k} ^{(e)}={\widetilde {\mathbf {T} }}^{{(e)}T}\,{\widetilde {\mathbf {k} }}^{(e)}\,{\widetilde {\mathbf {T} }}^{(e)}}$ (7.2)

### Given

The local stiffness matrix:

 ${\displaystyle \mathbf {{\hat {k}}^{(e)}} =k^{(e)}{\begin{bmatrix}1&-1\\-1&1\\\end{bmatrix}}}$ (7.3)
 ${\displaystyle \mathbf {T^{(e)}} ={\begin{bmatrix}l^{(e)}&m^{(e)}&0&0\\0&0&l^{(e)}&m^{(e)}\end{bmatrix}}}$ (7.4)
 ${\displaystyle \mathbf {T^{(e)^{T}}} ={\begin{bmatrix}l^{(e)}&0\\m^{(e)}&0\\0&l^{(e)}\\0&m^{(e)}\end{bmatrix}}}$ (7.5)
 ${\displaystyle {\widetilde {\mathbf {k} }}^{(e)}={\begin{bmatrix}1&0&-1&0\\0&0&0&0\\-1&0&1&0\\0&0&0&0\end{bmatrix}}}$ (7.6)
 ${\displaystyle {\widetilde {\mathbf {T} }}^{(e)}={\begin{bmatrix}l^{(e)}&m^{(e)}&0&0\\-m^{(e)}&l^{(e)}&0&0\\0&0&l^{(e)}&m^{(e)}\\0&0&-m^{(e)}&l^{(e)}\end{bmatrix}}}$ (7.7)
 ${\displaystyle {\widetilde {\mathbf {T} }}^{{(e)}T}={\begin{bmatrix}l^{(e)}&-m^{(e)}&0&0\\m^{(e)}&l^{(e)}&0&0\\0&0&l^{(e)}&-m^{(e)}\\0&0&m^{(e)}&l^{(e)}\end{bmatrix}}}$ (7.8)

### Solution

The element stiffness matrix ${\displaystyle \mathbf {{\hat {k}}^{(e)}} }$ is expressed in local coordinates. To transform the matrix to global coordinates: ${\displaystyle \mathbf {k^{(e)}} =\mathbf {T^{(e)-1}} \mathbf {{\hat {k}}^{(e)}} \mathbf {T^{(e)}} }$, where ${\displaystyle \mathbf {T^{(e)}} }$ is the transformation matrix.

It is shown that ${\displaystyle \mathbf {T^{(e)-1}} =\mathbf {T^{(e)T}} }$ because the interior vectors are orthogonal.

 ${\displaystyle \mathbf {T^{(e)}} ={\begin{bmatrix}l^{(e)}&m^{(e)}&0&0\\0&0&l^{(e)}&m^{(e)}\end{bmatrix}}}$ (7.9)

To prove orthogonality, the dot product of two vectors must equal zero. Using the dot product definition: The dot product of two vectors a = [a1, a2, ..., an] and b = [b1, b2, ..., bn] is defined as:[4]

 ${\displaystyle \mathbf {a} \cdot \mathbf {b} =\sum _{i=1}^{n}a_{i}b_{i}=a_{1}b_{1}+a_{2}b_{2}+\cdots +a_{n}b_{n}}$ (7.10)

The ${\displaystyle \mathbf {T^{(e)}} }$ matrix in illustrative form:

 ${\displaystyle \mathbf {T^{(e)}} ={\begin{bmatrix}a_{1}&a_{2}&a_{3}&a_{4}\\b_{1}&b_{2}&b_{3}&b_{4}\end{bmatrix}}}$ (7.11)
 Therefore, ${\displaystyle \mathbf {a} \cdot \mathbf {b} =\sum _{i=1}^{2}a_{i}b_{i}=a_{1}b_{1}+a_{2}b_{2}=l^{(e)}*0+m^{(e)}*0+0*l^{(e)}+0*m^{(e)}=0}$ (7.11)

Now we can say,

 ${\displaystyle \mathbf {k} ^{(e)}=\mathbf {T} ^{(e)T}\mathbf {\hat {k}} ^{(e)}\mathbf {T} ^{(e)}}$ (7.12)

${\displaystyle \mathbf {T} ^{(e)T}\mathbf {\hat {k}} ^{(e)}={\begin{bmatrix}l^{(e)}&0\\m^{(e)}&0\\0&l^{(e)}\\0&m^{(e)}\end{bmatrix}}k^{(e)}{\begin{bmatrix}1&-1\\-1&1\\\end{bmatrix}}}$
${\displaystyle \mathbf {T} ^{(e)T}\mathbf {\hat {k}} ^{(e)}=k^{(e)}{\begin{bmatrix}l^{(e)}*1+0*-1&l^{(e)}*-1+0*1\\m^{(e)}*1+0*-1&m^{(e)}*-1+0*1\\0*1+l^{(e)}*-1&0*-1+l^{(e)}*1\\0*1+m^{(e)}*-1&0*-1+m^{(e)}*1\end{bmatrix}}=k^{(e)}{\begin{bmatrix}l^{(e)}&-l^{(e)}\\m^{(e)}&-m^{(e)}\\-l^{(e)}&-l^{(e)}\\-m^{(e)}&m^{(e)}\end{bmatrix}}}$
${\displaystyle \mathbf {k} ^{(e)}=[\mathbf {T} ^{(e)T}\mathbf {\hat {k}} ^{(e)}]\mathbf {T} ^{(e)}=k^{(e)}{\begin{bmatrix}l^{(e)}&-l^{(e)}\\m^{(e)}&-m^{(e)}\\-l^{(e)}&-l^{(e)}\\-m^{(e)}&m^{(e)}\end{bmatrix}}\mathbf {T} ^{(e)}=k^{(e)}{\begin{bmatrix}l^{(e)}&-l^{(e)}\\m^{(e)}&-m^{(e)}\\-l^{(e)}&-l^{(e)}\\-m^{(e)}&m^{(e)}\end{bmatrix}}{\begin{bmatrix}l^{(e)}&m^{(e)}&0&0\\0&0&l^{(e)}&m^{(e)}\end{bmatrix}}}$

The element stiffness matrix in global coordinates is:

 ${\displaystyle \mathbf {k^{(e)}} =k^{(e)}{\begin{bmatrix}(l^{(e)})^{2}&l^{(e)}m^{(e)}&-(l^{(e)})^{2}&-l^{(e)}m^{(e)}\\l^{(e)}m^{(e)}&(m^{(e)})^{2}&-l^{(e)}m^{(e)}&-(m^{(e)})^{2}\\-(l^{(e)})^{2}&-l^{(e)}m^{(e)}&(l^{(e)})^{2}&l^{(e)}m^{(e)}\\-l^{(e)}m^{(e)}&-(m^{(e)})^{2}&l^{(e)}m^{(e)}&(m^{(e)})^{2}\end{bmatrix}}}$ (7.13)

Proving

 ${\displaystyle \mathbf {k} ^{(e)}={\widetilde {\mathbf {T} }}^{{(e)}T}\,{\widetilde {\mathbf {k} }}^{(e)}\,{\widetilde {\mathbf {T} }}^{(e)}}$ (7.14)

${\displaystyle {\widetilde {\mathbf {T} }}^{{(e)}T}{\widetilde {\mathbf {k} }}^{(e)}={\begin{bmatrix}l^{(e)}&-m^{(e)}&0&0\\m^{(e)}&l^{(e)}&0&0\\0&0&l^{(e)}&-m^{(e)}\\0&0&m^{(e)}&l^{(e)}\end{bmatrix}}k^{(e)}{\begin{bmatrix}1&0&-1&0\\0&0&0&0\\-1&0&1&0\\0&0&0&0\end{bmatrix}}}$
${\displaystyle {\widetilde {\mathbf {T} }}^{{(e)}T}{\widetilde {\mathbf {k} }}^{(e)}=k^{(e)}{\begin{bmatrix}l^{(e)}*1+-m^{(e)}*0+0*-1+0*0&l^{(e)}*0+-m^{(e)}*0+0*0+0*0&l^{(e)}*-1+-m^{(e)}*0+0*1+0*0&l^{(e)}*0+-m^{(e)}*0+0*0+0*0\\m^{(e)}*1+l^{(e)}*0+0*-1+0*0&m^{(e)}*0+l^{(e)}*0+0*0+0*0&m^{(e)}*-1+l^{(e)}*0+0*1+0*0&m^{(e)}*0+l^{(e)}*0+0*0+0*0\\0*1+0*0+l^{(e)}*-1+-m^{(e)}*0&0*0+0*0+l^{(e)}*0+-m^{(e)}*0&0*-1+0*0+l^{(e)}*1+-m^{(e)}*0&0*0+0*0+l^{(e)}*0+-m^{(e)}*0\\0*1+0*0+m^{(e)}*-1+l^{(e)}*0&0*0+0*0+m^{(e)}*0+l^{(e)}*0&0*-1+0*0+m^{(e)}*1+l^{(e)}*0&0*0+0*0+m^{(e)}*0+l^{(e)}*0\\\end{bmatrix}}}$
${\displaystyle k^{(e)}{\begin{bmatrix}l^{(e)}&0&-l^{(e)}&0\\m^{(e)}&0&-m^{(e)}&0\\l^{(e)}&0&l^{(e)}&0\\-m^{(e)}&0&m^{(e)}&0\end{bmatrix}}{\widetilde {\mathbf {T} }}^{(e)}=k^{(e)}{\begin{bmatrix}l^{(e)}&0&-l^{(e)}&0\\m^{(e)}&0&-m^{(e)}&0\\l^{(e)}&0&l^{(e)}&0\\-m^{(e)}&0&m^{(e)}&0\end{bmatrix}}{\begin{bmatrix}l^{(e)}&m^{(e)}&0&0\\-m^{(e)}&l^{(e)}&0&0\\0&0&l^{(e)}&m^{(e)}\\0&0&-m^{(e)}&l^{(e)}\end{bmatrix}}}$
${\displaystyle {\widetilde {\mathbf {T} }}^{{(e)}T}\,{\widetilde {\mathbf {k} }}^{(e)}\,{\widetilde {\mathbf {T} }}^{(e)}=k^{(e)}{\begin{bmatrix}(l^{(e)})^{2}&l^{(e)}m^{(e)}&-(l^{(e)})^{2}&-l^{(e)}m^{(e)}\\l^{(e)}m^{(e)}&(m^{(e)})^{2}&-l^{(e)}m^{(e)}&-(m^{(e)})^{2}\\-(l^{(e)})^{2}&-l^{(e)}m^{(e)}&(l^{(e)})^{2}&l^{(e)}m^{(e)}\\-l^{(e)}m^{(e)}&-(m^{(e)})^{2}&l^{(e)}m^{(e)}&(m^{(e)})^{2}\end{bmatrix}}}$

Once again, it is shown that:
${\displaystyle \mathbf {k^{(e)}} =k^{(e)}{\begin{bmatrix}(l^{(e)})^{2}&l^{(e)}m^{(e)}&-(l^{(e)})^{2}&-l^{(e)}m^{(e)}\\l^{(e)}m^{(e)}&(m^{(e)})^{2}&-l^{(e)}m^{(e)}&-(m^{(e)})^{2}\\-(l^{(e)})^{2}&-l^{(e)}m^{(e)}&(l^{(e)})^{2}&l^{(e)}m^{(e)}\\-l^{(e)}m^{(e)}&-(m^{(e)})^{2}&l^{(e)}m^{(e)}&(m^{(e)})^{2}\end{bmatrix}}}$

## Problem 8

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

### Find

The system natural circular frequency ${\displaystyle \omega }$ for the given problem, expression for the particular solution ${\displaystyle y_{p}(t)}$, the complete solution ${\displaystyle y(t)}$, the amplification factor A and then plot the homogeneous, particular and complete solution in separate figures and then on the same figure for comparison.

### Given

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

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

Excitation:

 ${\displaystyle f(t)=f_{0}\cos {\bar {\omega }}t}$ (8.2)
 ${\displaystyle f_{0}/m=r_{0}=1}$ (8.3)
 ${\displaystyle {\bar {\omega }}=0.9\omega }$ (8.4)

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

### Solution

Case: Underdamped

 ${\displaystyle \Delta =a^{2}-4b<0}$ (8.5)

${\displaystyle \alpha =-0.5,\beta =2}$

 ${\displaystyle (\lambda -(-0.5+2i))(\lambda -(-0.5-2i))=0}$ (8.6)
 ${\displaystyle \lambda ^{2}+1\lambda +4.25=0}$ (8.7)

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

Since,

 ${\displaystyle \omega ^{2}=b}$ (8.8)

${\displaystyle \omega =2.06}$

 ${\displaystyle \zeta ={\frac {a}{2\omega }}}$ (8.9)

${\displaystyle \zeta =0.243}$

Therefore L2-ODE-CC:

 ${\displaystyle y_{h}''+y_{h}'+2.25y_{h}=0}$ (8.10)
 ${\displaystyle y_{p}''+y_{p}'+4.25y_{p}=f_{0}\cos {\bar {\omega }}t}$ (8.11)

The particular solution is as follows:

 ${\displaystyle y_{p}(t)=de^{i({\bar {\omega }}t-\phi )}}$ (8.12)

Where,

 ${\displaystyle d={\frac {f_{0}/m}{\sqrt {(\omega ^{2}-{\bar {\omega }}^{2})^{2}+(2\zeta \omega {\bar {\omega }})^{2}}}}}$ (8.13)
 ${\displaystyle \tan \phi ={\frac {2\zeta \omega {\bar {\omega }}}{\omega ^{2}-{\bar {\omega }}^{2}}}}$ (8.14)

Therefore,

 ${\displaystyle y_{p}(t)={\frac {f_{0}/m}{\sqrt {(\omega ^{2}-{\bar {\omega }}^{2})^{2}+(2\zeta \omega {\bar {\omega }})^{2}}}}e^{i({\bar {\omega }}t-\arctan {\frac {2\zeta \omega {\bar {\omega }}}{\omega ^{2}-{\bar {\omega }}^{2}}})}}$ (8.15)

Calculating,

 ${\displaystyle y_{p}(t)=0.494e^{i(1.854t-1.161)}}$ (8.16)

Taking only the real part, ${\displaystyle e^{it}=\cos(t)+i\sin(t)}$

${\displaystyle y_{p}(t)=0.494\cos(1.854t-1.161)}$


To calculate the amplification factor A the following equations are used:

 ${\displaystyle \rho ={\frac {\bar {\omega }}{\omega }}}$ (8.17)
 ${\displaystyle A={\frac {1}{\left[(1-\rho ^{2})^{2}+(2\zeta \rho )^{2}\right]^{1/2}}}}$ (8.18)

${\displaystyle A=2.1}$
The homogeneous solution is as follows:

 ${\displaystyle y_{h}(t)=C_{1}e^{(-0.5t)}\cos(2t)+C_{2}e^{(-0.5t)}\sin(2t)}$ (8.19)

The complete solution:

 ${\displaystyle y(t)=y_{h}(t)+y_{p}(t)}$ (8.20)
 ${\displaystyle y(t)=C_{1}e^{(-0.5x)}\cos(2x)+C_{2}e^{(-0.5x)}\sin(2x)+0.494\cos(1.854t-1.161)}$ (8.21)
 ${\displaystyle y'(t)=-0.5C_{1}e^{(-0.5t)}\cos(2t)-2C_{1}e^{(-0.5t)}\sin(2t)-0.5C_{2}e^{(-0.5t)}\sin(2t)+2C_{2}e^{(-0.5t)}\cos(2t)-0.916\sin(1.854t-1.161)}$ (8.21)

Using the two initial conditions and solving equations the coefficients are as follows:
${\displaystyle C_{1}=0.803,C_{2}=-0.22}$

Therefore,

${\displaystyle y(t)=0.803e^{(-0.5t)}\cos(2t)-0.22e^{(-0.5t)}\sin(2t)+0.494\cos(1.854t-1.161)}$


## References

2. Kim,Sankar. Introduction to Finite Element Analysis and Design. 2008 Obtained from p.100-Exercise 10
3. S. Lipschutz, M. Lipson (2009). Linear Algebra (Schaum’s Outlines) (4th ed.). McGraw Hill. ISBN 978-0-07-154352-1.

## Contributing Team Members

Problem Assignments
Problem # Solved & Typed by Reviewed by
1 Kristin Howe All
2 Kevin Li All
3 Spencer Herran All
4 Joshua Plicque All
5 Brandon Wright All
6 Matthew Gidel All
7 Joshua Plicque All
8 Zeyn Kermani All

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