Chapter 30 The Key to Approaching This Problem is to Recast the PDE as a System of ODEs

CHAPTER 30 30.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 node dTi T − 2Ti + Ti +1 = k i −1 dt ∆x 2 By writing this equation for each node, the solution reduces to solving 4 simultaneous ODEs with Heun’s 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 shown t 0 0.1 0.2 • • • 3 6 9 12 x=0 100 100 100 x=2 0 2.04392 4.00518 x=4 0 0.02179 0.08402 x=6 0 0.01089 0.04267 x=8 0 1.02196 2.00259 x = 10 50 50 50 100 100 100 100 37.54054 53.24295 62.39033 68.71329 10.2745 24.66054 36.64937 46.03496 6.442321 17.46032 27.84901 36.5421 18.95732 27.92252 34.34692 39.53549 50 50 50 50 100 t=0 80 t=3 60 t=6 40 t=9 20 t = 12 0 0 5 10 30.2 Because we now have derivative boundary conditions, the boundary nodes must be simulated. For node 0, l +1 T0 l l l l = T0 + λ( T1 − 2T0 + T−1 ) (i) This introduces an exterior node into the solution at i = −1. The derivative boundary condition can be used to eliminate this node, dT dx = 0 T1 − T−1 2∆x which can be solved for T−1 = T1 − 2∆x dT0 dx which can be substituted into Eq. (i) to give l +1 T0  l dT l l l = T0 + λ  2T1 − 2T0 − 2 ∆x 0  dx      For our case, dT0/dx = 1 and ∆x = 2, and therefore T−1 = T1 + 4. This can be substituted into Eq. (i) to give, T0l +1 = T0l + λ ( 2T1l − 2T0l + 4) A similar analysis can be used to embed the zero derivative in the equation for the fifth node. The result is T5l +1 = T5l + λ( 2 T4l − 2 T5l ) Together with the equations for the interior nodes, the entire system can be iterated 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 results is also shown t x=0 x=2 x=4 x=6 x=8 x = 10 0 0.1 0.2 0.3 0.4 0.5 • • • 200 400 600 800 1000 50.0000 49.9165 49.8365 49.7597 49.6861 49.6153 50.0000 50.0000 49.9983 49.9949 49.9901 49.9840 50.0000 50.0000 50.0000 50.0000 49.9999 49.9997 50.0000 50.0000 50.0000 50.0000 49.9999 49.9997 50.0000 50.0000 49.9983 49.9949 49.9901 49.9840 50.0000 49.9165 49.8365 49.7597 49.6861 49.6153 30.00022 13.30043 -3.40115 -20.1055 -36.8103 31.80019 15.10041 -1.60115 -18.3055 -35.0103 33.20009 16.50035 -0.20115 -16.9055 -33.6103 34.19992 17.50028 0.798846 -15.9055 -32.6103 34.79981 18.10024 1.398847 -15.3055 -32.0103 34.99978 18.30023 1.598847 -15.1055 -31.8103 40 20 0 0 2 4 6 8 10 0 200 400 600 800 1000 -20 -40 Notice what’s happening. 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), t 0 0.1 0.2 x=0 x=2 x=4 x=6 x=8 100 0 0 0 0 100 2.0875 0 0 1.04375 100 4.087847 0.043577 0.021788 2.043923 x = 10 50 50 50 For ∆t = 0.05, it is t 0 0.05 0.1 0.15 0.2 x=0 x=2 x=4 x=6 x=8 100 0 0 0 0 100 1.04375 0 0 0.521875 100 2.065712 1.09E-02 5.45E-03 1.032856 100 3.066454 3.23E-02 0.016228 1.533227 100 4.046528 6.38E-02 3.22E-02 2.023265 x = 10 50 50 50 50 50 To assess the differences between 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 relative errors of the two Euler solutions. The results are summarized as x=0 x=2 x=4 x=6 x = 8 x = 10 100 4.006588 0.083044 0.042377 2.003302 50 Heun (h = 0.001) Euler (h = 0.1) Error relative to Heun 100 4.087847 0.043577 0.021788 2.043923 2.0% 47.5% 48.6% 2.0% 50 Euler (h = 0.05) Error relative to Heun 100 4.046528 0.063786 0.032229 2.023265 1.0% 23.2% 23.9% 1.0% 50 Notice, that as would be expected for Euler’s 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, Eq. (30.8) is first written for that node as − λT4l +1 + (1 + 2λ )T5l +1 − λT6l +1 = T5l (i) The value outside the system (i = 6) can be eliminated by writing the finite difference relationship for the derivative at node 5 as dT dx = 5 T6 − T4 2∆x which can be solved for T6 = T4 − 2∆x dT dx 5 For our case, dT/dx = 0, so T6 = T4 and Eq. (i) becomes l +1 −2 λT4 l +1 + (1 + 2 λ ) T5 l = T5 Thus, the simultaneous equations to be solved at the first step are −0.020875  1.04175  T11   −0.020875 1.04175  1 −0.020875   T2    T31  = −0.020875 1.04175 −0.020875   −0.020875 1.04175 −0.020875 T41      −0.04175 1.04175  T51   2.0875  0     0   0     0  which can be solved for  2 .004645   0.040186     0.000806   1.62 × 10−5     6.47 × 10 −7  For the second step, the right-hand side is modified to reflect these computed values of T at t = 0.1, −0.020875  1.04175  T11   4 .092145   −0.020875 1.04175   1   0.040186  −0.020875   T2      T31  =  0.000806  −0.020875 1.04175 −0.020875   −0.020875 1.04175 −0.020875 T41   1.62 × 10 −5        −0.04175 1.04175  T51  6.47 × 10−7  which can be solved for  3.930497   0117399  .    0.003127   7.83 × 10 −5    3.76 × 10 −6  30.5 The solution is identical to Example 30.3, but with 6 segments. Thus, the simultaneous equations to be solved at the first step are − 0.03006  2.06012  T11  6.012   −0.020875 2.06012 − 0.03006  1   0    T2      T31  =  0  − 0.03006 2.06012 − 0.03006   − 0.03006 2.06012 − 0.03006 T41   0        − 0.03006 2.06012  T51  3.006  which can be solved for 2.91890 0.04260     0.00093  0.02131    1.45945  For the second step, the right-hand side is modified to reflect these computed values of T at t = 0.1, − 0.03006  2.06012  T11  11.67559   −0.020875 2.06012 − 0.03006   1   0.17042    T2      T31  =  0.00373  − 0.03006 2.06012 − 0.03006   − 0.03006 2.06012 − 0.03006 T41   0.08524        − 0.03006 2.06012  T51   5.83780  which can be solved for 5.66986   0.16553    0 . 00543   0.08282     2.83493 The solution at t = 10 for this problem (n = 6) along with the results determined for n = 5, as in Example 30.3, are displayed in the following plot: 100 50 0 0 5 10 30.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 tridiagonal equations −0.0835  2.167  T1,1   6.2625       −0.0835 2.167 −0.0835T1,2  =  6.2625     −0.0835 2.167  T1, 3  18.7875 which can be solved for T1,1 = 3.018843 T1,2 = 3.345301 T1, 3 = 8.798722 In a similar fashion, tridiagonal equations can be developed and solved for T2,1 = 0.130591 T2, 2 = 0.370262 T2 ,3 = 6133184 . and T3,1 = 11017962 . T3,2 = 1.287655 T3,3 = 7.029137 For the second step to t = 10, Eq. (30.22) is applied to nodes (1,1), (2,1), and (3,1) to yield  2.167 −0.0835  T1,1  12.07537   −0.0835 2.167 −0.0835T  =  0.27029     1,2       −0.0835 2.167  T1, 3  4.060943 which can be solved for T1,1 = 5.5883 T1, 2 = 0.412884 T1,3 = 1889903 . Tridiagonal equations for the other rows can be developed and solved for T2,1 = 6.308761 T2, 2 = 0.902193 T2 ,3 = 2.430939 and T3,1 = 16.8241 T3,2 = 12.1614 T3,3 = 13.25121 Thus, the result at the end of the first step can be summarized as i=0 j=4 j=3 j=2 j=1 j=0 75 75 75 i=1 150 16.824 6.309 5.588 0 i=2 150 12.161 0.902 0.413 0 i=3 150 13.251 2.431 1.89 0 i=4 25 25 25 The computation can be repeated, and the results for t = 2000 s are below: i=0 j=4 j=3 j=2 j=1 j=0 75 75 75 i=1 150 98.214 70.089 44.643 0 i=2 150 97.768 62.5 33.482 0 i=3 150 80.357 48.661 26.786 0 i=4 25 25 25 30.7 Although this problem can be modeled with the finite-difference approach (see Sec. 32.1), the control-volume method provides a more straightforward way to handle the boundary conditions. The boundary fluxes and the reaction term can be used to develop the discrete form of the advection-diffusion equation for the interior volumes as ∆x dcil dt = −D cil − cil−1 ∆x +D cil+1 − cil ∆x +U cil + cil−1 2 −U cil+1 + cil 2 − k ∆xcil or dividing both sides by ∆x, dc il c l − 2 cil + cil−1 cil+1 + cil−1 = D i +1 − U − kcil dt 2 ∆x ∆x 2 which is precisely the 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, ∆x dc1l c l − c1l c l + c1l =D 2 + Uc0l − U 2 − k∆xc1l dt ∆x 2 or dividing both sides by ∆x, dc1l c l − c l 2c l − c 2l − c1l = D 2 21 + 0 − kc1l dt 2 ∆ x ∆x For the last boundary node, no diffusion is allowed through the exit pipe and advection out of the tank is again handled with a backward difference, ∆x dc nl dt = −D c nl − cnl −1 ∆x +U c nl + cnl −1 2 − Ucnl − k∆xc nl or dividing both sides by ∆x, dc nl cl − cl c l − cnl = − D n 2n −1 + U n −1 − kcnl dt 2 ∆ x ∆x By writing these equations for each equally-spaced volume, the PDE is transformed into a system of ODEs. Explicit methods like Euler’s method or other higher-order RK methods can then be used to solve the system. The results with and 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 are x c 0.5 1.5 2.5 3.5 14.2042 12.6506 11.2610 10.0385 The results are also plotted below: 4.5 8.9847 5.5 8.1025 6.5 7.3928 7.5 6.8583 8.5 6.5000 9.5 6.3201 16 12 8 4 0 0 2 4 6 8 10 30.8 Option Explicit Sub EulerPDE() Dim i As Integer, j As Integer, np As Integer, ns As Integer Dim Te(20) As Single, dTe(20) As Single, tpr(20) As Single, Tepr(20, 20) As Single Dim k As Single, dx As Single, L As Single, tc As Single, tf As Single Dim tp As Single, t As Single, tend As Single, h As Single L = 10 ns = 5 dx = 2 k = 0.835 Te(0) = 100 Te(5) = 50 tc = 0.1 tf = 1 tp = 0.1 np = 0 tpr(np) = t For i = 0 To ns Tepr(i, np) = Te(i) Next i Do 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 To ns Tepr(j, np) = Te(j) Next j If t >= tf Then Exit Do Loop Sheets("sheet1").Select Range("a4").Select For i = 0 To np ActiveCell.Value = tpr(i) For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = Tepr(j, i) Next j ActiveCell.Offset(1, -ns - 1).Select Next i End Sub Sub Derivs(Te, dTe, ns, dx, k) Dim j As Integer For j = 1 To ns - 1 dTe(j) = k * (Te(j - 1) - 2 * Te(j) + Te(j + 1)) / dx ^ 2 Next j End Sub 30.9 This program is set up to either use Dirichlet or gradient boundary conditions depending on the values of the parameters istrt and iend. Option Explicit Sub EulerPDE() Dim i As Integer, j As Integer, np As Integer, ns As Integer Dim istrt As Integer, iend As Integer Dim Te(20) As Single, dTe(20) As Single, tpr(200) As Single, Tepr(20, 200) As Single Dim k As Single, dx As Single, L As Single, tc As Single, tf As Single Dim tp As Single, t As Single, tend As Single, h As Single Dim dTedx(20) As Single L = 10 ns = 5 dx = 2 k = 0.835 dTedx(0) = 1 istrt = 0 dTedx(ns) = 0 iend = ns Te(0) = 50 Te(1) = 50 Te(2) = 50 Te(3) = 50 Te(4) = 50 Te(5) = 50 tc = 0.1 tf = 1000 tp = 200 np = 0 tpr(np) = t For i = 0 To ns Tepr(i, np) = Te(i) Next i Do tend = t + tp If tend > tf Then tend = tf h = tc Do If t + h > tend Then h = tend - t Call Derivs(Te(), dTe(), istrt, iend, ns, dx, k, dTedx()) For j = istrt To iend 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 To ns Tepr(j, np) = Te(j) Next j If t >= tf Then Exit Do Loop Sheets("sheet1").Select Range("a4").Select For i = 0 To np ActiveCell.Value = tpr(i) For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = Tepr(j, i) Next j ActiveCell.Offset(1, -ns - 1).Select Next i End Sub Sub Derivs(Te, dTe, istrt, iend, ns, dx, k, dTedx) Dim j As Integer If istrt = 0 Then dTe(0) = k * (2 * Te(1) - 2 * Te(0) - 2 * dx * dTedx(0)) / dx ^ 2 End If For j = 1 To ns - 1 dTe(j) = k * (Te(j - 1) - 2 * Te(j) + Te(j + 1)) / dx ^ 2 Next j If iend = ns Then dTe(ns) = k * (2 * Te(ns - 1) - 2 * Te(ns) + 2 * dx * dTedx(ns)) / dx ^ 2 End If End Sub 30.10 Option Explicit Sub SimpImplicit() Dim Dim Dim Dim np, ns, i, j, n Te(10), dTe(10), tpr(100), Tepr(10, 100), Tei As Single k, dx, L, tc, tf, tp, t, tend, h, lambda e(10), f(10), g(10), r(10), x(10), xrod L = 10# ns = 5 dx = L / ns k = 0.835 Te(0) = 100# Te(ns) = 50# Tei = 0 For i = 1 To ns - 1 Te(i) = Tei Next i t = 0 np = 0 tpr(np) = t For i = 0 To ns Tepr(i, np) = Te(i) Next i tc = 0.1 tp = 0.1 tf = 1 Do tend = t + tp If tend > tf Then tend = tf h = tc Do If t + h > tend Then h = tend - t lambda = k * h / dx ^ 2 f(1) = 1 + 2 * lambda g(1) = -lambda r(1) = Te(1) + lambda * Te(0) For j = 2 To ns - 2 e(j) = -lambda f(j) = 1 + 2 * lambda g(j) = -lambda r(j) = Te(j) Next j e(ns - 1) = -lambda f(ns - 1) = 1 + 2 * lambda r(ns - 1) = Te(ns - 1) + lambda * Te(ns) Call Tridiag(e(), f(), g(), r(), Te(), ns - 1) t = t + h If t >= tend Then Exit Do Loop np = np + 1 tpr(np) = t For j = 0 To ns Tepr(j, np) = Te(j) Next j If t >= tf Then Exit Do Loop Range("b5").Select xrod = 0 For j = 0 To ns ActiveCell.Value = xrod ActiveCell.Offset(0, 1).Select xrod = xrod + dx Next j Range("a6").Select For i = 0 To np ActiveCell.Value = "t = " & tpr(i) For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = Tepr(j, i) Next j ActiveCell.Offset(1, -ns - 1).Select Next i Range("a6").Select End Sub Sub Tridiag(e, f, g, r, x, n) Call Decomp(e, f, g, n) Call Substit(e, f, g, r, n, x) End Sub Sub Decomp(e, f, g, n) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k End Sub Sub Substit(e, f, g, r, n, x) Dim k As Integer For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub 30.11 Option Explicit Sub CrankNic() Dim Dim Dim Dim np, ns, i, j, n Te(10), dTe(10), tpr(100), Tepr(10, 100), Tei As Single k, dx, L, tc, tf, tp, t, tend, h, lambda e(10), f(10), g(10), r(10), x(10), xrod L = 10# ns = 5 dx = L / ns k = 0.835 Te(0) = 100# Te(5) = 50# Tei = 0 t = 0 np = 0 tpr(np) = t For i = 0 To ns Tepr(i, np) = Tei Next i tc = 0.1 tf = 10# tp = 1# Do tend = t + tp If tend > tf Then tend = tf h = tc Do If t + h > tend Then h = tend - t lambda = k * h / dx ^ 2 f(1) = 2 * (1 + lambda) g(1) = -lambda r(1) = lambda * Te(0) + 2 * (1 - lambda) * Te(1) + lambda * Te(2) r(1) = r(1) + lambda * Te(0) For j = 2 To ns - 2 e(j) = -lambda f(j) = 2 * (1 + lambda) g(j) = -lambda r(j) = lambda * Te(j - 1) + 2 * (1 - lambda) * Te(j) + lambda * Te (j + 1) Next j e(ns - 1) = -lambda f(ns - 1) = 2 * (1 + lambda) r(ns - 1) = lambda * Te(ns - 2) + 2 * (1 - lambda) * Te(ns - 1) + lambda * Te(ns) r(ns - 1) = r(ns - 1) + lambda * Te(ns) Call Tridiag(e(), f(), g(), r(), Te(), ns - 1) t = t + h If t >= tend Then Exit Do Loop np = np + 1 tpr(np) = t For j = 0 To ns Tepr(j, np) = Te(j) Next j If t >= tf Then Exit Do Loop Range("b5").Select xrod = 0 For j = 0 To ns ActiveCell.Value = xrod ActiveCell.Offset(0, 1).Select xrod = xrod + dx Next j Range("a6").Select For i = 0 To np ActiveCell.Value = "t = " & tpr(i) For j = 0 To ns ActiveCell.Offset(0, 1).Select ActiveCell.Value = Tepr(j, i) Next j ActiveCell.Offset(1, -ns - 1).Select Next i Range("a6").Select End Sub Sub Tridiag(e, f, g, r, x, n) Call Decomp(e, f, g, n) Call Substit(e, f, g, r, n, x) End Sub Sub Decomp(e, f, g, n) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k End Sub Sub Substit(e, f, g, r, n, x) Dim k As Integer For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub 30.12 Here is VBA code to solve this problem. The Excel output is also attached showing values for the first two steps along with selected snapshots of the solution as it evolves in time. Option Explicit Sub Dim Dim Dim Dim Dim Dim Dim Dim Dim Dim Dim ADI() np As Integer, i As Integer, j As Integer nx As Integer, ny As Integer Lx As Single, dx As Single Ly As Single, dy As Single Te(10, 10) As Single, dTe(10, 10) As Single tpr(100) As Single, Tepr(10, 10, 100) As Single, Tei As Single k As Single dt As Single, ti As Single, tf As Single, tp As Single t As Single, tend As Single, h As Single lamx As Single, lamy As Single e(10) As Single, f(10) As Single, g(10) As Single, r(10) As Single, Ted(10) As Single 'set computation parameters Lx = 40 nx = 4 dx = Lx / nx Ly = 40 ny = 4 dy = Ly / ny k = 0.835 dt = 10 tf = 500 ti = 0 tp = 10 Tei = 0 'set top boundary For i = 1 To nx - 1 Te(i, ny) = 100 Next i 'set bottom boundary For i = 1 To nx - 1 Te(i, 0) = 0 Next i 'set left boundary For j = 1 To ny - 1 Te(0, j) = 75 Next j 'set right boundary For j = 1 To ny - 1 Te(nx, j) = 50 Next j 'set corners for plot Te(0, 0) = (dy * Te(1, 0) + dx * Te(0, 1)) / (dy + dx) Te(nx, 0) = (dy * Te(nx - 1, 0) + dx * Te(nx, 1)) / (dy + dx) Te(0, ny) = (dy * Te(1, ny) + dx * Te(0, ny - 1)) / (dy + dx) Te(nx, ny) = (dy * Te(nx - 1, ny) + dx * Te(nx, ny - 1)) / (dy + dx) 'set interior For i = 1 To nx - 1 For j = 1 To ny - 1 Te(i, j) = Tei Next j Next i 'save initial values for output np = 0 t = ti tpr(np) = t For i = 0 To nx For j = 0 To ny Tepr(i, j, np) = Te(i, j) Next j Next i Do tend = t + tp If tend > tf Then tend = tf h = dt Do If t + h > tend Then h = tend - t 'Sweep y lamx = k * h / dx ^ 2 lamy = k * h / dy ^ 2 For i = 1 To nx - 1 f(1) = 2 * (1 + lamy) g(1) = -lamy r(1) = lamx * Te(i - 1, 1) + 2 * (1 - lamx) * Te(i, 1) + lamx * Te(i + 1, 1) _ + lamy * Te(i, 0) For j = 2 To ny - 2 e(j) = -lamy f(j) = 2 * (1 + lamy) g(j) = -lamy r(j) = lamx * Te(i - 1, j) + 2 * (1 - lamx) * Te(i, j) + lamx * Te(i + 1, j) Next j e(ny - 1) = -lamy f(ny - 1) = 2 * (1 + lamy) r(ny - 1) = lamx * Te(i - 1, ny - 1) + 2 * (1 - lamx) * Te(i, ny - 1) _ + lamx * Te(i + 1, ny - 1) + lamy * Te(i, nx) Call Tridiag(e(), f(), g(), r(), Ted(), nx - 1) For j = 1 To ny - 1 Te(i, j) = Ted(j) Next j Next i t = t + h / 2 'Sweep x For j = 1 To ny - 1 f(1) = 2 * (1 + lamx) g(1) = -lamx r(1) = lamy * Te(1, j - 1) + 2 * (1 - lamy) * Te(1, j) + lamy * Te(1, j + 1) _ + lamx * Te(0, j) For i = 2 To nx - 2 e(i) = -lamx f(i) = 2 * (1 + lamx) g(i) = -lamx r(i) = lamy * Te(i, j - 1) + 2 * (1 - lamy) * Te(i, j) + lamy * Te(i, j + 1) Next i e(nx - 1) = -lamx f(nx - 1) = 2 * (1 + lamx) r(nx - 1) = lamy * Te(nx - 1, j - 1) + 2 * (1 - lamy) * Te(nx - 1, j) _ + lamy * Te(nx - 1, j + 1) + lamx * Te(ny, j) Call Tridiag(e(), f(), g(), r(), Ted(), nx - 1) For i = 1 To nx - 1 Te(i, j) = Ted(i) Next i Next j t = t + h / 2 If t >= tend Then Exit Do Loop 'save values for output np = np + 1 tpr(np) = t For i = 0 To nx For j = 0 To ny Tepr(i, j, np) = Te(i, j) Next j Next i If t >= tf Then Exit Do Loop 'output results back to sheet Range("a5").Select Range("a5:e2005").ClearContents For k = 0 To np ActiveCell.Value = "t = " & tpr(k) ActiveCell.Offset(1, 0).Select For j = ny To 0 Step -1 For i = 0 To nx ActiveCell.Value = Tepr(i, j, k) ActiveCell.Offset(0, 1).Select Next i ActiveCell.Offset(1, -nx - 1).Select Next j ActiveCell.Offset(1, 0).Select Next k Range("a5").Select End Sub Sub Tridiag(e, f, g, r, x, n) Call Decomp(e, f, g, n) Call Substit(e, f, g, r, n, x) End Sub Sub Decomp(e, f, g, n) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k End Sub Sub Substit(e, f, g, r, n, x) Dim k As Integer For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub 30.13 MATLAB solution: %PDE Parabolic Problem - Heat conduction in a rod % u[xx]=u[t] % BC u(0,t)=0 u(1,t)=1 % IC u(x,0)=0 x<1 % % % % % % % % i=spatial index, from 1 to imax imax = no. of x points n=time index from 1 to nmax nmax = no. of time steps, Crank-Nicolson Formulation imax=61; nmax=60; % last time step = nmax+1 Constants dx=1/(imax-1); dx2=dx*dx; dt=dx2; % Setting dt to dx2 for good stability and results Independent space variable x=0:dx:1; Sizing matrices u=zeros(imax,nmax+1); t=zeros(1,nmax+1); a=zeros(1,imax); b=zeros(1,imax); c=zeros(1,imax); d=zeros(1,imax); ba=zeros(1,imax); ga=zeros(1,imax); up=zeros(1,imax); % Boundary Conditions u(1,1)=0; u(imax,1)=1; % % Time step loop n=1 represents 0 time, n+1 = next time step t(1)=0; for n=1:nmax t(n+1)=t(n)+dt; % Boundary conditions & Constants u(1,n+1)=0; u(imax,n+1)=1; dx2dt=dx2/dt; % coefficients b(2)=-2-2*dx2dt; c(2)=1; d(2)=(2-2*dx2dt)*u(2,n)-u(3,n); for i=3:imax-2 a(i)=1; b(i)=-2-2*dx2dt; c(i)=1; d(i)=-u(i-1,n)+(2-2*dx2dt)*u(i,n)-u(i+1,n); end a(imax-1)=1; b(imax-1)=-2-2*dx2dt; d(imax-1)=-u(imax-2,n)+(2-2*dx2dt)*u(imax-1,n)-2; % Solution by Thomas Algorithm ba(2)=b(2); ga(2)=d(2)/b(2); for i=3:imax-1 ba(i)=b(i)-a(i)*c(i-1)/ba(i-1); ga(i)=(d(i)-a(i)*ga(i-1))/ba(i); end % Back substitution step u(imax-1,n+1)=ga(imax-1); for i=imax-2:-1:2 u(i,n+1)=ga(i)-c(i)*u(i+1,n+1)/ba(i); end dt=1.1*dt; end % end of time step loop % Plot % Storing plot value of u as up, at every 5 time steps, np=5 %j=time index %i=space index np=5; for j=np:np:nmax for i=1:imax up(i)=u(i,j); end plot(x,up) hold on end grid title('u[xx]=u[t]; u(0,t)=0, Formulation') xlabel('x - ND Space') ylabel('u - ND Temperature') hold off u(1,t)=1, u(x,0)=0; Crank-Nicolson % Storing times for temp. profiles %These can be saved in a data file or examined in the command file tp=zeros(1,(nmax-1)/np); i=1; tp(1)=0; for k=np:np:nmax i=i+1; tp(i)=t(k); end tp gtext('n=60');gtext('n=50');gtext('n=40');gtext ('n=30'); gtext('n=20');gtext('n=10');gtext('n=1'); gtext('t=.766'); gtext('t=.1115');gtext('t=.0413');gtext('t=.014'); gtext('t=.0038');gtext('t=0') u[x x ]= u[t]; u(0,t)= 0, u(1,t)= 1, u(x ,0)= 0; Crank -N ic ols on F orm ulation 1 0.9 0.8 u - N D Tem perature 0.7 0.6 n= 60 t= .766 0.5 n= 50 0.4 n= 40 t= .1115 0.3 n= 30 t= .0413 0.2 n= 20 t= .014 0.1 n= 1 t= 0 n= 10 t= .0038 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x - ND S pac e tp = Columns 1 through 7 0 0.0013 0.0038 0.0078 0.0142 0.0246 0.0413 Columns 8 through 13 0.0682 0.1115 0.1813 0.2937 0.4746 0.7661 30.14 ∂2u ∂r 2 + 1 ∂u ∂u = r ∂r ∂t Substituting of second order correct Crank-Nicolson analogues 1  u i+1, n+1 − u i ,n +1 + ui −1, n+1 u i +1,n − u i ,n + u i−1,n  +   ∂r 2 2  ∆r 2 ∆r 2  ∂u 1  u i +1, n+1 − u i−1, n+1 u i+1,n − u i −1, n  =  +  ∂r 2  2∆r 2 ∆r  r = (i − 1)∆r ∂ 2u = ∂u u i , n+1 − u i, n = ∂t ∆t into the governing equation give the following finite difference equations:   1  ∆r 2   1   1  1 − 2(i − 1)  u i−1,n+1 +  − 2 − 2 ∆t u i ,n +1 + 1 + 2(i − 1) u i +1, n+1 =  − 1 + 2(i − 1)  u i−1,n           ∆r 2  1  − 2 − 2 u i ,n + − 1 +  u i+1, n ∆t  2(i − 1)    For the end points: x = 1 (i = R), substitute the value of uR = 1 into the above FD equation x = 0 (i = 1), set the FD analog to the first derivative = 0 1  u 2,n +1 − u 0, n+1 u 2,n − u 0,n  ∂u  +  ∂r  = 2  2 ∆r 2 ∆r   i =1    =0  Also substitute in i = 1 into the finite difference equation and algebraically eliminate u 0,n +1 + u 0,n from the two equations and get the FD equation at i = 1:   ∆r 2  ∆r 2  − 2 − 2 u1, n +1 + [ 2 ]u 2 ,n +1 = −  2 − 2 u1, n + [ − 2]u 2 , n ∆t  ∆t    %PDE Parabolic Problem - Heat conduction in the radial direction in a circular rod % u[rr]+(1/r)u[r]=u[t] 0