Matlab Project on Numerically Solving Differential Equations

Project 3: Numerically solving a differential equation

Due: July 17th by the beginning of class

10% of total grade

Background: In order to more effectively dissipate heat, fins can be added to a hot surface. In

the field of heat transfer a fin is known as an ‘extended surface’. We want to determine the time

required for the temperature in the fin to reach steady-state. Then we want to optimize the length

of the fin to maximize the amount of heat transfer per unit of material.

The fin will be 2 cm wide (into the page), 0.4 cm thick, and the length will vary. At the base, the

fin is attached to a hot surface that is maintained at 100 C and the air around the fin is 25 C. The

fin is made of aluminum and has a thermal conductivity of 240 W/m-K, a specific heat of 900

J/kg-K and a density of 2700 kg/m3. Assume that the convection coefficient on the surface of the

fin is a constant of 25 W/m2-K.

To analyze this fin we will divide the fin up into N segments along the length. We will assume

the temperature profile along the width and thickness does not change. In fact, we will assume

each segment of the fin has a constant temperature throughout that segment.

After dividing the fin into N even segments (control volumes), we will now apply the 1st law of

Thermodynamics to each segment of the fin. In rate form (energy per unit time),

�̇�𝑖𝑛 − �̇�𝑜𝑢𝑡 = ∆�̇�𝑠𝑦𝑠 (1)

where heat is conducted in to each segment on the left, heat is conducted out of each segment on

the right, and heat is also lost due to convection along the outer surfaces (top, bottom, and sides).

�̇�𝑖𝑛,𝑖 = �̇�𝑐𝑜𝑛𝑑,𝑙𝑒𝑓𝑡,𝑖 = −𝑘𝐴 𝑑𝑇

𝑑𝑥 (2)

Hot

surface

(T=100 C)

0.4 cm

Length

�̇�𝑜𝑢𝑡,𝑖 = �̇�𝑐𝑜𝑛𝑑,𝑟𝑖𝑔ℎ𝑡,𝑖 + �̇�𝑐𝑜𝑛𝑣 = −𝑘𝐴 𝑑𝑇

𝑑𝑥 + ℎ𝐴𝑠(𝑇𝑖 − 𝑇𝑎) (3)

∆�̇�𝑠𝑦𝑠,𝑖 = 𝑑𝑈𝑖

𝑑𝑡 = 𝑚𝑐

𝑑𝑇𝑖

𝑑𝑡 (4)

However, the fin segments by the base and at the tip have to be treated differently. The first

segment has conduction from the left adding heat, but the conduction occurs over half the normal

length (from the base to the center of the first element), slightly modifying the energy equation.

At the tip, the convection occurs over a larger area (𝐴𝑠 + 𝐴) because the tip also loses heat to convection instead of conduction into another segment.

At i = 1, �̇�𝑖𝑛,1 = �̇�𝑐𝑜𝑛𝑑,𝑙𝑒𝑓𝑡,1 = −𝑘𝐴 𝑑𝑇

𝑑𝑥 = −𝑘𝐴

𝑇1−𝑇𝑏𝑎𝑠𝑒

0.5𝑑𝑥 (5)

At i = N, �̇�𝑜𝑢𝑡,𝑁 = �̇�𝑐𝑜𝑛𝑣 = ℎ(𝐴 + 𝐴𝑠)(𝑇𝑖 − 𝑇𝑎) (6)

Combining the above equations at the non-end segments yields:

−𝑘𝐴 ( 𝑑𝑇

𝑑𝑥 )

𝑙𝑒𝑓𝑡 − (−𝑘𝐴 (

𝑑𝑇

𝑑𝑥 )

𝑟𝑖𝑔ℎ𝑡 + ℎ𝐴𝑠(𝑇𝑖 − 𝑇𝑎)) = 𝑚𝑐

𝑑𝑇𝑖

𝑑𝑡 (7)

Note that 𝑚 = 𝜌𝑉 = 𝜌𝐴𝑑𝑥 and 𝑑𝑥 = 𝐿 𝑁⁄ assuming N even segments. Substituting into

equation 7 yields:

−𝑘𝐴 ( 𝑑𝑇

𝑑𝑥 )

𝑙𝑒𝑓𝑡 − (−𝑘𝐴 (

𝑑𝑇

𝑑𝑥 )

𝑟𝑖𝑔ℎ𝑡 + ℎ𝐴𝑠(𝑇𝑖 − 𝑇𝑎)) = 𝜌𝑐𝐴𝑑𝑥

