Class Note for CHEE 502 with Professor Saez at UA
Class Note for CHEE 502 with Professor Saez at UA
Popular in Course
Popular in Department
This 23 page Class Notes was uploaded by an elite notetaker on Friday February 6, 2015. The Class Notes belongs to a course at University of Arizona taught by a professor in Fall. Since its upload, it has received 15 views.
Reviews for Class Note for CHEE 502 with Professor Saez at UA
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 02/06/15
MODEL 3 Heat Transfer in Circular Pipes Mathematical Aspect Numerical solution of PDEs initial value problems Heat transfer in pipes Consider laminar developed ow in a pipe of radius R The uid is at a uniform temperature To for zlt0 Figure l and at z0 it enters a region at which the temperature of the wall changes to TW T0 We would like to study how the temperature of the uid evolves with z To analyze this problem we will consider that the process is at steady state that viscous dissipation is negligible and that physical properties of the uid are uniform In addition we will assume that the temperature eld is symmetric in the angular direction of a cylindrical coordinate system so that we seek a solution for the temperature eld of the form TTrz TTW zgt0 velocity l Zlt0 I Z pro le Figure 1 A uid owing at uniform temperature To encounters a sudden change in wall temperature to a value TW at z0 First consider the expected shape for the radial temperature pro les along the pipe for the case TWgtT0 Figure 2 The pro le starts uniform at z0 l but the wall temperature increases suddenly At short distances 2 the uid at the core still will be at T0 but sharp temperature gradients occur near the wall Further downstream 3 even the uid at the center of the pipe will start to heat up and the temperature pro le will span the whole cross section of the pipe At long distances 4 the temperature will tend to approach the wall temperature everywhere TW 7 l I T TO To To iTW z0 i 3 l I 1 2 3 4 Figure 2 Evolution of radial temperature pro les along a pipe see text for explanation For steady incompressible developed ow the only nonzero component of the velocity vector is r2 vZ2ltvZgt 1 l where ltvzgt is the crosssectional averaged velocity related to the volumetric ow rate by ltvzgtQTER2 Under the conditions listed the thermal energy equation in cylindrical coordinates simpli es to 3T 1 2 3T BZT WP 32 kr ar r ar 322 2 where p and cp are density and constantpressure heat capacity respectively and k is the uids thermal conductivity The boundary conditions are 1 It will be assumed that the uid reaches the section at which the wall temperature changes 20 with uniform temperature TT0 20 3 2 For zgt0 TTW rR 4 3 The temperature is symmetric around the pipe aXis Or0 5 Each term in the PDE 2 can be identi ed with an energy transport mechanism the term on the lefthand side represents the aXial 2 transport of thermal energy by convection the rst term on the righthand side represents heat conduction in the radial direction note that this mechanism will be directly responsible for the temperature change in the uid since it will account for the heat exchange between the uid and the wall and the last term represents aXial heat conduction The relative magnitude of the two transport mechanisms in the aXial direction can be compared by performing an order of magnitude analysis Let L be a characteristic length in the z direction over which signi cant temperature changes AT occur We can state 3T AT pcpv20pcp ltvZ gtT 6 aZT AT k32 O k F 7 We will explore the conditions for which convection dominates over conduction in the axial direction This will happen when AT c ltv gt7gtgtk7 8 p p z L L Simplifying and rearranging leads to pop lt Vz gt L gtgt 1 convectlon dominates over conduction 1n 2 9 The dimensionless number that appears in this constraint compares the relative magnitude of convection and diffusion in the axial direction In general L is not known so a dimensionless number is defined using the pipe diameter known as characteristic length This dimensionless number is called the Peclet number pcpltVzgt2RltVZgt2R Fe 10 k a where Fldpcp is the thermal dispersivity The constraint 10 can be expressed as L Pei l l 1 2R If this constraint is satis ed we can neglect axial conduction in the thermal energy equation This is the most interesting problem from the point of view of convection since uid motion is responsible for energy transport along the pipe We will consider in our solution that constraint 1 l is satisfied The thermal energy equation 2 combined with the velocity profile 1 simplifies to 2 2ltvZ gtl1 Llal33ralj 12 K RZJBZ rBr Br Note that the absence of the axial conduction term has transformed this equation into a first order differential equation in 2 Consequently only one boundary condition is needed in 2 equation 3 The problem to solve consists of the partial differential equation 12 subject to boundary conditions 3 to 5 To reduce the number of parameters in the formulation we will express the problem in dimensionless form Dimensionless temperature and radial coordinate are defined by 87TT0 13 TW T0 1 7 14 T R Substituting into equation 12 we get after manipulations and use of equation 10 PeR1n2aa n 15 To eliminate all parameters from the PDE we select the following dimensionless aXial distance 2 7 16 E PeR to get 38 l 3 38 1 112 737617 17 BE T 3n 3n The boundary conditions 3 to 5 transform to 80 Efo 18 38 7 0 T10 19 3n 81 111 20 If a solution for this problem is attempted using separation of variables the problem in the T direction SturmLiouville problem does not have a known analytical solution What complicates the formulation is the variable coefficient lnz in the PDE 17 Hence a solution to this problem should be found numerically Numerical algorithms used to solve PDEs are strongly in uenced by the type of differential equation specifically the order of the derivatives with respect to the independent variables might affect the type of algorithm to be used The present problem is a boundaryvalue problem with respect to T since the PDE contains second order derivatives of the dependent variable with respect to T and the two boundary conditions equations 19 and 20 are imposed on two different extreme points of the domain However the problem is an initialvalue problem with respect to E since only first derivatives with respect to this variable are present in the PDE Initialvalue problems IVPs are typically easier to treat than boundaryvalue problems so that the IVP nature of the present problem will be the most important characteristic in establishing the numerical method to be used The approach used here to nd a numerical solution for the problem will be based on the general context of nite differences Introduction to nite difference discretization for PDEs The method of nite differences is based on de ning a certain number of points in the domain called nodes or discretization points that taken together form a discretization mesh or grid After these points have been de ned the following steps are applied 1 The derivatives that appear in the PDE are expressed in terms of values of the dependent variable at discretization points using approximate relations derived from truncated Taylor series 2 The differential equation is evaluated at the discretization points substituting the derivatives by the expressions mentioned in 1 This leads to a system of equations having as many equations as there are discretization points 3 The boundary conditions are used to represent boundary nodes 4 The system of equations is solved to nd approximate values of the dependent variable at all discretization points The method of nite differences was used in Model 1 to solve ODEs We will sta1t our analysis following some of the derivations presented before First consider a function of a single variable yyx and a discretization of the independent variable in an interval anSb formed by equidistant points de ned by xiaih 21 where b a h 22 Nl This scheme de nes Nl intervals and N2 discretization points with x0a xN1b and N interior points To nd an expression to approximate the rst derivative of y at point xi consider a thirdorder Taylor series for yx around xxi YX YX1 Y39X1X X1 lYquotX1X X02 1 1 2 23 EY 39X1X Xi3 ZyiVCX X04 where C is between x and x If we evaluate this equation at xxi1 and recall that h xiH xi we get 1 1 1 YX11 YX1 Y39X1hEYquotX1h2 Y 39Xih3 EYIVC1h4 24 with xiSCiniH We can solve this equation for the rst derivative to obtain 2 3 y39mW gywxi yquot39xi yw 25 The first of the last three terms of this equation is dominant for small h so that this equation can be written as yogi YXi1h YXi 001 26 Truncation of this equation leads to an approximation of the first derivative yogi YX1lYX1 27 h This is a firstorder forward difference formula because the derivative at a point depends on the value of the dependent variable at that point and at the next point To approximate the second derivative we can see that equation 24 will not be enough due to the presence of the first derivative Therefore we need an additional equation which will be equation 23 evaluated at xxi1xi h YXi 1 YXi Y39Xihyuxih2 Y39HXih3 yivCi 1h4 28 with xHSCHSxi Addition of equations 24 and 28 eliminates the first and third derivatives Solving the resulting equation for the second derivative yields yHXi W lyivCi1yiVCi 1lh2 29 Or yuXi YXil2y1fiYXi l Oh2 30 Truncation of this equation leads to an approximation of the second derivative yuXiz YX1l 2331YX1 l 31 This is a secondorder centered difference formula for the second derivative since it contains values of the dependent variable in the two adjacent points to the point where the derivative is evaluated A centered difference formula can be obtained for the first derivative too To do this we subtract equation 28 from equation 24 so that the second derivative is eliminated The result is h2 h3 y39oq W ywxo 2X 4y1VCi 1 YIVC11l 32 01 yogi YXi12 hYXi 1 0012 33 So that the centereddifference formula YXi1 YXi 1 34 mo x 2h has a truncation error of 0h2 An alternative formula for the first derivative can be obtained directly from equation 28 yogi YXi hYXi 1Oh 35 which leads to the approximation YXi YXi 1 36 Y39Xi h This is the firstorder backward difference formula for the first derivative The formulas derived in this section will be used in the nitedifference discretization of the problem given by equations 17 to 20 First we will look at a scheme that is based on the forwarddifference formula 27 for the first derivative The forward difference method The domain in which the PDE 17 will be solved is 031131 035 Since the domain in Ehas no defined limit we will pick a step size h in Ethat defines unlimited discretization points of the form J ih i012 37 On the other hand in the T direction we will use a step size g given by see equation 22 1 m 38 g so that discretization points in this direction are defined by nj gjj0l2Nl 39 A discretization point in the grid thus created is given by the coordinates nj i The discretization grid is shown in Figure 3 boundary nodes i0 i2 13 i4 Q G C quot C jNl p1pe wall 111 I I I b o o o c j3 l I I I I g I interior node I I I 039 o o o o j1 I I n 1 e e e e j0 pipe center 110 0310 a T Figure 3 Example of a nite difference discretization grid for the problem of heat transfer in a pipe with N3 Nodes for which j0 pipe center jNl pipe wall and i0 O are boundary nodes The rest are interior nodes To discretize the differential equation we need approximate expressions for derivatives with respect to E and 11 First consider the first derivative with respect to To approximate it we use the forwarddifference formula 27 33 1 1 gimp z W 40 a 3 where we have evaluated the derivative at a fixed radial position m From now on to designate the approximate nodal values of the dimensionless temperature we will use the following simplif1ed notation gianj ij 41 Equation 40 can be expressed as 8 1 armps 42 3 3E To discretize the righthand side of the differential equation 17 first we expand the derivative 1 E 38 1 38 328 if n7 7772 43 n 3n 3n n 3n 3n The second derivative is approximated by the centereddifference formula 31 328 5 121 Warniij 44 Since this formula has a truncation error of 0g2 we will use the centereddifference formula 34 to approximate the first derivative in equation 43 to keep the same level of error in all approximations 3 3n 5 11 inm l IJ 45 2g Now we use all these results equations 42 to 45 into the evaluation of PDE 17 at a nodal point 51m to get i1j 45g i 81j1 81j 1 ij1 Z g ij 1 h m 2g g2 l nJZ 46 Note that all the values of the dependent variable in this equation except one are evaluated at This means that this equation can be used to calculate the dimensionless temperature at EH1 from values at the previous step in axial position 54 which leads to an explicit numerical scheme similar to Euler39s method as seen in Model 1 Now we solve equation 46 for 91 to get A h l l A 2h A 817 77 11 1 gln g 2m 1 g21nJ2 1 i01j1 to N 47 h 1 1 A 72 quot GLJl g1nJ g 2m This equation applies to interior nodes only Boundary nodes are governed by the discretized form of the boundary conditions Consider rst boundary condition 20 This is the only non homogeneous condition and it will be therefore the condition that generates a nontrivial temperature pro le Application of this boundary condition yields 61M 1 io12 48 note that nN1l The other boundary condition in the radial direction equation 19 specifies a value for a radial derivative We will use the forwarddifference formula 27 to approximate the derivative to get 8 H 511 510 an 0 g 49 so that boundary condition 19 yields 30 on i012 50 Finally boundary condition 18 leads to oj0jltoN 51 An algorithm used to solve the problem with the equations stated above can be formulated as follows 1 The calculation starts at EEO i0 imposing Boundary condition 51 Boundary condition 48 for i0 At this point all values of the dependent variable are known at the line of nodes i0 Figure 3 2 For each line of nodes i apply the following equations sequentially until the desired final value for the E coordinate is reached Equation 47 starting at i0 Equation 48 Equation 50 When this algorithm was applied to solve the problem it was found that for certain values of g unphysical oscillatory solutions were obtained unless h was small enough For example Figure 4 shows how a clearly unstable oscillatory solution is obtained for h0000202 at n075 while a stable solution is found for a slightly lower value h0000200 It should be pointed out that unstable profiles are found for the latter h at other values of T An interesting aspect of the solution is that once h is low enough that the solution becomes stable further reductions in h do not change the solution appreciably which indicates that the exact solution has been found within the scale of the plot in Figure 4 10 09 7 075 08 7 005 07 7 06 7 05 h0000202 04 h0000200 03 r 02 r 01 0 001 002 003 004 005 006 E Figure 4 Calculated axial dimensionless temperature pro le at n075 for two values of h Further reduction of h below 0000200 yields the same solution at the scale of the gure Figure 5 shows calculated radial temperature profiles for relatively low values of Note that as seen in Figure 4 the centerline temperature EEO will continue to increase until the asymptotic pro le 851 for all T is obtained as i wo eventually the uid will reach the wall temperature Another important physical parameter of this problem is the heat ux at the wall since the surroundings must provide the thermal energy necessary to keep the wall temperature uniform that is if the wall is hotter than the uid The heat ux in the radial direction us given by Fourier s law 1r 3 52 To express this equation in dimensionless form we make use of equations 13 and 14 to get Rqr EB 53 kTwTo 3T1 The lefthand side of this equation is a dimensionless heat ux Let the value of this quantity at rR be ll R qurR 54 kTW T0 1 8 09 7 005 h00001 0 02 04 06 08 1 Tl Figure 5 Calculated radial dimensionless temperature pro les for various aXial positions Curves are shown at 001 intervals in Lower values of h would give the same solution at the scale of the plot Note how the pro les progressively penetrate into the uid until the centerline temperature 110 starts to rise Then from equation 53 we have 88 07 55 an 1 After the numerical solution for the dimensionless temperature has been found this dimensionless heat ux can be calculated by approximating the derivative using the backward difference formula 36 which yields 1N1 9m 21 9m g g 61 56 Note that the Nth node is adjacent to the pipe wall nNl g Figure 6 shows the dimensionless heat ux at the wall calculated from this equation using the numerical solution Note that the ux 12 decreases as E increases due to the reduction in heat transfer driving force as the temperature of the uid warms up The ux will asymptotically approach zero as i wo and the temperature of the uid becomes uniform and equal to the wall temperature 14 Oquot h00001 005 12 10 0 l l l l 0 0005 001 0015 002 0025 E Figure 6 Dimensionless heat ux at the wall as a function of axial position We have shown how the forwarddifference method can be used to find a numerical solution of a PDE that forms an IVP with respect to one of the independent variables E in this case However the instability of the solution requires relatively low values of h to be used to find a meaningful solution which involves a relatively high computational effort Important questions at this point would be how to find out if a specific algorithm may be unstable and how to formulate an algorithm to reduce or eliminate the instability problem These issues will be addressed in the following sections Stability analysis The instability observed in the solution presented above is typical of the forwarddifference method and other explicit methods Here we will explore how to determine if a method might lead to unstable solutions taking the forwarddifference method as a working example We will start with the most important equation of the algorithm equation 47 which can be rewritten as follows i1j Aj ij1Bj ij Cj ij 1 57 where A B and C are constants for a given j note that they depend on nj 13 Now suppose that in the calculation of 51 j we have incurred in an error with absolute value 5 constant for all j so that instead of calculating 51 we have calculated an erroneous value 51 j here we will purposely leave out the cause of the calculation error This means that iJ Cij 5forallj 58 Consequently when we apply equation 57 to determine 514 we will calculate erroneous values by doing 1Lj Aj ij1Bjij Cj ij 1 59 since the values at the previous step i are already erroneous Subtracting equation 57 which would involve the correct values for dimensionless temperature from this equation yields i1j i1j Aj 1j1 ij1 Bj C514 1j Cjij 1 ij 1 60 Now we take the absolute value of this equation and make use of the identities labcl lallbllcl 61 labl W 62 to write equation 60 as follows 811j HLj SlAjH LjH ij1 lBjH ij ij leH ij l 8ij 1 63 Using equation 58 we get l im 451le S lAJ l lBJ l lCJ D5 64 This equation gives an upper bound for the error that would be made in step il in terms of the error 5 introduced in step i We can see that if Aj Bj Cj S l the error at step il will be at most 5 so that successive calculations will either keep the error the same or will decrease it On the other hand if Aj Bj Cj gtl it is possible that the error will grow with successive calculation steps This could be a cause for instability since calculation errors are always present eg round off and they can grow uncontrollably In conclusion we can say that the method will be stable if 14 Aj Bj Cj 1 stability criterion 65 It is said that the forwarddifference method for the PDE studied here is conditionally stable since equation 65 represents a sufficient condition for stability Let us explore this condition in detail First comparing equations 47 and 57 we see that A 2L2 1 66 g1nj g 2n j EFF Z 2h 2 67 g l nj Cl ZLZ l 68 g1nj g 2m We can see that Aj and Cj are always positive Regarding Cj note that 1ggt12nj always since the minimum value that nj can have is fg for interior nodes This means that AJ39 Aj C j cj Using these facts and substituting equations 66 to 68 into equation 65 leads to after rearrangement 2h 2h 1 S 1 stability criterion 69 g2lnJ2 g2lnJ2 The absolute value of a number can never be less than the number itself This means that equation 69 can only be satisfied by the equality The fact that a number is equal to its absolute value means that the number is positive or zero We then conclude that equation 69 is satisfied only if 2h lt1 2 gza nj stability criterion 70 Since this equation should be valid for any value of m we will consider the least favorable case which corresponds to the value of nj closest to 1 for an interior node which is nN lg Figure 3 This turns equation 70 into 0lt1 2h stability criterion 71 g21 1 g2 which leads to 15 0 S l i stability criterion 72 g32g Since we expect that gltlt2 desirable to nd a solution close to the exact solution we can make the approximation 2 g z 2 and equation 72 becomes after manipulations h S g3 stability criterion 73 For example for the results presented in Figures 4 to 6 for g005 equation 73 states that the solution is stable if hS0000125 which is consistent with the results presented Next we explore a numerical method that is always stable for the problem under consideration The backward difference method In this method we approximate the first derivative with respect to E in equation 17 by using the backwarddifference formula 36 91 39 81 39 muVJ 74 3 BE The derivatives on the righthand side of equation 17 will be once again represented by the centered formulas 44 and 45 This leads to the following discretized form of the PDE 1j 91 1 i 8ij1 81j 1 8ij1 29g QM 1 h n j 2g g2 lnJZ 75 The main difference of this equation with the forwarddifference method is that this equation cannot be applied to explicitly calculate the dimensionless temperature in a stepbystep manner since the equation contains 3 different values of temperature at step i 51 ji j1i j1 Instead a linear system of equations is found if the equation is rearranged as follows h l l A 2h A quot7 GM 9m g nJ g m g nJ h 1 1 A A 7 77 818L g1 nJ2lg 2m I 1 J i12j1 toN 76 For a given i this is a system of N linear algebraic equations with N unknowns 51 j jl to N The equations corresponding to jl and N can be simplified further by using the boundary conditions in T as follows 1 For jl fg and equation 50 discretized form of boundary condition 19 applies This 16 transforms equation 76 to after manipulations A 311 A 1 77 2g21 g2l 1 1 2g21 g2 2 For jN fl g and equation 48 discretized form of boundary condition 20 applies 12 1 11 3h This transforms equation 76 to after manipulations 2 3 h A 2h A A h 3 g 81N 1 1 flaw 81 1N 3 78 2g ZgXlg g Zg 2g lg The resulting system of equations can be written in matrix form as follows b1 01 0 9U 11 212 b2 02 812 d2 0 213 b3 C3 813 d3 79 aN i bN 1 CN l iN 1 dN l 0 0 aN bN LN dN where the coefficients are given by b1 1 80 2g lg bj 1i2 j23N 81 g21 Tj aj L2 l i j23N 82 g1nj g 2m cjL2 li j12N1 83 g1nJ g 2m dj i1jj12N1 84 85 A h dN 839 1N 1 2g31 g 17 The matrix in equation 79 has as only nonzero coefficients those along the diagonal and two lines parallel and adjacent to the diagonal For this reason this is a tri diagonal matrix A linear system of equations with a tridiagonal matrix can be solved by the relatively simple Thomas algorithm The Thomas algorithm Consider the following linear system of equations identical to equation 79 b1 01 0 0 0 0 0 X1 11 a 2 b 2 c 2 0 0 0 0 X 2 d 2 0 b 0 0 0 d a 3 c X 3 86 O O O aN i bN 1 CN i XN l dN 1 0 aN bN XN dN O O O O This system can be manipulated by replacing a row ie equation by a different row that results from the linear combination of the row in question with any other rows in the system A simplification of the system of equations is obtained if all rows j for j23N are modified by using row jl to eliminate the jl column For example we can replace row 2 by the following operation row 2 9 row 2 a72rowl 87 b1 This will make the first coefficient of the second row equal to zero b1 c1 0 0 0 X1 d1 0 b azc c 0 0 d azd 7 X 7 2 b1 1 2 2 2 b1 1 88 0 a3 b3 C3 0 X3 d3 After this step the second equation has only two nonzero coefficients those corresponding to x and X3 This new equation can be used to eliminate the coefficient of x a3 in equation 88 from row 3 Sequential application of this procedure will lead to a matrix with zeroes in place of the aj39s of the original matrix Once this process reaches the last row after eliminating aN the last row will have only one nonzero coefficient corresponding to xN so that xN can be calculated directly Once xN is known xN1 can be calculated directly from row Nl since this row will have nonzero coefficients only for xN1 and xN Continuing this back substitution until row 1 is reached gives the solution to the problem This is the Thomas algorithm which can be formulated as follows Step 1 7 Elimination From j2 to N eliminate coefficient of Xj1 in equation j 18 bJebJ bJ cJ1 89 J l a diedj J dJ1 90 J l Step 2 7 Back substitution Find xN XN LN 91 bN FromjNl to 1 nd xj l Xjdj ijjl 92 J Use of the Thomas algorithm to solve equation 79 yields the nodal values of the dimensionless temperature for a given step i in the E direction Successive application for increasing values of i yield the nal numerical solution of the problem A comparison of the backward and forward difference methods leads to the following observations 1 Both methods have a truncation error of Ohg2 2 The forwarddifference method is explicit since values of the dependent variable at a given step in E are expressed in terms of the values at the previous step On the other hand the backwarddifference method is implicit since values of the dependent variable at a given step in imust be calculated by solving a system of equations 3 The forwarddifference method is conditionally stable whereas the backwarddifference method is unconditionally stable regardless of the discretization This last fact can be proven by performing an analysis similar to the one presented above for the forwarddifference method Improving stability by using an implicit method apparently leads to an increase in computational effort e g needing to solve a system of equations instead of doing a direct calculation However the fact that a stable solution can be found with a relatively coarse grid might make the backwarddifference method more computationally efficient in some cases As was the case for ODEs Model 1 it would be worth exploring higher order methods Specifically the fact that both methods presented above have errors of Oh in the axial direction might not be optimal for the numerical solution Next we consider a method that improves the error in the axial direction while keeping unconditional stability and introducing a minimum of extra computational effort The Crank Nicolson method Consider the firstorder forward difference formula used before equation 42 Going back to the original Taylor series equation 25 the exact expression is 19 3 3E 81 39 8139 2 ism ltmpomz 93 Similarly consider the rstorder backward difference formula 74 The exact expression based on the Taylor series 28 is 3 ij 81 1jh32 35 8 2 2 5mm h 2 Biz 1m0h 94 If this equation is evaluated at the next step il instead of i we have 3 811j 81jh32 35 ii1 j f Zg mnpmmz 95 We will take advantage of the fact that in the two approximations for the derivative equations 93 and 95 the second derivative term appears with different signs To see why this will be helpful consider the arithmetic average of equations 93 and 95 19 3 M 2a 1m BE 5115T1J h h 328 328 0112 96 Z EHIM E j nm Here we are interested in the order with respect to h of the terms involving second derivatives in this equation To find out we expand the function 328 3amp2 in a Taylor series in E at constant 111 around Ei and evaluate the series at Ei1 to get 828 828 838 Bij i lj iNjhg iamOh2 97 It is clear from this equation that 828 828 Eil j iam0h 98 so that the second term on the righthand side of equation 96 is of 0h2 and can be assimilated with other terms of the same order which simplifies equation 96 to 188 3 Eai iallj BE 81 39 8139 51lanj J J0h2 99 20 Next we evaluate the PDE 17 at the points Edmj and i1n j after discretizing the right hand side using equations 43 to 45 If we keep the error terms the resulting equations are 38 8 181 8 12881 7iinj 1 I 1 2 2 0g2 100 BE ZgnJlnj g lnj 8 8 28 8 a g lanj 1ll 131 1 1ll 2 1l2 1l 1Ogz 101 BE ZgnJlnj g lnj We can take the arithmetic average of the rst derivatives from these two equations and then substitute it into equation 99 The result is after manipulations eroalz h 4gnj1nj 102 811j1 81j1 2 ij 2 iLj 8Lj 1 81L11 0g2 2g20 nj2 Eliminating higher order terms in this equation leads to an approximation scheme with a local truncation error of Oh2g2 which improves upon the approximation of the forward and backwarddifference methods analyzed before The resulting equation is the basic approximation of the Crank Nicolson method The truncated form of equation 102 can be rearranged as follows note use of circum ex for approximated values aj i1j 1 l39bj HLj lquot3j i1j1 aj ij 1 kj 1j Cj 1j1j1N 103 where a h 2 i 104 2glnJ g 2m h b 1 105 1 gza nf cj 2 li 106 2g1nJ g 2m h kj 1 107 g21n 21 When equation 103 is applied to calculate values of the dimensionless temperature at step i1 its righthand side contains values evaluated at step i so that equation 103 constitutes a tri diagonal system of equations To build the tridiagonal matrix we use equation 103 for j2 to Nl and the top and bottom rows are modi ed to account for the boundary conditions as follows 1 Forj1 boundary condition 19 implies see equation 50 8410 9m 108 1 and 510 511 109 and equation 103 becomes 31 b1 i11 cl i12 a1 k0611 01613 110 2 For jN boundary condition 20 requires 811N1 81N1 1 111 and equation 103 becomes aN811N 1 bN i1N aN81N71 kN iN 20N 112 The resulting system of equations can be solved using Thomas algorithm Figure 7 shows a comparison of the error in the calculation of the dimensionless temperature at a point between the two implicit methods developed in this chapter as a function of h The difference in the dependence of the errors on h can be clearly seen At low h clearly eOh for the backwarddifference method and eOh2 for the CrankNicolson method It is worthwhile to point out that the forwarddifference method is unstable in the range of h values used in Figure 7 l0E00 10E01 10E02 10E03 10E04 10E05 10E06 10E07 10E08 10E09 10E10 10E05 Figure 7 Comparison of the error at 005 and 1105 as a function of h using g001 E El CrankNicolson O Backward Differences 10E04 10E03 10E02 10E01 22 Arbitrary lines of slopes l and 2 are shown The exact solution for the calculation of the error was taken as the CrankNicolson solution with h10397 and g10393
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'