Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 560492, 12 pages
Research Article

3D Unsteady Diffusion and Reaction-Diffusion with Singularities by GFEM with 27-Node Hexahedrons

Basic and Environmental Sciences Department, Engineering School of Lorena, University of São Paulo, 12602-810 Lorena, SP, Brazil

Received 16 April 2014; Revised 2 July 2014; Accepted 11 July 2014; Published 13 August 2014

Academic Editor: Weizhong Dai

Copyright © 2014 Estaner Claro Romão. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The Galerkin Finite Element Method (GFEM) with 8- and 27-node hexahedrons elements is used for solving diffusion and transient three-dimensional reaction-diffusion with singularities. Besides analyzing the results from the primary variable (temperature), the finite element approximations were used to find the derivative of the temperature in all three directions. This technique does not provide an order of accuracy compatible with the one found in the temperature solution; thereto, a calculation from the third order finite differences is proposed here, which provide the best results, as demonstrated by the first two applications proposed in this paper. Lastly, the presentation and the discussion of a real application with two cases of boundary conditions with singularities are proposed.

1. Introduction

The vast majority of physical problems are either ruled or represented by partial differential equations. Some mathematical methods are able to produce analytical solutions for physical problems, such as heat transfer [13], however for only a few and very simplified problems. Therefore, the numerical methods become essential for solving heat transfer problems. For decades, numerical methods have been used to solve such problems; among them the Finite Difference Method [4], the Finite Volume [5], and the Finite Element [68] are highlighted. Since the early 50s with [913], the Finite Element Method has been used with great success in various branches of engineering.

In [14] the authors present a Finite Element Method to solve transient problems governed by linear and nonlinear equations with dominant advective terms, owing to the many difficulties in numerical simulation of transient and permanent problems, with the dominant advective term present in the application of the classic Galerkin method with linear elements. The authors use different methods for solving such problems. The first is the generalized Galerkin method, which provides excellent results due to a correct relationship between the spatial and the temporal variations, expressed by the theory of characteristics. The authors also mention the time discretization methods, among them the Explicit Euler Method, which does not produce good results, when the problem is evaluated by unstructured finite-element meshes. However, when using methods based on the Taylor series on time, the Taylor-Galerkin method presents significant advantages; among them are the simple implementation and the third order precision in nonlinear problems. Therefore, the Least Squares Method employed contains the simplicity of the Taylor-Galerkin method and the unconditional stability of the method of characteristics; however, its accuracy is compromised for the Courant numbers greater than the unity.

It is presented in [15] the numerical solution for a three-dimensional transient heat conduction case, in which several temporal discretization methods are tested, with emphasis on the method where the value is admitted as 0, 0.5, 0.75, and 1, in addition to using the Least Squares Method (LSM). In the spatial discretization by the Finite Element Method, hexahedrons linear elements (8-nodes) and quadratic elements (20 nodes) are used in the mesh, in which meshes with 1000, 8000, 27000, 64000, and 125000 elements are compared. Important results are presented in this paper; however, the author compares his results neither with an analytical solution nor with the numerical results from other studies in the open literature.

Two years later, [16] presents a study aimed at solving the 2D transient convection-diffusion using the Cranck-Nicolson method and the high order Padé ADI method for the time discretization and a fourth order scheme for the space discretization. The author validates its formulation by applying it to solve a case of a Gaussian pulse transient convection-diffusion in a square domain . Compared to two other studies, which use the HOC-ADI [17] and the PR-ADI schemes [18], the author shows that his results are more suitable for the same refinements in space and in time.

The transient advection-diffusion-one-dimensional response is studied in [19]. The authors use the Crank-Nicolson Method for the time discretization, whereas for the spatial discretization some finite elements schemes are used, including the Streamline-Upwind Petrov-Galerkin (SUPG), the DWG Method (first proposed in [20]), and the Galerkin-Least Squares (GLS), and at the top of the list is the Link-Cutting Bubble (LCB) strategy. The authors present applications with advection-diffusion-reaction and 1D transient pollutants dispersion. The methods show good results for refined meshes, with particular emphasis on the LCB, which generally provides the best results.

In [21], the authors use a Finite Element Method with the semidiscrete element to numerically solve transient convection-diffusion-one-dimensional problems. The results are compared with the problem’s analytical solution and the finite difference streamline diffusion method. It provides excellent results, even with the most dominant convective problem. In this case, the authors test certain values for the conductivity term, leaving the fixed and the unitary as the constant following the convective term.

This paper proposes a study of the numerical solution of the transient three-dimensional diffusion-reaction by the Galerkin method, with 27-node hexahedrons. Also using 8-node hexahedrons, two numerical applications are presented and discussed, thus, emphasizing and validating the effectiveness of such a proposal. Lastly, an actual application is presented.