𝑑𝑇𝑖

𝑑𝑡 (8)

Dividing through by kAdx and rearranging yields:

( 𝑑𝑇

𝑑𝑥 ) 𝑟𝑖𝑔ℎ𝑡

−( 𝑑𝑇

𝑑𝑥 ) 𝑙𝑒𝑓𝑡

𝑑𝑥 −

ℎ𝐴𝑠

𝑘𝐴𝑑𝑥 (𝑇𝑖 − 𝑇𝑎) =

𝜌𝑐

𝑘

𝑑𝑇𝑖

𝑑𝑡 (9)

But notice that the first term is just the central difference approximation of 𝑑𝑇

𝑑𝑥 , so the first term

becomes:

( 𝑑𝑇

𝑑𝑥 ) 𝑟𝑖𝑔ℎ𝑡

−( 𝑑𝑇

𝑑𝑥 ) 𝑙𝑒𝑓𝑡

𝑑𝑥 ≈ (

𝑑2𝑇𝑖

𝑑𝑥2 ) (10)

Now to simplify, note that 𝐴𝑠 = 𝑃𝑑𝑥, where P is the perimeter of the cross section, substitute the

thermal diffusivity ∝= 𝑘

𝜌𝑐 , and let 𝑏 =

ℎ𝑃

𝑘𝐴 and equation 9 becomes:

𝑑2𝑇𝑖

𝑑𝑥2 − 𝑏(𝑇𝑖 − 𝑇𝑎) =

1

∝

𝑑𝑇𝑖

𝑑𝑡 (11)

Now we need to apply finite difference equations to the derivative terms. Notice that

temperature is both a function of position and time, so we will use two subscripts. The first

subscript will indicate the position and the second subscript the time. We will use FTCS, or

‘Forward Time Central Space’ in our finite difference equations. If you remember, the central

difference equation is second order accurate and the forward difference is first order accurate.

So we will have second order accuracy spatially and first order accuracy temporally. There are

two different approaches to handling the time dimension that will be discussed – implicit and

explicit methods.

Explicit method:

Applying the finite difference approximations, equation 11 becomes:

𝑇𝑖+1,𝑡−2𝑇𝑖,𝑡+𝑇𝑖−1,𝑡

(∆𝑥)2 − 𝑏(𝑇𝑖,𝑡 − 𝑇𝑎) =

1

∝

𝑇𝑖,𝑡+1−𝑇𝑖,𝑡

∆𝑡 (12)

Notice that the time subscript has a “+1” in only one location. This makes solving for the future

temperatures very straightforward and simple. Solve equation 12 for the one “+1” future

temperature:

𝑇𝑖,𝑡+1 = 𝑇𝑖,𝑡+∝ ∆𝑡 ( 𝑇𝑖+1,𝑡−2𝑇𝑖,𝑡+𝑇𝑖−1,𝑡

(∆𝑥)2 − 𝑏(𝑇𝑖,𝑡 − 𝑇𝑎)) (13)

Equation 13 provides a way to calculate the interior segment temperatures, but different

equations are needed at the boundaries:

at i=N

Modifying equation 8 to account for the differences in the final segment yields:

−𝑘𝐴 ( 𝑑𝑇

𝑑𝑥 )

𝑙𝑒𝑓𝑡 − (0 + ℎ(𝐴𝑠 + 𝐴)(𝑇𝑖 − 𝑇𝑎)) = 𝜌𝑐𝐴𝑑𝑥

𝑑𝑇𝑖

𝑑𝑡 (14)

Rearranging and applying the central difference equation on the spatial derivative yields:

1

∝

𝑑𝑇𝑁

𝑑𝑡 = −

1

𝑑𝑥 (

𝑑𝑇

𝑑𝑥 )

𝑙𝑒𝑓𝑡 − (

ℎ(𝑃𝑑𝑥+𝐴)

𝑘𝐴𝑑𝑥 (𝑇𝑁 − 𝑇𝑎)) = −

1

𝑑𝑥 (

𝑇𝑁−𝑇𝑁−1

𝑑𝑥 ) − ((𝑏 +

ℎ

𝑘𝑑𝑥 ) (𝑇𝑁 − 𝑇𝑎)) (15)

Using the explicit forward time method:

1

∝

𝑇𝑁,𝑡+1−𝑇𝑁,𝑡

∆𝑡 = −

1

