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

TheGalerkin Finite ElementMethod (GFEM)with 8and 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.


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][2][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][7][8] are highlighted.Since the early 50s with [9][10][11][12][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 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 (8-nodes) and quadratic elements (20 nodes) are used in the mesh, in which meshes with 1000, 8000, 27000, 64000, and 125000 elements 2 Mathematical Problems in Engineering 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 convectiondiffusion in a square domain [0, 2] × [0, 2].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  1 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 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.

Model Equation
Here, is introduced a numerical study of partial differential equation, which models the transient diffusive-reactive threedimensional phenomenon defined in the Θ = Ξ ⊗ Ω ⊂ R 3 , Ξ ⊂ R, Ω ⊂ R 3 domain, in which Ξ and Ω are limited and closed domains, described as where it is assumed that   ,   ,   = constants ̸ = 0,  = (, , , ),  = (, , ), and  = (, , ) with , , ,  ∈ R and boundary conditions of the first and second type, as well as the initial condition.

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, followed by a set of weight functions V  ( = 1, 2, . . .,  nodes ) and the definition of an inner product (, V  ).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 V  .This study will use the Galerkin method formulation, which considers the weight function to be similar to the interpolation function.

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 T 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 T 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) However, the boundary conditions are considered as being of the first and the second kinds, which mathematically can be written as follows: where Γ  ∪ Γ  = Γ and Γ  ∩ Γ  = 0, Γ 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 T by the function T is as follows: Substituting ( 14) in (13) we have with ,  = 1, . . .,  Nodes .
Since the last term on the right side of (15) carries the contour information and T,  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

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  2 (average error in all domains) as proposed in [23]: . In this equation,  nost is the total number of nodes in the mesh, and   = | (num)  −  (an)  |, where  (num) is the result from the numerical solution and  (an) 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..In this application will be considered the solution of the equation

Application 1: Diffusion Equation
In a domain of unit cube, its analytical solution will be considered as In Tables 1 and 2, it is evident that, for Δ = 0.01 to use 8or 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 Δ = 0001 and Δ = 0.0001, 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  = 1,2, or, 3, with  1 = ,  2 = , and  3 = , 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  nodes and  nos 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.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.

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 Δ = 0.01, 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 Δ = 0.001 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 (, , ) = (1, 1, 1).It may be noted that for the norm  2 (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 Δ = 0.0001, since in Δ = 0001 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.

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,  = 0.72 W/mK, whereas for the same temperature we have regular brick:  = 1920 kg/m 3 and   = 835 J/kg⋅K; cement mortar:  = 1860 kg/m 3 and   = 780 J/kg⋅K.
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 ( = 0.12 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 0 ≤  ≤ 0.12 m, with  =  = 0, toward the  = 0 to  = 0.12).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  = 0.1,  = 0.5, and 0 ≤  ≤ 1.0 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 ( = 0.12) toward the inner face ( = 0).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  = 0.5.Thus, in Figure 3, it can be noted that the maximum temperature of this line rotates around  = 14.3 ∘ C.
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  = 0.5 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  = 0.5 (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  = 0.5.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  = 0.5 is around  = 14.3 ∘ C. 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.(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 (, ) = (0.12, 1).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.

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.
in which  Nodes is the number of nodes in each element,   are the values of the functions in the element nodes, and 1] and Φ  ,  = 1, . . ., 6, are defined from the form

Table 1 :
Norm  ∞ from the error made by the GFEM  at  = 0.1.

Table 2 :
Norm  2 from the error made in  by the GFEM in  = 0.1.

Table 3 :
Norm  ∞ from the error made at   ≅   ≅   by the GFEM/FEM in  = 0.1.

Table 4 :
Norm  ∞ from the error made in   ≅   ≅   by the GFEM/FDM in t = 0.1.

Table 5 :
Norm  2 from the error made in   ≅   ≅   by the GFEM/FEM in  = 0.1.

Table 6 :
Norm  2 from the error made in   ≅   ≅   by the GFEM/FDM in  = 0.1.

Table 7 :
Norm  ∞ from the error made in  by the GFEM in  = 0.1.

Table 8 :
Norm  2 from the error made in  by the GFEM in  = 0.1.

Table 9 :
Norm  ∞ from the error made in   ≅   ≅   in  = 0.1.

Table 10 :
Norm  2 from the error made in   ≅   ≅   in  = 0.1.