2. Model Equation

Here, is introduced a numerical study of partial differential equation, which models the transient diffusive-reactive three-dimensional phenomenon defined in the , , domain, in which and are limited and closed domains, described as where it is assumed that , , , and with and boundary conditions of the first and second type, as well as the initial condition.

3. Method of Weighted Residuals

The goal here is to use the Method of Weighted Residuals to obtain an approximate solution to the differential equation (1).

By introducing the following set of functions, for each element, in the form, in which is the number of nodes in each element, are the values of the functions in the element nodes, and are the interpolation functions. Since (2) is an approximate solution, by substituting (2) in (1), it does not completely satisfy the governing differential equation. Thus define a residual as

followed by a set of weight functions and the definition of an inner product . Determine that this inner product is equal to zero; that is,

It is the equivalent of forcing the approximation error of the differential equation, on the average, to be equal to zero. There are several ways of choosing the weight functions . This study will use the Galerkin method formulation, which considers the weight function to be similar to the interpolation function.

4. Formulation by the Galerkin Method

Equation (1) can be written as

In this work, the time discretization equation (5) is realized with the Cranck-Nicolson method [8] approximating the function by the function as follows:

Following what is defined in (4) and in (6), it yields the following integral:

For the extension integration (only for analysis of this extension will be replaced by , without loss of generality),

Integration by parts will be used as defined on page 153, of [6]; thus, the section described in (8) can be written as

However, the boundary conditions are considered as being of the first and the second kinds, which mathematically can be written as follows: where and , represent the boundary, , , and represent the directions cosines, is the heat transfer coefficient, is the environmental temperature, and is the flux concentration at the boundary. Now, by applying (11) to (9), we will have

Then substituting (12) in (7) yields

Now, the function approximation of by the function is as follows:

Substituting (14) in (13) we have with .

Since the last term on the right side of (15) carries the contour information and is already known from the previous step in the time “”, we have the following system: in which

The numerical integrations are performed using the Gauss-Legendre [8], such that the coordinates transformation from global to local will be required (,  ,  ), where (17) and (19) are as follows: where the Jacobian is defined as ,    and  ,  , are defined from the form in which

Further details of the above equations can be found in [22, pages 43 and 44].

5. Numerical Applications

Here we propose the numerical solution of (1) by solving the linear system (16) which represents it. For this, a computer code was developed in FORTRAN language, and the numerical results were compared with the analytical solution (applications 1 and 2) through the norm (maximum error committed) and (average error in all domains) as proposed in [23]: . In this equation, is the total number of nodes in the mesh, and , where is the result from the numerical solution and is the result form the analytical solution, respectively.

Importantly, the construction of the meshes from the following applications, , , and , is the edges of the regular hexahedron in the , , and directions, respectively. It will be considered in the 1 and 2 applications.

5.1. Application 1: Diffusion Equation

In this application will be considered the solution of the equation

In a domain of unit cube, its analytical solution will be considered as

In Tables 1 and 2, it is evident that, for to use 8- or 27-node hexahedrons, there is no significant difference in their numerical results for the solution; however, even with the mesh’s spatial refinement, the results are not satisfactory. In turn, with and , using 27-node hexahedrons the results become very satisfactory, where, for the adopted meshes, the accuracy varies from approximately 10−3 to 10−5, whereas for the 8-node hexahedron the same is not true. Importantly for the engineering applications, for example, in the thermal area, there is not much point reviewing a precision of several decimal places, thus becoming justifiable the use of elements with at most 27 nodes (quadratic) in this work.

Table 1: Norm from the error made by the GFEM at .
Table 2: Norm from the error made in by the GFEM in .

In other studies regarding the GFEM use, the following technique was used for calculating the derivatives in the three spatial directions for [24, 25]: where , or, 3, with ,  , and , with unsatisfactory results, as it can also be verified in Tables 3, 4, 5, and 6. This work presents an alternative to the calculation of these derivatives. The application of the third order Finite Differences Method is as follows [26]: where and are the number of nodes in the regular mesh in the and directions, respectively. The three previous expressions represent the calculation of the three derivatives, using the third order forward finite differences method. Since the mesh proposed in this paper consists of regular hexahedrons elements, there are points in the mesh, where one or more of these three expressions cannot be applied; in these cases, the third order backward finite differences were used as follows: where other derivatives are similar.

Table 3: Norm from the error made at by the GFEM/FEM in .
Table 4: Norm from the error made in by the GFEM/FDM in t = 0.1.
Table 5: Norm from the error made in by the GFEM/FEM in .
Table 6: Norm from the error made in by the GFEM/FDM in .

