Experiment 2

Objective
Finite Element Formulation and analysis of one dimensional problem by developing computer program.

Software used
Matlab 2018

One dimensional Bar Problem

A stepped bar is subjected to an axial load of P = 2000 N as shown in figure. If the material of the
stepped bar is steel with a Young’s modulus of 207 × 109 Pa, find the stresses developed in each of the
steps.

Step 1- Discretization

The domain is subdivided into three elements and four nodes as shown in the figure.

Step 2- Element stiffness matrices

The two element stiffness matrices k1, k2 and k3 are obtained by making calls to the MATLAB function
. linear triangle element stiffness. Each matrix has size 2×2.

1 −1 1 −1
k1 = 5.175×1010 × [ ] k2 = 2.07×1010 × [ ]
−1 1 −1 1
1 −1
k3 = 5.175×109 × [ ]
−1 1

Step 3- Assembling the global stiffness Matrix

Since the structure has four nodes, the size of the global matrix is 4×4. So by assembling k1, k2 and k3
. we get K.

5.175 − 5.175 0 0
−5.175 7.245 −2.07 0
K = 1010 ×[ ]
0 −2.07 2.5875 −0.5175
0 0 −0.5175 0.5175

Step 4- Applying the Boundary conditions

The boundary conditions for this problem are given as:

𝐹3 = 2000𝑁,
 So using these boundary conditions:

5.175 − 5.175 0 0 𝑢1 F1
−5.175 7.245 −2.07 0 𝑢2 F2
1010 ×[ ] ×[𝑢 ] = [ ]
0 −2.07 2.5875 −0.5175 3 2000
0 0 −0.5175 0.5175 𝑢4 F4

Step 5- Solving the equations

Solution of the system equation is performed by Gaussian elimination (with MATLAB).

Step 6- Post processing

In this step, we obtain stresses and strains in both element using MATLAB.

MATLAB functions used

oneD_ElementStiffness(k)
where k is this function calculates the element stiffness matrix for each linear line element where k is element
stiffness which is calculated by EA/L where E is modulus of elasticity, A is the area of cross section, and L is
the length of element. It returns the 2×2 element stiffness matrix k.

oneD_Assemble (K,k,i,j)
this function assembles the element stiffness matrix k of the linear line element joining nodes i and j into the
global stiffness matrix K. it returns the n×n global stiffness matrix K every time an element is assembled.

oneD_ElementStiffness(k)
function y = oneD_ElementStiffness(k)
%element_stiffness i.s this function call the stiffness matrix of 2*2 size
y=[k -k;-k k];
end

oneD_Assemble (K,k,i,j)
function y = oneD_Assemble(K,k,i,j )
% Assemble-this function assemble the element stiffness matrix k
% of element with nodes i and j into the global stiffness
% matrix K
%This function call K after matrix k is assembled
K(i,i)=K(i,i) +k(1,1);
K(i,j)=K(i,j) +k(1,2);
K(j,i)=K(j,i) +k(2,1);
K(j,j)=K(j,j) +k(2,2);
y=K;
end

Generalized MATLAB program for a one dimensional bar problem:

6/23/21 6:18 PM E:\M.TECH\SEM...\Axial_Bar_LinearElement.m 1 of 4

clc; clear; close all

warning('off');
format short g
disp(' ')
Prob_no= input('Question.No:- 1.19','s');
%% Start from here---
% %(NOTE: this programme is only for simplex type element)
%
%% LINEAR BAR ELEMENT PROBLEM*
%%
fprintf('\b [LINEAR BAR ELEMENT PROBLEM]\n\n' )
%% Input Parameter
disp(' ')
disp('GIVEN PARAMETERS ARE:')
disp('---------------------')
m=input('Total no of elements (m) = ');
n=m+1; % for axial bar problem and linear element n=m+1;
dof=1:n; % total node numbering
K=zeros(n); U=zeros(n,1);F=zeros(n,1);
Ke=zeros(1,m); i=1:m; j=2:m+1; E=zeros(1,m); A=zeros(1,m); L=zeros(1,m);
for m=1:m
fprintf('FOR ELEMENT-(%i) => E%i,A%i,L%i = ',m,m,m,m);
EAL=input('','s');
EAL=str2num(EAL);
E(m)=EAL(1);
A(m)=EAL(2);
L(m)=EAL(3);
Ke(m)=E(m)*A(m)/L(m);
end

%% Boundary condition
% using dialoge box , creative way
% I have created a function of dialoge box to access all boundary condition
[dof1,U,dof2,F]=Boundary_condn_dlg(U,F);
% dof1=fixed node or degree of restriction due to given displacement (that nodes where
displacement are given)
% dof2= that node where force value are given
% dof3= free node or that nodes where displacement value are not given
dof3=setxor(1:n,dof1);
for a=dof1
fprintf('\nDisplacement at node-%i (U%i)= %g',a,a,U(a))
end
for a=dof2
fprintf('\nExternal force at node-%i (F%i)= %g',a,a,F(a))
end
%% Solution on the basis of given Input
fprintf('\n');

%% *STEP-1: DISCRETIZATION*

fprintf('STEP-1: DISCRETIZATION \n')

6/23/21 6:18 PM E:\M.TECH\SEM...\Axial_Bar_LinearElement.m 2 of 4

disp('~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nTotal no of node = %i',n);
fprintf('\nTotal no of element = %i\n',m );
elements=1:m;
nodes=[i;j];
T1=table(elements',nodes',L',E',A',Ke','VariableNames',
{'Element','nodes','Length','Young_modulus','Area','Stiffness'});
disp(' ');
disp(T1)

