- Email: [email protected]
Efficient discretization schemes for Graetz-type problems
Computers and Chemical Engineering 23 (1999) 377—383
Efficient discretization schemes for Graetz-type problems E. Bruce Nauman*, Joseph Savoca The Howard P. Isermann Department of Chemical Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, USA Received 22 August 1998; revised 30 September 1998; accepted 30 September 1998
Abstract Uniform grid spacing for the finite difference solution of partial differential equations is computationally inefficient. Two nonuniform discretization schemes are presented and applied to Graetz-type problems. Results are compared to those for a uniform discretization scheme and shown to be able to decrease the amount of computational work required for the same accuracy. 1999 Elsevier Science Ltd. All rights reserved.
1. Motivation and problem statement Solving thousands of ordinary differential equations (ODEs) is now routine but solving thousands of partial differential equations (PDEs) still taxes most computational facilities. Our interest in solving a large set of PDEs stems from a problem in polymer reaction engineering, a reactive—convective—diffusion problem involving many polymeric species of different chain lengths (Kim and Nauman, 1997). In this problem, the diffusion coefficient depends upon the chain length of the particular species, and a rigorous but computationally intensive approach is to solve thousands of PDEs, one for each species. One approach to decreasing the computational effort for such problems is the use of a more nearly optimum grid spacing. Uniform grid spacings for the finite difference solution of PDEs are easily implemented but computationally inefficient. Small step sizes may be necessary to meet stability requirements, and the total computational work can be enormous. The general problem under consideration is a form of the convective diffusion equation, *c v "D c#Source term, X *z
orthogonal to z, c is a scalar such as temperature or concentration which is conserved in the absence of a source term, and D is the diffusivity. A property of real flow systems is that v "0 at solid boundaries such as X tube walls. This fact leads to numerical difficulties in most solution techniques. A classic example of Eq. (1) is the Graetz problem governing heat transfer to a fluid in laminar flow in a circular tube
1*c *c r *c 2uN 1! "D # r *r *r R *z
(2)
subject to the boundary conditions
*c "0. (3) *r A closely related problem governs the concentration of a first-order reactant, c(0, r)"0,
c(z, R)"1,
r *c 1 *c *c 2uN 1! "D # !kc r *r *r R *z
(4)
subject to the boundary conditions (1)
where v is the velocity component in the axial or z direcX tion, is the Laplacian in one or two dimensions
* Corresponding author. Tel.: 001 518 276 6726; fax: 001 518 276 4030; e-mail: [email protected]
*c *c "0, "0. (5) *r *r 0 Analytical solutions are available for these systems (Huang et al., 1984), and they are commonly used to test computation methods. We use them here for that purpose. Generalization to more complex situations will be considered subsequently.
c(0, r)"1,
0098-1354/99/$ — see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 9 8 ) 0 0 2 8 0 - 4
378
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Most solution techniques for equations such as (2) and (4) begin with discretization in the radial direction. Examples include the well-known method of lines and other techniques such as the control volume formulation or Band (J) (Newman, 1968) and semianalytical techniques (Haran and White, 1996). Regardless of method, the choice is almost invariably a uniform radial spacing with *r"1/N where N is the number of radial increments. The major purpose of this paper is to show that such uniform grids may be inefficient from a computational viewpoint. We illustrate this point using the method of lines, although the general conclusion that uniform grids are suboptimal is expected to apply to any solution technique using radial discretization. An example of nonuniform axial spacing is also explored.
2. Variable radial discretization Second-order finite difference approximations with unequal increments are *c !*r c (*r !*r )c G+ G> G\ # G> G G *r *r *r (*r #*r )*r G G> G G G> *r c G G> # (*r #*r )*r G G> G> and
(6)
*c 2c 2c G+ G\ G ! *r (*r #*r )*r *r *r G G> G G G> 2c G> # , (7) (*r #*r )*r G > G> where *r "r !r , i"1, 2,2, N. Substitution into G G G\ Eq. (2) gives
dc (!*r #2r )c (*r !*r !2r )c G"D G> G G\# G> G G G dz r *r (*r #*r ) r *r *r G G G G> G G G> (*r #2r )c r G G G> # 2uN 1! G . (8) r (*r #*r )*r R G G G> G> A special form applies at the centerline, i"0, r "0 dc !4c 4c "D # (2uN ). (9) dz *r *r The centerline boundary condition has been applied in the derivation of Eq. (9). The wall boundary condition is just c "1 for the Graetz problem of Eqs. (2) and (3). , The wall boundary condition for the reactive case of Eqs. (4) and (5) is
!*r c (*r #*r )c , ,\ , ,\ ,\ . c " # , (*r #2*r )*r (*r #2*r )*r ,\ , ,\ ,\ , ,\
(10)
This result is obtained by using a second-order approximation for *c/*r"0 at the wall. The reactive case also includes a source term, !kc . G The set of ODEs defined by Eqs. (8) and (9) can be rewritten as dc G"A c #B c #C c #Source term evaluated G G\ G G G G> dz at c , 1(i(N!1, G dc "B c #C c #Source term evaluated at c , (11) dz where the A , B , C , i"0, 1, 2, N!1, form a tridiagG G G onal matrix of coefficients. The set is stiff in the sense that the eigenvalues of the matrix differ greatly in magnitude. Some stiffness is inherent in the method of lines and arises even in the case of diffusion in a flat plate geometry with a flat velocity profile. The analytical solutions to convective diffusion problems have an infinite, semi-periodic set of eigenvalues, and ODE approximations to the analytical solutions must become increasingly stiff as the number of radial increments is increased. Indeed, stiffness is the price of an accurate solution and limits the allowable axial step size, *z, for any explicit integration scheme. However, substantial latitude exists for minimizing the decrease in allowable *z. The diagonal coefficients are !2D B " , *ruN D(*r !*r !2r )/r *r *r G\ G G G G G>, B" G 2uN (1!r/R) G
0(i(N!1. (12)
These diagonal coefficients can all be made equal by a judicious choice of the r . A fast, simple algorithm for G doing this is given in Table 1. Table 1 Algorithm for determining r for equal diagonal coefficients G FOR m"0 TO n v(m)"(1!(m/n) 2) NEXT r(0)"0 r(1)"1 FOR nn"1 TO 16 c"r(1) 2/4 FOR m"1 TO n!1 d¸"r(m)!r(m!1) dH"(d¸#2*r(m))/(1#v(m)*r(m)*d¸/c) r(m #1) " r(m) # dH NEXT m FOR m"1 TO n r(m)"r(m)/r(n) v(m)"(1!r(m) 2) NEXT NEXT nn
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
379
Table 2 Nonuniform radial increments for a parabolic velocity profile and N"8 Nonuniform increments
Uniform increments
Increment number, i
Radius, r /R G
Maximum step size, *z /¸
Radius, r /R G
Maximum step size, *z /¸
0 1 2 3 4 5 6 7 8
0.0 0.138322 0.222604 0.334594 0.432771 0.550024 0.667027 0.811314 1.0
0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 na
0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0
0.078125 0.153809 0.146484 0.134277 0.117188 0.095215 0.068359 0.036621 na
Equalization of the B dramatically lowers the stiffness G of Eq. (11). When this set of equations is solved by Euler’s method, equalization of the B has the directly tangible G benefit of giving the same value for the maximum permissible *z at each radial point. With Euler’s method, the maximum stable step size is given by v *z " X .
B G
(13)
The allowable *z generally varies with radial position. To solve the entire set of Eq. (11), the most restrictive of the maximum step sizes must be used. With uniform *r, this most restrictive size occurs at the radial position nearest the wall, R-*r, and is small since v is small near X the wall. Table 2 shows results of the algorithm for the case of N"8. The results for *z /¸ apply to the parabolic
velocity profile, v "2uN (1!r/R). For comparison, X Table 2 also shows *z /L for the parabolic profile with
uniform radial increments. The nonuniform radial increments give a permissible axial step size of 0.096 compared to 0.037 for the uniform case. Accuracy of the algorithm in Table 1 depends on the limit of the nn loop. As shown, this limit is 16, and the resulting radial increments are accurate to about eight decimal places. Note that the discretization technique does not require highly accurate values of the r to be effective. G 3. Variable axial discretization Another approach to computational efficiency is to use a uniform radial grid with a nonuniform axial grid. After the radial discretization, each radial point was assigned a constant axial step size which was an integer multiple of the most restrictive step size. The use of different step sizes for a group of stiff ODEs has been used before, for
example, to solve kinetics problems with widely varying time scales (Martinez and Drozdowicz, 1989) but required an approximation of the governing ODEs. The radial point closest to the tube wall has the most stringent stability requirement as in Eq. (5). This stability requirement was used to set the axial step size for this radial point. This axial increment was used as the basis step size for all of the other radial points. The inner radial points were assigned axial step sizes which were integral multiples of the basis step size, depending upon the ratio of the axial velocity at that point to the axial velocity at the near wall point. To perform the integration, each radial point is marched ahead by Euler’s method by its assigned axial step size. Each radial point is checked and the one with an unknown temperature at the smallest axial distance has a new value calculated for it further down the tube. If a point on either side of the point to be integrated has already been integrated further down the tube, then linear interpolation is used on the last two calculated temperatures at the neighboring point to estimate the temperature at the current axial position. This method reduces the computation necessary since much larger axial steps are used at the radial positions where the solutions change on slower time scales. The three discretization methods, equal radial increments, equal diagonal coefficients and variable axial increments were used to solve the Graetz and reactive problems for varying number of radial increments. The values of c used for comparison with the analytical solutions were the bulk or mixing cup average as calculated using the trapezoidal rule and the centerline value. The wall value was also calculated for the reactive case. The exact solutions were calculated from an infinite series of confluent hypergeometric functions as formulated by Huang et al. (1984). The series gave identical values to eight decimal places for 10 and 13 terms so convergence was assumed. The dimensionless group governing the solution is DtM /R where tM "¸/uN is the mean residence time. The
380
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
reactive case also includes the dimensionless rate constant, ktM . Figs. 1 and 2 show errors in the bulk and centerline values as calculated for the Graetz problem with DtM /R"0.1. Other values produced similar results. The abscissas are the number of function evaluations; that is, the number of times a new value for c(r, z) is calculated. With each doubling of the number of radial increments,
the number of function evaluations increased by a factor of 16 (2) for the uniform grid case, a factor of 8 (2) for the nonuniform radial grid and a factor of 9 for the nonuniform axial grid case. Actual running times were roughly proportional to the number of function evaluations, although the variable axial grid case did require approximately an order of magnitude more time per evaluation compared to the other two methods.
Fig. 1. Error in bulk temperature for Graetz problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
Fig. 2. Error in centerline temperature for Graetz problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Figs. 3—5 show the bulk, centerline and wall errors for the reactive case with DtM /R"1.0 and ktM "1.0. Although there are some anomalous results for crude grids, the general trends seem clear. The nonuniform axial grid outperformed the uniform grid in all cases. The nonuniform radial grid outperformed the uniform grid for the Graetz case but was worse for the reactive case. Where differences exist, they are likely to be large. For
381
highly accurate results, the nonuniform axial grid outperformed the uniform grid by two orders of magnitude. The focus of this note is on obtaining high precision results. The comparison of methods is intended to reflect their relative merits in the limit of small grid sizes and consequently very many functional evaluations. It may well be that conclusions regarding computational efficiency are different if low levels of precision are
Fig. 3. Error in bulk concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
Fig. 4. Error in centerline concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
382
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Fig. 5. Error in wall concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
acceptable. For example, the uniform grid is superior at low precision levels for the cases illustrated in Figs. 1 and 2. Fig. 2 even shows the anomalous result that the precision worsens in one instance when the grid size is reduced. Again, the conclusions apply to what appears to be the asymptotic behavior of the various methods. Although the majority of our results used Euler’s method to solve the ODEs of Eq. (8), limited tests were done using Runge—Kutta integration. We specifically used the NAG subroutine D02PCF with orders 2—3 and 4—5. It was found that the accuracy of the solution was limited by the radial discretization so that the more sophisticated axial integration techniques gave marginal improvements at the same number of radial increments and were less efficient on the basis of function evaluations.
4. Generalizations and future work A judicious choice of discretization scheme has been shown to significantly decrease the computational work required for solving Graetz-type problems. Indeed, the computational effort needed to achieve a given accuracy can be reduced by several orders of magnitude. It is not yet possible to predict a priori what may be the best method for a particular problem, and this is clearly a direction for future work. The ideas presented in this paper can be generalized in several ways. The radial discretization scheme can be applied to solution techniques other than the method of lines. It can also be applied to other physical problems.
The application to non-Newtonian velocity profiles is straightforward. The method is directly suited to simultaneous equations involving the diffusion of several chemical species as well as heat. The same choice of radial increments will give equal values of the diagonal coefficients for all the species. Turning to versions of Eq. (1) where c depends on two spatial variables, the diagonal coefficients will depend on discretization in two directions, but the choices can still be made to set all the diagonal coefficients equal. Thus the method can be applied to flow in noncircular ducts. The approach of nonuniform axial steps as used here can also be applied to a variety of physical problems: non-Newtonian flow, simultaneous equations and noncircular ducts. Its application to solution techniques other than the method of lines with Euler integration is complicated, but may be worthwhile for parallel computing since the algorithm minimizes the need to exchange information at different radial locations.
Nomenclature c D i k ¸ N r R tM
temperature or concentration diffusivity increment number reaction rate constant tube length number of radial increments radial coordinate tube radius average residence time
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
uM v X z
average axial velocity axial velocity axial coordinate
References Haran, B. S., & White, R. E. (1996). A semianalytical technique for solving nonlinear partial differential equations. Comp. Appl. Engng. Ed, 4, 229—240.
383
Huang, C., Matlosz, M., Pan, W., & Snyder, Jr., W. (1984). Heat transfer to a laminar flow fluid in a circular tube. AIChE J., 30, 833—835. Kim, D., & Nauman, E. B. (1997). Nonterminating polymerizations in continuous flow systems. IECR, 36, 1088—1094. Martinez, E. C., & Drozdowicz, B. (1989). Multitime-scale approach to real-time simulation of stiff dynamic systems. Comp. Chem. Engng, 13, 767—778. Newman, J. (1968). Numerical solution of coupled ordinary differential equations. IECF, 7, 514—517.
Efficient discretization schemes for Graetz-type problems E. Bruce Nauman*, Joseph Savoca The Howard P. Isermann Department of Chemical Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, USA Received 22 August 1998; revised 30 September 1998; accepted 30 September 1998
Abstract Uniform grid spacing for the finite difference solution of partial differential equations is computationally inefficient. Two nonuniform discretization schemes are presented and applied to Graetz-type problems. Results are compared to those for a uniform discretization scheme and shown to be able to decrease the amount of computational work required for the same accuracy. 1999 Elsevier Science Ltd. All rights reserved.
1. Motivation and problem statement Solving thousands of ordinary differential equations (ODEs) is now routine but solving thousands of partial differential equations (PDEs) still taxes most computational facilities. Our interest in solving a large set of PDEs stems from a problem in polymer reaction engineering, a reactive—convective—diffusion problem involving many polymeric species of different chain lengths (Kim and Nauman, 1997). In this problem, the diffusion coefficient depends upon the chain length of the particular species, and a rigorous but computationally intensive approach is to solve thousands of PDEs, one for each species. One approach to decreasing the computational effort for such problems is the use of a more nearly optimum grid spacing. Uniform grid spacings for the finite difference solution of PDEs are easily implemented but computationally inefficient. Small step sizes may be necessary to meet stability requirements, and the total computational work can be enormous. The general problem under consideration is a form of the convective diffusion equation, *c v "D c#Source term, X *z
orthogonal to z, c is a scalar such as temperature or concentration which is conserved in the absence of a source term, and D is the diffusivity. A property of real flow systems is that v "0 at solid boundaries such as X tube walls. This fact leads to numerical difficulties in most solution techniques. A classic example of Eq. (1) is the Graetz problem governing heat transfer to a fluid in laminar flow in a circular tube
1*c *c r *c 2uN 1! "D # r *r *r R *z
(2)
subject to the boundary conditions
*c "0. (3) *r A closely related problem governs the concentration of a first-order reactant, c(0, r)"0,
c(z, R)"1,
r *c 1 *c *c 2uN 1! "D # !kc r *r *r R *z
(4)
subject to the boundary conditions (1)
where v is the velocity component in the axial or z direcX tion, is the Laplacian in one or two dimensions
* Corresponding author. Tel.: 001 518 276 6726; fax: 001 518 276 4030; e-mail: [email protected]
*c *c "0, "0. (5) *r *r 0 Analytical solutions are available for these systems (Huang et al., 1984), and they are commonly used to test computation methods. We use them here for that purpose. Generalization to more complex situations will be considered subsequently.
c(0, r)"1,
0098-1354/99/$ — see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 9 8 ) 0 0 2 8 0 - 4
378
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Most solution techniques for equations such as (2) and (4) begin with discretization in the radial direction. Examples include the well-known method of lines and other techniques such as the control volume formulation or Band (J) (Newman, 1968) and semianalytical techniques (Haran and White, 1996). Regardless of method, the choice is almost invariably a uniform radial spacing with *r"1/N where N is the number of radial increments. The major purpose of this paper is to show that such uniform grids may be inefficient from a computational viewpoint. We illustrate this point using the method of lines, although the general conclusion that uniform grids are suboptimal is expected to apply to any solution technique using radial discretization. An example of nonuniform axial spacing is also explored.
2. Variable radial discretization Second-order finite difference approximations with unequal increments are *c !*r c (*r !*r )c G+ G> G\ # G> G G *r *r *r (*r #*r )*r G G> G G G> *r c G G> # (*r #*r )*r G G> G> and
(6)
*c 2c 2c G+ G\ G ! *r (*r #*r )*r *r *r G G> G G G> 2c G> # , (7) (*r #*r )*r G > G> where *r "r !r , i"1, 2,2, N. Substitution into G G G\ Eq. (2) gives
dc (!*r #2r )c (*r !*r !2r )c G"D G> G G\# G> G G G dz r *r (*r #*r ) r *r *r G G G G> G G G> (*r #2r )c r G G G> # 2uN 1! G . (8) r (*r #*r )*r R G G G> G> A special form applies at the centerline, i"0, r "0 dc !4c 4c "D # (2uN ). (9) dz *r *r The centerline boundary condition has been applied in the derivation of Eq. (9). The wall boundary condition is just c "1 for the Graetz problem of Eqs. (2) and (3). , The wall boundary condition for the reactive case of Eqs. (4) and (5) is
!*r c (*r #*r )c , ,\ , ,\ ,\ . c " # , (*r #2*r )*r (*r #2*r )*r ,\ , ,\ ,\ , ,\
(10)
This result is obtained by using a second-order approximation for *c/*r"0 at the wall. The reactive case also includes a source term, !kc . G The set of ODEs defined by Eqs. (8) and (9) can be rewritten as dc G"A c #B c #C c #Source term evaluated G G\ G G G G> dz at c , 1(i(N!1, G dc "B c #C c #Source term evaluated at c , (11) dz where the A , B , C , i"0, 1, 2, N!1, form a tridiagG G G onal matrix of coefficients. The set is stiff in the sense that the eigenvalues of the matrix differ greatly in magnitude. Some stiffness is inherent in the method of lines and arises even in the case of diffusion in a flat plate geometry with a flat velocity profile. The analytical solutions to convective diffusion problems have an infinite, semi-periodic set of eigenvalues, and ODE approximations to the analytical solutions must become increasingly stiff as the number of radial increments is increased. Indeed, stiffness is the price of an accurate solution and limits the allowable axial step size, *z, for any explicit integration scheme. However, substantial latitude exists for minimizing the decrease in allowable *z. The diagonal coefficients are !2D B " , *ruN D(*r !*r !2r )/r *r *r G\ G G G G G>, B" G 2uN (1!r/R) G
0(i(N!1. (12)
These diagonal coefficients can all be made equal by a judicious choice of the r . A fast, simple algorithm for G doing this is given in Table 1. Table 1 Algorithm for determining r for equal diagonal coefficients G FOR m"0 TO n v(m)"(1!(m/n) 2) NEXT r(0)"0 r(1)"1 FOR nn"1 TO 16 c"r(1) 2/4 FOR m"1 TO n!1 d¸"r(m)!r(m!1) dH"(d¸#2*r(m))/(1#v(m)*r(m)*d¸/c) r(m #1) " r(m) # dH NEXT m FOR m"1 TO n r(m)"r(m)/r(n) v(m)"(1!r(m) 2) NEXT NEXT nn
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
379
Table 2 Nonuniform radial increments for a parabolic velocity profile and N"8 Nonuniform increments
Uniform increments
Increment number, i
Radius, r /R G
Maximum step size, *z /¸
Radius, r /R G
Maximum step size, *z /¸
0 1 2 3 4 5 6 7 8
0.0 0.138322 0.222604 0.334594 0.432771 0.550024 0.667027 0.811314 1.0
0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 0.095664 na
0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0
0.078125 0.153809 0.146484 0.134277 0.117188 0.095215 0.068359 0.036621 na
Equalization of the B dramatically lowers the stiffness G of Eq. (11). When this set of equations is solved by Euler’s method, equalization of the B has the directly tangible G benefit of giving the same value for the maximum permissible *z at each radial point. With Euler’s method, the maximum stable step size is given by v *z " X .
B G
(13)
The allowable *z generally varies with radial position. To solve the entire set of Eq. (11), the most restrictive of the maximum step sizes must be used. With uniform *r, this most restrictive size occurs at the radial position nearest the wall, R-*r, and is small since v is small near X the wall. Table 2 shows results of the algorithm for the case of N"8. The results for *z /¸ apply to the parabolic
velocity profile, v "2uN (1!r/R). For comparison, X Table 2 also shows *z /L for the parabolic profile with
uniform radial increments. The nonuniform radial increments give a permissible axial step size of 0.096 compared to 0.037 for the uniform case. Accuracy of the algorithm in Table 1 depends on the limit of the nn loop. As shown, this limit is 16, and the resulting radial increments are accurate to about eight decimal places. Note that the discretization technique does not require highly accurate values of the r to be effective. G 3. Variable axial discretization Another approach to computational efficiency is to use a uniform radial grid with a nonuniform axial grid. After the radial discretization, each radial point was assigned a constant axial step size which was an integer multiple of the most restrictive step size. The use of different step sizes for a group of stiff ODEs has been used before, for
example, to solve kinetics problems with widely varying time scales (Martinez and Drozdowicz, 1989) but required an approximation of the governing ODEs. The radial point closest to the tube wall has the most stringent stability requirement as in Eq. (5). This stability requirement was used to set the axial step size for this radial point. This axial increment was used as the basis step size for all of the other radial points. The inner radial points were assigned axial step sizes which were integral multiples of the basis step size, depending upon the ratio of the axial velocity at that point to the axial velocity at the near wall point. To perform the integration, each radial point is marched ahead by Euler’s method by its assigned axial step size. Each radial point is checked and the one with an unknown temperature at the smallest axial distance has a new value calculated for it further down the tube. If a point on either side of the point to be integrated has already been integrated further down the tube, then linear interpolation is used on the last two calculated temperatures at the neighboring point to estimate the temperature at the current axial position. This method reduces the computation necessary since much larger axial steps are used at the radial positions where the solutions change on slower time scales. The three discretization methods, equal radial increments, equal diagonal coefficients and variable axial increments were used to solve the Graetz and reactive problems for varying number of radial increments. The values of c used for comparison with the analytical solutions were the bulk or mixing cup average as calculated using the trapezoidal rule and the centerline value. The wall value was also calculated for the reactive case. The exact solutions were calculated from an infinite series of confluent hypergeometric functions as formulated by Huang et al. (1984). The series gave identical values to eight decimal places for 10 and 13 terms so convergence was assumed. The dimensionless group governing the solution is DtM /R where tM "¸/uN is the mean residence time. The
380
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
reactive case also includes the dimensionless rate constant, ktM . Figs. 1 and 2 show errors in the bulk and centerline values as calculated for the Graetz problem with DtM /R"0.1. Other values produced similar results. The abscissas are the number of function evaluations; that is, the number of times a new value for c(r, z) is calculated. With each doubling of the number of radial increments,
the number of function evaluations increased by a factor of 16 (2) for the uniform grid case, a factor of 8 (2) for the nonuniform radial grid and a factor of 9 for the nonuniform axial grid case. Actual running times were roughly proportional to the number of function evaluations, although the variable axial grid case did require approximately an order of magnitude more time per evaluation compared to the other two methods.
Fig. 1. Error in bulk temperature for Graetz problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
Fig. 2. Error in centerline temperature for Graetz problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Figs. 3—5 show the bulk, centerline and wall errors for the reactive case with DtM /R"1.0 and ktM "1.0. Although there are some anomalous results for crude grids, the general trends seem clear. The nonuniform axial grid outperformed the uniform grid in all cases. The nonuniform radial grid outperformed the uniform grid for the Graetz case but was worse for the reactive case. Where differences exist, they are likely to be large. For
381
highly accurate results, the nonuniform axial grid outperformed the uniform grid by two orders of magnitude. The focus of this note is on obtaining high precision results. The comparison of methods is intended to reflect their relative merits in the limit of small grid sizes and consequently very many functional evaluations. It may well be that conclusions regarding computational efficiency are different if low levels of precision are
Fig. 3. Error in bulk concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
Fig. 4. Error in centerline concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
382
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
Fig. 5. Error in wall concentration for reactive problem. Uniform grid (- -), nonuniform radial grid ()), nonuniform axial grid (— —).
acceptable. For example, the uniform grid is superior at low precision levels for the cases illustrated in Figs. 1 and 2. Fig. 2 even shows the anomalous result that the precision worsens in one instance when the grid size is reduced. Again, the conclusions apply to what appears to be the asymptotic behavior of the various methods. Although the majority of our results used Euler’s method to solve the ODEs of Eq. (8), limited tests were done using Runge—Kutta integration. We specifically used the NAG subroutine D02PCF with orders 2—3 and 4—5. It was found that the accuracy of the solution was limited by the radial discretization so that the more sophisticated axial integration techniques gave marginal improvements at the same number of radial increments and were less efficient on the basis of function evaluations.
4. Generalizations and future work A judicious choice of discretization scheme has been shown to significantly decrease the computational work required for solving Graetz-type problems. Indeed, the computational effort needed to achieve a given accuracy can be reduced by several orders of magnitude. It is not yet possible to predict a priori what may be the best method for a particular problem, and this is clearly a direction for future work. The ideas presented in this paper can be generalized in several ways. The radial discretization scheme can be applied to solution techniques other than the method of lines. It can also be applied to other physical problems.
The application to non-Newtonian velocity profiles is straightforward. The method is directly suited to simultaneous equations involving the diffusion of several chemical species as well as heat. The same choice of radial increments will give equal values of the diagonal coefficients for all the species. Turning to versions of Eq. (1) where c depends on two spatial variables, the diagonal coefficients will depend on discretization in two directions, but the choices can still be made to set all the diagonal coefficients equal. Thus the method can be applied to flow in noncircular ducts. The approach of nonuniform axial steps as used here can also be applied to a variety of physical problems: non-Newtonian flow, simultaneous equations and noncircular ducts. Its application to solution techniques other than the method of lines with Euler integration is complicated, but may be worthwhile for parallel computing since the algorithm minimizes the need to exchange information at different radial locations.
Nomenclature c D i k ¸ N r R tM
temperature or concentration diffusivity increment number reaction rate constant tube length number of radial increments radial coordinate tube radius average residence time
E.B. Nauman, J. Savoca / Computers and Chemical Engineering 23 (1999) 377—383
uM v X z
average axial velocity axial velocity axial coordinate
References Haran, B. S., & White, R. E. (1996). A semianalytical technique for solving nonlinear partial differential equations. Comp. Appl. Engng. Ed, 4, 229—240.
383
Huang, C., Matlosz, M., Pan, W., & Snyder, Jr., W. (1984). Heat transfer to a laminar flow fluid in a circular tube. AIChE J., 30, 833—835. Kim, D., & Nauman, E. B. (1997). Nonterminating polymerizations in continuous flow systems. IECR, 36, 1088—1094. Martinez, E. C., & Drozdowicz, B. (1989). Multitime-scale approach to real-time simulation of stiff dynamic systems. Comp. Chem. Engng, 13, 767—778. Newman, J. (1968). Numerical solution of coupled ordinary differential equations. IECF, 7, 514—517.