Tables 36 show that this proposal to use the GFEM/FDM (Galerkin Finite Element Method/Finite Difference Method) for the calculation of the first derivatives generally shows an order (or two) of accuracy higher than using FEM for the same calculation.

5.2. Application 2: Reaction-Diffusion

For this case, the following equation is taken as model: In a domain of unit cube, its analytical solution will be considered as

Unlike the application 1, generally, the results for this reaction-diffusion case are significantly better; for example, for a , we have that for the 8-node hexahedrons the accuracy is around 10−4 (three orders of accuracy higher than in application 1), while, for the 27-node hexahedrons, the accuracy is around 10−6 (five orders of precision higher than in application 1), for such, just compare Tables 1 and 7. Tables 7, 8, 9, and 10 show poor results for the and 8-node hexahedrons, thus, demonstrating that, in this case, a greater refinement in time and 8-node hexahedrons are not a good combination. It is noteworthy that, after a thorough analysis of the results, it was noted that the results are bad in the norm . With the analysis of the largest error made in all the mesh, where it generally occurred at the internal node closest to the vertex . It may be noted that for the norm (average error in all the mesh) the numerical error is not so significant. However, for this same mesh, poor results were found in the solutions of the first derivatives of temperature, something expected, since their calculations depend directly on the results obtained from the calculation of temperature. Unlike the application 1, no analyses are performed here for the results of , since in and 27-node hexahedrons a numerical accuracy on the order of 10−8 can already be found, which is extremely suitable for application in engineering problems of heat transfer.

Table 7: Norm from the error made in by the GFEM in .
Table 8: Norm from the error made in by the GFEM in .
Table 9: Norm from the error made in in .
Table 10: Norm from the error made in in .
5.3. Application 3

In this application, we assume a situation of heat transfer by conduction in a regular brick wall, 10 cm thick with cement mortar (1 cm thickness). According to Table A.3 from [27], for both the regular brick and the cement mortar at a temperature of 300 K,  W/mK, whereas for the same temperature we haveregular brick:  kg/m3 and  J/kg·K;cement mortar:  kg/m3 and  J/kg·K.

For this application we will use a weighted mean to determine the value of the specific mass () and of the specific heat () as follows:

By the heat conduction equation described in (2.13), page 43 of [27], it follows that  W/mK and . Therefore, the governing equation of this application will be this way:

Thus, we assume a rectangular wall with thickness of 0.12 m ( m), width 1.00 m ( m), and a height of 1.00 m ( m).

Let us suppose a challenging situation where the initial temperature of the wall is 12°C (early in the morning), and because of the sun’s presence, suddenly the outside temperature rises to 20°C; that is, from the first step of time, the conditions will be of 12°C for all the walls, except for the outer wall ( m).

There are some challenges in this application because both the domain and its initial geometric boundaries present a unique region. Figure 1 shows the example.

Figure 1: Outline of boundary conditions near the vertex .

Figure 1 clearly shows that the situation highlighted several possible nodes with a temperature of 12°C and, suddenly, a 20°C temperature (e.g., in line  m, with , toward the to ). This situation can lead to extremely high temperature gradients in this region, particularly in the direction to the lower refining adopted.

Figures 2 to 6 show some variations of temporal and spatial refinements. For this analysis the line where  , , and was adopted as reference, because of its closeness to the singularities. This analysis begins by varying the values of (Figures 2 to 4). Figure 2 shows the mean temperature on that line, where the proposed refinement shows altogether that the average becomes constant at the beginning of the curves. Conversely, in Figure 3 the maximum temperature in this line is studied. It is noteworthy that, for the proposed problem, its physical and geometric characteristics lead us to conclude that the three-dimensional temperature profile is equivalent to a paraboloid entering into the domain through the outer face () toward the inner face (). Therefore, even with our analysis being held in a line of this three-dimensional domain, one can easily conclude that the maximum temperature of this line will be at . Thus, in Figure 3, it can be noted that the maximum temperature of this line rotates around .

Figure 2: Mean temperature for , , and (mesh: , , , 8 nodes), (mesh: , , , 27 nodes), , : number of the nodes in all the mesh.
Figure 3: Temperature in for , ,  and (mesh: , , , 8 nodes) (mesh: , , , 27 nodes), .
Figure 4: Maximum temperature for , ,  and (mesh: , , , 8 nodes) (mesh: , , , 27 nodes), .

