Covalidation of Integral Transforms and Method of Lines in Nonlinear Convection-Diffusion with Mathematica

Covalidation of Integral Transforms and Method of Lines in Nonlinear Convection-Diffusion with Mathematica

Leonardo S. de B. Alves

Renato M. Cotta

Mikhail D. Mikhailov

Laboratório de Transmissão e Tecnologia do Calor - LTTC

Engenharia Mecânica - EE/COPPE/UFRJ

Universidade Federal do Rio de Janeiro

Caixa. Postal 68503, Cidade Universitária

21945-970 Rio de Janeiro, RJ. Brazil

[e-mail], [e-mail], [e-mail]

The Mathematica system (version 4.0) is employed in the solution of nonlinear difusion and convection-difusion problems, formulated as transient one-dimensional partial diferential equations with potential dependent equation coefficients. The Generalized Integral Transform Technique (GITT) is first implemented for the hybrid numerical-analytical solution of such classes of problems, through the symbolic integral transformation and elimination of the space variable, followed by the utilization of the built-in Mathematica function NDSolve for handling the resulting transformed ODE system. This approach ofers an error-controlled final numerical solution, through the simultaneous control of local errors in this reliable ODE's solver and of the proposed eigenfunction expansion truncation order. For covalidation purposes, the same built-in function NDSolve is employed in the direct solution of these partial diferential equations, as made possible by the algorithms implemented in Mathematica (versions 3.0 and up), based on application of the method of lines. Various numerical experiments are performed and relative merits of each approach are critically pointed out.

Keywords: Integral Transforms, symbolic computation, mathematica system, convection-diffusion

Introduction

Nonlinear diffusion and convection-diffusion problems are among the most usual mathematical formulations in modern engineering and physical sciences, and different solution methodologies have been proposed along the last few decades for this class of problems. Besides the more classical purely numerical approaches for these nonlinear partial differential equations, such as the best known finite differences and finite elements methods, a number of hybrid numerical-analytical procedures were also advanced within recent years. One such hybrid approach, originated within the heat and fluid flow application area, is known as the Generalized Integral Transform Technique (GITT) as an extension to classical integral transforms for linear transformable problems, and is well documented in a number of books and review articles [1, 2, 3, 4, 5, 6]. This computational approach, based on the Classical Integral Transform Method [1], offers the advantages of automatic global error control, such as in a purely analytical solution, and only mild increase in overall computational costs for increasing number of independent variables (spatial coordinates), since the numerical tasks in any problem are always concentrated in one single independent variable (time, for instance).

More recently, this approach has been progressively implemented within a mixed symbolic-numerical environment [3, 4, 7], based on the Mathematica platform [8], allowing for a marked reduction on formulae derivation and analysis effort, while providing integrated algorithms for numerical computations and graphical representation. A variety of notebooks have been constructed which automatically reproduce previously obtained numerical results for different classes of problems in heat and fluid flow, originally achieved through conventional algorithmic programming languages and lengthy hand-derived analytical expressions. In addition, a number of extensions have been proposed along the last few years, broadening the scope of application of this PDEs solution methodology coupled to this hybrid computational implementation strategy. The evident success of theses initiatives points towards emphasizing the more advanced use of these analysis tools in engineering problems.

On the other hand, the Mathematica system itself [8], provides a built-in function NDSolve which besides being an excellent ordinary differential systems solver, itself employed within the GITT computational methodology, allows for the straightforward numerical solution of certain classes of partial differential equations, based on the implementation of the method of lines for the time variable and a variable-order discrete scheme for the spatial coordinate. As a PDE solver, the Mathematica function NDSolve offers to a certain extent an unbeatable simplicity in use, coupled to the recognized robustness of the well-known method of lines, providing a very attractive automatic tool for analysis in certain classes of diffusion and convection-diffusion problems.

In the present covalidation exercise, it is our intention to briefly examine the accuracy expectations and limitations of these two solution alternatives for nonlinear one-dimensional diffusion and convection-diffusion problems, especially aimed at directing future extensions to multi-dimensional situations. First, the NDSolve function is employed as an ODE solver in combination to the integral transformation process for elimination of the space coordinate, and second the same built-in function is directly employed in the numerical solution of the proposed partial differential formulations. Representative examples are selected, extracted from the available literature on the GITT approach [9, 10], and related respectively to heat conduction with nonlinear thermal conductivity and to a nonlinear Burgers equation formulation.

Nomenclature

