COLLEGE OF ENGINEERING
DEPARTMENT OF MECHANICAL AND INDUSTRIAL ENGINEERING
GENG 300 / Numerical Methods
Semester
Spring 2018
Course Project
Final Report
Student Information
Amre Dardor 201508721
Mohamed Ali 201509408
Mohammed Alharbi 201407508
Mohammed Al-Harami 201405072
Section
L03
Date
27/05/2018
Coordinator Information
Dr. Sayed El-Sayed
Score
Contents
TOC o “1-3” h z u Problem 1 PAGEREF _Toc515179260 h 31.Problem Statement: PAGEREF _Toc515179261 h 32.Introduction: PAGEREF _Toc515179262 h 33.Methodology PAGEREF _Toc515179263 h 43.1Model Equations PAGEREF _Toc515179264 h 43.2MATLAB Solution Method PAGEREF _Toc515179265 h 63.3Excel Solution Method PAGEREF _Toc515179266 h 64.Solution PAGEREF _Toc515179267 h 74.
1MATLAB Solution PAGEREF _Toc515179268 h 74.2Excel Solution PAGEREF _Toc515179269 h 95.Interpretation PAGEREF _Toc515179270 h 10Problem 2 PAGEREF _Toc515179271 h 111.Problem Statement: PAGEREF _Toc515179272 h 112.Introduction: PAGEREF _Toc515179273 h 113.Methodology PAGEREF _Toc515179274 h 113.1Model Equations PAGEREF _Toc515179275 h 113.2MATLAB Solution Method PAGEREF _Toc515179276 h 123.3Excel Solution Method PAGEREF _Toc515179277 h 134.Solution PAGEREF _Toc515179278 h 134.1MATLAB Solution PAGEREF _Toc515179279 h 134.2Excel Solution PAGEREF _Toc515179280 h 155.Interpretation PAGEREF _Toc515179281 h 17
Problem 1Problem Statement:1303655167386000A long bar of rectangular cross section, 0.4m X 0.6m on a side and having a thermal conductivity of 1.
5 W/m.K, is subjected to the boundary conditions shown below. Two of the sides are maintained at a uniform temperature of 200 C. One of the sides is adiabatic, and the remaining side is subjected to a convection process with Tinf = 30 C and h = 50 W/m2.K. Using an appropriate numerical technique with a grid spacing of 0.1m, determine the temperature distribution in the bar.
Introduction:The principle used to find the temperature at different points with the bar is by dividing the bar into two halves (upper and lower) since the temperature distribution is symmetrical about the horizontal. Then the remaining half of the bar will be divided into a grid with 0.1m X 0.1m squares which gives 20 nodes in each half but only 15 are unknown since the edge has a uniform known temperature of 200 C. 15 equations can be found to determine the temperature of the 15 nodes by using Nodal Finite Difference Equations used in heat transfer calculations.
center285115Line of Symmetry
Line of Symmetry
MethodologyModel Equations5638801000760In order to find the temperature, a model of 15 equations was found using the finite difference method in heat transfer. Each equation depends on the position of the node and follows this convention:
Using the above models, the following equations were found:
Node no. Equation
1 -4T1+2T2+2T4=02 -4T2+T1+T3+2T5=03 -4T3+200+2T6+T2=04 -4T4+T1+2T5+T7=05 -4T5+T2+T6+T8+T4=06 -4T6+T5+T3+200+T9=07 -4T7+T4+2T8+T10=08 -4T8+T7+T5+T9+T11=09 -4T9+T8+T6+200+T12=010 -4T110+T7+2T11+T13=011 -4T11+T10+T8+T12+T14=012 -4T12+T11+T9+200+T15=013 2T10+T14+6.6666?30-10.6666T13=014 2T11+T13+T15+6.666?30-23.333+2T14=015 2T12+T14+200+6.666?30-23.333+2T15=0center93789500Which can be integrated into the following matrix [A] [T] = [C] with [T] being the unknown temperatures.AT=CMATLAB Solution MethodIn order to solve this system of equations using MATLAB, Gauss elimination with partial pivoting was used. Gauss elimination uses the pivot or the diagonal element in each column in order to eliminate the elements below it. This is called forward elimination. After this is done, backward substitution is performed which finds the values of the variables from the bottom to the top of the matrix. This results in the final solution vector T. Partial pivoting is an addition to Na?ve Gauss elimination in that it reorders the rows to put the largest element in the pivot position in order to obtain better results and avoid issues associated with ill-conditioned systems.
Excel Solution MethodThe solution in excel used built-in excel functions. First the system of equations was input as matrix A and vector C. Then the inverse of matrix A was found using the MINVERSE function that takes an array as input. Then the inverse of A was multiplied by C using MMULT function. The result was vector T which is the temperatures at each node.
SolutionMATLAB SolutionCode
A=[ -4 2 0 2 0 0 0 0 0 0 0 0 0 0 0;
1 -4 1 0 2 0 0 0 0 0 0 0 0 0 0;
0 1 -4 0 0 -2 0 0 0 0 0 0 0 0 0;
1 0 0 -4 2 0 1 0 0 0 0 0 0 0 0;
0 1 0 1 -4 1 0 1 0 0 0 0 0 0 0;
0 0 1 0 1 -4 0 0 1 0 0 0 0 0 0;
0 0 0 1 0 0 -4 2 0 1 0 0 0 0 0;
0 0 0 0 1 0 1 -4 1 0 1 0 0 0 0;
0 0 0 0 0 1 0 1 -4 0 0 1 0 0 0;
0 0 0 0 0 0 1 0 0 -4 2 0 1 0 0;
0 0 0 0 0 0 0 1 0 -1 -4 1 0 1 0;
0 0 0 0 0 0 0 0 1 0 1 -4 0 0 1;
0 0 0 0 0 0 0 0 0 2 0 0 -10.66 2 0;
0 0 0 0 0 0 0 0 0 0 2 0 1 -10.66 1;
0 0 0 0 0 0 0 0 0 0 0 2 0 1 -10.66;
] ;
C=[0 0 -200 0 0 -200 0 0 -200 0 0 -200 -200 -200 -400]’;
[m,n]=size(A);
if m~=n, error(‘Matrix A must be square’); end
nb=n+1;
Aug=[A C];
% forward elimination
for k = 1:n-1
% partial pivoting
[big,i]=max(abs(Aug(k:n,k)));
ipr=i+k-1;
if ipr~=k
Aug([k,ipr],:)=Aug([ipr,k],:);
end
for i = k+1:n factor=Aug(i,k)/Aug(k,k);
Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
end
end
% back substitution
T=zeros(n,1);
T(n)=Aug(n,nb)/Aug(n,n);
for i = n-1:-1:1
T(i)=(Aug(i,nb)-Aug(i,i+1:n)*T(i+1:n))/Aug(i,i);
end
format short
disp(T)
Results
Temp
65.9792
58.9538
13.0532
73.0046
78.3915
103.3706
69.2561
78.2370
122.0377
47.5458
43.2626
106.5433
34.4019
35.8161
60.8727
Excel Solutioncenter18097500
center287528000
21418553125470center000
InterpretationThe main target of this problem is to find the temperature distribution at different locations on a bar. It was found that the results from MATLAB and Excel were almost identical. Which means that the two techniques that have been used (Gauss Elimination and Inversion) produce the same results for this particular system. However, the solutions using excel was easier and faster to complete. On the other hand, MATLAB allows greater flexibility as this scrip that was used could be converted into a function to solve most similar systems even with different sizes.
Problem 2Problem Statement:63817537655500Determine the force in each member of the truss shown
Introduction:In this problem the method of joints will be used to solve the problem, which is by calculating the forces on every joint in this structure. Furthermore, in this question the Newtonian static force equations will be used which are Fx=0 and Fy=0, then there will be 15 unknowns and 15 equations, and the matrix of this problem will be ready to construct.
MethodologyModel EquationsMA=0Fx=0Fy=018Gy=1044FBc+FBAcos53.8=0FAC+FBAsin53.8+Ay=0FAH-FBAcos53.8=0-FAC-FCHsin53.8=0-FBc+FCD+FCHcos53.8=0FDH+FEHsin53.8+sin53.8=60-FAH+FGH+FEHcos53.8-FCHcos53.8=0FDH=0FEF-FEHcos53.8-FDE=0-FGE-FEHsin53.8=0-FFGcos53.8-FEF=0-FFGsin58.3=24-FGH+FFGcos53.8=0FGE+FFGsin53.8+Gy=0MATLAB Solution MethodThe inverse matrix method will be used to solve this matrix problem by using MATLAB. For matrices, normal division method is inapplicable like the normal equations. So, this method is using the multiplication of the original matrix inverse for both sides of the equation, in order to get the final solution.
Matrix original form
Ax=b . (equation 1)
Where A and b are matrices. So,
x?bAOn the other hand, multiplying (equation 1) by matrix (A) inverse
A-1Ax=A-1bfinal answer form
x=A-1bExcel Solution MethodThe method used of the solution in excel is multiplying the inverse of the main matrix by the vector B which is the values of all equations, then the values of the forces will occur. Thus, to do these steps in excel, built in functions have been used which are MINVERSE to find the inverse of the main matrix, and MMULT to multiply the inverse of the matrix by the vector B to find values of the forces.
SolutionMATLAB SolutionCode
a=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 18;
-0.9238 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
-0.38289 0 1 0 1 0 0 0 0 0 0 0 0 0 0;
0.92379 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 -1 0 0.38289 0 0 0 0 0 0 0 0;
0 -1 0 0 0 1 -0.92379 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 -0.38289 -0.38289 1 0 0 0 0 0 0;
0 0 0 -1 0 0 0.92379 -0.92379 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0.92379 0 -1 0 1 0 0 0;
0 0 0 0 0 0 0 0.38289 0 0 -1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 -1 0.92379 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0.38289 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -0.92379 -1 0;
0 0 0 0 0 0 0 0 0 0 1 0 -0.38289 0 1;];
b=[1044;0;0;0;0;0;60;0;0;0;0;0;24;0;0];
forces=inv(a)*b
Results
Forces
41.7875
38.6033
42.0000
-38.6028
-26.0000
-24.1264
-67.9046
-88.7983
0
-24.1268
-34.0000
57.9043
62.6812
-57.9043
58.0000
Excel SolutionAx = B
InterpretationThe main issue of this problem is to calculate the forces on each node in the truss. The result for both EXCEL and MATLAB were very near to each other, which means that the technique that been used produced same answer for this problem. The reason behind this difference in results could be due to rounding and the way that each software (excel and MATLAB) calculate the results internally.