From this reasoning, we analyzed the maximum temperature (Figure 4) in this line; it seemed logical after the above reasoning; however, owing to high gradients in the singularities, for the less refined meshes, the maximum temperatures found were higher than those in point , especially in the singularity near . Figure 4 also shows that, with the refinement, the maximum temperature tends the temperature at the point , which should actually be.

Then the mesh is refined in the direction (toward the highest gradients), and the other refinements are fixed (Figure 5). With this proposal we note that the average temperature, the maximum temperature, and the temperature at become stagnant for the first refinements; that is, for the proposed situation, the refinement in the direction brings no better results. Furthermore, the maximum temperature is generally higher than the temperature at (the actual maximum temperature in the line), thus, demonstrating the influence of the singularities in the numerical results. Another way of analyzing this situation is to think that with a smaller more mesh points will have a situation where the results should come out suddenly, from 12°C to 20°C, that is, the lower the the greater the values of in this region.

Figure 5: Temperature for , , and (mesh: , , 8 nodes) (mesh: , , 27 nodes), , .
Figure 6: Temperature for , , and (mesh: , , and , 8 nodes) (mesh: , , and , 27 nodes), , .

Lastly, the temporal refinement is a variable , as shown in Figure 6. In this situation, we note that the temporal refinement is not enough to reduce the maximum temperature down to the temperature value at . However, it is important to emphasize that the use of 27-node hexahedrons presents smaller errors in the maximum temperature compared to the 8-node hexahedrons in the less refined meshes; notwithstanding, their results were comparable after the refinement. Significantly, the average temperature becomes constant with the temporal refinement and again the temperature at is around . Furthermore, it is expected that, with the lower , a larger amount of steps in time existed, and thus, the significant errors in the singularity areas are transmitted to the remaining proposed three-dimensional domain. Figure 7 shows a temperature profile and its derivative in . The numerical oscillations presented in occur near the singularities.

Figure 7: Profile of (left) and (right) for the plane with for a mesh with 27-node hexahedrons, , , , and in .

To demonstrate a more realistic engineering practice, let us assume, for the same problem again, that the initial condition is 12°C, but the following will be taken as boundary condition.(i)In (inner wall), (assuming that the internal coolant is able to maintain a constant temperature).(ii)In  m (outer wall) the expression   (°C/s) is taken.(iii)The other sides will be regarded as geometric symmetry, that is, equivalent to the insulation.

Figure 8 shows the temperature profile in a plane, at the 900 s and 2700 s instants of time. In the second case, a slight numerical oscillation was noted in the region of () = (). It is important to mention that the mesh adopted was not as expressively refined. Furthermore, the numerical results of the temperature variation in the direction (Figure 9) present oscillations with a little more expression. To this end, it led to the construction of a somewhat more refined mesh to analyze whether or not these oscillations existed. These results are shown in Figure 10, and it is noted that, at the 2700 s moment of time, in both the temperature solution and its variation toward , numerical oscillations no longer occur, thus, demonstrating the efficiency of GFEM with 27-node hexahedrons with a mesh not significantly refined, since it contains only 7225 nodes.

Figure 8: Temperature profile (°C) in the plane with in left and right (mesh: , , , , 27 nodes).
Figure 9: profile (°C/m) in the plane with in left and right (mesh: , , , , 27 nodes).
Figure 10: Profile of (left) and (right) for the plane with for a mesh with , , , and in .

6. Conclusions

The Galerkin Finite Element Method (GFEM) showed significant efficiency in the numerical solution of the problems presented, with special attention to the situations where the mesh was successfully built with 27-node hexahedrons. The presentation of a realistic and a challenging application demonstrated that even in relatively complex situations, in relation to the difficulty of the numerical solution, the method proposed in this work had a good performance. Importantly, the use of FDM in this study was restricted to structured grids, with an objective this author the use of FDM for the unstructured grids to calculate derivatives.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.


