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

Authors: Estaner Claro Romão
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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, thenumerical 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 theCranck-Nicolson method and thehigh 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 [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 P1 to numerically solve transient convection-diffusion-one-dimensional problems. The results are compared with the problem’s analytical solution and thefinite 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Θ=Ξ⊗Ω⊂R3, Ξ⊂R, Ω⊂R3 domain, in which Ξ and Ω are limited and closed domains, described as
(1)ψ∂T∂t+kx∂2T∂x2+ky∂2T∂y2+kz∂2T∂z2+BT=0,
where it is assumed that kx,ky,kz=constants≠0, T=T(x,y,z,t), B=B(x,y,z), and ψ=ψ(x,y,z) with x,y,z,t∈R 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 theMethod 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,(2)T≈T-=∑i=1NnodesTiNi,
in which NNodes is the number of nodes in each element, Ti are the values of the functions in the element nodes, and Ni 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 R as
(3)R≈ψ∂T-∂t+kx∂2T-∂x2+ky∂2T-∂y2+kz∂2T-∂z2+BT-,followed by a set of weight functionsvi(i=1,2,…,Nnodes) and the definition of an inner product (R,vi). Determine that this inner product is equal to zero; that is,
(4)∫ΩRvidΩ=0.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 functionsvi. 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
(5)ψ∂T∂t=-kx∂2T∂x2-ky∂2T∂y2-kz∂2T∂z2-BT.In this work, the time discretization equation (5) is realized with the Cranck-Nicolson method [8] approximating the T function by the T~ function as follows:
(6)ψT~s+1-T~sΔt=α·[-kx∂2T~∂x2-ky∂2T~∂y2-kz∂2T~∂z2-BT~]s+1+(1-α)·[-kx∂2T~∂x2-ky∂2T~∂y2-kz∂2T~∂z2-BT~]s⟹{ψT~Δt+α·[kx∂2T~∂x2+ky∂2T~∂y2+kz∂2T~∂z2+BT~]}s+1+{-ψT~Δt+(1-α)·[kx∂2T~∂x2+ky∂2T~∂y2+kz∂2T~∂z2+BT~]}s=0.Following what is defined in (4) and in (6), it yields the following integral:
(7)∫Ω{ψT~Δt+α·[kx∂2T~∂x2+ky∂2T~∂y2+kz∂2T~∂z2+BT~]}s+1NjdΩ+∫Ω{[kx∂2T~∂x2+ky∂2T~∂y2+kz∂2T~∂z2+BT~]-ψT~Δt+(1-α)·[kx∂2T~∂x2+ky∂2T~∂y2+kz∂2T~∂z2+BT~]}sNjdΩ=0.For the extension integration (only for analysis of this extensionT~ will be replaced by T, without loss of generality),
(8)∫Ω[kx∂2T∂x2+ky∂2T∂y2+kz∂2T∂z2]·NjdΩ.Integration by parts will be used as defined on page 153, of [6]; thus, the section described in (8) can be written as
(9)∫Ω[kx∂2T∂x2+ky∂2T∂y2+kz∂2T∂z2]·NjdΩ=-∫Ω[kx∂Ni∂x∂T∂x+ky∂Ni∂y∂T∂y+kz∂Ni∂z∂T∂z]dΩ+∫ΓqNikx∂T∂xldΓq+∫ΓqNiky∂T∂ymdΓq+∫ΓqNikz∂T∂zndΓq.However, the boundary conditions are considered as being of the first and the second kinds, which mathematically can be written as follows:(10)T=TbinΓb,(11)kx∂T∂xl+ky∂T∂ym+kz∂T∂zn+q+h(T-Ta)=0inΓq,
where Γb∪Γq=Γ and Γb∩Γq=0, Γ represent the boundary, l, m, and n represent the directions cosines, h is the heat transfer coefficient, Ta is the environmental temperature,and q is the flux concentration at the boundary. Now, by applying (11) to (9), we will have
(12)∫Ω[kx∂2T∂x2+ky∂2T∂y2+kz∂2T∂z2]·NjdΩ=-∫Ω[kx∂Ni∂x∂T∂x+ky∂Ni∂y∂T∂y+kz∂Ni∂z∂T∂z]dΩ-∫ΓqNj(q+h(T-Ta))dΓq.Then substituting (12) in (7) yields
(13)∫Ω{α·[-kx∂T~∂x∂Nj∂x-ky∂T~∂y∂Nj∂y-kz∂T~∂z∂Nj∂z]}s+1dΩ+∫Ω{ψT~Δt+α·[BT~]}s+1NjdΩ+∫Ω{[-kx∂T~∂x∂Nj∂x-ky∂T~∂y∂Nj∂y-kz∂T~∂z∂Nj∂z](1-α)·[-kx∂T~∂x∂Nj∂x-ky∂T~∂y∂Nj∂y-kz∂T~∂z∂Nj∂z]}sdΩ+∫Ω{-ψT~Δt+α·[BT~]}sNjdΩ=∫ΓqNj(q+h(T-Ta))dΓq.Now, the function approximation ofT~ by the function T^e is as follows:
(14)T~≅T^e=∑i=1NNodesNiT^ie.Substituting (14) in (13) we have
(15)∫Ωe{ψΔtNiNj+α·[-kx∂Ni∂x∂Nj∂x-ky∂Ni∂y∂Nj∂y-kz∂Ni∂z∂Nj∂z+BNiNj]}dΩ·T^is+1,e+∫Ωe{-ψΔtNiNj+(1-α)·[-kx∂Ni∂x∂Nj∂x-ky∂Ni∂y∂Nj∂y-kz∂Ni∂z∂Nj∂z+BNiNj]}dΩ·T^is,e=∫ΓqNj(q+h(T-Ta))dΓq,
with i,j=1,…,NNodes.Since the last term on the right side of (15) carries the contour information and T^is,e is already known from the previous step in the time “s”, we have the following system:
(16)[K]{Tie,s+1}={F}-[H]{Tie,s}
in which
(17)Kij=-∫Ωeα·kx∂Ni∂x∂Nj∂xdΩ-∫Ωeα·ky∂Ni∂y∂Nj∂ydΩ-∫Ωeα·kz∂Ni∂z∂Nj∂zdΩ+∫Ωe(ψΔt+αB)NiNjdΩ,(18)Fi=∫ΓqNj(q+h(T-Ta))dΓq,(19)Hij=∫Ωe(1-α)·kx∂Ni∂x∂Nj∂xdΩ+∫Ωe(1-α)·ky∂Ni∂y∂Nj∂ydΩ+∫Ωe(1-α)·kz∂Ni∂z∂Nj∂zdΩ-∫Ωe(-ψΔt+(1-α)B)NiNjdΩ.The numerical integrations are performed using the Gauss-Legendre [8], such that the coordinates transformation from global to local will be required (x→ξ,  y→η,  z→ζ), where (17) and (19) are as follows:
(20)Kij=-∫Ω-eα·kxdet⁡(J)Φ1Φ4dΩ--∫Ωeα·kydet⁡(J)Φ2Φ5dΩ--∫Ω-eα·kzdet⁡(J)Φ3Φ6dΩ-+∫Ω-e(ψΔt+αB)NiNjdΩ-,Hij=∫Ω-e(1-α)·kxdet⁡(J)Φ1Φ4dΩ-+∫Ω-e(1-α)·kydet⁡(J)Φ2Φ5dΩ-+∫Ω-e(1-α)·kzdet⁡(J)Φ3Φ6dΩ--∫Ω-e(-ψΔt+(1-α)B)NiNjdΩ-,
where the Jacobian is defined as
(21)det⁡(J)=∂x∂ξ(∂y∂η∂z∂ζ-∂z∂η∂y∂ζ)+∂x∂η(∂z∂ξ∂y∂ζ-∂y∂ξ∂z∂ζ)+∂x∂ζ(∂y∂ξ∂z∂η-∂z∂ξ∂y∂η)≠0,dΩ-=dζdηdξ,  Ω-=[-1,1]×[-1,1]×[-1,1]  and  Φi,  i=1,…,6, are defined from the form
(22)Φ1=A1∂Nie∂ξ+A2∂Nie∂η+A3∂Nie∂ζ,Φ2=A4∂Nie∂ξ+A5∂Nie∂η+A6∂Nie∂ζ,Φ3=A7∂Nie∂ξ+A8∂Nie∂η+A9∂Nie∂ζ,Φ4=A1∂Nje∂ξ+A2∂Nje∂η+A3∂Nje∂ζ,Φ5=A4∂Nje∂ξ+A5∂Nje∂η+A6∂Nje∂ζ,Φ6=A7∂Nje∂ξ+A8∂Nje∂η+A9∂Nje∂ζ,
in which
(23)A1=∂y∂η∂z∂ζ-∂z∂η∂y∂ζ,A2=∂x∂ζ∂z∂η-∂x∂η∂z∂ζ,A3=∂x∂η∂y∂ζ-∂x∂ζ∂y∂η,A4=∂z∂ξ∂y∂ζ-∂y∂ξ∂z∂ζ,A5=∂x∂ξ∂z∂ζ-∂x∂ζ∂z∂ξ,A6=∂y∂ξ∂x∂ζ-∂y∂ζ∂x∂ξ,A7=∂y∂ξ∂z∂η-∂z∂ξ∂y∂η,A8=∂x∂η∂z∂ξ-∂z∂η∂x∂ξ,A9=∂x∂ξ∂y∂η-∂x∂η∂y∂ξ.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 L∞ (maximum error committed) and L2 (average error in all domains) as proposed in [23]: ∥e∥=[(∑i=1Nnostei2)/Nnost]1/2. In this equation, Nnost is the total number of nodes in the mesh, and ei=|T(num)i-T(an)i|, where T(num) is the result from the numerical solution and T(an) is the result form the analytical solution, respectively.Importantly, the construction of the meshes from the following applications,Δx, Δy, and Δz, is the edges of the regular hexahedron in the x, y, and z directions, respectively. It will be considered h=Δx=Δy=Δz in the 1 and 2 applications.


### 5.1. Application 1: Diffusion Equation

In this application will be considered the solution of the equation(24)∂2T∂x2+∂2T∂y2+∂2T∂z2=∂T∂t.In a domain of unit cube, its analytical solution will be considered as(25)T(x,y,z,t)=100e-3π2tsin⁡(π·x)sin⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂x=100π·e-3π2tcos⁡⁡(π·x)sin⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂y=100π·e-3π2tsin⁡(π·x)cos⁡⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂z=100π·e-3π2tsin⁡(π·x)sin⁡(π·y)cos⁡⁡(π·z).In Tables1 and 2, it is evident that, for Δt=0.01 to use 8- or 27-node hexahedrons, there is no significant difference in their numerical results for the T(x,y,z) solution; however, even with the mesh’s spatial refinement, the results are not satisfactory. In turn, with Δt=0001 and Δt=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.Table 1
NormL∞ from the error made by the GFEM T at t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
4.53
E
-
01
1.14
E
-
01
3.42
E
-
01
2.59
E
-
03
3.41
E
-
01
1.60
E
-
03
1/8
3.06
E
-
01
1.13
E
-
01
1.95
E
-
01
1.58
E
-
03
1.94
E
-
01
5.42
E
-
04
1/10
2.37
E
-
01
1.12
E
-
01
1.26
E
-
01
1.30
E
-
03
1.24
E
-
01
2.34
E
-
04
1/16
1.61
E
-
01
1.12
E
-
01
5.01
E
-
02
1.14
E
-
03
1.18
E
-
01
4.60
E
-
05Table 2
NormL2 from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.27
E
-
01
3.57
E
-
02
9.61
E
-
02
8.19
E
-
04
9.58
E
-
02
4.72
E
-
04
1/8
9.08
E
-
02
3.64
E
-
02
5.78
E
-
02
5.14
E
-
04
5.75
E
-
02
1.56
E
-
04
1/10
7.27
E
-
02
3.69
E
-
02
3.86
E
-
02
4.31
E
-
04
3.82
E
-
02
6.76
E
-
05
1/16
5.20
E
-
02
3.79
E
-
02
1.62
E
-
02
3.87
E
-
04
1.58
E
-
02
1.38
E
-
05In other studies regarding the GFEM use, the following technique was used for calculating the derivatives in the three spatial directions forT(x,y,z) [24, 25]:
(26)∂T∂xk≈∂T^e∂xk=∑i=1Nnodes∂Ni∂xkT^ie,
where k=1,2, or, 3, with x1=x,  x2=y, and x3=z, 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]:
(27)∂Ti∂x≅-11Ti+18Ti+1-9Ti+2+2Ti+36Δx,∂Ti∂y≅-11Ti+18Ti+Nnodesx-9Ti+2*Nnodesx+2Ti+3*Nnodesx6Δy,∂Ti∂z≅(-11Ti+18Ti+Nnodesx*NNosy-9Ti+2*Nnodesx*NNosy+2Ti+3*Nnodesx*NNosy)×(6Δy)-1,
where Nnodesx and Nnosy are the number of nodes in the regular mesh in the x and y directions, respectively. The three previous expressions represent the calculation of the T three derivatives, using the third orderforward 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 orderbackward finite differences were used as follows:
(28)∂Ti∂x≅11Ti-18Ti-1+9Ti-2-2Ti-36Δx,
where other derivatives are similar.Table 3
NormL∞ from the error made at Tx≅Ty≅Tz by the GFEM/FEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
4.33
E
-
00
5.22
E
-
01
4.24
E
-
00
3.54
E
-
01
4.24
E
-
00
3.57
E
-
01
1/8
3.25
E
-
00
4.49
E
-
01
3.19
E
-
00
2.01
E
-
01
3.18
E
-
00
2.04
E
-
01
1/10
2.60
E
-
00
4.14
E
-
01
2.55
E
-
00
1.28
E
-
01
2.55
E
-
00
1.31
E
-
01
1/16
1.65
E
-
00
3.76
E
-
01
1.59
E
-
00
4.84
E
-
02
1.59
E
-
00
5.19
E
-
02Table 4
NormL∞ from the error made in Tx≅Ty≅Tz by the GFEM/FDM in t = 0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.12
E
-
00
3.33
E
-
01
7.66
E
-
01
7.75
E
-
02
7.62
E
-
01
7.84
E
-
02
1/8
8.59
E
-
01
3.45
E
-
01
5.08
E
-
01
3.28
E
-
02
5.04
E
-
01
3.35
E
-
02
1/10
7.01
E
-
01
3.49
E
-
01
3.51
E
-
01
1.70
E
-
02
3.47
E
-
01
1.73
E
-
02
1/16
4.99
E
-
01
3.51
E
-
01
1.50
E
-
01
5.24
E
-
03
5.16
E
-
01
4.25
E
-
03Table 5
NormL2 from the error made in Tx≅Ty≅Tz by the GFEM/FEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.26
E
-
00
1.20
E
-
01
1.23
E
-
00
1.01
E
-
01
1.23
E
-
00
1.02
E
-
01
1/8
9.86
E
-
01
1.11
E
-
01
9.59
E
-
01
5.74
E
-
02
9.59
E
-
01
5.79
E
-
02
1/10
8.13
E
-
01
1.12
E
-
01
7.88
E
-
01
3.67
E
-
02
7.88
E
-
01
3.72
E
-
02
1/16
5.39
E
-
01
1.18
E
-
01
5.16
E
-
01
1.42
E
-
02
1.47
E
-
01
1.45
E
-
02Table 6
NormL2 from the error made in Tx≅Ty≅Tz by the GFEM/FDM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
3.21
E
-
01
1.11
E
-
01
2.15
E
-
01
2.31
E
-
02
2.14
E
-
01
2.36
E
-
02
1/8
2.67
E
-
01
1.18
E
-
01
1.55
E
-
01
9.78
E
-
03
1.54
E
-
01
1.00
E
-
02
1/10
2.28
E
-
01
1.20
E
-
01
1.13
E
-
01
5.08
E
-
03
1.12
E
-
01
5.22
E
-
03
1/16
1.70
E
-
01
1.22
E
-
01
5.12
E
-
02
1.67
E
-
03
5.00
E
-
02
1.30
E
-
03Tables3–6 show that this proposal to use the GFEM/FDM (Galerkin Finite Element Method/Finite Difference Method) for the calculation of the T 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:(29)2∂2T∂x2+2∂2T∂y2+2∂2T∂z2-T=∂T∂t.
In a domain of unit cube, its analytical solution will be considered as
(30)T(x,y,z,t)=et(ex+ey+ez),∂T∂x=etex,∂T∂y=etey,∂T∂z=etez.Unlike the application 1, generally, the results for this reaction-diffusion case are significantly better; for example, for aΔt=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 Δt=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 L∞. With the analysis of the largest error made in all the mesh, where it generally occurred at the internal node closest to the vertex (x,y,z)=(1,1,1). It may be noted that for the norm L2 (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 Δt=0.0001, since in Δt=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.Table 7
NormL∞ from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
1/6
7.26
E
-
04
1.99
E
-
06
7.16
E
-
01
1.86
E
-
06
1/8
3.99
E
-
04
1.40
E
-
06
1.66
E
-
01
6.02
E
-
07
1/10
2.52
E
-
04
1.29
E
-
06
9.36
E
-
02
2.49
E
-
07
1/16
9.75
E
-
05
1.25
E
-
06
3.82
E
-
01
3.84
E
-
08Table 8
NormL2 from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
1/6
2.60
E
-
04
5.82
E
-
07
5.49
E
-
02
5.49
E
-
07
1/8
1.51
E
-
04
4.84
E
-
07
8.81
E
-
03
1.84
E
-
07
1/10
9.91
E
-
05
5.04
E
-
07
4.12
E
-
03
7.87
E
-
08
1/16
4.00
E
-
05
5.35
E
-
07
8.10
E
-
03
1.38
E
-
08Table 9
NormL∞ from the error made in Tx≅Ty≅Tz in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
1/6
2.36
E
-
01
2.85
E
-
03
6.53
E
-
03
4.06
E
-
04
4.53
E
-
00
1.24
E
+
01
6.53
E
-
03
4.23
E
-
04
1/8
1.80
E
-
01
1.42
E
-
03
3.73
E
-
03
1.75
E
-
04
1.54
E
-
00
2.75
E
-
00
3.73
E
-
03
1.85
E
-
04
1/10
1.45
E
-
01
9.24
E
-
04
2.41
E
-
03
9.38
E
-
05
7.91
E
-
01
2.69
E
-
00
2.41
E
-
03
9.69
E
-
05
1/16
9.19
E
-
02
3.87
E
-
04
9.61
E
-
04
2.83
E
-
05
6.02
E
-
00
1.76
E
+
01
9.55
E
-
04
2.44
E
-
05Table 10
NormL2 from the error made in Tx≅Ty≅Tz in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
1/6
1.60
E
-
01
1.93
E
-
03
3.60
E
-
03
2.96
E
-
04
4.72
E
-
01
1.15
E
-
00
3.59
E
-
03
2.96
E
-
04
1/8
1.21
E
-
01
8.49
E
-
04
2.02
E
-
03
1.25
E
-
04
1.54
E
-
01
2.49
E
-
01
2.02
E
-
03
1.25
E
-
04
1/10
9.71
E
-
02
4.65
E
-
04
1.29
E
-
03
6.45
E
-
05
1.10
E
-
01
1.38
E
-
01
1.29
E
-
03
6.40
E
-
05
1/16
6.10
E
-
02
1.52
E
-
04
5.07
E
-
04
1.59
E
-
05
1.80
E
-
01
4.40
E
-
01
5.07
E
-
04
1.55
E
-
05


### 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, k=0.72 W/mK, whereas for the same temperature we haveregularbrick: ρ=1920 kg/m3 and cp=835 J/kg·K;cement mortar: ρ=1860 kg/m3 and cp=780 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 (cp) as follows:
(31)ρm=0.02×1,860+0.1×1.9200.12=1,910Kg/m3,cpm=0.02×780+0.1×8350.12≅825.83J/kg·K.By the heat conduction equation described in (2.13), page 43 of [27], it follows that kx=ky=kz=0.72 W/mK and ψ=-ρcp=-1,910×825.83≅-1,577,341.67J/m3·K. Therefore, the governing equation of this application will be this way:
(32)0.72(∂2T∂x2+∂2T∂y2+∂2T∂z2)=1,577,341.67∂T∂t.Thus, we assume a rectangular wall with thickness of 0.12 m (0≤x≤0.12 m), width 1.00 m (0≤z≤1.00 m), and a height of 1.00 m (0≤y≤1.00 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 (x=0.12 m).There are some challenges in this application because both the domain and its initial geometric boundaries present a unique region. Figure1 shows the example.Figure 1
Outline of boundary conditions near the vertex(x,y,z)=(0.12,0,0).Figure1 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≤x≤0.12 m, with y=z=0, toward the x=0 to x=0.12). This situation can lead to extremely high temperature gradients in this region, particularly in the x direction to the lower refining adopted.Figures2 to 6 show some variations of temporal and spatial refinements. For this analysis the line where  x=0.1, y=0.5, and 0≤z≤1.0 was adopted as reference, because of its closeness to the singularities. This analysis begins by varying the values of Δz (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 (x=0.12) toward the inner face (x=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 z=0.5. Thus, in Figure 3, it can be noted that the maximum temperature of this line rotates around T=14.3°C.Figure 2
Mean temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes), (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s, NNost: number of the nodes in all the mesh.Figure 3
Temperature inz=0.5 for x=0.1, y=0.5,  and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s.Figure 4
Maximum temperature forx=0.1, y=0.5,  and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s.From this reasoning, we analyzed the maximum temperature (Figure4) 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 z=0.5, especially in the singularity near z=1.0. Figure 4 also shows that, with the refinement, the maximum temperature tends the temperature at the point z=0.5, which should actually be.Then the mesh is refined in thex 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 z=0.5 become stagnant for the first refinements; that is, for the proposed situation, the refinement in the x direction brings no better results. Furthermore, the maximum temperature is generally higher than the temperature at z=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 Δx 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 Δx the greater the values of ∂T/∂x in this region.Figure 5
Temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δy=Δz=1/16, 0.12/6≤Δx≤0.12/60, 8 nodes) (mesh: Δy=Δz=1/8, 0.12/3≤Δx≤0.12/30, 27 nodes), Δt=1/20, t=400s.Figure 6
Temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, and Δz=1/16, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, and Δz=1/8, 27 nodes), 2.5s≤Δt≤20s, t=400s.Lastly, the temporal refinement is a variableΔt, 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 z=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 z=0.5 is around T=14.3°C. Furthermore, it is expected that, with the lower Δt, 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 x. The numerical oscillations presented in Tx occur near the singularities.Figure 7
Profile ofT (left) and Tx (right) for the plane XZ with y=0.5 for a mesh with 27-node hexahedrons, Δx=0.12/6, Δy=1/4, Δz=1/40, and Δt=20s in t=400s.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)
Inx=0 (inner wall), T=12°C (assuming that the internal coolant is able to maintain a constant temperature).(ii)
Inx=0.12 m (outer wall) the expression T(x,y,z,t)=12°C+(t/180)  (°C/s) is taken.(iii)
The other sides will be regarded as geometric symmetry, that is, equivalent to the insulation.Figure8 shows the temperature profile in a XZ 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 (x,z) = (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 x 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 x, 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 planeXZ with y=0.5 in t=900s left and t=2700s right (mesh: Δx=0.12/6, Δy=1/4, Δz=1/4, Δt=30s, 27 nodes).Figure 9
T
x profile (°C/m) in the plane XZ with y=0.5 in t=900s left and t=2700s right (mesh: Δx=0.12/6, Δy=1/4, Δz=1/4, Δt=30s, 27 nodes).Figure 10
Profile ofT (left) and Tx (right) for the plane XZ with y=0.5 for a mesh with Δx=0.12/12, Δy=1/8, Δz=1/8, and Δt=30s in t=2700s.


## 5.1. Application 1: Diffusion Equation

In this application will be considered the solution of the equation(24)∂2T∂x2+∂2T∂y2+∂2T∂z2=∂T∂t.In a domain of unit cube, its analytical solution will be considered as(25)T(x,y,z,t)=100e-3π2tsin⁡(π·x)sin⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂x=100π·e-3π2tcos⁡⁡(π·x)sin⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂y=100π·e-3π2tsin⁡(π·x)cos⁡⁡(π·y)sin⁡(π·z),∂T(x,y,z,t)∂z=100π·e-3π2tsin⁡(π·x)sin⁡(π·y)cos⁡⁡(π·z).In Tables1 and 2, it is evident that, for Δt=0.01 to use 8- or 27-node hexahedrons, there is no significant difference in their numerical results for the T(x,y,z) solution; however, even with the mesh’s spatial refinement, the results are not satisfactory. In turn, with Δt=0001 and Δt=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.Table 1
NormL∞ from the error made by the GFEM T at t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
4.53
E
-
01
1.14
E
-
01
3.42
E
-
01
2.59
E
-
03
3.41
E
-
01
1.60
E
-
03
1/8
3.06
E
-
01
1.13
E
-
01
1.95
E
-
01
1.58
E
-
03
1.94
E
-
01
5.42
E
-
04
1/10
2.37
E
-
01
1.12
E
-
01
1.26
E
-
01
1.30
E
-
03
1.24
E
-
01
2.34
E
-
04
1/16
1.61
E
-
01
1.12
E
-
01
5.01
E
-
02
1.14
E
-
03
1.18
E
-
01
4.60
E
-
05Table 2
NormL2 from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.27
E
-
01
3.57
E
-
02
9.61
E
-
02
8.19
E
-
04
9.58
E
-
02
4.72
E
-
04
1/8
9.08
E
-
02
3.64
E
-
02
5.78
E
-
02
5.14
E
-
04
5.75
E
-
02
1.56
E
-
04
1/10
7.27
E
-
02
3.69
E
-
02
3.86
E
-
02
4.31
E
-
04
3.82
E
-
02
6.76
E
-
05
1/16
5.20
E
-
02
3.79
E
-
02
1.62
E
-
02
3.87
E
-
04
1.58
E
-
02
1.38
E
-
05In other studies regarding the GFEM use, the following technique was used for calculating the derivatives in the three spatial directions forT(x,y,z) [24, 25]:
(26)∂T∂xk≈∂T^e∂xk=∑i=1Nnodes∂Ni∂xkT^ie,
where k=1,2, or, 3, with x1=x,  x2=y, and x3=z, 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]:
(27)∂Ti∂x≅-11Ti+18Ti+1-9Ti+2+2Ti+36Δx,∂Ti∂y≅-11Ti+18Ti+Nnodesx-9Ti+2*Nnodesx+2Ti+3*Nnodesx6Δy,∂Ti∂z≅(-11Ti+18Ti+Nnodesx*NNosy-9Ti+2*Nnodesx*NNosy+2Ti+3*Nnodesx*NNosy)×(6Δy)-1,
where Nnodesx and Nnosy are the number of nodes in the regular mesh in the x and y directions, respectively. The three previous expressions represent the calculation of the T three derivatives, using the third orderforward 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 orderbackward finite differences were used as follows:
(28)∂Ti∂x≅11Ti-18Ti-1+9Ti-2-2Ti-36Δx,
where other derivatives are similar.Table 3
NormL∞ from the error made at Tx≅Ty≅Tz by the GFEM/FEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
4.33
E
-
00
5.22
E
-
01
4.24
E
-
00
3.54
E
-
01
4.24
E
-
00
3.57
E
-
01
1/8
3.25
E
-
00
4.49
E
-
01
3.19
E
-
00
2.01
E
-
01
3.18
E
-
00
2.04
E
-
01
1/10
2.60
E
-
00
4.14
E
-
01
2.55
E
-
00
1.28
E
-
01
2.55
E
-
00
1.31
E
-
01
1/16
1.65
E
-
00
3.76
E
-
01
1.59
E
-
00
4.84
E
-
02
1.59
E
-
00
5.19
E
-
02Table 4
NormL∞ from the error made in Tx≅Ty≅Tz by the GFEM/FDM in t = 0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.12
E
-
00
3.33
E
-
01
7.66
E
-
01
7.75
E
-
02
7.62
E
-
01
7.84
E
-
02
1/8
8.59
E
-
01
3.45
E
-
01
5.08
E
-
01
3.28
E
-
02
5.04
E
-
01
3.35
E
-
02
1/10
7.01
E
-
01
3.49
E
-
01
3.51
E
-
01
1.70
E
-
02
3.47
E
-
01
1.73
E
-
02
1/16
4.99
E
-
01
3.51
E
-
01
1.50
E
-
01
5.24
E
-
03
5.16
E
-
01
4.25
E
-
03Table 5
NormL2 from the error made in Tx≅Ty≅Tz by the GFEM/FEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
1.26
E
-
00
1.20
E
-
01
1.23
E
-
00
1.01
E
-
01
1.23
E
-
00
1.02
E
-
01
1/8
9.86
E
-
01
1.11
E
-
01
9.59
E
-
01
5.74
E
-
02
9.59
E
-
01
5.79
E
-
02
1/10
8.13
E
-
01
1.12
E
-
01
7.88
E
-
01
3.67
E
-
02
7.88
E
-
01
3.72
E
-
02
1/16
5.39
E
-
01
1.18
E
-
01
5.16
E
-
01
1.42
E
-
02
1.47
E
-
01
1.45
E
-
02Table 6
NormL2 from the error made in Tx≅Ty≅Tz by the GFEM/FDM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.0001
8 nodes
27 nodes
8 nodes
27 nodes
8 nodes
27 nodes
1/6
3.21
E
-
01
1.11
E
-
01
2.15
E
-
01
2.31
E
-
02
2.14
E
-
01
2.36
E
-
02
1/8
2.67
E
-
01
1.18
E
-
01
1.55
E
-
01
9.78
E
-
03
1.54
E
-
01
1.00
E
-
02
1/10
2.28
E
-
01
1.20
E
-
01
1.13
E
-
01
5.08
E
-
03
1.12
E
-
01
5.22
E
-
03
1/16
1.70
E
-
01
1.22
E
-
01
5.12
E
-
02
1.67
E
-
03
5.00
E
-
02
1.30
E
-
03Tables3–6 show that this proposal to use the GFEM/FDM (Galerkin Finite Element Method/Finite Difference Method) for the calculation of the T 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:(29)2∂2T∂x2+2∂2T∂y2+2∂2T∂z2-T=∂T∂t.
In a domain of unit cube, its analytical solution will be considered as
(30)T(x,y,z,t)=et(ex+ey+ez),∂T∂x=etex,∂T∂y=etey,∂T∂z=etez.Unlike the application 1, generally, the results for this reaction-diffusion case are significantly better; for example, for aΔt=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 Δt=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 L∞. With the analysis of the largest error made in all the mesh, where it generally occurred at the internal node closest to the vertex (x,y,z)=(1,1,1). It may be noted that for the norm L2 (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 Δt=0.0001, since in Δt=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.Table 7
NormL∞ from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
1/6
7.26
E
-
04
1.99
E
-
06
7.16
E
-
01
1.86
E
-
06
1/8
3.99
E
-
04
1.40
E
-
06
1.66
E
-
01
6.02
E
-
07
1/10
2.52
E
-
04
1.29
E
-
06
9.36
E
-
02
2.49
E
-
07
1/16
9.75
E
-
05
1.25
E
-
06
3.82
E
-
01
3.84
E
-
08Table 8
NormL2 from the error made in T by the GFEM in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
1/6
2.60
E
-
04
5.82
E
-
07
5.49
E
-
02
5.49
E
-
07
1/8
1.51
E
-
04
4.84
E
-
07
8.81
E
-
03
1.84
E
-
07
1/10
9.91
E
-
05
5.04
E
-
07
4.12
E
-
03
7.87
E
-
08
1/16
4.00
E
-
05
5.35
E
-
07
8.10
E
-
03
1.38
E
-
08Table 9
NormL∞ from the error made in Tx≅Ty≅Tz in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
1/6
2.36
E
-
01
2.85
E
-
03
6.53
E
-
03
4.06
E
-
04
4.53
E
-
00
1.24
E
+
01
6.53
E
-
03
4.23
E
-
04
1/8
1.80
E
-
01
1.42
E
-
03
3.73
E
-
03
1.75
E
-
04
1.54
E
-
00
2.75
E
-
00
3.73
E
-
03
1.85
E
-
04
1/10
1.45
E
-
01
9.24
E
-
04
2.41
E
-
03
9.38
E
-
05
7.91
E
-
01
2.69
E
-
00
2.41
E
-
03
9.69
E
-
05
1/16
9.19
E
-
02
3.87
E
-
04
9.61
E
-
04
2.83
E
-
05
6.02
E
-
00
1.76
E
+
01
9.55
E
-
04
2.44
E
-
05Table 10
NormL2 from the error made in Tx≅Ty≅Tz in t=0.1.
h
Δ
t
=
0.01
Δ
t
=
0.001
8 nodes
27 nodes
8 nodes
27 nodes
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
GFEM/FEM
GFEM/FDM
1/6
1.60
E
-
01
1.93
E
-
03
3.60
E
-
03
2.96
E
-
04
4.72
E
-
01
1.15
E
-
00
3.59
E
-
03
2.96
E
-
04
1/8
1.21
E
-
01
8.49
E
-
04
2.02
E
-
03
1.25
E
-
04
1.54
E
-
01
2.49
E
-
01
2.02
E
-
03
1.25
E
-
04
1/10
9.71
E
-
02
4.65
E
-
04
1.29
E
-
03
6.45
E
-
05
1.10
E
-
01
1.38
E
-
01
1.29
E
-
03
6.40
E
-
05
1/16
6.10
E
-
02
1.52
E
-
04
5.07
E
-
04
1.59
E
-
05
1.80
E
-
01
4.40
E
-
01
5.07
E
-
04
1.55
E
-
05


## 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, k=0.72 W/mK, whereas for the same temperature we haveregularbrick: ρ=1920 kg/m3 and cp=835 J/kg·K;cement mortar: ρ=1860 kg/m3 and cp=780 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 (cp) as follows:
(31)ρm=0.02×1,860+0.1×1.9200.12=1,910Kg/m3,cpm=0.02×780+0.1×8350.12≅825.83J/kg·K.By the heat conduction equation described in (2.13), page 43 of [27], it follows that kx=ky=kz=0.72 W/mK and ψ=-ρcp=-1,910×825.83≅-1,577,341.67J/m3·K. Therefore, the governing equation of this application will be this way:
(32)0.72(∂2T∂x2+∂2T∂y2+∂2T∂z2)=1,577,341.67∂T∂t.Thus, we assume a rectangular wall with thickness of 0.12 m (0≤x≤0.12 m), width 1.00 m (0≤z≤1.00 m), and a height of 1.00 m (0≤y≤1.00 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 (x=0.12 m).There are some challenges in this application because both the domain and its initial geometric boundaries present a unique region. Figure1 shows the example.Figure 1
Outline of boundary conditions near the vertex(x,y,z)=(0.12,0,0).Figure1 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≤x≤0.12 m, with y=z=0, toward the x=0 to x=0.12). This situation can lead to extremely high temperature gradients in this region, particularly in the x direction to the lower refining adopted.Figures2 to 6 show some variations of temporal and spatial refinements. For this analysis the line where  x=0.1, y=0.5, and 0≤z≤1.0 was adopted as reference, because of its closeness to the singularities. This analysis begins by varying the values of Δz (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 (x=0.12) toward the inner face (x=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 z=0.5. Thus, in Figure 3, it can be noted that the maximum temperature of this line rotates around T=14.3°C.Figure 2
Mean temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes), (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s, NNost: number of the nodes in all the mesh.Figure 3
Temperature inz=0.5 for x=0.1, y=0.5,  and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s.Figure 4
Maximum temperature forx=0.1, y=0.5,  and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, 1/6≤Δz≤1/120, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, 1/3≤Δz≤1/60, 27 nodes), t=400s.From this reasoning, we analyzed the maximum temperature (Figure4) 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 z=0.5, especially in the singularity near z=1.0. Figure 4 also shows that, with the refinement, the maximum temperature tends the temperature at the point z=0.5, which should actually be.Then the mesh is refined in thex 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 z=0.5 become stagnant for the first refinements; that is, for the proposed situation, the refinement in the x direction brings no better results. Furthermore, the maximum temperature is generally higher than the temperature at z=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 Δx 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 Δx the greater the values of ∂T/∂x in this region.Figure 5
Temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δy=Δz=1/16, 0.12/6≤Δx≤0.12/60, 8 nodes) (mesh: Δy=Δz=1/8, 0.12/3≤Δx≤0.12/30, 27 nodes), Δt=1/20, t=400s.Figure 6
Temperature forx=0.1, y=0.5, and 0≤z≤1.0 (mesh: Δx=0.12/12, Δy=1/8, and Δz=1/16, 8 nodes) (mesh: Δx=0.12/6, Δy=1/4, and Δz=1/8, 27 nodes), 2.5s≤Δt≤20s, t=400s.Lastly, the temporal refinement is a variableΔt, 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 z=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 z=0.5 is around T=14.3°C. Furthermore, it is expected that, with the lower Δt, 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 x. The numerical oscillations presented in Tx occur near the singularities.Figure 7
Profile ofT (left) and Tx (right) for the plane XZ with y=0.5 for a mesh with 27-node hexahedrons, Δx=0.12/6, Δy=1/4, Δz=1/40, and Δt=20s in t=400s.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)
Inx=0 (inner wall), T=12°C (assuming that the internal coolant is able to maintain a constant temperature).(ii)
Inx=0.12 m (outer wall) the expression T(x,y,z,t)=12°C+(t/180)  (°C/s) is taken.(iii)
The other sides will be regarded as geometric symmetry, that is, equivalent to the insulation.Figure8 shows the temperature profile in a XZ 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 (x,z) = (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 x 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 x, 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 planeXZ with y=0.5 in t=900s left and t=2700s right (mesh: Δx=0.12/6, Δy=1/4, Δz=1/4, Δt=30s, 27 nodes).Figure 9
T
x profile (°C/m) in the plane XZ with y=0.5 in t=900s left and t=2700s right (mesh: Δx=0.12/6, Δy=1/4, Δz=1/4, Δt=30s, 27 nodes).Figure 10
Profile ofT (left) and Tx (right) for the plane XZ with y=0.5 for a mesh with Δx=0.12/12, Δy=1/8, Δz=1/8, and Δt=30s in t=2700s.


## 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.



---

*Source: 560492-2014-08-13.xml*