MATLAB

ME 508, Spring 2017, Mechanical Engineering, University of Kansas Page 1 of 4

ME 508: Project 1 Due Date and Time: Friday March 3rd 2017, beginning of the lecture (12:30 PM).

A thermos which is made out of three concentric cylinders is shown in Figure 1. The innermost cylinder, C0, is at a temperature of T0, the middle cylinder, C1, is at a temperature of T1, and the outermost cylinder, C2, is at a temperature of T2. It is exposed to the outside atmospheric air which is at a temperature of T3.

Figure 1: Thermos with concentric cylinders. Ref. 1

Cylinder C0 and C1 are separated by vacuum and the only means of heat transfer between them is by thermal radiation, and the heat transfer rate per unit area (heat flux) is q01. There is some air between Cylinder C1 and C2, and thus along with radiation there is also some convection heat transfer as well, and the heat flux between them is q12. Finally, the heat flux between the outermost cylinder C2 and the surrounding air is by convection and the heat flux is given by q23. Heat fluxes q01, q12, and q23, are the rates once the thermos has reached steady sate (so temperature does not change with respect to time). Thus, the following relation holds:

01 12 23q q q= = (1)

The expression for the above heat transfer rates (assuming 1D heat transfer in radial direction only) in terms of the temperatures and material coefficients is given below in equations

( )4 4 201 0 1 Wq T T mσε= − (2) ( ) 212 1 24 Wq T T m= − (3)

ME 508, Spring 2017, Mechanical Engineering, University of Kansas Page 2 of 4

( ) 4

3 223 2 31.3 Wq T T m

= − (4)

Where the temperature T0, T1, T2, and T3 in equations (2), (3), and (4) are in kelvin. σ = 5.67x10- 8 W/m2K4 is the Stefan-Boltzmann constant, and ε is the emissivity coefficient, and the product σε = 10-9 W/m2K4.

Given that the temperatures of the innermost cylinder is at T0 = 6000C and the surrounding atmospheric air temperature, T3=200C, calculate the value of the temperature of the inner middle and the outer cylinder, T1 and T2, respectively.

You will write a MATLAB program to solve for the values of the unknown temperatures T1 and T2, by using Generalized Newton-Raphson method (GN-R method) to solve the two system of non-linear algebraic equations that you must derive by making use of equations (1), (2), (3), and (4).

1. First obtain a system of two non-linear algebraic equations in T1 and T2. Then write a MATLAB code to solve for T1 and T2 using GN-R method.

2. Determine a suitable starting solution for T1 and T2. 3. Convergence criteria: The norm of the residual ||R|| has to be less than or equal to 10 – 12 . 4. After the iteration has converged, your program should print the final values of T1 and T2

on the command window, in 0C. 5. Also display the number of iterations it took to converge (this excludes 0th iteration), ||R||

for the final iteration in scientific notation, and approximate relative percent errors of T1 and T2 (𝜀𝜀𝑎𝑎,𝑇𝑇1 and 𝜀𝜀𝑎𝑎,𝑇𝑇2, respectively), on the command window.

6. Plot 𝜀𝜀𝑎𝑎,𝑇𝑇1 and 𝜀𝜀𝑎𝑎,𝑇𝑇2 as a function of the iteration number in a separate figure in your MATLAB code. NOTE: While plotting 𝜀𝜀𝑎𝑎, only plot from the 1st iteration and onwards. 𝜀𝜀𝑎𝑎 for the 0th iteration does not exist, so don’t include that while plotting it as a function of iteration number.

7. Plot ||R||, as a function of the iteration number in another figure. Here||R|| exists right from the 0th iteration till the last iteration.

8. For both the error and norm plots, make sure the plot has markers at the data points (could be a square or circle or anything that’s presentable), along with connecting lines. The y-axis should be in log scale.

9. From your two system of non-linear algebraic equations you derived in step (1), express T2 as a function of T1 for both the equations. Then, plot T2 from these two equations for a range of T1 (from a value of air temperature T3 to a value of innermost cylinder temperature T0) on the same plot so that you will have two curves for these two equations. The units of the temperatures on this figure must be in 0C. Look at the coordinates of the point of intersection of these two curves and see if it is close to the solution of T1 and T2 you obtained from GN-R method. Comment on this.

ME 508, Spring 2017, Mechanical Engineering, University of Kansas Page 3 of 4

Please remember to label your plot axes, legend, and title! Also have the grids on all the plots.

10. Hint: For starting solution, please note that since T0>T3, T1>T2. You could pick a value of T1 that is close to T0 and a value of T2 closer to T3, and run your program and check for convergence.

Submission Format: Write a report that contains a title page, executive summary, introduction, numerical analysis procedure, results and discussion, conclusion, and appendix. Must also include table of contents, lift of figures, etc.

• Title page must include title, group number, names, date of submission, location, school, etc.

• Executive summary is an abstract of your report, after the title page but before table of contents, list of figures, etc.

• The problem statement, the heat transfer equations given above, a note about GN-R method, etc. should all be in the introduction.

• Numerical Analysis and Procedure must include your system of two-nonlinear algebraic equations, procedure to solve this using GN-R method, basic algorithm, Jacobean matrix [J], determination of initial solutions required for GN-R method, convergence criteria and tolerance, reformulating the two equations such that both are expressed with T2 as a function of T1, etc.

• Results and Discussion section must contain your final solution of T1 and T2, all the plots, any other results that has been asked in the problem statement, etc. Plots and results must be discussed. Also for all the iterations, tabulate the following values: iteration number, T1, T2, f(x,y), g(x,y), 𝜀𝜀𝑎𝑎,𝑇𝑇1, 𝜀𝜀𝑎𝑎,𝑇𝑇2, ||R||.Where, f(x,y) is the residual error of the first equation, g(x,y) is the residual error of the second equation. The following is an example of what I expect in the contents of your table.

Iteration Number

T1 (0C)

T2 (0C)

f(x,y) g(x,y) εa,x (%)

εa,y (%)

||R||

0 - NA -

- NA -

1 . . .

nth iteration

ME 508, Spring 2017, Mechanical Engineering, University of Kansas Page 4 of 4

• Your MATLAB code (and pseudo code if you include that as well) must be in the appendix.

• Both the hard and soft copy of the report must be submitted by the due date and time. The soft copy must be in a pdf format. Both the report and your MATLAB code(s) must be inside a folder and this folder must be zipped and uploaded on to the Blackboard by any one member of your team in the assignments page link for Project 1.

• The hard copy must be printed (you may print it in black and white) on only one side, and stapled at the top left hand corner. Do not fold the report in half.

• MATLAB: Your MATLAB code must be modularized. So your program will have a main program and function(s) that the main program calls. One such function that is necessary is a function that takes a vector of unknown of the current iteration ( {x}i ) as the input argument, and calculates and returns the vector of functions ( {f}i ) and the Jacobean matrix ([J]i ) to the calling program. You may have more functions as you deem is necessary.

• Your code must be well commented. Every function must have a comment on the top that lists the input and output arguments of that function!

Additional Notes: * Enclose your final solutions in a box. * You may NOT use readily available codes online like MATLAB fileexchange, stack exchange, etc. The code will be checked automatically by Blackboard for such submissions. * You may NOT use MATLAB inbuilt functions to automatically solve for the non-linear equations (such as fsolve, fzero, etc.), except perhaps to check your final answers. * You may use MATLAB symbolic toolbox to calculate the derivatives in the [J], but your report will still have expression for the [J] in the section on procedure.

Reference:

1. Chapra, S., Canale, R., “Numerical Methods for Engeineers,” 7th Edition, McGraw-Hill, April 2009.