Project 1 Numerical Methods I: MIE334 H1S
Department of Mechanical and Industrial Engineering University of Toronto
Winter 2020
Project 1
MIE334: Numerical Methods I
Due Date: Tuesday, March 3rd, 2020, 11:59 pm
Problem Statement:
In heat transfer, fins are used to increase the effective surface area of a heat exchanger and, thus, enhance the rate of convective heat transfer. As an example, fins are extensively used in electronic circuit boards to dissipate the generated heat by transistors to the environment.
Figure 1. A finned heatsink used to cool an integrated circuit.
Here, we are dealing with a single fin with triangular cross-section attached to a base object
as shown in Figure 2. A constant heat flux of 𝑞′′ = 3.2 𝑊/𝑐𝑚2 enters the fin through its 𝑏
base. Also, the fin is being cooled with a stream of air having a temperature of 𝑇 = 300 𝐾. ∞
The convective heat transfer coefficient between the fin and cooling air is assumed to be h𝑏 = 18 𝑊/𝑚2𝐾 at the fin base and decrease linearly towards the tip of the fin where h𝑡 = 6 𝑊/𝑚2𝐾. The fin is made of aluminum for which the properties―density, specific heat capacity, conductivity, and melting point―are listed in Table 1 along with the geometrical dimensions of the fin. The fin and the base are designed in such a way that 𝑤 is much longer than 𝐿 (i.e., 𝑤 ≫ 𝐿).
Figure 2. Schematics of the triangular fin.
Table 1. Geometrical dimensions of the fin along with properties of aluminum. 𝑏 (𝑐𝑚) 𝐿 (𝑐𝑚) 𝜌 (𝑘𝑔/𝑚3) 𝑐𝑝 (𝐽/𝑘𝑔 𝐾) 𝑘 (𝑊/𝑚𝐾) 𝑀𝑒𝑙𝑡𝑖𝑛𝑔 𝑃𝑜𝑖𝑛𝑡 (𝐾)
1 7 2720 871 202 930
1
We are trying to obtain the temperature profile inside the fin as a function of 𝑥. For this purpose, a control volume approach can be used which is explained in the following section.
Solution Approach
As mentioned above, the width 𝑤 of the fin is much bigger than its length 𝐿. Therefore, temperature variations along the width of the fin can be neglected and 𝑤 can be set to unity (i.e., 𝑤 = 1). Also, since the fin is relatively thin (due to its thickness being smaller than its length) and the conductivity of aluminum is high, temperature gradients in the 𝑦 direction can also be neglected without losing much accuracy. Based on these assumptions, the problem can be solved one-dimensionally only in the 𝑥 direction. To further simplify the problem, radiation from the fin to the surrounding is neglected and the fin is considered to be in a steady state; meaning that the temperature profile does not change with time.
Figure 3. Schematics of the solution domain divided into 𝑛 control volumes. The length of all the control volumes in the 𝑥-direction is ∆𝑥 except the first and last ones for which the length is ∆𝑥.
2
As seen in Figure 3, the fin is divided into 𝑛 control volumes. The first and the last control volumes include the boundary conditions. Due to the triangular shape of the fin, the cross- section of the fin changes along the x-direction which should be carefully accounted for in the solution procedure.
Figure 4. Geometrical parameters of a control volume (left) and heat fluxes entering/exiting the control volume (right). Note that 𝐴𝑐 is NOT equal to ∆𝑥.
𝑞′′ 1 𝑖−2
𝑞′′ 1 𝑖+2
𝑞′′ 𝑐,𝑖
𝑞′′ 𝑐,𝑖
2
For the control volume at node 𝑖, the energy equation can be written by balancing the energies entering and exiting the boundaries of the control volume as follows:
∑(𝑞̇ −𝑞̇ )=0 → 𝑞′′ ×𝐴 1−𝑞′′ ×𝐴 1−2×𝑞′′ ×𝐴 =0
𝑖𝑛 𝑜𝑢𝑡 𝑖−1 𝑖− 𝑖+1 𝑖+ 𝑐,𝑖 𝑐 (1)
2222
where 𝑞′′ 1 and 𝑞′′ 1 show the rate of heat conduction per unit area entering/exiting the left (𝑖 − 1) 𝑖−2 𝑖+2 2
and right (𝑖 + 1) boundary surfaces of the control volume 𝑖, respectively, and 𝑞′′ shows the rate of 2 𝑐,𝑖
heat convection exiting from the upper and lower boundaries to the environment per unit area. All the fluxes entering the control volume are taken to be positive (adding to the energy content) while the exiting fluxes are considered to be negative (reducing the energy content). From heat transfer you know that the rate of heat conduction per unit area (flux) can be obtained using the Fourier’s law:
𝑞′′ = −𝑘𝑑𝑇 𝑑𝑥
(2)
Whereas, the convective heat transfer per unit area can be obtained using the Newton’s cooling law:
𝑞′′ = h(𝑇 − 𝑇 ) (3) 𝑐∞
Using these relation, Equation 1 can be rewritten as:
𝑖−2 𝑑𝑥 𝑖−1 𝑖+2 𝑑𝑥 𝑖+1 22
𝑑𝑥
expansion―at the left and right boundaries of the control volume―as you have learnt in
−𝑘𝐴 1(𝑑𝑇) +𝑘𝐴 1(𝑑𝑇) −2×h(𝑥)×𝐴 (𝑇 −𝑇 )=0
𝑐𝑖∞ (4) In the above relation the two derivatives 𝑑𝑇 can be approximated using the Taylor’s series
this course. For (𝑑𝑇)
𝑑𝑥 𝑖−1
2
a 2nd order accurate central approximation can be developed as follows: (5)
𝑑𝑇 𝑇−𝑇 ()=𝑖 𝑖−1
𝑑𝑥 𝑖−1 ∆𝑥 2
where ∆𝑥 is the distance between the centers of two adjacent control volumes. A similar relation can also be developed for (𝑑𝑇) 1:
𝑑𝑇 𝑇−𝑇 ()=𝑖+1 𝑖
𝑑𝑥 𝑖+1 ∆𝑥 2
(6)
𝑑𝑥 𝑖+2
Consequently, Equation 4 will be transformed into a system of equations which should be solved to obtain the values of temperature at every control volume.
3
Project Requirements: Online Report Submissions:
You should upload your report in the designated place on Quercus. You have to submit a single *.pdf file containing all parts.
1. First, using the Taylor’s series expansion show why the numerical approximation developed for (𝑑𝑇) 1 is central and 2nd order accurate.
𝑑𝑥 𝑖−2
Hint: Notice that (𝑑𝑇) 1 is the temperature derivative at point (𝑖 − 1) NOT at (𝑖). 𝑑𝑥 𝑖−2 2
Submission: Provide the complete details of your solution for this part along with some discussion to support your answer.
2. Substitute the derived relations for 𝑑𝑇 (Equations 5 and 6) back into Equation 4, factor 𝑑𝑥
out the terms containing 𝑇, and rearrange the obtained relation to reach the following form:
𝑒𝑇 +𝑓𝑇+𝑔𝑇 =𝑟 (7) 𝑖 𝑖−1 𝑖 𝑖 𝑖 𝑖+1 𝑖
which clearly results in a tridiagonal system of equations.
Submission: Show all the relations you obtain for 𝑒 , 𝑓 , 𝑔 , and 𝑟 . In each case include
the complete details of the derivation procedure.
3. Equation7canonlybeusedfortheintermediatecontrolvolumes(i.e.for𝑛=2,3,…,𝑛− 1). For the first and last control volumes the situation is slightly different. Using the figure below, develop an expression for the first control volume in the same manner as in Equation 1 by starting with energy balancing. In this procedure, you should finally find relations for 𝑓 , 𝑔 , 𝑒 , and 𝑟 .
111 1
𝑞′′ 𝑏
𝑞′′ 𝑐,1
𝑖𝑖𝑖 𝑖
Now, follow a similar procedure to find the relations for 𝑒 , 𝑓 , 𝑔 , and 𝑟 . 𝑛𝑛𝑛𝑛
𝑞′′ 𝑐,1
4
𝑞′′ 1→2
Submission:Showtheallrelationsyouobtainedfor𝑓,𝑔 ,𝑒 ,and𝑟 aswellas𝑒 ,𝑓, 1111𝑛𝑛
𝑔 , and 𝑟 . In each case include the complete details of the derivation procedure. 𝑛𝑛
4. In the final step you are required to develop a program in MATLAB. Using the relations you obtained for 𝑒 , 𝑓 , and 𝑔 (and the boundary conditions), form the tridiagonal matrix
𝑖𝑖𝑖
of coefficients. Also form the matrix of constants using the formula obtained for 𝑟. 𝑖
Subsequently, employ the Tridiagonal Matrix Algorithm (also known as the Thomas Algorithm) to solve the system of equations and obtain the values of 𝑇 for all the control
volumes.
Submission: Plot the temperature variations against 𝑥 for different number of control volumes (𝑛 = 5, 10, 20, 40, and 80) all in the same figure. The figures should be complete in terms of axis labels, units, legends, title, line style, etc. Comment on your results. How does, in your opinion, increasing the number of control volumes affect your results?
Also, report the values of temperature at the tip and base of the fin for different values of 𝑛. Which one of these values you think represents the true value? Explain your answer.
5. Assumethatthestreamofcoolingairissuddenlyturnedoff.Thiscausestheconvective heat transfer coefficient to fall down to 2 𝑊/𝑚2𝐾 on the whole surface of the fin. What would happen in this situation?
Submission: Plot the temperature variation resulted in this case together with the one you obtained in Part 4 in the same figure. Use only 𝑛 = 40 for both cases. Write down your observations based on the results.
Online Code Submissions:
You should upload your m-files in the designated place on Quercus. You have to submit
a single *.zip file containing all of your m-files. The name of the *.zip file should be: lastname_studentID_project1.zip
For example: einstein_9999999999_project1.zip
To test your program, the TAs will first declare 𝒏, an array containing the values of 𝑻𝒃, 𝑻∞, 𝒉𝒃, 𝒉𝒕, 𝒌 and another array containing 𝑳 and 𝒃. Then they will call your program in the Command Window. For example:
>> n = 20;
>> properties = [32000; 300; 18; 6; 202];
>> dimensions = [0.07; 0.01];
>> einstein_9999999999_project1(properties, dimensions, n)
You should define the name of your main function to make possible such an action in MATLAB. Once called, you program should give out an array containing the values of 𝑻, an array containing the values of 𝒙, and a figure where the temperature is plotted against 𝒙. It is strongly recommended that you check your program in the same
𝑖
5
manner to avoid any problems. It is your responsibility to make sure that your code works properly without any further adjustments by the TAs.
Remember not to capitalize your last name at any point. The commented header at the beginning of each file which describes your function must be completed.
Do NOT write the whole code in one function. Break your code into 4 functions. Calculate the values of 𝑔′ and 𝑟′ in 2 separate functions and develop the tridiagonal matrix solver in
another function as you did previously for your Assignment 2. The name of all your functions should be in the following form:
lastname_studentID_functionname
where the functionname can be chosen by you. The 4th function should be your main
function which will be called in the Command Window.
𝑖𝑖
6