numerical analysis chapra 5ESM30.docx

1、numerical analysis chapra 5ESM30CHAPTER 3030.1 The key to approaching this problem is to recast the PDE as a system of ODEs. Thus, by substituting the finite-difference approximation for the spatial derivative, we arrive at the following general equation for each nodeBy writing this equation for eac

2、h node, the solution reduces to solving 4 simultaneous ODEs with Heuns method. The results for the first two steps along with some later selected values are tabulated below. In addition, a plot similar to Fig. 30.4, is also showntx = 0x = 2x = 4x = 6x = 8x = 1001000000500.11002.0439230.0217880.01089

3、41.021962500.21004.0051780.0840220.0426722.00259350310037.5405410.274496.44231918.9573250610053.2429424.6605217.460327.9225150910062.3903236.6493727.8490134.34692501210068.7133146.0349836.5421339.53555030.2 Because we now have derivative boundary conditions, the boundary nodes must be simulated. For

4、 node 0, (i)This introduces an exterior node into the solution at i = 1. The derivative boundary condition can be used to eliminate this node,which can be solved forwhich can be substituted into Eq. (i) to giveFor our case, dT0/dx = 1 and x = 2, and therefore T1 = T1 4. This can be substituted into

5、Eq. (i) to give,A similar analysis can be used to embed the zero derivative in the equation for the nth node, (ii)This introduces an exterior node into the solution at n + 1. The derivative boundary condition can be used to eliminate this node,which can be solved forwhich can be substituted into Eq.

6、 (ii) to give For our case, n = 5 and dTn/dx = 0, and thereforeTogether with the equations for the interior nodes, the entire system can be solved with a step of 0.1 s. The results for some of the early steps along with some later selected values are tabulated below. In addition, a plot of the later

7、 results is also showntx = 0x = 2x = 4x = 6x = 8x = 10025.000025.000025.000025.000025.000025.00000.124.916525.000025.000025.000025.000025.00000.224.836524.998325.000025.000025.000025.00000.324.759724.994925.000025.000025.000025.00000.424.686124.990124.999925.000025.000025.00000.524.615324.984024.999

8、725.000025.000025.00002005.0000816.8000748.2000599.2000489.80004210.00004400-11.6988-9.89883-8.49883-7.49882-6.89881-6.69881600-28.4008-26.6008-25.2008-24.2008-23.6007-23.4007800-45.1056-43.3056-41.9056-40.9056-40.3056-40.10561000-61.8104-60.0104-58.6104-57.6104-57.0104-56.8104Notice whats happening

9、. The rod never reaches a steady state, because of the heat loss at the left end (unit gradient) and the insulated condition (zero gradient) at the right.30.3 The solution for t = 0.1 is (as computed in Example 30.1),tx = 0x = 2x = 4x = 6x = 8x = 1001000000500.11002.0875001.04375500.21004.0878470.04

10、35770.0217882.04392350For t = 0.05, it istx = 0x = 2x = 4x = 6x = 8x = 1001000.0000000.0000000.0000000.000000500.051001.0437500.0000000.0000000.521875500.11002.0657120.0108940.0054471.032856500.151003.0664540.0322840.0162281.533227500.21004.0465280.0637860.0322292.02326550To assess the differences b

11、etween the results, we performed the simulation a third time using a more accurate approach (the Heun method) with a much smaller step size (t = 0.001). It was assumed that this more refined approach would yield a prediction close to true solution. These values could then be used to assess the relat

12、ive errors of the two Euler solutions. The results are summarized asx = 0x = 2x = 4x = 6x = 8x = 10Heun (h = 0.001)1004.0065880.0830440.0423772.00330250Euler (h = 0.1)1004.0878470.0435770.0217882.04392350Error relative to Heun2.0%47.5%48.6%2.0%Euler (h = 0.05)1004.0465280.0637860.0322292.02326550Err

13、or relative to Heun1.0%23.2%23.9%1.0%Notice, that as would be expected for Eulers method, halving the step size approximately halves the global relative error.30.4 The approach described in Example 30.2 must be modified to account for the zero derivative at the right hand node (i = 5). To do this, E

14、q. (30.8) is first written for that node as (i)The value outside the system (i = 6) can be eliminated by writing the finite difference relationship for the derivative at node 5 aswhich can be solved forFor our case, dT/dx = 0, so T6 = T4 and Eq. (i) becomesThus, the simultaneous equations to be solv

15、ed at the first step arewhich can be solved forFor the second step, the right-hand side is modified to reflect these computed values of T at t = 0.1,which can be solved for30.5 The solution is identical to Example 30.3, but with 9 interior nodes. Thus, the simultaneous equations to be solved at the

