# 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

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

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.

Assignment Help Features
Assignment Help Services
Calculator