The National Council of Scientific Development and Technology, CNPq, Brazil, (Proc. 408250/2013-5) and FAPESP (Proc. 2014/06679-8) supported the present work.


  1. V. S. Arpaci, Conduction Heat Transfer, Addison-Wesley, Reading, Mass, USA, 1966.
  2. A. Bejan, Transferência de Calo, Editora Edgard Blücher Ltda, 1996.
  3. H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Claredon Press, Oxford, UK, 2nd edition, 1986.
  4. G. D. Smith, Numerical Solution of Partial Differential Equations, Oxford Mathematical Handbooks, New York, NY. USA, 1965. View at MathSciNet
  5. T. J. Chung, Computational Fluid Dynamics, Cambridge University Press, Cambridge, UK, 2002. View at Publisher · View at Google Scholar · View at MathSciNet
  6. R. W. Lewis, P. Niyhiarasu, and K. N. Seetharamu, Fundamentals of the Finite Element Method for Heat and Fluid Flow, John Wiley & Sons, 2004.
  7. J. Donea and A. Huerta, Finite Element Methods for Flow Problems, John Wiley & Sons, 2003.
  8. J. N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill, 2nd edition, 1993.
  9. M. J. Turner, R. W. Clough, H. C. Martin, and L. P. Topp, “Stiffness and deflection analysis of complex structures,” Journal of the Aeronautical Sciences, vol. 23, no. 9, pp. 805–823, 1956. View at Google Scholar
  10. R. W. Clough, “The finite element method in plane stress analysis,” in Proceedings of the 2nd Conference on Electronic Computation, pp. 345–378, Pittsburg, Pa, USA, 1960.
  11. J. H. Argyris, Recent Advances in Matrix Methods of Structural Analysis, vol. 4 of Progress in Aeronautical Science, Pergamon Press, Elmsford, NY, USA, 1963.
  12. O. C. Zienkiewicz and Y. K. Cheung, “Finite element in the solution of field problems,” The Engineer, p.p. 507–510, 1965.
  13. J. T. Oden and L. C. Wellford Jr., “Analysis of flow of viscous fluids by the finite element method,” AIAA Journal, vol. 10, no. 12, pp. 1590–1599, 1972. View at Publisher · View at Google Scholar · View at Scopus
  14. J. Donea and L. Quartapelle, “An introduction to finite element methods for transient advection problems,” Computer Methods in Applied Mechanics and Engineering, vol. 95, no. 2, pp. 169–203, 1992. View at Publisher · View at Google Scholar · View at Scopus
  15. M. R. Vujicic and S. G. R. Brown, “Iterative solvers in the finite element solution of transient heat conduction,” FME Transactions, vol. 32, pp. 61–68, 2004. View at Google Scholar
  16. D. You, “A high-order Padé method for unsteady convection-diffusion equations,” Journal of Computational Physics, vol. 214, no. 1, pp. 1–11, 2006. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. S. Karaa and J. Zhang, “High order ADI method for solving unsteady convection-diffusion problems,” Journal of Computational Physics, vol. 198, no. 1, pp. 1–9, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  18. D. W. Peaceman and H. H. Rachford Jr., “The numerical solution of parabolic and elliptic differential equations,” Journal of the Society for Industrial and Applied Mathematics, vol. 3, no. 1, pp. 28–41, 1959. View at Google Scholar
  19. M. I. Asensio, B. Ayuso, and G. Sangalli, “Coupling stabilized finite element methods with finite difference time integration for advection-diffusion-reaction problems,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 35-36, pp. 3475–3491, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  20. J. Douglas. Jr. and J. P. Wang, “An absolutely stabilized finite element method for the Stokes problem,” Mathematics of Computation, vol. 52, pp. 495–508, 1989. View at Publisher · View at Google Scholar
  21. Z. Si, X. Feng, and A. Abduwali, “The semi-discrete streamline diffusion finite element method for time-dependented convection-diffusion problems,” Applied Mathematics and Computation, vol. 202, no. 2, pp. 771–779, 2008. View at Publisher · View at Google Scholar · View at Scopus
  22. G. Dhatt and G. Touzot, The Finite Element Method Displayed, John Wiley & Sons, 1984.
  23. M. Zlámal, “Superconvergence and reduced integration in the finite element method,” Mathematics of Computation, vol. 32, pp. 663–685, 1978. View at Publisher · View at Google Scholar · View at MathSciNet
  24. E. C. Romão, M. D. Campos, and L. F. M. Moura, “Application of the Galerkin and least-squares finite element methods in the solution of 3D Poisson and Helmholtz equations,” Computers & Mathematics with Applications, vol. 62, no. 11, pp. 4288–4299, 2011. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  25. E. C. Romão and L. F. M. de Moura, “Galerkin and least squares methods to solve a 3D convection-diffusion-reaction equation with variable coefficients,” Numerical Heat Transfer, Part A: Applications, vol. 61, no. 9, pp. 669–698, 2012. View at Publisher · View at Google Scholar · View at Scopus
  26. E. C. Romão, J. C. Z. Aguilar, M. D. de Campos, and L. F. M. de Moura, “Central difference method of O(x6) in solution of the CDR equation with variable coefficients and robin condition,” International Journal of Applied Mathematics, vol. 25, no. 1, pp. 139–153, 2012. View at Google Scholar
  27. F. P. Incropera and D. P. DeWitt, Fundamentos de Transferência de Calor e Massa, Quinta Edição, Editora LTC, 2003, (Portuguese).