a, b = constant coefficients

nq = otal number of terms in filtered potential expansion

= dimensionless transformed potential

x = dimensionless spatial coordinate

= dimensionless transformed temperature

Subscripts

K = dimensionless coefficient in diffusion term

t = dimensionless time

T_f = dimensionless filtering potential

Greek Symbols

l_i = temperature eigenvalue

f = filtering potential

nT = total number of terms in potential expansion

T = dimensionless potential

U = dimensionless convection coefficient

q = dimensionless filtered temperature

Y_i = dimensionless eigenfunction

i;j ;k = expansion indices

Nonlinear Heat Conduction Problem

Statement of the Problem

The transient nonlinear one-dimensional conduction problem to be here considered is defined in dimensionless form by,

where thermal conductivity has, for instance, a linear dependence on temperature,

The temperature boundary and initial conditions used are given by,

Transformation Pairs

According to the formalism in the integral transform approach, the transformation pair to be used during the integral transformation process depends upon the auxiliary eigenvalue problem previously chosen, which is given by,

which has the following normalized solution,

So, the integral transformation can be defined as,

and the inversion formula as,

Transforming the PDE to an ODE System

The integral transformation of the partial differential equation (1) in the x direction gives,

Operating the linear diffusion term above by parts and using the boundary information from both eigenfunction and temperature potentials, we obtain,

Using equations (10) and (7), equation (9) becomes,

Using the temperature inversion formula (8) in equation (11) we find,

Transforming the temperature initial condition (4), we obtain,

After solving analytically the integrals within equations (12) and (13) using the Mathematica built-in function Integrate, another built-in function named NDSolve can now be applied to solve this transformed ordinary differential initial value system. Then, the inversion formula (8) can be used to recover the original temperature potential.

Nonlinear Convection-Diffusion Problem

Statement of the Problem

The transient nonlinear one-dimensional convection-diffusion problem governed by Burger's equation is defined in dimensionless form by,

where velocity has, for instance, a linear dependence on the potential,

The boundary and initial conditions used are given by,

Filtering

Before we begin to solve problem (14- 17), a filter is used to remove the non-homogeneous boundary condition from Eq. (16), and thus enhance the convergence behavior of the integral transform solution. The transformation introduced by the filter is,

where T_f (x) is the solution of the linear steady-state version of equation (14), given by,

Equation (19) is subjected to the following boundary conditions,

The analytical solution of the filter problem (19, 20) can be found using Mathematica built-in function DSolve and is given by,

Applying the filter (18) to equations (14 - 17) and using equation (19) to simplify the problem, we obtain the new homogeneous partial differential system that is going to be subjected to the integral transformation technique. Thus, the partial differential equation (18) becomes,

and is now subjected to the new boundary conditions,

and to the new initial condition,

Transformation Pairs

Once again we choose the auxiliary eigenvalue problem in order to define the inverse/transform pair,

which has the following normalized solution,

So, the integral transformation can be defined as,

and the inversion formula as,

Transforming the PDE to an ODE System

The integral transformation of the partial differential equation (22) in the x direction gives,

Operating the linear diffusion term above by parts and using the boundary information from both eigenfunction and temperature potentials, we obtain,

Using equations (30) and (27), equation (29) becomes,

Using the temperature inversion formula (28) in equation (31) we find,

Transforming the temperature initial condition (24), we obtain,

Once again, after analytically solving the integrals within equations (32) and (33) using Mathematica built-in function Integrate, the built-in function NDSolve can be applied to solve this transformed ordinary differential initial value system. Then, the inversion formula (28) is used with transformation (18) to recover the original temperature potential shown below,

Nonlinear Heat Conduction Results

Constant Thermal Conductivity (Benchmark)

The first step of our analysis is to compare the results of the linear transient heat conduction problem obtained from both the Generalized Integral Transform Technique (exact solution for this linear case) and NDSolve. To do so, we just simplify equation (1) by letting b = 0. Equation (1), which is nonlinear, has been solved throughout this work and its associated notebook by using the GITT, making it possible to have a more general solution, so we just need to let b = 0 at eq. (12). Then, we obtain the following system of ordinary differential equations (ODE) for the transformed temperature,

where i is the index of the equation.

The solution of this system, shown below, can be explicitly obtained from the built-in function DSolve,

Using the inversion formula (8), we find,

In tables 1 and 2 below, we show results for the temperature profile for values of t = 0.1 and 0.5, respectively. In both tables we have a = 1, and in each line we show the solution of the series expansion for increasing number of terms, as presented in the first column.

