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

d₃= 2 cm, d₄= -1 cm,d₅= -3 cm, d₆= 5 cm

Figure 1.1: KS 2008 p.81 Example 2.5

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. ($d_{3}=0$ and $d_{4}=0$)

Figure 2.1: Spring-Mass-Damper System

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

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

(2.1)

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

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

(2.2)

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

(2.3)

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

(2.4)

Solution

Free Body Diagrams

FBD of Mass 1:

Figure 2.2: FBD of Mass 1

FBD of Mass 2:

Figure 2.3: 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:

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

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

$\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 $\mathbf {M}$ is the mass matrix, $\mathbf {C}$ is the damping matrix, and $\mathbf {K}$ is the stiffness matrix and $d(t)$ and $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:
$\mathbf {d} ={\begin{Bmatrix}d_{1}\\d_{2}\end{Bmatrix}}$
and

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

$\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${\dot {d}}_{1}$ and ${\dot {d}}_{2}$. So we now have:

$\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 $d_{1}$ and $d_{2}$. Which is:

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

$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

$m_{1}=3,m_{2}=2$
$c_{1}=1/2,c_{2}=1/4,c_{3}=1/3$
$k_{1}=10,k_{2}=20,k_{3}=15$
$F_{1}(t)=0,F_{2}(t)=0$
$d_{1}(0)=-1,d_{2}(0)=2$
$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, $\lambda$, and eigenvectors, $x$, of the generalized eigenvalue problem:

$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
%-------------------------------------------------------------
 
% LAST MODIFIED: H Carlsson  1993-09-21
% 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, $\omega _{1}$ , using the formula:

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

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

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

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

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

$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
%-------------------------------------------------------------
 
% LAST MODIFIED: G Sandberg  1994-01-29
% 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, $d_{1}(t)$ and $d_{2}(t)$, over the time interval are shown below:

$d_{1}(t)$

$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 $R_{L}$ and $R_{R}$

Given

Figure for Exercise 10 on pg. 100

$E=100GPa,F=10,000N,$ area of cross sections of the three portions shown are, starting from the left,$A_{1}=10^{-4}m^{2},A_{2}=2*10^{-4}m^{2},$ and$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.

Annotated Nodes and Elements

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:

Graph of Bar System with Displacements

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

Members are made of Aluminum

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

Plots the Elements without Displacement

Shows Elements with Dsiplacement

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: F₁₃ = F₂₇ = 10,000 N

Load Case A: F₁₄ = F₂₈ = 10,000 N

Load Case C: F₁₃ = 10,000 N and F₂₇ = -10,000 N

Given

Figure 6.1: 25 member truss with L=0.3 m

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

Figure 6.2: Deflection under load case A using MATLAB

Load Case B:

           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

Figure 6.3: Deflection under load case B using MATLAB

Load Case C:

           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

Figure 6.4: Deflection under load case C using MATLAB

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

Figure 6.5: Deflection under load case A using CALFEM

For Load Case B:

           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

Figure 6.6: Deflection under load case B using CALFEM

For Load Case C:

           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

Figure 6.7: Deflection under load case C using CALFEM

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: $\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: $\mathbf {k} ^{(e)}={\widetilde {\mathbf {T} }}^{{(e)}T}\,{\widetilde {\mathbf {k} }}^{(e)}\,{\widetilde {\mathbf {T} }}^{(e)}$

(7.2)

Given

The local stiffness matrix:

$\mathbf {{\hat {k}}^{(e)}} =k^{(e)}{\begin{bmatrix}1&-1\\-1&1\\\end{bmatrix}}$

(7.3)

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

(7.4)

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

(7.5)

${\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)

${\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)

${\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 $\mathbf {{\hat {k}}^{(e)}}$ is expressed in local coordinates. To transform the matrix to global coordinates: $\mathbf {k^{(e)}} =\mathbf {T^{(e)-1}} \mathbf {{\hat {k}}^{(e)}} \mathbf {T^{(e)}}$, where $\mathbf {T^{(e)}}$ is the transformation matrix.

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

$\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 = [a₁, a₂, ..., a_n] and b = [b₁, b₂, ..., b_n] is defined as:

$\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 $\mathbf {T^{(e)}}$ matrix in illustrative form:

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

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

(7.12)

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

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

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

(7.14)

${\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}}$
${\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}}$
$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}}$
${\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:
$\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 $\omega$ for the given problem, expression for the particular solution $y_{p}(t)$, the complete solution $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:

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

(8.1)

Excitation:

$f(t)=f_{0}\cos {\bar {\omega }}t$

(8.2)

$f_{0}/m=r_{0}=1$

(8.3)

${\bar {\omega }}=0.9\omega$

(8.4)

Initial Conditions:
$y(0)=1,y'(0)=0$

Solution

Case: Underdamped

$\Delta =a^{2}-4b<0$

(8.5)

$\alpha =-0.5,\beta =2$

$(\lambda -(-0.5+2i))(\lambda -(-0.5-2i))=0$

(8.6)

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

(8.7)

$a=1,b=4.25$

Since,

$\omega ^{2}=b$

(8.8)

$\omega =2.06$

$\zeta ={\frac {a}{2\omega }}$

(8.9)

$\zeta =0.243$

Therefore L2-ODE-CC:

$y_{h}''+y_{h}'+2.25y_{h}=0$

(8.10)

$y_{p}''+y_{p}'+4.25y_{p}=f_{0}\cos {\bar {\omega }}t$

(8.11)

The particular solution is as follows:

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

(8.12)

Where,

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

(8.13)

$\tan \phi ={\frac {2\zeta \omega {\bar {\omega }}}{\omega ^{2}-{\bar {\omega }}^{2}}}$

(8.14)

Therefore,

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

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

(8.16)

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

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

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

$\rho ={\frac {\bar {\omega }}{\omega }}$

(8.17)

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

(8.18)

$A=2.1$
The homogeneous solution is as follows:

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

(8.19)

The complete solution:

$y(t)=y_{h}(t)+y_{p}(t)$

(8.20)

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

$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:
$C_{1}=0.803,C_{2}=-0.22$

Therefore,

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

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.