∆𝑥 (

𝑇𝑁,𝑡−𝑇𝑁−1,𝑡

∆𝑥 ) − ((𝑏 +

ℎ

𝑘∆𝑥 ) (𝑇𝑁,𝑡 − 𝑇𝑎)) (16)

Rearranging to solve for the “+1” time term, the final segment is calculated using the following:

𝑇𝑁,𝑡+1 = 𝑇𝑁,𝑡+∝ ∆𝑡 (− 1

∆𝑥 (

𝑇𝑁,𝑡−𝑇𝑁−1,𝑡

∆𝑥 ) − ((𝑏 +

ℎ

𝑘∆𝑥 ) (𝑇𝑁,𝑡 − 𝑇𝑎))) (17)

at i=1,

Applying the explicit forward time and central space approximations to equation 8 to account for

the differences in the initial segment yields:

−𝑘𝐴 𝑇1,𝑡−𝑇𝑏𝑎𝑠𝑒

0.5𝑑𝑥 − (−𝑘𝐴

𝑇2,𝑡−𝑇1,𝑡

𝑑𝑥 + ℎ𝐴𝑠(𝑇1,𝑡 − 𝑇𝑎)) = 𝜌𝑐𝐴𝑑𝑥

𝑇1,𝑡+1−𝑇1,𝑡

𝑑𝑡 (18)

Rearranging, and solving for the “+1” time term, the initial segment is calculated using the

following:

𝑇1,𝑡+1 = 𝑇1,𝑡+∝ ∆𝑡 ( 1

(∆𝑥)2 (𝑇2,𝑡 − 3𝑇1,𝑡 + 2𝑇𝑏𝑎𝑠𝑒) −

ℎ𝑃

𝑘𝐴 (𝑇1,𝑡 − 𝑇𝑎)) (19a)

𝑇1,𝑡+1 = 𝑇1,𝑡 + (𝛽(𝑇2,𝑡 − 3𝑇1,𝑡 + 2𝑇𝑏𝑎𝑠𝑒) − 𝑏 ∝ ∆𝑡(𝑇1,𝑡 − 𝑇𝑎)) (19b)

Where (∝ ∆𝑡

(∆𝑥)2 ) = 𝛽 is defined for algebraic simplicity (and for something else on the

next page).

These equations can easily be implemented, provided that temperature values are known at t=0

(initial conditions). We can assume the entire fin starts at 25 C at t=0. An example calculation is

given below:

Use equation 19b for the first element:

Where ∝ = 240𝑊 𝑚𝐾⁄

2700 𝑘𝑔

𝑚3 ⁄ ∗900

𝐽 𝑘𝑔𝐾⁄

= 0.000099 𝑚 2

𝑠⁄

𝑎𝑛𝑑 𝑏 = (25 𝑊

𝑚2𝐾 ⁄ ) (2 ∗ .02𝑚 + 2 ∗ .004𝑚)

(240 𝑊 𝑚𝐾⁄ ) (. 02𝑚 ∗ .004𝑚)

= 62.5

𝑚2

𝑎𝑛𝑑 𝛽 = ∝ ∆𝑡

(∆𝑥)2 = 0.000099 𝑚

2

𝑠⁄ ∗ 0.01𝑠

( 0.5𝑚 40

) 2 = 0.006321

𝑇1,0+1 = 𝑇1,0 + (𝛽(𝑇2,0 − 3𝑇1,0 + 2𝑇𝑏𝑎𝑠𝑒) − 𝑏 ∝ ∆𝑡(𝑇1,0 − 𝑇𝑎))

= 25

+ (0.006321 ∗ (25 − 3 ∗ 25 + 2 ∗ 100) − 62.5 ∗ .000099 ∗ .01𝑠(25 − 25))

= 25.9482 𝐶

Then, using equation 13 at i=2,

𝑇2,0+1 = 𝑇2,0+∝ ∆𝑡 ( 𝑇2+1,0 − 2𝑇2,0 + 𝑇2−1,0

(∆𝑥)2 − 𝑏(𝑇2,0 − 𝑇𝑎))

= 𝛽(𝑇2+1,0 − 2𝑇2,0 + 𝑇2−1,0) − 𝑏 ∝ ∆𝑡(𝑇2,0 − 𝑇𝑎)

= 0.006321(25 − 2 ∗ 25 + 25) − 𝑏 ∝ ∆𝑡(25 − 25) = 25 𝐶

So after 1 time step of 0.01s, the temperature at the first segment is about 1 degree higher and the

temperature at the second (and 3rd and 4th and so on) segment is unchanged. Once new values

for all segments have been obtained, then a new time step can begin back at i=1.

