Abstract
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 [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 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.
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.
Tables 3–6 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.
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 have regular 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 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 .
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.
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.
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.
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.
Acknowledgment
The National Council of Scientific Development and Technology, CNPq, Brazil, (Proc. 408250/2013-5) and FAPESP (Proc. 2014/06679-8) supported the present work.