To obtain the results from NDSolve, we must perform some adaptation. The initial and boundary conditions must be consistent, as required by the method of lines implemented in this built-in function. So, we have to use a function at the boundary that has a value of 1 at t = 0 and that falls to zero as fast as possible for t > 0. This way, eqs. (3) turn into,

Another built-in function, named Timing, can be employed to evaluate the performance of the default use of NDSolve to compute the solution shown in table 3.

As one can see from the two sets of results above, from the integral transform solution and the built-in function, some accuracy improvement is to be sought in the numerical built-in solution. Some possible paths of improvement, according to the options available and to the suggestions in [8], are now numerically investigated. The accuracy from NDSolve results may be increased by using WorkingPrecision ® 20 and AccuracyGoal ® ¥ . The PrecisionGoal option is automatically changed once the WorkingPrecision has been increased, while we keep the AccuracyGoal to an infinite value, not to deal with absolute errors stopping criteria within the computation. These results are shown in table 4.

The accuracy of the NDSolve results may be further increased for comparison with the above results by using WorkingPrecision ® 24. The results are presented in table 5.

We can notice that increasing the overall precision used in computation we can obtain better results with NDSolve, which are closer to the GITT benchmark results, and that these results are apparently better for higher values of t. However, it is noticeable some saturation on the gain when just the working precision is further increased. Thus, we now attempt to increase the finite difference formula order employed in the discretization of the space variable, through the option DifferenceOrder. The default value is { 4, Automatic }, i.e., 4th order finite differences for the space variable (assuming that the domain for the x variable has been specified first in the call to NDSolve), and automatic scheme order in the ODE system solution in the time variable. Here we attempt to resolve the problem with a 6th order finite differences scheme for the x variable (table 6).

We have now reached a very good agreement between the analytical and the numerical solutions to within five significant digits, at increased computational effort for the numerical built-in function, once the related parameter options were modified from their default values to match the exact solution. Before proceeding towards the comparison of the two approaches for more involved nonlinear formulations, we consider now reducing the working precision, once the use of higher differentiation formulae for the space variable improved the overall accuracy quite markedly, aimed at achieving a reduced computational cost. The results are shown in table 7.

Finally one more attempt is made to alter the finite differences formulae, within the same working precision, which results in essentially the same results as above, but with a considerable increase in computational cost, with respect to the sixth order finite differences, as presented in table 8.

Nonlinear Thermal Conductivity

Now we are going to compare the results for the nonlinear case, using the full solution generated by the Generalized Integral Transform Technique. Again, we are going to use a = 1, but will generate results for b = 0.1, 1 and 10.

We generate the system of ODEs for the transformed temperature and solve it with NDSolve.

Fixing b = 0.1, we can solve the problem for different truncation orders of the infinite system, and generate results for t = 0:1 and 0:5 (tables 9 and 10) in the same way as done before:

Using here the same modification of the boundary conditions (38), NDSolve gives (table 11),

Using now the same parameters that best reproduced the benchmark case, we present the results on table 12 with the solution for both values of t used above.

Some improvement is observed, though not as noticeably as in the linear situation above. Thus, an attempt is made to increase also the working precision, towards a better agreement between the two proposed solutions (table 13), which finally provides a more evident agreement to the fourth significant digit with the GITT solution above obtained, at the price of increasing the computational effort on the method of lines solution.

Fixing b = 1.0, and thus increasing the nonlinearity importance, we can solve the problem for different truncation orders of the system, and again generate two tables (¹⁴ and ¹⁵) in the same way we had done before:

Once again using boundary conditions (38), NDSolve gives (table 16),

An agreement around two significant digits is achieved against the GITT solution, with the use of the default parameters in the NDSolve function, and this solution is now attempted to be improved through the same parameters modifications as employed before (table 17).

which finally provides a more evident agreement to the third significant digit with the GITT solution above obtained, at the price of increasing the computational effort on the NDSolve direct solution.

Fixing b = 10, for an even more nonlinear problem, we can solve the problem for different truncation orders of the ODE system, and generate two tables (¹⁸ and ¹⁹) as done before,

And again, using boundary conditions (38), NDSolve yields the results shown in table 20,

Now, for this even more nonlinear situation, the default usage of NDSolve leads to final results already disturbed at the second significant digit in different positions within the medium, and again we attempt to improve this behavior by modifying the default parameters:

which at least confirms more closely the agreement to the second significant digit with the GITT solution above obtained, at the price of substantial increase in the computational effort on the method of lines direct solution.

To remove still any concern about the accuracy of the GITT solution itself, we rerun this most nonlinear situation under more strict precision requirements within the NDSolve scheme for ODE's, as employed within the GITT solution, reconstructing with further precision the convergence tables above for lower truncation orders (tables 22 and 23).

The above new tables, with slight changes on the last two digits, confirm the controlled accuracy within the GITT solution here proposed. We can see that the parametric modifications implemented for the linear case in the direct NDSolve solution do not produce a proportional improvement in accuracy as for situations with increasing nonlinearity in the problem formulation, even for higher values of t, when compared with the error controlled solution obtained through the GITT. It is also quite noticeable from all the above tables that low truncation order solutions obtained through the eigenfunction expansion approach may provide reasonably accurate results for most engineering purposes, in fact comparable to those obtained through the method of lines after adjustment of the options parameters for better accuracy. In this sense, the adaptive procedure of truncation order reduction proposed in [2], can be employed in the construction of an automatic global error-controlling procedure, that works towards an user prescribed accuracy requirement at reduced computational cost.

Nonlinear Convection-Diffusion Results

Purely Diffusive Problem (Benchmark)

The first step of our next analysis is to compare the results of the linear transient heat conduction problem obtained from both the Classical Integral Transform Technique (exact solution) and NDSolve. To do so, we have to simplify equation (14) or (22) by letting a = b = 0. However, eq. (14), which is non-linear, has been solved throughout this work using the GITT instead, making it possible to have a more general solution, and here just let a = b = 0 at eq. (32), to obtain the following system of ODEs for the transformed temperature,

where i is the index of the equation.

The solution of this system, shown below, can be obtained from the built-in function DSolve,

Using the inversion formula (28) with transformation (18) to recover the original temperature potential we find,

Numerical results are here reported for K = 1 and presented in the tables below, for the tem- perature profiles at t = 0.1 and 0.5, respectively. In both tables 24 and 25, we show in each line the solution of the series expansion for increasing number of terms, which is shown in the first column.

To obtain the results from NDSolve (method of lines), we must perform some adaptation again. The initial and boundary conditions must be consistent, so we use the same exponential function used before, turning boundary conditions (16) into,

The results obtained from solving equation (14) subjected to conditions (42) and (17) using NDSolve are then presented in table 26, with the computational time needed by the built-in function,

Increasing the accuracy of NDSolve by using WorkingPrecision = 24 and AccuracyGoal = ¥ , we obtain results that have a better agreement with the Benchmark ones, as shown in table 27:

Nonlinear Convection-Diffusion Problem

We now compare the results for the nonlinear case, using the full solution generated by GITT, as in the previous section for the nonlinear heat conduction problem. Again, we take K = 1, a = 1 and use b = 0.1, 1.0 and 10.

Tables 28 and 29 present results for t = 0.1 and 0.5 for b = 0.1 obtained from GITT in the same way as in the previous case.

The direct solution of the nonlinear convection-diffusion problem using NDSolve shown in table 30 below differs in the fourth significant digit from the one obtained through the GITT. However, by setting WorkingPrecision = 24, AccuracyGoal = ¥ and DifferenceOrder = 6, the difference between these results drops to the sixth significant digit, as one may see in table 31, but with the price of having higher computational cost in the latter case.

The strength of the nonlinear term is now increased by letting b = 1, and the results obtained through the GITT are presented in tables 32 and 33 for t = 0.1 and 0.5.

Once again we obtain the direct solution of the nonlinear PDE using NDSolve, and show the results in table 34. One may notice that still there is a difference in the fourth significant digit against the results obtained through GITT. However, due to the higher nonlinearity of this case, we need to set WorkingPrecision = 24, AccuracyGoal = ¥ and DifferenceOrder = 8 in order to decrease this difference to the sixth significant digit, as one may see in table 35.

Finally, we make the problem highly nonlinear by setting b = 10, and the new results obtained through GITT are presented in tables 36 and 37 for t = 0.1 and 0.5. In this case, the GITT solution still needs more terms in order to achieve full convergence at the fifth significant digit.