%% *STEP-2: ELEMENTS MATRIX*

fprintf('\nSTEP-2: ELEMENTS MATRIX \n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
for m=1:m
fprintf('\nFOR ELEMENT%i:\n',m);
Kb=oneD_ElementStiffness(Ke(m));
fprintf('K%i=\n',m)
disp(Kb);
K=oneD_Assemble(K,Kb,i(m),j(m));
end
%% *STEP-3: ASSEMBLY OF ELEMENTS MATRIX*

fprintf('\nSTEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
disp('K = ')
disp(K);
%% *STEP-4: APPLYING BOUNDARY CONDITIONS*
%u=sym('U',[n,1]);
%f=sym(F);
fprintf('\nSTEP-4: APPLYING BOUNDARY CONDITIONS\n');
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nGiven :')
if dof1~=0
for a=dof1
fprintf('\nU%i = %g ',a,U(a));
% u(a)=U(a);
% f(a)=sym(sprintf('R%i',a));
end
end
for a=dof2
fprintf('\nF%i = %g',dof2,F(dof2));
end
% fprintf('\nAfter putting the boundary values\n\n')
% u=string(u);u=char(u);
% f=string(f);f=char(f);
% ST1=table(K,u,f,'VariableNames',{'K','U','F'});
% disp(ST1);
% fprintf('\n\nNow delete the row and column crossponding to fixed node\n\n')
% ST2=table(K(dof3,dof3),u(dof3,:),f(dof3,:),'VariableNames',{'K','U','F'});
% disp(ST2);
6/23/21 6:18 PM E:\M.TECH\SEM...\Axial_Bar_LinearElement.m 3 of 4

%% *STEP-5: SOLUTION ( FIND OUT U=? FROM K*U=F )*

fprintf('\n\nSTEP-5: SOLUTION \n');
disp('~~~~~~~~~~~~~~~~~~~~~~~')
disp('Now solve the Equation K*U=F')
S=K;
if dof1~=0
for a=dof1
F=F-K(:,a)*U(a);
K(:,a)=0;
K(a,:)=0;
K(a,a)=1;
F(a)=U(a);
end
end
if size(F,1)~=n && size(U,1)~=n
F=F';
end
U=K\F;
F=S*U;

for dof3=dof3
fprintf('\nU%i= %1.3g ',dof3,U(dof3))
end

if dof1~=0
for a=dof1
fprintf('\nR%i= %1.3g ',a,F(a))
end
end
%% *STEP-6: POST PROCESING*

fprintf('\n\nSTEP-5: POST PROCESING (i.s Element wise solution)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
strain=zeros(1,m);
element_load=zeros(1,m);
stress=zeros(1,m);
for m=1:m
fprintf('\nFOR ELEMENT-%i and node-(%i,%i)',m,i(m),j(m));
j1=i(m);j2=j(m);
strain(m)=(U(j2)-U(j1))/L(m);
stress(m)=E(m)*strain(m);
element_load(m)=Ke(m)*(U(j2)-U(j1));
if element_load(m)>0
X=' (tensile force)';
elseif element_load(m)<0
X=' (compressive force)';
else
X=' (no force)';
end
fprintf('\n (1) displacement: U%i=%1.3g and U%i=%1.3g',i(m),U(i(m)),j(m),U(j(m)))
fprintf('\n (2) strain (e%i)={(U%i-U%i)/L } = %1.3g',m,j(m),i(m),strain(m))
6/23/21 6:18 PM E:\M.TECH\SEM...\Axial_Bar_LinearElement.m 4 of 4

fprintf('\n (3) element load (f%i)={K%i*(U%i-U%i)}= %1.3g ',m,m,j(m),i(m),

element_load(m)),disp(X);
fprintf(' (4) stress=E%i*e%i=%1.3g \n',m,m,stress(m))

end

%% * CONCLUSION :--*

fprintf('\n\nSTEP-6: CONCLUSION :-- \n');

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~')
T2=table([1:n]',U,F,'VariableNames',{'Node','Displacement','Force'});
T3=table(elements', nodes', strain', element_load', stress' ,'VariableNames',
{'Element','nodes','strain','element_load','stress'});
disp('1) Nodes wise Result:')
disp('----------------------')
disp(T2)
disp('2) Elements wise Result:')
disp('------------------------')
disp(T3)
fprintf('\nNOTE: (+ve)element load means tensile in nature')
fprintf('\n and (-ve)element load means compressive in nature\n')
disp('===================================================================')
MATLAB Command Window Page 1

Question.No:- 1.19 [LINEAR BAR ELEMENT PROBLEM]

GIVEN PARAMETERS ARE:

---------------------
Total no of elements (m) = 3
FOR ELEMENT-(1) => E1,A1,L1 = [207*10^9,5,20]
FOR ELEMENT-(2) => E2,A2,L2 = [207*10^9,3,30]
FOR ELEMENT-(3) => E3,A3,L3 = [207*10^9,1,40]

Boundary condition
------------------------
Displacement at node-1 (U1)= 0
Displacement at node-4 (U4)= 0
External force at node-3 (F3)= 2000
STEP-1: DISCRETIZATION
~~~~~~~~~~~~~~~~~~~~~~~~~

Total no of node = 4
Total no of element = 3

Element nodes Length Young_modulus Area Stiffness

_______ ______ ______ _____________ ____ _________

1 1 2 20 2.07e+11 5 5.175e+10
2 2 3 30 2.07e+11 3 2.07e+10
3 3 4 40 2.07e+11 1 5.175e+09

STEP-2: ELEMENTS MATRIX

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT1:
K1=
5.175e+10 -5.175e+10
-5.175e+10 5.175e+10

FOR ELEMENT2:
K2=
2.07e+10 -2.07e+10
-2.07e+10 2.07e+10

FOR ELEMENT3:
K3=
5.175e+09 -5.175e+09
-5.175e+09 5.175e+09
MATLAB Command Window Page 2

STEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
K =
5.175e+10 -5.175e+10 0 0
-5.175e+10 7.245e+10 -2.07e+10 0
0 -2.07e+10 2.5875e+10 -5.175e+09
0 0 -5.175e+09 5.175e+09

STEP-4: APPLYING BOUNDARY CONDITIONS

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Given :
U1 = 0
U4 = 0
F3 = 2000

STEP-5: SOLUTION
~~~~~~~~~~~~~~~~~~~~~~~
Now solve the Equation K*U=F

U2= 2.86e-08
U3= 1e-07
R1= -1.48e+03
R4= -519

STEP-5: POST PROCESING (i.s Element wise solution)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT-1 and node-(1,2)

(1) displacement: U1=0 and U2=2.86e-08
(2) strain (e1)={(U2-U1)/L } = 1.43e-09
(3) element load (f1)={K1*(U2-U1)}= 1.48e+03 (tensile force)
(4) stress=E1*e1=296

FOR ELEMENT-2 and node-(2,3)

(1) displacement: U2=2.86e-08 and U3=1e-07
(2) strain (e2)={(U3-U2)/L } = 2.39e-09
(3) element load (f2)={K2*(U3-U2)}= 1.48e+03 (tensile force)
(4) stress=E2*e2=494

FOR ELEMENT-3 and node-(3,4)

(1) displacement: U3=1e-07 and U4=0
(2) strain (e3)={(U4-U3)/L } = -2.5e-09
(3) element load (f3)={K3*(U4-U3)}= -519 (compressive force)
(4) stress=E3*e3=-519

STEP-6: CONCLUSION :--

~~~~~~~~~~~~~~~~~~~~~~~~~~~
1) Nodes wise Result:
MATLAB Command Window Page 3

----------------------
Node Displacement Force
____ ____________ __________

1 0 -1481.5
2 2.8628e-08 1.0346e-13
3 1.002e-07 2000
4 0 -518.52

2) Elements wise Result:

------------------------
Element nodes strain element_load stress
_______ ______ ___________ ____________ _______

1 1 2 1.4314e-09 1481.5 296.3

2 2 3 2.3856e-09 1481.5 493.83
3 3 4 -2.5049e-09 -518.52 -518.52

NOTE: (+ve)element load means tensile in nature

and (-ve)element load means compressive in nature
===================================================================
>>
Generalized MATLAB program for a one dimensional spring problem:

Q.1 Five springs, having stiffnesses k1 = 105 N/m, k2 = 2 × 105 N/m, k3 = 3 × 105 N/m, k4 = 4 × 105 N/m,
and k5 = 5 × 105 N/m are connected as a parallel-series system and is subjected to a load P = 1000 N at
node 4 as shown in Figure 1. Determine the displacements of nodes 2, 3, and 4 using the finite
element method. State the assumptions made in your solution. (Example 1.6, Chapter 1)

Figure 1
6/23/21 6:28 PM E:\M.TECH\SEM 2\...\Spring_LinearElement.m 1 of 4

clc; clear; close all

warning('off');
format short g
disp('===================================================================')
Prob_no= input('Question.No:- ','s');
%% Start from here---
% %(NOTE: this programme is only for simplex type element)
%
%% SPRING ELEMENT PROBLEM*
%%
fprintf('\b [SPRING ELEMENT PROBLEM]\n\n' )
%% Input Parameter
disp('===================================================================')
disp('GIVEN PARAMETERS ARE:')
disp('===================================================================')
n=input('Total no of nodes (n) = ');
m=input('Total no of elements (m) = ');
dof=1:n; % total node numbering
K=zeros(n); U=zeros(n,1);F=zeros(n,1);
Ke=zeros(1,m); i=zeros(1,m); j=zeros(1,m);
for m=1:m
fprintf('FOR ELEMENT-(%i) => i,j,Ke%i = ',m,m);
IJKe=input('','s');
IJKe=str2num(IJKe);
i(m)=IJKe(1);
j(m)=IJKe(2);
Ke(m)=IJKe(3);
end
%% Boundary condition
% using dialoge box , creative way , looking good
% I have created a function of dialoge box to access all boundary condition
[dof1,U,dof2,F]=Boundary_condn_dlg(U,F);
% dof1=fixed node or degree of restriction due to given displacement (that nodes where
displacement are given)
% dof2= that node where force value are given
% dof3= free node or that nodes where displacement value are not given
dof3=setxor(1:n,dof1);
for a=dof1
fprintf('\nDisplacement at node-%i (U%i)= %g',a,a,U(a))
end
for a=dof2
fprintf('\nExternal force at node-%i (F%i)= %g',a,a,F(a))
end
%% Solution on the basis of given Input
fprintf('\n');
disp('===================================================================')
disp('SOLUTION OF THE PROBLEM STEP WISE STEP:')
disp('===================================================================')
%% *STEP-1: DISCRITIZATION*
6/23/21 6:28 PM E:\M.TECH\SEM 2\...\Spring_LinearElement.m 2 of 4

fprintf('STEP-1: DISCRITIZATION \n')

disp('-----------------------')
fprintf('\nTotal no of node = %i',n);
fprintf('\nTotal no of element = %i\n',m );
elements=1:m;
nodes=[i;j];
T1=table(elements',nodes',Ke','VariableNames',{'Element','nodes','Stiffness'});
disp(' ');
disp(T1)
%% *STEP-2: ELEMENTS MATRIX*

fprintf('\nSTEP-2: ELEMENTS MATRIX \n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
for m=1:m
fprintf('\nFOR ELEMENT%i:\n',m);
Kb=oneD_ElementStiffness(Ke(m));
fprintf('K%i=\n',m)
disp(Kb);
K=oneD_Assemble(K,Kb,i(m),j(m));
end
%% *STEP-3: ASSEMBLY OF ELEMENTS MATRIX*

fprintf('\nSTEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
disp('K = ')
disp(K);
%% *STEP-4: APPLYING BOUNDARY CONDITIONS*
%u=sym('U',[n,1]);
%f=sym(F);
fprintf('\nSTEP-4: APPLYING BOUNDARY CONDITIONS\n');
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nGiven :')
if dof1~=0
for a=dof1
fprintf('\nU%i = %g ',a,U(a));
%u(a)=U(a);
%f(a)=sym(sprintf('R%i',a));
end
end
for a=dof2
fprintf('\nF%i = %g',dof2,F(dof2));
end
% fprintf('\nAfter putting the boundary values\n\n')
% u=string(u);u=char(u);
% f=string(f);f=char(f);
% ST1=table(K,u,f,'VariableNames',{'K','U','F'});
% disp(ST1);
% fprintf('\n\nNow delete the row and column crossponding to fixed node\n\n')
% ST2=table(K(dof3,dof3),u(dof3,:),f(dof3,:),'VariableNames',{'K','U','F'});
% disp(ST2);
6/23/21 6:28 PM E:\M.TECH\SEM 2\...\Spring_LinearElement.m 3 of 4

%% *STEP-5: SOLUTION ( FIND OUT U=? FROM K*U=F )*

fprintf('\n\nSTEP-5: SOLUTION \n');
disp('~~~~~~~~~~~~~~~~~~~~~~~')
disp('Now solve the Equation K*U=F')
S=K;
if dof1~=0
for a=dof1
F=F-K(:,a)*U(a);
K(:,a)=0;
K(a,:)=0;
K(a,a)=1;
F(a)=U(a);
end
end
if size(F,1)~=n && size(U,1)~=n
F=F';
end
U=K\F;
F=S*U;

for dof3=dof3
fprintf('\nU%i= %1.3g ',dof3,U(dof3))
end

if dof1~=0
for a=dof1
fprintf('\nR%i= %1.3g ',a,F(a))
end
end
%% *STEP-6: POST PROCESING*

fprintf('\n\nSTEP-6: POST PROCESING (i.s element wise solution)\n')

disp('--------------------------------------------------------')
element_load=zeros(1,m);
for m=1:m
fprintf('\nFOR ELEMENT-%i and node-(%i,%i)',m,i(m),j(m));
j1=i(m);j2=j(m);
element_load(m)=Ke(m)*(U(j2)-U(j1));
if element_load(m)>0
X=' (tensile force)';
elseif element_load(m)<0
X=' (compressive force)';
else
X=' (no force)';
end
fprintf('\n (1) displacement: U%i=%1.3g and U%i=%1.3g',i(m),U(i(m)),j(m),U(j(m)))
fprintf('\n (2) element load (f%i)={K%i*(U%i-U%i)}= %1.3g ',m,m,j(m),i(m),
element_load(m)),disp(X);

end
6/23/21 6:28 PM E:\M.TECH\SEM 2\...\Spring_LinearElement.m 4 of 4

%% * CONCLUSION :--*
fprintf('\n\nSTEP-7: CONCLUSION :-- \n');
disp('--------------------------~~')
T2=table([1:n]',U,F,'VariableNames',{'Node','Displacement','Force'});
T3=table(elements', nodes', element_load' ,'VariableNames',
{'Element','nodes','element_load'});
disp('(1) Nodes wise Result:')
disp('----------------------')
disp(T2)
disp('(2) Elements wise Result:')
disp('-------------------------')
disp(T3)
fprintf('\nNOTE: (+ve)element load means tensile in nature')
fprintf('\n and (-ve)element load means compressive in nature\n')
disp('===================================================================')
MATLAB Command Window Page 1

===================================================================
Question.No:- Example 1.6 [SPRING ELEMENT PROBLEM]

===================================================================
GIVEN PARAMETERS ARE:
===================================================================
Total no of nodes (n) = 4
Total no of elements (m) = 5
FOR ELEMENT-(1) => i,j,Ke1 = [1,2,10^5]
FOR ELEMENT-(2) => i,j,Ke2 = [1,2,2*10^5]
FOR ELEMENT-(3) => i,j,Ke3 = [2,3,3*10^5]
FOR ELEMENT-(4) => i,j,Ke4 = [3,4,4*10^5]
FOR ELEMENT-(5) => i,j,Ke5 = [3,4,5*10^5]

Boundary condition
------------------------
Displacement at node-1 (U1)= 0
Displacement at node-4 (U4)= 0
External force at node-4 (F4)= 1000
===================================================================
SOLUTION OF THE PROBLEM STEP WISE STEP:
===================================================================
STEP-1: DISCRITIZATION
-----------------------

Total no of node = 4
Total no of element = 5

Element nodes Stiffness

_______ ______ _________

1 1 2 1e+05
2 1 2 2e+05
3 2 3 3e+05
4 3 4 4e+05
5 3 4 5e+05

STEP-2: ELEMENTS MATRIX

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT1:
K1=
100000 -100000
-100000 100000

FOR ELEMENT2:
K2=
200000 -200000
-200000 200000
MATLAB Command Window Page 2

FOR ELEMENT3:
K3=
300000 -300000
-300000 300000

FOR ELEMENT4:
K4=
400000 -400000
-400000 400000

FOR ELEMENT5:
K5=
500000 -500000
-500000 500000

STEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
K =
300000 -300000 0 0
-300000 600000 -300000 0
0 -300000 1200000 -900000
0 0 -900000 900000

STEP-4: APPLYING BOUNDARY CONDITIONS

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Given :
U1 = 0
U4 = 0
F4 = 1000

STEP-5: SOLUTION
~~~~~~~~~~~~~~~~~~~~~~~
Now solve the Equation K*U=F

U2= 0
U3= 0
R1= 0
R4= 0

STEP-6: POST PROCESING (i.s element wise solution)

--------------------------------------------------------

FOR ELEMENT-1 and node-(1,2)

(1) displacement: U1=0 and U2=0
MATLAB Command Window Page 3

(2) element load (f1)={K1*(U2-U1)}= 0 (no force)

FOR ELEMENT-2 and node-(1,2)

(1) displacement: U1=0 and U2=0
(2) element load (f2)={K2*(U2-U1)}= 0 (no force)

FOR ELEMENT-3 and node-(2,3)

(1) displacement: U2=0 and U3=0
(2) element load (f3)={K3*(U3-U2)}= 0 (no force)

FOR ELEMENT-4 and node-(3,4)

(1) displacement: U3=0 and U4=0
(2) element load (f4)={K4*(U4-U3)}= 0 (no force)

FOR ELEMENT-5 and node-(3,4)

(1) displacement: U3=0 and U4=0
(2) element load (f5)={K5*(U4-U3)}= 0 (no force)

STEP-7: CONCLUSION :--

--------------------------~~
(1) Nodes wise Result:
----------------------
Node Displacement Force
____ ____________ _____

1 0 0
2 0 0
3 0 0
4 0 0

(2) Elements wise Result:

-------------------------
Element nodes element_load
_______ ______ ____________

1 1 2 0
2 1 2 0
3 2 3 0
4 3 4 0
5 3 4 0

NOTE: (+ve)element load means tensile in nature

and (-ve)element load means compressive in nature
===================================================================
>>
Generalized MATLAB program for a one dimensional Heat Transfer problem:
6/23/21 6:05 PM E:\M.TECH...\Heat_Transfer_LinearElement.m 1 of 4

clc; clear; close all

warning('off');
format short g
disp(' ')
fprintf('Example 1.9 ');
fprintf('\b HEAT TRANSFER PROBLEM]\n\n' )

disp('GIVEN PARAMETERS ARE:')

disp('---------------------')
m=input('Total no of elements (m) = ');
n=m+1;
dof=1:n;
K=zeros(n); T=zeros(n,1);Q=zeros(n,1);
Ke=zeros(1,m); i=1:m; j=2:m+1; c=zeros(1,m); A=zeros(1,m); L=zeros(1,m);
for m=1:m
fprintf('FOR ELEMENT-(%i) => c%i,A%i,L%i = ',m,m,m,m);
cAL=input('','s');
cAL=str2num(cAL);
c(m)=cAL(1);
A(m)=cAL(2);
L(m)=cAL(3);
Ke(m)=c(m)*A(m)/L(m);
end
%% Boundary Conditions

[dof1,T,dof2,Q]=Boundary_condn_dlg(T,Q);
dof3=setxor(1:n,dof1);
for a=dof1
fprintf('\nTemperature at node-%i (T%i)= %g',a,a,T(a));
end
for a=dof2
fprintf('\nHeat_flow at node-%i (F%i)= %g',a,a,Q(a))
end

fprintf('\n');

%% DISCRETIZATION

fprintf('STEP-1: DISCRITIZATION \n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nTotal no of node = %i',n);
fprintf('\nTotal no of element = %i\n',m );
elements=1:m;
nodes=[i;j];
T1=table(elements',nodes',L',c',A',Ke','VariableNames',
{'Element','nodes','Length','conductivity','Area','Stiffness'});
disp(' ');
disp(T1)

%% ELEMENTS MATRIX
6/23/21 6:05 PM E:\M.TECH...\Heat_Transfer_LinearElement.m 2 of 4

fprintf('\nSTEP-2: ELEMENTS MATRIX \n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
for m=1:m
fprintf('\nFOR ELEMENT%i:\n',m);
Kb=oneD_ElementStiffness(Ke(m));
fprintf('K%i=\n',m)
disp(Kb);
K=oneD_Assemble(K,Kb,i(m),j(m));
end
%% ASSEMBLY OF ELEMENTS MATRIX

fprintf('\nSTEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
disp('K = ')
disp(K);
%% APPLYING BOUNDARY CONDITIONS
t=sym('T',[n,1]);
q=sym(Q);
fprintf('\nSTEP-4: APPLYING BOUNDARY CONDITIONS\n');
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nGiven :')
if dof1~=0
for a=dof1
fprintf('\nT%i = %g ',a,T(a));
t(a)=T(a);
q(a)=sym(sprintf('R%i',a));
end
end
for a=dof2
fprintf('\nQ%i = %g',dof2,Q(dof2));
end
fprintf('\nAfter putting the boundary values\n\n')
t=string(t);t=char(t);
q=string(q);q=char(q);
ST1=table(K,t,q,'VariableNames',{'K','T','Q'});
disp(ST1);
fprintf('\n\nNow delete the row and column crossponding to fixed node\n\n')
ST2=table(K(dof3,dof3),t(dof3,:),q(dof3,:),'VariableNames',{'K','T','Q'});
disp(ST2);
%% SOLUTION ( FIND OUT T=? FROM K*T=Q )
fprintf('\n\nSTEP-5: SOLUTION \n');
disp('~~~~~~~~~~~~~~~~~~~~~~~')
disp('Now solve the Equation K*T=Q')
S=K;
if dof1~=0
for a=dof1
Q=Q-K(:,a)*T(a);
K(:,a)=0;
K(a,:)=0;
K(a,a)=1;
Q(a)=T(a);
6/23/21 6:05 PM E:\M.TECH...\Heat_Transfer_LinearElement.m 3 of 4

end
end
if size(Q,1)~=n && size(T,1)~=n
Q=Q';
end
T=K\Q;
Q=S*T;

for dof3=dof3
fprintf('\nT%i= %1.3g ',dof3,T(dof3))
end

if dof1~=0
for a=dof1
fprintf('\nR%i= %1.3g ',a,Q(a))
end
end
%% POST PROCESING

fprintf('\n\nSTEP-5: POST PROCESING (i.s Element wise solution)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
strain=zeros(1,m);
Heat_flow=zeros(1,m);
stress=zeros(1,m);
for m=1:m
fprintf('\nFOR ELEMENT-%i and node-(%i,%i)',m,i(m),j(m));
j1=i(m);j2=j(m);
strain(m)=(T(j2)-T(j1))/L(m);
stress(m)=c(m)*strain(m);
Heat_flow(m)=Ke(m)*(T(j2)-T(j1));
if Heat_flow(m)>0
X=' (positive heat flow)';
elseif Heat_flow(m)<0
X=' (negative heat flow)';
else
X=' (no heat flow)';
end
fprintf('\n (1) temperature: T%i=%1.3g and T%i=%1.3g',i(m),T(i(m)),j(m),T(j(m)));
fprintf('\n (2) strain (e%i)={(T%i-T%i)/L } = %1.3g',m,j(m),i(m),strain(m));
fprintf('\n (3) Heat_flow (Q%i)={K%i*(T%i-T%i)}= %1.3g ',m,m,j(m),i(m),Heat_flow(m)),
disp(X);
fprintf(' (4) stress=c%i*e%i=%1.3g \n',m,m,stress(m))

end

%% CONCLUSION :

fprintf('\n\nSTEP-6: CONCLUSION :-- \n');

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~')
T2=table((1:n)',T,Q,'VariableNames',{'Node','Temperature','Heat_flow'});
T3=table(elements', nodes', strain', Heat_flow', stress' ,'VariableNames',
6/23/21 6:05 PM E:\M.TECH...\Heat_Transfer_LinearElement.m 4 of 4

{'Element','nodes','strain','Heat_flow','stress'});
disp('1) Nodes wise Result:')
disp('----------------------')
disp(T2)
disp('2) Elements wise Result:')
disp('------------------------')
disp(T3)

disp(' ')
MATLAB Command Window Page 1

Example 1.9 HEAT TRANSFER PROBLEM]

GIVEN PARAMETERS ARE:

---------------------
Total no of elements (m) = 2
FOR ELEMENT-(1) => c1,A1,L1 = [0.25,1,3]
FOR ELEMENT-(2) => c2,A2,L2 = [0.08,1,6]

Boundary condition
------------------------
Temperature at node-3 (T3)= 30
Heat_flow at node-1 (F1)= -0.24
STEP-1: DISCRITIZATION
~~~~~~~~~~~~~~~~~~~~~~~~~

Total no of node = 3
Total no of element = 2

Element nodes Length conductivity Area Stiffness

_______ ______ ______ ____________ ____ _________

1 1 2 3 0.25 1 0.083333
2 2 3 6 0.08 1 0.013333

STEP-2: ELEMENTS MATRIX

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT1:
K1=
0.083333 -0.083333
-0.083333 0.083333

FOR ELEMENT2:
K2=
0.013333 -0.013333
-0.013333 0.013333

STEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
K =
0.083333 -0.083333 0
-0.083333 0.096667 -0.013333
0 -0.013333 0.013333

STEP-4: APPLYING BOUNDARY CONDITIONS

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
MATLAB Command Window Page 2

Given :
T3 = 30
Q1 = -0.24
After putting the boundary values

K T Q
___________________________________ __ _____

0.083333 -0.083333 0 T1 -6/25

-0.083333 0.096667 -0.013333 T2 0
0 -0.013333 0.013333 30 R3

Now delete the row and column crossponding to fixed node

K T Q
______________________ __ _____

0.083333 -0.083333 T1 -6/25

-0.083333 0.096667 T2 0

STEP-5: SOLUTION
~~~~~~~~~~~~~~~~~~~~~~~
Now solve the Equation K*T=Q

T1= 9.12
T2= 12
R3= 0.24

STEP-5: POST PROCESING (i.s Element wise solution)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT-1 and node-(1,2)

(1) temperature: T1=9.12 and T2=12
(2) strain (e1)={(T2-T1)/L } = 0.96
(3) Heat_flow (Q1)={K1*(T2-T1)}= 0.24 (positive heat flow)
(4) stress=c1*e1=0.24

FOR ELEMENT-2 and node-(2,3)

(1) temperature: T2=12 and T3=30
(2) strain (e2)={(T3-T2)/L } = 3
(3) Heat_flow (Q2)={K2*(T3-T2)}= 0.24 (positive heat flow)
(4) stress=c2*e2=0.24

STEP-6: CONCLUSION :--

~~~~~~~~~~~~~~~~~~~~~~~~~~~
MATLAB Command Window Page 3

1) Nodes wise Result:

----------------------
Node Temperature Heat_flow
____ ___________ __________

1 9.12 -0.24
2 12 1.0755e-16
3 30 0.24

2) Elements wise Result:

------------------------
Element nodes strain Heat_flow stress
_______ ______ ______ _________ ______

1 1 2 0.96 0.24 0.24

2 2 3 3 0.24 0.24

>>
Generalized MATLAB program for a one dimensional Fluid flow problem:
6/23/21 5:37 PM E:\M.TECH\SEM 2\F...\Fluid_LinearElement.m 1 of 4

clc; clear; close all

warning('off');
format short g

fprintf('Example 1.7 ');

fprintf('\b [Fluid ELEMENT PROBLEM]\n\n' )

%% Given Parameters

disp('GIVEN PARAMETERS ARE:')

disp('---------------------')
n=input('(i) Total number of nodes = ');
m=n-1;
fprintf('(ii) Total number of elements = %i\n',m);
dof=1:n;
K=zeros(n); P=zeros(n,1);Q=zeros(n,1); D=zeros(1,m);
Ke=zeros(1,m); i=zeros(1,m); j=zeros(1,m);
dofn=1;
DOF=n*dofn;
fprintf('(iii) Degree of Freedom / node = %d\n',dofn);
fprintf('(iv) Degree of freedom of problem is %d\n',DOF);
u=input('(v) Viscosity of fluid is : ','s');
u=str2num(u);
fprintf('(vi) Diameter of pipes are : \n D1,D2,D3,.... - ');
D=input('','s');
D=str2num(D);
fprintf('(vii) Length of pipes are : \n L1,L2,L3,.... - ');
L=input('','s');
L=str2num(L);

C=(1/u)*((D.^4)./L);
Ke=C.*(pi/128);
disp(' ')
disp('(viii)Element Stiffness')
for m=1:m

fprintf(' Ke%i = ',m);

disp(Ke(m));
end
for m=1:m
fprintf('FOR ELEMENT-(%i) => i,j,Ke%i = ',m,m);
IJKe=input('','s');
IJKe=str2num(IJKe);
i(m)=IJKe(1);
j(m)=IJKe(2);
Ke(m)=IJKe(3);
end
%% Boundary Conditions

[dof1,P,dof2,Q]=Boundary_condn_dlg(P,Q);
6/23/21 5:37 PM E:\M.TECH\SEM 2\F...\Fluid_LinearElement.m 2 of 4

dof3=setxor(1:n,dof1);
for a=dof1
fprintf('\nPressure at node-%i (P%i)= %g',a,a,P(a))
end
for a=dof2
fprintf('\nDischarge at node-%i (Q%i)= %g',a,a,Q(a))
end
%% Solution

fprintf('\n');

fprintf('STEP-1: DISCRETIZATION \n')

disp('-----------------------')
fprintf('\nTotal no of node = %i',n);
fprintf('\nTotal no of element = %i\n',m );
elements=1:m;
nodes=[i;j];
T1=table(elements',nodes',Ke','VariableNames',{'Element','nodes','Stiffness'});
disp(' ');
disp(T1)
%% Element Matrix

fprintf('\nSTEP-2: ELEMENTS MATRIX \n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
for m=1:m
fprintf('\nFOR ELEMENT%i:\n',m);
Kb=oneD_ElementStiffness(Ke(m));
fprintf('K%i=\n',m)
disp(Kb);
K=oneD_Assemble(K,Kb,i(m),j(m));
end
%% Assembly of Element Matrix

fprintf('\nSTEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)\n')

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
disp('K = ')
disp(K);
%% Applying Boundary Conditions

p=sym('P',[n,1]);
q=sym(Q);
fprintf('\nSTEP-4: APPLYING BOUNDARY CONDITIONS\n');
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
fprintf('\nGiven :')
if dof1~=0
for a=dof1
fprintf('\nP%i = %g ',a,P(a));
p(a)=P(a);
q(a)=sym(sprintf('R%i',a));
end
6/23/21 5:37 PM E:\M.TECH\SEM 2\F...\Fluid_LinearElement.m 3 of 4

end
for a=dof2
fprintf('\nQ%i = %g',dof2,Q(dof2));
end
fprintf('\n\nAfter putting the boundary values\n\n')
p=string(p);p=char(p);
q=string(q);q=char(q);
ST1=table(K,p,q,'VariableNames',{'K','P','Q'});
disp(ST1);
fprintf('\n\nNow delete the row and column crossponding to fixed node\n\n')
ST2=table(K(dof3,dof3),p(dof3,:),q(dof3,:),'VariableNames',{'K','P','Q'});
disp(ST2);
%% SOLUTION ( FIND OUT P=? FROM K*P=Q )

fprintf('\n\nSTEP-5: SOLUTION \n');

disp('~~~~~~~~~~~~~~~~~~~~~~~')
disp('Now solve the Equation K*P=Q')
S=K;
if dof1~=0
for a=dof1
Q=Q-K(:,a)*P(a);
K(:,a)=0;
K(a,:)=0;
K(a,a)=1;
Q(a)=P(a);
end
end
if size(Q,1)~=n && size(P,1)~=n
Q=Q';
end
P=K\Q;
Q=S*P;

for dof3=dof3
fprintf('\nP%i= %1.3g ',dof3,P(dof3))
end

if dof1~=0
for a=dof1
fprintf('\nR%i= %1.3g ',a,Q(a))
end
end
%% POST PROCESING

fprintf('\n\nSTEP-6: POST PROCESING (i.s element wise solution)\n')

disp('--------------------------------------------------------')
Discharge=zeros(1,m);
for m=1:m
fprintf('\nFOR ELEMENT-%i and node-(%i,%i)',m,i(m),j(m));
j1=i(m);j2=j(m);
Discharge(m)=Ke(m)*(P(j1)-P(j2));
6/23/21 5:37 PM E:\M.TECH\SEM 2\F...\Fluid_LinearElement.m 4 of 4

if Discharge(m)>0
X=' (positive flow)';
elseif Discharge(m)<0
X=' (negative flow)';
else
X=' (no discharge)';
end
fprintf('\n (1) pressure: P%i=%1.3g and P%i=%1.3g',i(m),P(i(m)),j(m),P(j(m)))
fprintf('\n (2) discharge:(q%i)={K%i*(P%i-P%i)}= %1.3g ',m,m,j(m),i(m),Discharge(m)),
disp(X);

end
MATLAB Command Window Page 1

Example 1.7 [Fluid ELEMENT PROBLEM]

GIVEN PARAMETERS ARE:

---------------------
(i) Total number of nodes = 5
(ii) Total number of elements = 4
(iii) Degree of Freedom / node = 1
(iv) Degree of freedom of problem is 5
(v) Viscosity of fluid is : 2*10^-6
(vi) Diameter of pipes are :
D1,D2,D3,.... - [3,2.5,2,1.5]
(vii) Length of pipes are :
L1,L2,L3,.... - [50000,60000,70000,40000]

(viii)Element Stiffness
Ke1 = 19.88

Ke2 = 7.9895

Ke3 = 2.805

Ke4 = 1.5532

FOR ELEMENT-(1) => i,j,Ke1 = [1,2,19.88]

FOR ELEMENT-(2) => i,j,Ke2 = [2,3,7.9895]
FOR ELEMENT-(3) => i,j,Ke3 = [2,4,2.805]
FOR ELEMENT-(4) => i,j,Ke4 = [2,5,1.5532]

Boundary condition
------------------------
Pressure at node-1 (P1)= 30
Pressure at node-3 (P3)= 22
Pressure at node-4 (P4)= 20
Pressure at node-5 (P5)= 25
Discharge at node-2 (Q2)= 0
STEP-1: DISCRETIZATION
-----------------------

Total no of node = 5
Total no of element = 4

Element nodes Stiffness

_______ ______ _________

1 1 2 19.88
2 2 3 7.9895
3 2 4 2.805
4 2 5 1.5532

STEP-2: ELEMENTS MATRIX

MATLAB Command Window Page 2

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FOR ELEMENT1:
K1=
19.88 -19.88
-19.88 19.88

FOR ELEMENT2:
K2=
7.9895 -7.9895
-7.9895 7.9895

FOR ELEMENT3:
K3=
2.805 -2.805
-2.805 2.805

FOR ELEMENT4:
K4=
1.5532 -1.5532
-1.5532 1.5532

STEP-3: ASSEMBLY OF ELEMENTS MATRICES (K)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
K =
19.88 -19.88 0 0 0
-19.88 32.228 -7.9895 -2.805 -1.5532
0 -7.9895 7.9895 0 0
0 -2.805 0 2.805 0
0 -1.5532 0 0 1.5532

STEP-4: APPLYING BOUNDARY CONDITIONS

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Given :
P1 = 30
P3 = 22
P4 = 20
P5 = 25
Q2 = 0

After putting the boundary values

K P Q
__________________________________________________ __ __
MATLAB Command Window Page 3

19.88 -19.88 0 0 0 30 R1
-19.88 32.228 -7.9895 -2.805 -1.5532 P2 0
0 -7.9895 7.9895 0 0 22 R3
0 -2.805 0 2.805 0 20 R4
0 -1.5532 0 0 1.5532 25 R5

Now delete the row and column crossponding to fixed node

K P Q
______ __ __

32.228 P2 0

STEP-5: SOLUTION
~~~~~~~~~~~~~~~~~~~~~~~
Now solve the Equation K*P=Q

P2= 26.9
R1= 61.5
R3= -39.2
R4= -19.4
R5= -2.96

STEP-6: POST PROCESING (i.s element wise solution)

--------------------------------------------------------

FOR ELEMENT-1 and node-(1,2)

(1) pressure: P1=30 and P2=26.9
(2) discharge:(q1)={K1*(P2-P1)}= 61.5 (positive flow)

FOR ELEMENT-2 and node-(2,3)

(1) pressure: P2=26.9 and P3=22
(2) discharge:(q2)={K2*(P3-P2)}= 39.2 (positive flow)

FOR ELEMENT-3 and node-(2,4)

(1) pressure: P2=26.9 and P4=20
(2) discharge:(q3)={K3*(P4-P2)}= 19.4 (positive flow)

FOR ELEMENT-4 and node-(2,5)

(1) pressure: P2=26.9 and P5=25
(2) discharge:(q4)={K4*(P5-P2)}= 2.96 (positive flow)
>>