TRƯỜNG ĐẠI HỌC KHOA HỌC VÀ CÔNG NGHỆ HÀ NỘI
UNIVERSITÈ DES SCIENCES ET DES TECHNOLOGIES DE HANOI
UNIVERSITY OF SCIENCE AND TECHNOLOGY OF HANOI
FINITE ELEMENT MODELING METHODS
REPORT
Student: NGÔ VĂN ĐƯỢC
ID: USTHBI6-039
DEPARTMENT OF SPACE AND AERONAUTICS
[PROBLEMS] – 1D
Compute the displacement along the bar with A=1 in2, E=10*10^6 psi, L=30 in., and P=1000 lb.
Here A, E, and L are generic symbols used for cross-sectional area, modulus of elasticity, and length,
respectively. P is the force at the center of the bar. We obtain the differential equation governing the
linear-elastic bar behavior as
d du
AE 0
dx dx
where u is the axial displacement function along the element in the x direction, and A and E will be
assumed constant over the whole length of the bar in our derivations to follow. And the displacements at
two ends are zero.
[Analytical solution]
Px
From x=0 to center: u
2 AE
Px PL
From center to ending point: u
2 AE 2 AE
1. CHOOSE ELEMENT TYPE AND GENERATE THE MESH
Element type: Linear
To find out the solution with arbitrary number of elements, suppose that the mesh is uniform.
The number of Elements is: Ne (arbitary)
So the number of Nodes is: Nd=Ne + 1
%GENERATE THE DISTANCE MATRIX TO PLOT THE FIGURE OF "DISPLACEMENT VS
DISTANCE"
Dist=zeros(1,Nd); %Distance matrix has the same size as Displacement matrix
for i=1:Nd
Dist(i)=(10/Ne)*(i-1);
end;
h=Dist(2)-Dist(1); %The distance between 2 ends of an element.
2. DERIVE THE GALERKIN’S INTERGRAL TO FIND THE ELEMENT
EQUATION
I now derive Galerkin’s integral to formulate the bar element stiffness equations for 1 arbitrary element.
I begin with the basic differential equation, without distributed load.
d du
AE 0 (1)
dx dx
where A and E are constants.
Applying Galerkin’s criterion with the weighting function N k , I have
xj
d du
dx AE dx N
xi
k dx 0 with k = i, j (2)
i: is the beginning point of the element .
j: is the ending point of the element.
Applying integration by parts udv uv vdu to solve Equation (2)
Letting
dNk
u Nk du dx
dx
d du du
dv AE dx v AE
dx dx dx
I find that Equation (2) becomes
xj xj
du du dN k
N k AE
dx xi xi
AE
dx dx
dx 0 (3)
Because u N d , we have
du dNi dN j
ui uj (4)
dx dx dx
xj x x xi
With N i Nj , Equation (4) becomes
x j xi x j xi
du 1 1 ui
( h x j xi )
dx h h u j
(5)
Using Equation (5) in Equation (3), we then express Equation (3) as
1 ui
xj x
1 du
j
dN
AE k h dx N k AE
h u j
(6)
xi dx dx xi
Equation (6) is really two equations (one for N k N i and one for Nk N j ).
First, using the weighting function N k N i , we have
1 ui
xj x
1 du
j
dN
AE i h dx Ni AE (7)
xi dx h u j dx xi
d xj x 1 1
Substituting for dNi / dx , we obtain
dx x j xi x j xi h
1 ui
xj
1 1
AE dx fix
h h h u j
(8)
xi
xj
du
where fix Ni AE
dx xi
Evaluating Equation (8)
We have
1 ui
xj
1 1
AE dx
xi h h h u j
xj
1 1 1
AE ui u j dx
xi h h h
xj
AE 1 1
ui u j x
h h h xi
AE 1 1
ui u j x j xi
h h h
AE 1 1
ui u j h
h h h
AE
h
ui u j
Therefore, Equation (8) is equivalent to
AE
h
ui u j fix (9)
Similarly, using
Nk N j
d x xi 1 1
dN j / dx , we obtain
dx x j xi x j xi h
1 ui
xj x
1 1 du
j
AE dx N j AE
xi h h h u j dx xi (10)
Simplifying Equation (10) yields
AE
h
ui u j f jx (11)
xj
du
where f jx N j AE
dx xi
So for 1 arbitary element, we have 2 equations: (9) and (11)
AE
h
ui u j fix
AE
h
ui u j f jx
Or in the matrix form
AE 1 1 ui fix
h 1 1 u j f jx
And that’s the result of Galerkin’s integral.
3. ASSEMBLE THE GLOBAL STIFFNESS MATRIX K AND GLOBAL
LOAD VECTOR F, IMPOSE BOUDARY CONDITIONS
Symbolically, we write
Ne Ne
F f (e)
and
K k (e)
e 1 e 1
It is easier to illustrate by choosing a certain number of element. Here I choose 6 elements
We have the Total equation (6 elements so we have 7 nodes, or 7 equations)
1 1 0 0 0 0 0 u1 F1x
1 1 1 1
0 0 0 0 u2 F2 x
0 1 1 1 1 0 0 0 u3 F3 x
AE
0 0 1 1 1 1 0 0 u4 F4 x
h
0 0 0 1 1 1 1 0 u5 F5 x
0 0 0 0 1 1 1 1 u6 F6 x
0
0 0 0 0 1 1 u7 F7 x
We have the initial conditions: u1 = u7 = 0
F4 x P 1000 (F at center = P)
Other Fs = 0 (except for F1x and F7 x )
So the matrices become
1 1 0 0 0 0 0 0 F1x
1 1 1 1
0 0 0 0 u2 0
0 1 1 1 1 0 0 0 u3 0
AE
0 0 1 1 1 1 0 0 u4 1000
h
0 0 0 1 1 1 1 0 u5 0
0 0 0 0 1 1 1 1 6u 0
0
0 0 0 0 1 1 0 F7 x
Now we can see that we have 5 variables (u2 to u6) and 5 equations. That’s enough!
1 1 0 0 0 0 0 0 F1x
1 1 1 1
0 0 0 0 u2 0
0 1 1 1 1 0 0 0 u3 0
AE
0 0 1 1 1 1 0 0 u4 1000
h
0 0 0 1 1 1 1 0 u5 0
0 0 0 0 1 1 1 1 6u 0
0
0 0 0 0 1 1 0 F7 x
But precisely, the matrices are:
0
1 1 1 1 0 0 0 0 u2 0
0 1 1 1 1 0 0
0 u3 0
AE
0 0 1 1 1 1 0 0 u4 1000
h
0 0 0 1 1 1 1 0 u5 0
0 0 0 0 1 1 1 1 57 u6 0 51
0 71
We now can solve this matrix equation by MATLAB
How to actually do it in MATLAB:
Firstly, I will generate the full stiffness matrix
1 1 0 0 0 0 0
1 1 1 1 0 0 0 0
0 1 1 1 1 0 0 0
AE 0 0 1 1 1 1 0 0
0 0 0 1 1 1 1 0
0 0 0 0 1 1 1 1
0 1 1
0 0 0 0
%GENERATE THE GLOBAL STIFFNESS MATRIX K(Nd,Nd)
Main_cen = 2*ones(1,Nd); %Vector of values at the main diagonal
Other_cen=-ones(1,Nd-1); %vector of values at the upper and lower diagonal
a=diag(Main_cen); %The matrix contents the main diagonal values
b=diag(Other_cen,-1); %The matrix contents the lower diagonal values
c=diag(Other_cen,+1); %The matrix contents the upper diagonal values
K=a+b+c; %Sum of three matrices to form the matrix K
K(1,1)=1; %Assign the "first value" to 1
K(Nd,Nd)=1; %Assign the "end value" to 1
Secondly, I delete the first and last rows. And we have
1 1 1 1 0 0 0 0
0 1 1 1 1 0 0 0
AE
0 0 1 1 1 1 0 0
h
0 0 0 1 1 1 1 0
0 0 0 0 1 1 1 1 57
%Very IMPORTANT
New_K=K([2:Nd-1],:)*(A*E/L); %Cut the first and last rows of K
%This matrix is for solving U in MATLAB
4. SOLVE TO FIND DISPLACEMENT VALUES
And now we can solve to find vector U by just one command in MATLAB.
U = New_K\F;
5. RESULTS AND POST-PROCESSING
When Ne=1000 elements
1) We can see that the graph is consistent with the analytical solution: a linear line y=Ax+B
2) Along the first-half distance of the bar, the displacement function has positive slope. And long
the remaining-half distance of the bar, the displacement function has negative slope. This agrees
that the force is at the center.
3) Actually I plot both Finite solution and Analytical solution, but they appear only the red color of
Finite solution. The reason is that 2 lines are very likely the same.
plot(Dist,Ana_total,'blue',Dist,U,'red')
4) And the error is very low.
6. CODE REFERENCES
%FINITE ELEMENT METHODS
%1D BAR WITH LINEAR ELEMENTS
%WRITTEN BY: LM.DC
%DEPARTMENT OF SPACE AND AERONAUTICS
%NOTE: YOU CAN CHOOSE ARBITRARY NUMBER OF ELEMENTS
clc
clear all
close all
disp('The displacement values is in matrix U')
%PHYSICAL CONDITIONS
A=1; %Cross-sectional area
E=10*10^6; %Modulus of elasticity
L=30; %The length of the bar
P=1000; %Force applied at the center of the bar
%U: The matrix of displacement values
%F: The global matrix of forces
%K: The global stiffness matrix
%CHOOSE THE NUMBER OF ELEMENTS
Ne=100 %The number of elements
Nd=Ne+1; %The number of nodes
%GENERATE THE DISTANCE MATRIX TO PLOT THE FIGURE OF "DISPLACEMENT VS
DISTANCE"
Dist=zeros(1,Nd); %Distance matrix has the same size as Displacement matrix
for i=1:Nd
Dist(i)=(L/Ne)*(i-1);
end;
h=Dist(2)-Dist(1); %The distance between 2 ends of an element.
%GENERATE THE GLOBAL STIFFNESS MATRIX K(Nd,Nd)
Main_cen = 2*ones(1,Nd); %Vector of values at the main diagonal
Other_cen=-ones(1,Nd-1); %vector of values at the upper and lower diagonal
a=diag(Main_cen); %The matrix contents the main diagonal values
b=diag(Other_cen,-1); %The matrix contents the lower diagonal values
c=diag(Other_cen,+1); %The matrix contents the upper diagonal values
K=a+b+c; %Sum of three matrices to form the matrix K
K(1,1)=1; %Assign the "first value" to 1
K(Nd,Nd)=1; %Assign the "end value" to 1
%Very IMPORTANT
New_K=K([2:Nd-1],:)*(A*E/h); %Cut the first and last rows of K
%This matrix is for solving U in MATLAB
%GENERATE THE GLOBAL FORCE MATRIX
F=zeros(Nd-2,1); %Cut the first and last rows of F to solve U in MATLAB
F(int64(Nd-2)/2,1)=P ;%Assign FORCE at the center of the bar
%int64: return to integer number from double number
%With odd number of elements: force at left-center
%With even number of elements: force at exact-center
%THE MAIN EQUATION: SOLVING THE VALUES OF DISPLACEMENT U
U=New_K\F;
%PLOT THE FIGURE OF "DISPLACEMENT VS DISTANCE"
p=plot(Dist,U,'blue');
xlabel('Distance (inch)');
ylabel('Displacement (inch)');
title('The displacement of the bar')
%THIS BLOCK AIMS TO COMPARE THE ANALYTICAL SOLUTION VS FINITE SOLUTION.
Dist_1=Dist(:,[1:int64(Nd/2)]); %First part of the bar
Dist_2=Dist(:,[int64(Nd/2)+1:end]); %Last part of thr bar
Ana_u1=Dist_1*(P/(2*A*E));%Ana slope %First Ana solution of the bar
Ana_u2=Dist_2*(-P/(2*A*E))+P*L/(2*A*E); %Last Ana solution of the bar
Ana_total=[Ana_u1,Ana_u2]; %Combine first and last parts to form
%Ana total solution
figure;
plot(Dist,Ana_total,'blue',Dist,U,'red')%Plot Ana line and Finite line
%together
xlabel('Distance (inch)');
ylabel('Displacement (inch)');
title('Comparision between Finite solution and Analytical Solution')
%THIS BLOCK AIMS TO FIND THE ERROR BETWEEN ANALYTICAL AND FINITE VALUES
Error=abs(([[Ana_total]'-U]./[Ana_total]'))*100; %transpose Ana_total
%to get the same dimension
figure;
plot(Dist,Error); %Plot Error vs Distance
xlabel('Distance (inch)');
ylabel('Error (%)');
title('Error along the bar distance')
%THIS BLOCK FUNTIONS TO MEASURE COMPUTING TIME
x = rand(100);
tic;
x([3,9],:) = [];
toc; % Elapsed time is 0.000230 seconds.
%END.