16、first step arewhich can be solved forFor the second step, the right-hand side is modified to reflect these computed values of T at t = 0.1,which can be solved for30.6 Using the approach followed in Example 30.5, Eq. (30.20) is applied to nodes (1,1), (1,2), and (1,3) to yield the following tridiagon

17、al equationswhich can be solved forIn a similar fashion, tridiagonal equations can be developed and solved forandFor the second step to t = 10, Eq. (30.22) is applied to nodes (1,1), (2,1), and (3,1) to yieldwhich can be solved forTridiagonal equations for the other rows can be developed and solved

18、forandThus, the result at the end of the first step can be summarized asi = 0i = 1i = 2i = 3i = 4j = 49012012012085j = 36013.462939.82028812.8694450j = 2605.0507560.8156854.2928150j = 1604.473950.4162053.73783350j = 03000025The computation can be repeated and the results for t = 2000 s are shown bel

19、ow:i = 0i = 1i = 2i = 3i = 4j = 49012012012085j = 36080.7142983.8392977.1428650j = 26059.0178657.554.7321450j = 16037.8571432.4107134.2857150j = 0300002530.7 Although this problem can be modeled with the finite-difference approach (see Sec. 32.1), the control-volume method provides a more straightfo

20、rward way to handle the boundary conditions. Thus, the tank is idealized as a series of control volumes:The boundary fluxes and the reaction term can be used to develop the discrete form of the advection-diffusion equation for the interior volumes asor dividing both sides by x,which is precisely the

21、 form that would have resulted by substituting centered finite difference approximations into the advection-diffusion equation.For the first boundary node, no diffusion is allowed up the entrance pipe and advection is handled with a backward difference,or dividing both sides by x,For the last bounda

22、ry node, no diffusion is allowed through the exit pipe and advection out of the tank is again handled with a backward difference,or dividing both sides by x,By writing these equations for each equally-spaced volume, the PDE is transformed into a system of ODEs. Explicit methods like Eulers method or

23、 other higher-order RK methods can then be used to solve the system.The results with an initial condition that the reactor has zero concentration with an inflow concentration of 100 (using Euler with a step size of 0.005) for t = 100 arex0.51.52.53.54.55.56.57.58.59.5c42.032041.512841.050940.646340.

24、298940.008739.776039.600839.483639.4248A plot of the results is shown below:30.8 Here is a VBA program that implements the explicit method. It is set up to duplicate Example 30.1.Option ExplicitSub Explicit()Dim i As Integer, j As Integer, np As Integer, ns As IntegerDim Te(20) As Double, dTe(20) As

25、 Double, tpr(20) As Double, Tepr(20, 20) As DoubleDim k As Double, dx As Double, L As Double, tc As Double, tf As DoubleDim tp As Double, t As Double, tend As Double, h As Double, tol As DoubleDim x As Doubletol = 0.000001L = 10ns = 5dx = L / nsk = 0.835Te(0) = 100Te(5) = 50tc = 0.1tf = 12tp = 3np =

26、 0tpr(np) = tFor i = 0 To ns Tepr(i, np) = Te(i)Next iDo tend = t + tp If tend tf Then tend = tf h = tc Do If t + h tend Then h = tend - t Call Derivs(Te, dTe, ns, dx, k) For j = 1 To ns - 1 Te(j) = Te(j) + dTe(j) * h Next j t = t + h If t = tend Then Exit Do Loop np = np + 1 tpr(np) = t For j = 0 T

27、o ns Tepr(j, np) = Te(j) Next j If t + tol = tf Then Exit DoLoopSheets(sheet1).SelectRange(a4:bb5000).ClearContentsRange(a4).SelectActiveCell.Value = timex = 0For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = x = & x x = x + dxNext jRange(a5).SelectFor i = 0 To np ActiveCell.Value =

28、tpr(i) For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = Tepr(j, i) Next j ActiveCell.Offset(1, -ns - 1).SelectNext iRange(a5).SelectEnd SubSub Derivs(Te, dTe, ns, dx, k)Dim j As IntegerFor j = 1 To ns - 1 dTe(j) = k * (Te(j - 1) - 2 * Te(j) + Te(j + 1) / dx 2Next jEnd SubWhen the pr

29、ogram is run, the result is30.9 This VBA program is set up to either use Dirichlet or gradient boundary conditions depending on the values of the parameters istrt and iend. It is set up to solve the first few steps of Prob. 30.2.Option ExplicitSub EulerPDE()Dim i As Integer, j As Integer, np As Integer, ns As IntegerDim is