Research Article  Open Access
3D Unsteady Diffusion and ReactionDiffusion with Singularities by GFEM with 27Node Hexahedrons
Abstract
The Galerkin Finite Element Method (GFEM) with 8 and 27node hexahedrons elements is used for solving diffusion and transient threedimensional reactiondiffusion 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 [1–3], 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 [6–8] are highlighted. Since the early 50s with [9–13], 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 finiteelement meshes. However, when using methods based on the Taylor series on time, the TaylorGalerkin 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 TaylorGalerkin 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 threedimensional 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 (8nodes) 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 convectiondiffusion using the CranckNicolson 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 convectiondiffusion in a square domain . Compared to two other studies, which use the HOCADI [17] and the PRADI schemes [18], the author shows that his results are more suitable for the same refinements in space and in time.
The transient advectiondiffusiononedimensional response is studied in [19]. The authors use the CrankNicolson Method for the time discretization, whereas for the spatial discretization some finite elements schemes are used, including the StreamlineUpwind PetrovGalerkin (SUPG), the DWG Method (first proposed in [20]), and the GalerkinLeast Squares (GLS), and at the top of the list is the LinkCutting Bubble (LCB) strategy. The authors present applications with advectiondiffusionreaction 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 convectiondiffusiononedimensional 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 threedimensional diffusionreaction by the Galerkin method, with 27node hexahedrons. Also using 8node 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 diffusivereactive threedimensional 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 CranckNicolson 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 GaussLegendre [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 27node 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 27node hexahedrons the results become very satisfactory, where, for the adopted meshes, the accuracy varies from approximately 10^{−3} to 10^{−5}, whereas for the 8node 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.


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.