Warning: this method can become unstable if the time step is too large or the dx is too small.

In order to maintain stability, the following criteria must be met:

(∝ ∆𝑡

(∆𝑥)2 ) = 𝛽 < 0.5

Implicit method:

The implicit method uses similar equations to the explicit method. For the central grid points

(everything except the first and last), equation 12 becomes the following for the implicit method:

𝑇𝑖+1,𝑡+1−2𝑇𝑖,𝑡+1+𝑇𝑖−1,𝑡+1

(∆𝑥)2 − 𝑏(𝑇𝑖,𝑡+1 − 𝑇𝑎) =

1

∝

𝑇𝑖,𝑡+1−𝑇𝑖,𝑡

∆𝑡 (20)

Notice now that the “future” temperature values appear in several places, not just one.

Rearranging equation 20 to put all “+1” terms on the LHS,

𝑇𝑖+1,𝑡+1−2𝑇𝑖,𝑡+1+𝑇𝑖−1,𝑡+1

(∆𝑥)2 − 𝑏(𝑇𝑖,𝑡+1) −

1

∝

𝑇𝑖,𝑡+1

∆𝑡 = −𝑏𝑇𝑎 −

𝑇𝑖,𝑡

∝∆𝑡 (21a)

(−𝛽)𝑇𝑖+1,𝑡+1 + (2𝛽 + 𝑏 ∝ ∆𝑡 + 1)𝑇𝑖,𝑡+1 + (−𝛽)𝑇𝑖−1,𝑡+1 = 𝑏 ∝ ∆𝑡𝑇𝑎 + 𝑇𝑖,𝑡 (21b)

Equation 21b has 3 unknowns so it cannot directly be solved for the future temperature values.

However, when the equation is applied to each interior segment, there will be N unknowns and

N-2 equations. The final two equations come from the boundary conditions at i = N and i = 1.

From equation 16, for i = N,

1

∝

𝑇𝑁,𝑡+1 − 𝑇𝑁,𝑡 ∆𝑡

= − 1

∆𝑥 ( 𝑇𝑁,𝑡+1 − 𝑇𝑁−1,𝑡+1

∆𝑥 ) − ((𝑏 +

ℎ

𝑘∆𝑥 ) (𝑇𝑁,𝑡+1 − 𝑇𝑎))

Rearranging to put the “+1” terms on the LHS

1

∝

𝑇𝑁,𝑡+1

∆𝑡 + (

𝑇𝑁,𝑡+1−𝑇𝑁−1,𝑡+1

(∆𝑥)2 ) + (𝑏 +

ℎ

𝑘∆𝑥 ) 𝑇𝑁,𝑡+1 =

𝑇𝑁,𝑡

∝∆𝑡 + (𝑏 +

ℎ

𝑘∆𝑥 ) 𝑇𝑎 (22a)

(𝛽 + 1+∝ ∆𝑡 (𝑏 + ℎ

𝑘∆𝑥 )) 𝑇𝑁,𝑡+1 − (𝛽)𝑇𝑁−1,𝑡+1 = +𝑇𝑁,𝑡+∝ ∆𝑡 (𝑏 +

ℎ

𝑘∆𝑥 ) 𝑇𝑎 (22b)

From equation 18, for i = 1,

−𝑘𝐴 𝑇1,𝑡+1−𝑇𝑏𝑎𝑠𝑒

0.5𝑑𝑥 − (−𝑘𝐴

𝑇2,𝑡+1−𝑇1,𝑡+1

𝑑𝑥 + ℎ𝐴𝑠(𝑇1,𝑡+1 − 𝑇𝑎)) = 𝜌𝑐𝐴𝑑𝑥

𝑇1,𝑡+1−𝑇1,𝑡

𝑑𝑡 (23a)

Rearranging to put the “+1” terms on the LHS

(1 + 3𝛽 + 𝑏 ∝ ∆𝑡)𝑇1,𝑡+1 − (𝛽)𝑇2,𝑡+1 = 𝑇1,𝑡 + 2𝛽𝑇𝑏𝑎𝑠𝑒 + 𝑏 ∝ ∆𝑡𝑇𝑎 (23b)

With N equations and N unknowns, we can now use a matrix to solve the set of simultaneous

linear equations to get the temperature values at the new time step! For example, with N = 5, the

first equation comes equation 23b. At t=0, let all temperatures be 25 C.

The second equation comes from letting i = 2 in equation 21b.

