The global stiffness matrix Assignment
Question 1
- Part a: The global stiffness matrix
- Part b: Displacements of Nodes 2 and 3
- Part c: The reaction force at nodes 1 and 4
- Part d: The force in the spring 2
clc
clear
Given
K_1=100;
K_2=200;
K_3=100;
P=500;
u_1=0;
u_4=u_1;
Part a: The global stiffness matrix
Number of nodes is 4 and number of elements is 3 Global stifffness matrix, [K] is given by: [K]=[K_1]+[K_2]+[K_3] where [K_1] the local stiffness matrix of element 1 [K_2] the local stiffness matrix of element 2 and [K_3] the local stiffness matrix of element 3
K_1m=[100 -100; -100 100]; %local stiffness matrix of element 1
K_2m=[200 -200; -200 200]; %local stiffness matrix of element 2
K_3m=[100 -100; -100 100]; %local stiffness matrix of element 3
%defining the global stiffness matrix
K=zeros(4, 4);
%updating the elements of the global stiffness matrix using the local
%stiffness matrices
K(1:2, 1:2)=K(1:2, 1:2)+K_1m;
K(2:3, 2:3)=K(2:3, 2:3)+K_2m;
K(3:4, 3:4)=K(3:4, 3:4)+K_3m;
K
K =
100.0000e+000 -100.0000e+000 0.0000e+000 0.0000e+000
-100.0000e+000 300.0000e+000 -200.0000e+000 0.0000e+000
0.0000e+000 -200.0000e+000 300.0000e+000 -100.0000e+000
0.0000e+000 0.0000e+000 -100.0000e+000 100.0000e+000
Part b: Displacements of Nodes 2 and 3
From equilibrium equation [K][d]=[F] where [d] is the displacement vector [F] the force vector
%We know that:
% d_m=[0; u_2; u_3; 0] since u_1=u_4=0
% F=[0; 0; 500; 0];
% eqns=K*d_m==F;
syms u_2 u_3
eqn1=K(2,2)*u_2+K(2,3)*u_3==500;
eqn2=K(3,2)*u_2+K(3,3)*u_3==0;
eqns=[eqn1,eqn2];
vars=[u_2 u_3];
[d2,d3]=solve(eqns,vars)%d2 is the displacement at node 2 and d3 is the displacement at node 3
d2 =
3
d3 =
2
Part c: The reaction force at nodes 1 and 4
%Formula for reaction
% R=[K][d]-F
u_1=0; u_2=2; u_3=3;u_4=0;
d=[u_1; u_2; u_3; u_4]; %displacement vector
F=[0; 0; 500; 0]; %force vector
R=K*d-F; %Finding the rection matrix
R_1=R(1,1) %reaction at node 1
R_2=R(4,1) %reaction at node 1
R_1 =
-200.0000e+000
R_2 =
-300.0000e+000
Part D: The force in the spring 2
F=k_2(u_3-u_2)
F_2=K_2*(d3-d2)
F_2 =
-200
Question 2
- Part a
- Part b: Determining u(x) using Rayleigh ritz method.
- Part c: Plots
clc
clear
Part a
E=30*10^6; %psi
A=2; %in^2
syms x a_0 a_1 a_2
%from the qiven equation:
% u(x)=a_0+a_1*x+a_2*x^2
% we know that at x=0, u(x)=0, which gives
% a_0=0;
a_0=0;
%we also know that at x=60, u(x)=0, which gives
% a_1=-6a_2
%a_1=-6*a_2;
% and u(x)=-6*a_2*x+a_2*x^2;
% u(x)=a_2*(x^2-6x)
u=a_2*(x^2-6*x);
% The potential energy (II) is given by say P_e
% P_e=0.5*int(((E*A)*(diff(u,x))^2),x,0,60)-int((T_1*u*A*dx), x, 0,30)-int((T_2*u*4*A*dx), x, 30,60)
T_1(x)=10*x;
T_2(x)=300*x/x; %T(x)=300 lb/in
P_e=0.5*(int(((E*A)*(diff(u,x))^2),x,0,60)-int((T_1(x)*u*A), x, 0,30)-int((T_2(x)*u*4*A), x, 30,60));
Potential_Energy=P_e
Potential_Energy =
7408800000000*a_2^2 - 67365000*a_2
Part b: Determining u(x) using Rayleigh ritz method.
a_2=solve(diff(P_e,a_2)==0,a_2) %solving for a_2
%using a_2 in u(x) expression
u_2=a_2*(x^2-6*x)
sigma_xx=E*diff(u_2,x)
a_2 =
499/109760000
u_2 =
(499*x^2)/109760000 - (1497*x)/54880000
sigma_xx =
(187125*x)/686 - 561375/686
Part c: Plots
x=1:0.5:60;
%plot of u versus x
figure (1)
u_3=499/109760000*(x.^2-6*x);
plot(x,u_3)
xlabel('x')
ylabel('u(x)')
title('plot of u versus x')
%plot of sigma_xx versus x
figure (2)
sig_xx=(187125*x)/686 - 561375/686;
plot(x,sig_xx)
xlabel('x')
ylabel('\sigma')
title('plot of \sigma versus x')
Question 3
- Using force displacement method
- Nodal elongations
- Element stress calculations
- Calculating the support reactions
clc
clear
Using force displacement method
%Elongation is given
% Elong=PL/AE
% Total elongation is given by:
% Elong=(PL/AE)_1+(PL/AE)_2+(PL/AE)_3
L1=0.15;A1=250*10^-6;P1=900*10^3;
L2=0.15;A2=250*10^-6;P2=600*10^3;
L3=0.2;A3=400*10^-6; P3=600*10^3;
L4=0.2;A4=400*10^-6; P4=0;
E=200*10^9; %N/m^2
Elong=(P1*L1)/(A1*E)+(P2*L2)/(A2*E)+(P3*L3)/(A3*E);
disp(Elong)
%It can be seen that the obtained elongation is larger than 3.5*10^-3 gap
% Hence there is a compression force at ends less than the deformation of
% the bars
%We know that total elongation is 3.5 mm
syms P5
Elong1=((P1-P5)*L1)/(A1*E)+((P2-P5)*L2)/(A2*E)+((P3-P5)*L3)/(A3*E)+((P4-P5)*L4)/(A4*E);
soln=solve(Elong1==3.5*10^-3);
%Reaction at unfixed end P5
P5=vpa(soln, 5)
6.0000e-003
P5 =
227277.0
Nodal elongations
u1=vpa(((P1-P5)*L1)/(A1*E),3)%node 1
u2=vpa(((P2-P5)*L2)/(A2*E),3) %node 2
u3=vpa(((P3-P5)*L3)/(A3*E),3) %node 3
u4=vpa(((P4-P5)*L4)/(A4*E),3) %node 4
u5=3.5*10^-3 %From boundary conditions
u1 =
0.00202
u2 =
0.00112
u3 =
9.32e-4
u4 =
-5.68e-4
u5 =
3.5000e-003
Element stress calculations
S1=vpa((P1-P5)/A1, 5) %Element 1
S2=vpa((P2-P5)/A2, 5) %Element 2
S3=vpa((P3-P5)/A3, 5) %Element 3
S4=vpa((P4-P5)/A4, 5) %Element 4
S1 =
2.6909e+9
S2 =
1.4909e+9
S3 =
9.3182e+8
S4 =
-5.6818e+8
Calculating the support reactions
P_end=P5
P_start=vpa((600+300)*10^3-P5, 5)
P_end =
227277.0
P_start =
672733.0
Question 4
- Obtaining the element stiffness matrices
- Obtaining the global stiffness matrix
- Obtaining the nodal displacements and Support reactions from the Force displacement equation
- Element stresses
clc;
clear;
There are 4 elements and 5 nodes
n=4; Nodes=5;
Obtaining the element stiffness matrices
%Element 1 and 2 have the same length and area. Since they are made of the
%same material, the stiffness matrices are equal:
E=200*10^9; %N/m^2
A1=250;A2=250;A3=400;A4=400; %area
L1=150; L2=150;L3=200; L4=200; %length
% element stiffness matrices
M=[1 -1; -1 1];
K1=((A1*E)/L1).*M %element 1 stiffness matrix
K2=((A2*E)/L2).*M %element 2 stiffness matrix
K3=((A3*E)/L3).*M %element 3 stiffness matrix
K4=((A4*E)/L4).*M %element 4 stiffness matrix
K1 =
333.3333e+009 -333.3333e+009
-333.3333e+009 333.3333e+009
K2 =
333.3333e+009 -333.3333e+009
-333.3333e+009 333.3333e+009
K3 =
400.0000e+009 -400.0000e+009
-400.0000e+009 400.0000e+009
K4 =
400.0000e+009 -400.0000e+009
-400.0000e+009 400.0000e+009
Obtaining the global stiffness matrix
K=zeros(5, 5);
%updating the elements of the global stiffness matrix using the local
%stiffness matrices
K(1:2, 1:2)=K(1:2, 1:2)+K1;
K(2:3, 2:3)=K(2:3, 2:3)+K2;
K(3:4, 3:4)=K(3:4, 3:4)+K3;
K(4:5, 4:5)=K(4:5, 4:5)+K4;
disp("The Global stiffness matrix is")
disp(K)
The Global stiffness matrix is
Columns 1 through 4
333.3333e+009 -333.3333e+009 0.0000e+000 0.0000e+000
-333.3333e+009 666.6667e+009 -333.3333e+009 0.0000e+000
0.0000e+000 -333.3333e+009 733.3333e+009 -400.0000e+009
0.0000e+000 0.0000e+000 -400.0000e+009 800.0000e+009
0.0000e+000 0.0000e+000 0.0000e+000 -400.0000e+009
Column 5
0.0000e+000
0.0000e+000
0.0000e+000
-400.0000e+009
400.0000e+009
Obtaining the nodal displacements and Support reactions from the Force displacement equation
Boundary conditions u_1=0; u_5=3.5; The known forces: F2=300 kN F3=0 kN F4=600 kN The support reactions are: F1 and F5 which are unknowns From the force displacement equation F=K*u
syms u2 u3 u4 F1 F5
u=[0 u2 u3 u4 3.5*10^-3]';
F=[F1 900*10^3 600*10^3 600*10^3 F5]';
soln=solve(F==K*u);
format short eng;
% Nodal displacemnets
u2=vpa(soln.u2, 3) %m
u3=vpa(soln.u3, 3) %m
u4=vpa(soln.u4, 3) %m
%support reactions
F1=vpa(soln.F1, 5) %N
F5=vpa(soln.F5, 5) %N
u2 =
9.58e-4
u3 =
0.00191
u4 =
0.00271
F1 =
-3.1925e+8
F5 =
3.1715e+8
Element stresses
syms s1 s2 s3 s4
F2=900*10^3; F3=600*10^3; F4=600*10^3;
FS=[F1 F2 F3 F4]';
Stress=[s1 s2 s3 s4]';
A_m=[1/A1 0 0 0; 0 1/A2 0 0; 0 0 1/A3 0; 0 0 0 1/A4];
sol=solve(Stress==A_m*FS);
s1=vpa(sol.s1, 5) %Pa
s2=vpa(sol.s2, 5) %Pa
s3=vpa(sol.s3, 5) %Pa
s4=vpa(sol.s4, 5) %Pa
s1 =
-1.277e+6
s2 =
3600.0
s3 =
1500.0
s4 =
1500.0
Works Cited
Allaire, P.E., et al. “Simplex Finite Element Analysis of Viscous Incompressible Flow with Penalty Function Formulation.” Finite Elements in Analysis and Design, vol. 1, no. 1, Apr. 2015, pp. 71–88, 10.1016/0168-874x(85)90009-5. Accessed 25 July 2020.
Felton, Lewis P, and Richard B Nelson. Matrix Structural Analysis. New York, J. Wiley, 2018.
Hughes, Thomas J R. The Finite Element Method : Linear Static and Dynamic Finite Element Analysis. New York, Dover Publication, Inc, 2017.