Table 38 presents the results obtained from direct solution of equation (14) subjected to the boundary conditions (42) and (17), and one may see that differences between these results and those obtained through GITT appear on the second significant digit near the boundary at x = 1. Convergence is expected to be worse in this region due to the larger temperature gradients. In order to narrow this difference between the two solutions, we had to set WorkingPrecision = 28 and DifferenceOrder = 8 (table 39), and even though the difference dropped to the sixth significant digit almost everywhere, it only reached the fourth digit at x = 0.9.

The higher values of the computational time required to obtain the solutions were a natural consequence of increasing the order of magnitude of the nonlinear terms, and as expected, nonlinear convection terms require greater computational effort than the correspondent nonlinear diffusion terms. However, both solutions obtained, using the ODE solver within NDSolve and the PDE solver also within NDSolve, could keep track of accuracy by properly setting their respective options, despite the amount of time needed.

Two kinds of boundary conditions were employed in this work, prescribed temperature and insulating condition, for each problem analyzed. A problem with a heat flux condition would be treated in the same way as the prescribed temperature one, since the difficulties that might be added are essentially the same, slower convergence due to a non-homogeneous boundary condition. A proper filter is in general sufficient to accelerate convergence to within practical usage limits.

Conclusions

The Mathematica system was employed in the mixed symbolic-numerical solution of nonlinear one-dimensional diffusion and convection-diffusion problems by following two alternative paths. The first one involves the symbolic integral transform elimination of the spatial coordinates combined with the numerical solution of the resulting ODEs system through the built-in function NDSolve. The second one adopts the method of lines discrete implementation directly available in the same built-in function NDSolve for one-dimensional partial differential systems. This covalidation effort indicates that the direct built-in integration is a very simple to use alternative, providing numerical results with suficient accuracy for several engineering purposes, even for increasing importance of nonlinear terms, once the appropriate choices of the built-in option parameters are investigated by the user, for each specific problem. On the other hand, the integral transform approach presents itself as an alternative methodology for both low and high accuracy requirements, coupled to the very robust ODEs system solver available in the Mathematica package. Also, future extensions of these ideas may be combined for multi-dimensional situations, when the integral transform methodology may be employed to reduce the original multi-variable equation to a one-dimensional partial differential system, which could then be handled in a straightforward way by the NDSolve built-in function. New developments in the integral transform technique recently made available, such as reordering schemes to improve convergence, have also been allowing for a more straightforward application of GITT in two or three dimensional media.

Acknowledgements

The authors would like to acknowledge the financial support provided by FAPERJ, CNPq, PRONEX and CAPES, all of them sponsoring agencies in Brazil. The support of Wolfram Research, Inc., USA, is also sincerely acknowledged, through the Mathematica Technical Center of COPPE/UFRJ.

Mikhailov, M. D. and Cotta, R. M., "Integral transform method with mathematica", in preparation.

Article received: January, 2000. Technical Editor: Átila P. S. Freire

• Mikhailov, M. D. and Ozisik, M. N., "Unified Analysis and Solutions of Heat and Mass Diffusion", John Wiley & Sons, New York, 1984.
• Cotta, R. M., "Integral Transforms in Computational Heat and Fluid Flow", CRC Press, Boca Raton, FL, 1993.
• Cotta, R. M. and Mikhailov, M. D., "Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computation", John Wiley & Sons, England, 1997.
• Cotta, R. M. (ed.), "The Integral Transform Method in Thermal and Fluids Science and Engineering", Begell House, Inc., New York, 1998.
• Cotta, R. M., "Benchmark results in computational heat and fluid flow: - integral transform method", International Journal of Heat and Mass Transfer, Vol. 37, pp. 381-394, 1994, Invited Paper.
• Cotta, R. M., "Integral transforms in transient convection:- benchmarks and engineering simulations", In ICHMT International Symposium on Transient Convective Heat Transfer, pp. 433-453, Turkey, August 1996. Invited Keynote Lecture.
• Wolfram, S., "The Mathematica Book", 4th ed., Wolfram Media, Cambridge, 1999.
• Serfaty, R. and Cotta, R. M., "Integral transform solutions of diffusion problems with nonlinear equation coeficients", International Communications in Heat and Mass Transfer, Vol. 17, No. 6, pp. 851-864, 1990.
• Serfaty, R. and Cotta, R. M., "Hybrid analysis of transient nonlinear convection-diffusion problems", International Journal of Numerical Methods in Heat and Fluid Flow, Vol. 2, pp. 55-62, 1992.

Publication Dates

• Publication in this collection
  19 Aug 2002
• Date of issue
  2001

History

• Received
  Jan 2000