(𝛽)𝑇2+1,0+1 + (−2𝛽 − 𝑏 ∝ ∆𝑡 − 1)𝑇2,0+1 + (𝛽)𝑇2−1,0+1 = −𝑏 ∝ ∆𝑡𝑇𝑎 − 𝑇2,0

When i = 3 and 4, equation 21b is applied. When i = 5 (i = N), equation 22b is applied.

[ (1 + 3𝛽 + 𝑏 ∝ ∆𝑡) −𝛽 0

−𝛽 (1 + 2𝛽 + 𝑏 ∝ ∆𝑡) −𝛽

0 −𝛽 (1 + 2𝛽 + 𝑏 ∝ ∆𝑡)

0 0

−𝛽

0 0 0

0 0 −𝛽 (1 + 2𝛽 + 𝑏 ∝ ∆𝑡) −𝛽

0 0 0 −𝛽 (1 + 𝛽+∝ ∆𝑡 (𝑏 + ℎ

𝑘∆𝑥 ))

]

[ 𝑇1,𝑡+1 𝑇2,𝑡+1 𝑇3,𝑡+1 𝑇4,𝑡+1 𝑇5,𝑡+1]

=

[ 𝑇1,0 + (2𝛽)𝑇𝑏𝑎𝑠𝑒 + 𝑏 ∝ ∆𝑡𝑇𝑎

𝑏 ∝ ∆𝑡𝑇𝑎 + 𝑇2,0 𝑏 ∝ ∆𝑡𝑇𝑎 + 𝑇3,0 𝑏 ∝ ∆𝑡𝑇𝑎 + 𝑇4,0

𝑇5,0+∝ ∆𝑡 (𝑏 + ℎ

𝑘∆𝑥 ) 𝑇𝑎 ]

A tri-diagonal matrix is produced. This could be solved using Gaussian Elimination or LU

decomposition, but it can also be solved using the inverse function or something similar in

Matlab. It will need to be solved for each time step, so it will need to be solved many, many

times. Efficiency is important.

Notice that the coefficient matrix will not change as the temperature values change. But the

RHS vector will change as those values are dependent on the temperatures from the previous

time step, so the RHS vector will need to be generated fresh on each time step.

Aside from being more physically realistic, the implicit method is inherently stable, so there is

not the same restriction on the dt or dx values.

Determining total heat dissipated by the fin at steady state:

To determine the total heat dissipated we will add up the convection off of every exterior surface

of every segment.

Using Newton’s Law of cooling:

�̇�𝑡𝑜𝑡𝑎𝑙 = ∑ ℎ𝑃Δ𝑥(𝑇𝑖 − 𝑇𝑎) 𝑁 𝑖=1 (22)

Procedure (use Matlab to do the following):

1. (60 points) Start with N=40 segments (you can experiment with this later, but do not go

much below this number) and a length of 0.5m. Solve for the temperature as a function

of time using the explicit method and determine the approximate time required to reach

steady-state. This can be determined by finding the absolute relative error from one time

step to the next for each segment, finding the maximum, and quitting when the maximum

change is less than, say 0.01% (provided your time step is not too small) for any segment.

Generate and display 3 graphs, one after 5 time steps, one roughly in the middle, and one

at steady state. Also, provide the approximate time to reach steady state.

2. (20 points) Use at least 40 segments and a length of 0.5 m and solve for the temperature

as a function of time using the implicit method. Generate and display 3 graphs, one after

5 time steps, one roughly in the middle, and one at steady state. Also, provide the

approximate time to reach steady state.

3. (10 points) Using either method, find the total heat being dissipated by the fin at steady

state.

4. (10 points) Using either method, vary the length and find the optimum value using the

following information: each watt of power dissipated by the fin creates a profit of $1.68

and each kg of aluminum, machined/rolled/formed into the fin has a cost of $3.27. Find

the length at which adding additional material ceases to become cost effective. Find the

optimum length accurate to the millimeter. You may use guess and check or something

more sophisticated.

5. (10 bonus points) Create an animation of your fin temperature profile as a function of

time. Save your animation as an avi file and come show me in my office, along with the

supporting code.

6. Explain all of your steps using comments in Matlab.

What to submit (all printed, stapled, and ready by start of class on the due date):

 Signed affidavit sheet

 Published mfile in html format

o Follow the sample project format including cell formatting, published html

format, commenting, etc

http://www.eng.usf.edu/~kaw/class/EML3041/homework/sample_experimental.html

 Upload your mfile to Canvas. If you created additional function mfiles, upload those as

well.