I have a project in Computational Methods class to solve differential equations numerically using matlab.
The project is attached.
Matlab Project on Numerically Solving Differential Equations
Project 3 : Numerically s olving a differential equation Due: July 17 th 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 know n 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 optim ize the length of the fin to maximize the amount of heat transfer per unit of material. The fin will be 2 cm wid e (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 specif ic heat of 900 J/kg -K and a density of 2700 kg/m 3. Assume that the convection coefficient on the surface of the fin is a constant of 25 W/m 2-K. To analyze this fin we will divide the fin up into N segments along the length. We will assume the te mperature 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 1 st 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 modif ying the energy equation. At the tip, th e 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 s ubscript 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, e qu ation 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 s traightforward 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 )− ((+ ℎ )(− )) (1 5) 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 ) Re arranging, and solving for the “+1 ” time term, the initial segment is calculated using the following: 1,+1= 1,+∝ ∆( 1 (∆)2(2,− 31,+ 2 )− ℎ (1,− )) (19a) 1,+1= 1,+ ((2,− 31,+ 2 )− ∝ ∆(1,− )) (19b) W here (∝ ∆ (∆)2)= is defined for algebraic simplicity (and for something else on the next page). These equations can easily be implemented, provided that temperature va lues are known at t=0 (initial conditions). We can assume the e ntire fin starts at 25 C at t=0 . An ex ample calculation is given below : Use equation 19b for the first element: W here ∝ = 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− 31,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− 22,0+ 2−1,0 (∆)2 − (2,0− )) = (2+1,0− 22,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 3 rd and 4 th 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 metho d. 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 ju st one. Rearranging equation 20 to put all “+1” terms on the LHS, +1,+1−2,+1+−1,+1 (∆)2 − (,+1)− 1 ∝ ,+1 ∆ = − − , ∝∆ (21 a) (−)+1,+1+ (2 + ∝ ∆+ 1),+1+ (−)−1,+1= ∝ ∆+ , (21 b) Equation 21 b has 3 unknowns so it cannot directly be solved for the future temperature values. However, when the equation is ap plied to each int erior segment, there will be N unknowns and N-2 equations. The final two equation s come from the boundary condition s 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= , ∝∆+ (+ ℎ ∆) (22 a) ( + 1+∝ ∆(+ ℎ ∆)),+1− ()−1,+1= +,+∝ ∆(+ ℎ ∆) (22 b) 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 + ∝ ∆ (23 b) 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 firs t equation comes equation 23b . At t=0, let all temperatures be 25 C . The second equation comes from letting i = 2 in equation 21 b. ()2+1,0+1+ (−2 − ∝ ∆− 1)2,0+1+ ()2−1,0+1= − ∝ ∆− 2,0 When i = 3 and 4 , equa tion 21 b is applied. When i = 5 (i = N) , equation 2 2b 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 i s 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 fre sh 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 ex terior surface of every segment. Using Newton ’s Law of cooling : ̇ = ∑ ℎΔ(− ) =1 (22) Procedure (use Matlab to do the following) : 1. (6 0 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 s egment, 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. (1 0 points) Using either method, vary the len gth and find the optimum value using the following information: each wat t of power dissipated by the fin creates a profit of $1 .68 and each kg of aluminum, machined/rolled/form ed into the fin has a cost of $3 .27. Find the length at which adding additiona l 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. Sa ve your animation as a n 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/E ML3041/homework/sample_experimental.html Upload your mfile to Canvas . If you created additional function mfiles, upload those as well.
