Generalized finite element method with time-independent enrichment functions for 3D transient heat diffusion problems

International Journal of Heat and Mass Transfer 149 (2020) 118969

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Generalized finite element method with time-independent enrichment functions for 3D transient heat diffusion problems Muhammad Iqbal a,⇑, Kashif Masood a, Aiman Aljuhni b, Afaq Ahmad c a

Creative Engineering & Management Solutions, Deans Centre Peshawar, Pakistan Khatib & Alami, Andalus Street Suhaili Business Centre, P.O. Box 9330, Jeddah 21413, Saudi Arabia c Department of Civil Engineering, University of Engineering and Technology Taxila, Pakistan b

a r t i c l e

i n f o

Article history: Received 20 August 2019 Received in revised form 3 October 2019 Accepted 26 October 2019 Available online xxxx Keywords: FEM GFEM Enrichment functions Transient diffusion problems

a b s t r a c t An enriched generalized finite element method (GFEM) is proposed for the efficient solution of threedimensional (3D) transient heat diffusion problems. The proposed GFEM formulation is shown to produce results with better accuracy and less degrees of freedoms (DOFs) as compared to the standard FEM with linear basis functions. For GFEM, the solution space is enriched with an approximate solution describing the heat diffusion decay. Multiple enrichment functions that mimic the solution behaviour are used to capture the high temperature gradients. The enrichment functions are independent of time; the spatial and time varying decay of the solution is embedded in the formulation of the enrichment functions. The resultant formulation significantly reduces the computational cost as compared to the enrichment functions with only spatial approximations. In the current formulation, the system matrix is assembled only at the first time step and retained for subsequent time steps. Different numerical examples in 3D spatial domains are considered to analyze the performance of the proposed method. Depending on the nature of the problem, two different enrichment functions; exponential enrichment functions and hyperbolic enrichment functions are used in the computations. It is shown that the proposed approach is more simple and efficient than the conventional h-refinement to increase the accuracy of the finite element method. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction For decades, the standard finite element method has been successfully used for the accurate and robust solution of 2D steady-state and transient heat diffusion equations. Although the computer architecture and design have become much more advanced in recent times, the solution of most practical 3D problems using classical finite element methods is still a complex and computationally extensive task. For an acceptable accuracy, an unreasonably large number of DOFs are required to solve 3D problems [1], which makes the simulation very difficult and challenging. This difficulty is further augmented with the transient nature of the problems. The solution of excessively large number of equations at every time step makes the solution very time consuming. Some representative application areas where the solution of time-dependent diffusion problems pose difficulties include heat transfer through thermal radiation [2], optical tomography [3] and the cooling process of molten glass [4], among the others. Another difficulty arises ⇑ Corresponding author. E-mail address: [email protected] (M. Iqbal). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118969 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

when diffusion equations need to be solved in combination with problems governed by hyperbolic PDEs [5]. The presence of steep fronts and boundary layers in these coupled parabolic-hyperbolic problems pose severe numerical complications [6,5]. The standard finite element with piecewise continuous polynomial interpolation functions cannot resolve these steep fronts unless a substantially refined mesh grid, or higher order polynomial basis functions are used [7]. In the past two decades, domain based methods with field enrichments have been developed to overcome such difficulties. The main idea of the field enrichments consist to include a priori information about the problem to be solved by introducing specifically designed functions into the approximation space. Such techniques include the Partition of Unity Finite Element Method (PUFEM) [8,9], the Ultra-Weak Variational Formulation (UWVF) [10], the Discontinuous Enrichment Method (DEM) [7,11], the Generalized Finite Element Method (GFEM) [12–14] and the Xtended Finite Element Method (XFEM) [15] among the others. Other methods formulated to improve the accuracy of FEM include nodebased smoothed point interpolation method (NS-PIM) [16,17], edge-based smoothed point interpolation method (ES-PIM)

2

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

[18,19], face-based smoothed finite element method (FS-FEM) [20] and the numerical manifold method (NMM) [21]. The PUFEM and Partition of Unity (PU) based GFEM and XFEM, have been successfully used for the solution of complex practical problems. The core idea of PU was originally proposed by Babuška et al. [22]. In classical FEM the approximation is made with polynomial trial functions, but in PUFEM the formulation of the trial functions incorporates known information about the solution of the problem under consideration. The Partition of Unity method (PUM) developed by Melenk and Babuška [8] allows using any PU with the local enrichment functions. A particular instant of the PUM is GFEM, where the conventional ‘‘hat-functions” of FEM serves as the PU [23]. Strouboulis et al. [12,24] introduced GFEM as a combination of standard finite elements and the enriched partition of unity finite element method. They multiplied the special functions of the partition of unity method with the standard vertex shape-functions which were then used in combination with the existing finite element basis to make an augmented conforming finite element space. The literature reveals that to date enriched FEM is mainly explored for fracture mechanics [25] and wave problems with high frequencies [26]. Munts et al. [27] are among the early researchers who used the enriched FEM for heat transfer problems. They used PUFEM to solve 1D convection–diffusion problems in using two different enrichment functions; exponential enrichment function based on approximate analytical solution and polynomial enrichment function as in [28]. They show that the choice of the enrichment function has effect on the results of PUFEM. O’Hara et al. [29] investigated the application of the GFEM with global-local enrichments (GFEMgl) for steady-state heat transfer problems with solutions that exhibit highly localized sharp thermal gradients. They defined a local problem on a fine-mesh to formulate enrichments which were then used to enrich the coarse mesh of the global problem. This methodology eliminated the need of a priori knowledge about the problem in hand. The authors extended their work to time dependent problems exhibiting localized thermal gradients [30,31]. Soghrati et al. used the interface-enriched FEM to capture the variation of temperature and it gradients in 2D [32] and 3D problems [33]. Van der Meer et al. [34] analysed transient geothermal problems using the enriched FEM. Their algorithm considered the use of enrichment functions that were updated at every time step and the solution at each time step was calculated using these updated enrichment functions. Cosimo et al. [35] also used timedependent enrichment functions to model the solidification process. Similar enrichment functions were previously proposed in [36] for two-phase flow problems. In a recent contribution, Zeller et al. [37] used XFEM with parametrized enrichment functions for one-dimensional transient problems. A notable feature of the above mentioned developments is that the basic formulation of the finite elements is time-dependent for the implementation of these approximation functions. Computing the solution for huge number of time steps where the system matrix has to be evaluated at every time step, makes the analysis process cumbersome and time consuming. Given that the problem is in 3D domain and time dependent, computing the solution progress over highly refined mesh grids and for large number of time steps can lead the procedure to become computationally expensive and impractical, even with advanced computing facilities. This current work is motivated by providing an efficient solution for 3D transient heat diffusion problems by the introduction of time-independent global enrichment functions in the framework of GFEM as in [38,39]. The time variation of temperature field is embedded in the definition of the enrichment functions. This eliminates the need for updating enrichment functions at each time step, as a result the system matrix is assembled only at the first time step and retained for the subsequent time steps. A

remarkable reduction in the computational cost is reported in the results. The rest of the paper is arranged as follows. Section 2 states the formulation of the considered boundary value problem with prescribed initial and boundary conditions along with its weak formulation. The numerical approximation of the problem by standard FEM and GFEM is detailed in Section 3. The main computations showing the efficiency of the GFEM are presented in Section 4. Section 5 contains some concluding remarks. 2. Boundary value problem and weak formulation Let X
R3 be a three-dimensional open, bounded domain with polygonal boundary C and time interval 0; T, we are interested to solve the following transient heat diffusion problem

@u  kr2 u ¼ f ðt; xÞ; @t

ðt; xÞ 20; T  X;

ð1Þ

here x ¼ ðx; y; zÞT are the spatial coordinates, t is the time variable, k > 0 is the diffusion coefficient and f ðt; xÞ describes the effect of internal sources/sinks. An initial condition is considered as

ðxÞ 2 X;

uðt ¼ 0; xÞ ¼ u0 ðxÞ;

ð2Þ

where u0 ðxÞ is a prescribed initial field. The aforementioned equations are subjected to the Robin boundary condition

@u þ au ¼ gðt; xÞ; @n

ðt; xÞ 20; T  C;

ð3Þ

where n is the outward unit normal on the outer boundary C and gðt; xÞ is a given boundary function, and a P 0 is the convection heat coefficient on C. To solve Eq. (1) numerically, we choose an implicit Euler scheme for time discretization. The time interval ½0; T is divided into quasi-distributed subintervals ½t n ; tnþ1  of size dt ¼ t nþ1  tn for n ¼ 0; 1; . . . N t . The time derivative in (1) is approximated by a difference quotient.

unþ1  un  kr2 unþ1 ¼ f ðt nþ1 ; xÞ: dt

ð4Þ

This can be rearranged as

 r2 unþ1 þ

1 nþ1 1 u ¼ ðdtf ðt nþ1 ; xÞ þ un Þ; kdt kdt

ð5Þ

or

r2 unþ1 þ ku

nþ1

¼ F nþ1 :

ð6Þ

with the definition of F nþ1 and k as

F nþ1 ¼ kðdtf ðtnþ1 ; xÞ þ un Þ;

k¼

1 : kdt

To approximate Eq. (6) numerically, we multiply the equation with a test function W, and then integrate it over X

Z 

X

W r2 unþ1 dX þ

Z

Wku X

nþ1

Z

dX ¼

X

WF nþ1 dX 8W

2 H1 ðXÞ:

ð7Þ

where H ðXÞ is the conventional Sobolev space. Now suppose we have a function S, such that 1

S ¼ W runþ1

8 9 nþ1 W @u@x > > > > < = nþ1 ¼ W @u@y > > > : @unþ1 > ; W @z

We find the divergence of the function S,

ð8Þ

3

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

r  S ¼ rW  runþ1 þ W r2 unþ1

ð9Þ

unþ1 ðxÞ ¼

The divergence theorem states that

Z

X

r  SdX ¼

Z

C

S  ndC

ð10Þ

X

rW  runþ1 þ W r2 unþ1 dX ¼

C

X

W r2 unþ1 dX ¼

Z C

W runþ1  ndC

W runþ1  ndC 

Z X

ð11Þ

rW  runþ1 dX

finally, substituting (12) into (7) results

Z



ZC

¼

X

W runþ1  ndC þ

Z

X

rW  runþ1 dX þ

Z

ð12Þ

Z Wku

nþ1

X

dX ð13Þ

ðrW  runþ1 þ Wku X Z ¼ WF nþ1 dX :

nþ1

Z ÞdX 

C

Z nþ1 ðrW  runþ1 þ Wku ÞdX þ W aunþ1 dC X C Z Z ¼ WF nþ1 dX þ Wg nþ1 dC : C

ð14Þ

ð15Þ

We aim to find an approximate numerical solution unþ1 of the weak form (15), using the standard FEM and GFEM. 3. FEM and GFEM formulations To solve the weak form (15) with finite elements, the domain X is discretized into finite number of elements. We assume that X is a polygon and generate a quasi-uniform partition X
Xh of N e nonoverlapping element domains Xe .

Xh ¼

Ne [

Xe with

e¼1

Ne \

Xe ¼ £

e¼1

The element edges are denoted by E. For standard FEM, the solution in each element is approximated by sum of the combination of nodal values and piecewise linear shape functions, i.e. we look for unþ1 of the form

unþ1 ðxÞ ¼

M X Nj ðxÞunþ1 j

ð16Þ

j¼1

where unþ1 are the unknown nodal values, N j are the piecewise linj ear shape functions and M refers to the total number of nodes in the element. In GFEM the solution space is enriched with basis functions that incorporates a priori knowledge about the problem in hand. This strategy provides better approximation as compared to conventional piecewise linear basis functions. For GFEM solution, (16) is modified and nodal values unþ1 are taken as combination j of the enrichment functions, i.e.

but in the FEM is used to calculate the nodal solution values unþ1 j Taking W ¼ Pr results in the GFEM discretization of the weak form (15) such that: Find unþ1 of the form (17) such that u0 ¼ u0 and for all r ¼ 1; . . . ; MQ

Z

nþ1

X

and with the boundary condition (3), the weak formulation of the heat diffusion problem is concluded as: Find a solution unþ1 on the domain X such that u0 ¼ u0 and for all test functions W on X and all n 2 N such that

X

where Gq ðxÞ are the enrichment functions and Aqj are their amplitudes at nodal points of the domain. Q denotes the total number of enrichment functions. In (17) the finite element shape functions Nj ðxÞ; j ¼ 1; . . . M, constitutes the partition of unity, i.e., PM j¼1 N j ðxÞ ¼ 1 for all x in X covered by the FE mesh. To solve (16),

ðrPr  runþ1 þ Pr ku ÞdX þ X Z Z ¼ Pr F nþ1 dX þ Pr g nþ1 dC:

W runþ1  ndC

X

Z

ð18Þ

(17) the system is solved for new unknowns Aqj instead of unþ1 . j

WF nþ1 dX;

or

Q X Aqj Gq ðxÞ q¼1

then rearranging gives

Z

with

unþ1 ¼ j Z

ð17Þ

j¼1 q¼1

substituting (8) and (9) into (10)

Z

Q M X X Aqj N j ðxÞGq ðxÞ

C

Z

C

Pr aunþ1 dC ð19Þ

It is worth mentioning that the functions Pr are defined in the global co-ordinates but modulated locally by multiplying with local shape functions Nj . The selection of enrichment functions is dictated by the physical behaviour of the problem at hand. For the example problems considered in this article, two different enrichment functions are used depending on the nature of the problem. Formulation of both enrichment functions is detailed in the respective example problems where they are used. 4. Numerical experiments This section investigates the performance of the proposed GFEM for three dimensional transient heat diffusion problem defined by (1)–(3). We use the classical FEM and GFEM approaches to solve the weak formulation (15) and (19), respectively. To assess the performance of the proposed GFEM algorithm against the standard FEM approach, three different example problems are considered. To compute the solution numerically, we used an 8 noded Hexahedral mesh with piecewise linear shape functions. To evaluate the numerical integrals, the Gauss Legendre quadrature is used. It is ensured that no integration errors are introduced due to the number of quadrature points. On the basis of a preliminary study, 20 integration points in each direction for GFEM whereas for FEM with piecewise linear basis functions only 2 integration points in each direction are used to evaluate the integrals. A direct solver is used to solve the resulting system of equations. All the computations are performed on IntelÒ XeonÒ ES-1620 CPU @ 3.50 GHz processor speed with 32 GB installed RAM using self-written FORTRAN codes with standard libraries. The codes are not parallel and only consider the default optimization of the computer. All the quantities in the numerical examples are measured in SI units. For the reasons of simplification, all the quantities are mentioned without units in the discussion. 4.1. Example Problem 1: Approximation of problem with a known analytical solution To asses the accuracy of the method, we take a problem with a known analytical solution as our first example problem. We consider a diffusion problem in 3D domain defined by X ¼ ½0; 23 . For

4

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

the proposed problem, the source term f ðt; xÞ, the boundary term gðt; xÞ and the initial condition u0 ðxÞ are selected in such a way that the exact solution is defined by c c

c

c c

c

kt

Uðx; tÞ ¼ x ð2  xÞ y ð2  yÞ z ð2  zÞ ð1  e

Þ

ð20Þ

where t is the time variable, x ¼ ðx; y; zÞT are the spatial coordinates and k > 0 is the heat diffusion coefficient.The parameter c controls the steepness of the thermal gradient such that higher value of c results in sharper thermal gradient and vice versa. To produce sharp thermal gradients, we use c ¼ 20 in our computations. In order to quantify the error in FEM and GFEM solutions, we calculate the relative L2 error (e2 ) defined as

e2 ¼

jju  UjjL2 ðXÞ jjUjjL2 ðXÞ

 100

ð21Þ

where u and U are the approximate numerical and exact solutions, respectively. To calculate the GFEM solution, the approximation space is enriched by suitable enrichment functions. The literature reveals instances where enrichments are applied locally [29–31], but for our proposed GFEM formulation, we use the following global functions of varying standard deviations.

Gq ðxÞ ¼

e

R q

Rc q

 e ð C Þ ; Rc q 1  eð C Þ 0 C

q ¼ 1; 2; . . . ; Q :

ð22Þ

Here R0 :¼ jx  xc j is the distance from the centre of enrichment function xc to the point x. Rc and C are constants which control the shape of Gq , and Q is the total number of enrichment functions. Similar enrichment functions are used for the PUFEM solution of 2D problems in [38]. With the definition of enrichment functions, the GFEM approximation defined by (17) can now be written as nþ1

u

Q M X X e ðxÞ ¼ Aqj Nj ðxÞ j¼1 q¼1

R q

Rc q

 eð C Þ Rc q 1  eð C Þ 0 C

ð23Þ

For simplicity the product of Gq ðxÞ and the polynomial shape function N j ðxÞ is taken as a new shape function Pn , i.e.,

Pn ¼ Nj ðxÞ

e

R q

Rc q

 eð C Þ Rc q 1  eð C Þ 0 C

ð24Þ

The derivatives of Pn are given as

 8 @Nj 9 > > > < @x > = @Nj @P n ¼  R q @y @y > ð c Þ > > > > > 1e C : @Pn > ; : @Nj > ; @z 8 @z 9 > R q < ðx  xc Þ > = 0 ðq2Þ q  C 1 e R0 Nj ðy  yc Þ Rc q Cq > > : ; 1eð C Þ ðz  zc Þ 8 9 @P n > > > < @x > =

  q R e



0 C



e

q

ðRCc Þ

ð25Þ

For numerical computations, we use the parameter a ¼ 1 and the time step value Dt ¼ 0:001 with a total solution time of t ¼ 0:1. The heat diffusion coefficient k is selected to be 0:01 for a first set of results. Another set of results is presented with k ¼ 0:1. In order to get a converged solution, we use uniform hrefinements for FEM, while for GFEM both h- and q-refinements are considered. To solve the problem with FEM, we start with a coarse mesh of 1000 elements with a total of 1331 degrees of freedom (DOFs), which results in an L2 error of 8.04% at t ¼ 0:1. The mesh is then refined gradually to improve the quality of results. Fig. 1 shows the meshes used for the FEM computations. For GFEM, the computations are started with a very coarse mesh of 64 elements in the 3D cube and the solution space is enriched with different number of enrichment functions, Q = 2,3. . .6. The mesh is then refined to a total of 216 elements in the whole domain with 6 elements in each direction. Again the same enrichment functions Q = 2,3. . .6, are considered. A third refinement of 512 elements with Q = 2,3. . .5 is also considered to improve the results further. In the coming discussion an analysis with 64 elements will be referred to as GFEM1, and with 216

Fig. 1. 3D meshes used for FEM computations.

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

and 512 elements as GFEM2 and GFEM3, respectively. Fig. 2 shows the meshes used for GFEM computations. Fig. 3 compares the L2 errors for both FEM and GFEM solutions at different time intervals with h-refinements for FEM and qrefinements for GFEM. In the figure, the DOFs are plotted on the abscissa while on the ordinate the L2 error is plotted on a logarithmic scale. Comparison of both the computations show that the GFEM solution converges much faster than FEM and gives much better L2 error with less DOFs. As shown in Fig. 3(a), in case of FEM a maximum accuracy of 0.05% is achieved with 226981 DOFs for k ¼ 0:01. With GFEM we get a minimum L2 error of 0.17% with only 3645 DOFs. In case of FEM a comparable accuracy of 0.15% is achieved with 68921 DOFs. For FEM, this accuracy is achieved with a very fine mesh of 64000 elements, while in case of GFEM we use a very coarse mesh of only 512 elements with Q ¼ 5. For k ¼ 0:1, the comparison of the GFEM and FEM results confirm the above observation. With FEM, the maximum achievable accuracy is 0.22% with a very fine mesh of 216000 elements and 226981 total DOFs (TOTDOF). GFEM gives a comparable accuracy of 0.19% with only 3645 TOTDOF. The TOTDOF in the case of GFEM is less than 2% of that used for FEM, resulting in a saving of more than 98% in the total DOFs. To process 216000 elements and build the linear system at the first time step, FEM solution took 228589s (63 h, 29 min and 49 s), and another 5818s (1 h, 36 min and 58 s) to solve the resulting system of equations. The corresponding times for GFEM to process 512 elements with 5 enrichment functions are 883.5s (14 min and 43.5 s) and 1.21s to solve the system of equations. Also, to update the right hand side of the system of equations and to recalculate the solution at subsequent time steps, FEM took on average 340s at every step. For GFEM the corresponding time is only 42s. As iterated previously, the enrichment functions are time-independent, therefore the system matrix is assembled only at the first time step and retained for the subsequent time steps. As a result the time saving at the subsequent time steps is also obvious in the case of GFEM. This is emphasized that despite the fact that GFEM uses high number of integration points, yet the very low number of elements with GFEM makes the time needed to assemble the system matrix much smaller than that for the FEM. The CPU time for the solution of GFEM system is only a small fraction of the FEM. An important observation made from Fig. 3 reveals that in case of FEM, we get higher errors for k ¼ 0:1 as compared to k ¼ 0:01 for the same computation time. This is because the rate of diffusion is increased by increasing the value of k, and the heat now travels to more distant areas of the domain for the same computation time. In case of k ¼ 0:01, minimum L2 error of 0.05% is achieved with 226981 TOTDOF, which is case of k ¼ 0:1 is 0.22%. In comparison the diffusion of heat to the larger area of the domain does not affect the L2 error considerably for the GFEM solution. The global nature

5

of enrichment functions ensure to efficiently capture the solution in the whole domain and a very small change in the solution accuracy is observed by changing the value of k from 0.01 to 0:1. Another comparison using h-refinements for both FEM and GFEM results is given in Fig. 4. GFEM2Enr denotes the GFEM results using 2 enrichment functions and with three different mesh densities of 64, 216 and 512 elements. Similarly GFEM3Enr, GFEM4Enr, GFEM5Enr and GFEM6Enr are the results of the GFEM computations with 3, 4, 5 and 6 enrichment functions, respectively. Comparison of the GFEM results show that the use of higher number of enrichment functions lead to better accuracy with less DOFs. For example, GFEM4Enr gives an L2 error of 1.39% with 1372 TOTDOF. With GFEM6Enr a better accuracy of 1.06% is achieved with 750 TOTDOF, resulting in a saving of 622 DOFs. Similarly with GFEM6Enr and with 2058 TOTDOF, the L2 error is 0.43% and in case of GFEM3Enr the same error is achieved with 2916 TOTDOF. This suggests to use higher number of enrichment functions in the GFEM rather than refined meshes, but one problem that limits the number of enrichment functions is the conditioning of the system matrix. The conditioning of the system matrix deteriorates with higher numbers of enrichment functions. In their work, Malek et al. [40] also reported this problem. The same was observed in [41] as well. For visualization of the results, Fig. 5 displays the temperature profiles at different simulation times for the exact, FEM and GFEM solutions using k ¼ 0:1. For GFEM simulations, the presented plots are obtained with 216 elements using 4 enrichment functions. For FEM the results are presented using 216000 elements. Both FEM and GFEM results show temperature trends similar to the exact solution but in the case of GFEM the total DOFs used to get the solution is only a small fraction of that used to obtain the FEM solution.

4.2. Example Problem 2: Approximation of problem with a cubic domain The second example considers a transient heat diffusion problem in a 3D domain X ¼ ½0; 13 with a heat source in the central part. A cross section through the middle of the domain is shown in Fig. 6. For the central part, x 2 ½0:4; 0:63 , the source dissipates heat at a constant rate of f ¼ 2000 and is zero in the rest of the domain. The total simulation time is taken to be t ¼ 0:4 with a time step value of Dt ¼ 0:001. The source dissipates heat for half of the simulation time, i.e. from t ¼ 0 to t ¼ 0:2, and is then switched off and the medium is allowed to cool down for the remaining half of the simulation. The initial temperature uo of the domain and the boundary source g, of expressions (2) and (3), are set to be 300. The convection heat transfer coefficient on the boundaries is taken

Fig. 2. 3D meshes used for the GFEM computations.

6

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

Fig. 3. L2 error (e2 ) for FEM and GFEM solutions with uniform h-refinements for FEM and q-refinements for GFEM using k ¼ 0:01 (top row) and k ¼ 0:1 (bottom row).

Fig. 4. L2 error (e2 ) for FEM and GFEM solutions with h-refinements using k ¼ 0:01 (top row) and k ¼ 0:1 (bottom row).

to be a ¼ 1, and the diffusion heat transfer coefficient for the medium is assumed to be k ¼ 0:01. Another more challenging case with k ¼ 0:001 is also considered to show the usefulness of the proposed GFEM approach. Example Problem 2 compares the results for a problem with no known analytical solution, using GFEM and the standard FEM. For GFEM we use a coarse mesh of 125 elements with 5 enrichment functions, whereas the standard FEM solution is obtained with 1000 elements. The mesh densities are chosen such that both GFEM and FEM have comparable DOFs. For GFEM the TOTDOF is

1080, whereas for FEM the TOTDOF is 1331. In order to calculate the errors of the GFEM and FEM solutions, a reference FEM solution (FEMR) on a very fine mesh of 216000 elements is used. Having FEMR as a reference solution, the accuracy of GFEM is compared against the FEM solution. Fig. 7 shows the temperature distribution for the reference solution FEMR as well as for FEM and GFEM. Results are presented at three different simulation times; t ¼ 0:1; 0:2 and 0:4, that is half of the time the source is on, when the source is switched off and another 0:2 after the source is switched off. The results captured

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

7

Fig. 5. Temperature distribution in the middle of the domain for exact (top row), FEM (centre row) and GFEM solutions (bottom row) at t ¼ 0:05 (left column) and t ¼ 0:1 (right column). The simulations are obtained using k ¼ 0:1.

Fig. 6. x-section of the 3D domain for Example Problem 2 with a heat source in the centre.

with GFEM show the same temperature profiles as obtained by the fine mesh reference solution FEMR, whereas the FEM solution shows very coarse profiles. This ascertains that GFEM with a very coarse mesh was able to capture the same solution dynamics as

those obtained on a very fine mesh solution FEMR. GFEM uses only 1081 DOFs to capture the sharp thermal gradients whereas the FEMR solution is obtained with 226981 TOTDOF. This significant reduction in the TOTDOF clearly shows the advantage of GFEM on the standard FEM. The FEM solution was not able to capture the solution accurately despite using higher number of DOFs than GFEM. The figures show that the temperature of the medium increases for half of the simulation time when the source is dissipating heat. The medium starts to cool down when the heat source is switched off. The temperature is maximum in the centre of the domain because of the heat source and it moves towards the boundaries where temperature is lower. To quantify the accuracy of the GFEM method, we calculate the relative L2 error (e2 ) of GFEM and FEM solutions against the reference solution FEMR. Fig. 8 shows the solution along a cross section located in the middle of the domain for both values of k. Results shown in Fig. 8(a) are obtained using k ¼ 0:01 and are presented at the selected simulation times. For FEM solution, the maximum error is observed at the centre of the domain. For half of the simulation time when the source is dissipating heat, the maximum L2 norm error in the centre of the domain is more than 20%. For GFEM, the error in the centre of the domain is less than 1%. In the case of GFEM, the maximum error of around 7% is observed only along the boundaries of the heat source. This error can be minimized by tailoring the enrichment functions according to the shape of the heat source. In the case of FEM, the error around the boundaries of the heat source is more than 15%. For later simulation times when the heat source is turned off and the medium is cooling down, the error decreases significantly for both FEM and

8

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

Fig. 7. Temperature distribution in the middle of domain for FEMR (left column), FEM (centre column) and GFEM solution (right column) at t ¼ 0:1 (top row), t ¼ 0:2 (middle row) and t ¼ 0:4 (bottom row). The simulations are obtained using k ¼ 0:01.

Fig. 8. L2 error (e2 ) for FEM and GFEM solutions calculated along a centre line at different simulation times.

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

GFEM. A similar set of results is presented in Fig. 8(b) using k ¼ 0:001. The local errors for this case are relatively higher both for FEM and GFEM solutions. In the previous case, the maximum error in the middle of the domain for FEM solution at t ¼ 0:2 was around 20% which rises to a value as high as 50% for k ¼ 0:001. For GFEM, we see an increase from 1% to 10%. Higher errors are observed in the central part even at later simulation times but far from the heat source the error still remains very low, even lower than the previous case. For the previous case, at t ¼ 0:4 the error dropped to 5% for the FEM solution which in the present case is still around 40%. For GFEM this value was around 1% in the previous case but in the present case the maximum error is around 10%. To further analyse the results of Example Problem 2, in Fig. 9 we plot the time evolution of the maximum temperature and the difference in the maximum temperatures. From start to the end of simulation, the variation of maximum temperature is plotted in Fig. 9(a) for FEMR, FEM and GFEM. Fig. 9(b) shows the difference of maximum temperatures between the reference FEMR and FEM solutions and the FEMR and GFEM solutions. Fig. 9(a) shows clearly that the time evolution of maximum temperature for the coarse mesh GFEM solution is almost identical to the very fine mesh FEMR solution. Both FEMR and GFEM solutions have very close values for the whole duration of the simulation. On the other hand, the FEM results are far away from the reference FEMR solution for most of the simulation time. Only at later times, when the medium is cooling down and the solution tends towards the steady state condi-

9

tion, the FEM results become closer to the FEMR solution. Fig. 9 (b) also shows that there is a significant difference in temperatures between the FEM and FEMR solutions while for GFEM the difference in maximum temperatures is very small for the whole simulation time. In the case of FEM, the maximum difference is about 130, while in the case of GFEM it stays well below 10. For FEM, the accuracy deteriorates further when using k ¼ 0:001 as shown in Fig. 10. As observed in Fig. 9(a), FEM results tend to become closer to the FEMR solution in the cooling stage and almost identical at the end of the simulation but in Fig. 10(a), FEM results are relatively far from the FEMR results even during the cooling stage. The same trend is observed in Fig. 10(b) where the maximum difference in temperature at t ¼ 0:2 reaches a high value of 400 and stays around 300 even at the end of the simulation. In the case of GFEM, the temperatures are nearly identical to the fine mesh FEMR results, ascertaining the advantage of GFEM. 4.3. Example Problem 3: Approximation of problem with an annular domain The third example simulates the transient heat transfer in a 3D annular domain with inner diameter di ¼ 1:0 and outer diameter do ¼ 1:5, as shown in Fig. 11. The depth of the domain in the zdirection is taken to be 0:3. This example presents one of the representative problems where the domain has curved edges, and the errors from the discretization of the geometry will have a contribution to the overall error of the solution. Mohamed [42] suggested to

Fig. 9. Time evolution of (a) maximum temperatures (b) differences in max temperatures. Results are obtained using k ¼ 0:01.

Fig. 10. Time evolution of (a) maximum temperatures (b) differences in max temperatures. Results are obtained using k ¼ 0:001.

10

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

Fig. 11. Annular domain for Example Problem 3.

use elements that describe the exact shape of the curved domain. It is shown that this approach reduces the overall error of the solution. In this current study, the same approach is used and elements are defined such that the geometry is described exactly. The element having curved boundaries are defined by the polar coordinates. It should be noted that if we don’t discretize the geometry exactly, then a relatively larger number of elements would be required to represent the geometry accurately. For more complex geometries, for example industrial problems, a relatively fine mesh will of-course offset some of the advantages of GFEM, but in the case of FEM even finer mesh grids would be required to accurately represent the geometry. To capture the sharp gradients with coarse meshes in complex shape geometries, the Iso-geometrics GFEM approach would be an alternate option to exactly represent the geometry [40]. For this example problem, a boundary source at the inner surface of the body is selected which dissipates heat with g ¼ 2000. The initial temperature of the domain is set at uo ¼ 300 and the heat source within the domain is set to zero. The time step value is fixed at Dt ¼ 0:0001, and the diffusion heat coefficient for the medium is taken as k ¼ 0:01. The convection heat transfer coefficient at the outer surface is set to be a ¼ 1. The same value of a is used for the top and bottom surfaces of the domain. To capture the high thermal gradients at the boundaries of the domain, hyperbolic tangent basis functions are used to enrich the

Fig. 12. Variation of enrichment function Gtanh along the radial direction with q different values of Hq .

solution space. These enrichment functions are selected due to the type of temperature behaviour of the considered problem. Mohamed et al. [39] used hyperbolic enrichment functions to capture high thermal gradients in domains which were localized in an otherwise uniform domain. For a circular edge with r ¼ re , the enrichment functions in polar coordinates (r; h) are defined as

Gtanh ¼ C 1 þ C 2 tanh q

  r  re ; Hq

q ¼ 1; 2; . . . ; Q :

ð26Þ

where C 1 and C 2 are the control parameters defining the amplitude . The global derivatives of the enrichment funcof the function Gtanh q tion are given as

   @Gtanh x 2 r  re q 1  tanh ¼ C2 rHq @x Hq

Fig. 13. Meshes used for the computations of Example Problem 3.

ð27Þ

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

11

Fig. 14. Temperature variation along the radial direction of the domain at different time intervals for FEMr, FEMf, FEMc and GFEM solutions.

Fig. 15. Time evolution of (a) maximum temperatures (b) difference in maximum temperatures.

   @Gtanh y 2 r  re q 1  tanh ¼ C2 rHq @y Hq

ð28Þ

   @Gtanh z 2 r  re q 1  tanh ¼ C2 rHq @z Hq

ð29Þ

For the considered problem, r e ¼ ri is the inner radius of the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi domain, while r ¼ x2 þ y2 is the radius at a particular point in the domain. For computations, the values of C 1 and C 2 are taken as C 1 ¼ 1:0 and C 2 ¼ 1:0. The enrichments are taken to vary only in the x and y directions and are constant in the z-direction. This is because in the selected problem the temperature varies only in the radial direction and is constant along the thickness of the domain. This results in

@Gtanh q @z

¼ 0. The steepness of the enrichment function

depends on the value of the parameter Hq . Fig. 12 illustrates Gtanh q the variation of the enrichment function along the radial direction of the domain for different values of Hq . It is evident that the use of smaller values of Hq leads to functions with high gradients while higher values of Hq produce slowly varying enrichment functions. The functions with high gradients can recover solutions at early

time steps in the vicinity of a thermal shock while enrichment functions that vary slowly can be useful at later time steps when the temperatures become more uniform. To compare the accuracy of GFEM with the conventional FEM, two different FEM meshes and one GFEM mesh are considered for the computations. A coarse FEM mesh termed as FEMc is selected such that it has comparable DOFs to the enriched GFEM model and a fine mesh FEMf which has about ten times more DOFs. Since the exact analytical solution is not known for the considered problem, a reference FEM solution (FEMr) is obtained using a very fine mesh of 30000 elements with 33150 DOFs. Fig. 13 shows the meshes used for the computations. For GFEM computations a total of 4 enrichment functions are used resulting in total 32 DOFs. In case of FEMc the total DOFs are 36 whereas FEMf solution is obtained with total 360 DOFs. To compare the results obtained with FEM and GFEM solutions, Fig. 14 shows the temperature variation along the radial direction of the domain. The temperatures obtained with FEM and GFEM are compared with the reference FEMr solution. For all the shown simulation times, GFEM produced very similar results to those produced by FEMr. Both FEMr and GFEM results are almost identical. FEMc, which has similar DOFs

12

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969

to GFEM, produces very large errors. The FEMf although having ten times more DOFs, still leads to considerable differences in the results at the earlier time steps due to the higher temperature gradients. The diffusion nature of the problem makes the solution smoother at later time steps and the difference between FEMr and FEMf solutions start to decrease. Also from the graphs, it is evident that the differences in the temperature are more prominent at the inner side of the domain because of the localized large temperature gradient. The temperate field is more uniform in the outer region of the domain and smaller errors were expected by all methods. The results are further assessed in Fig. 15 by plotting the maximum domain temperatures and the difference in the maximum temperatures between the reference FEMr solution and FEMc, FEMf and GFEM solutions. Fig. 15(a) shows the variation of the maximum domain temperature for the whole simulation time while in Fig. 15(b) is shown the absolute differences in maximum temperatures between the FEMr and FEMc, FEMr and FEMf, and FEMr and GFEM solutions. From Fig. 15(a) it is clear that the temperatures produced by the coarse mesh GFEM solution are almost identical to the fine mesh FEMr solution. The FEMc solution with similar DOFs to the GFEM, produces results that differ significantly from the FEMr results. Even the FEMf results have lower accuracy although they are produced using a mesh grid with ten time more DOFs than the GFEM. This clearly indicates the superiority of the GFEM approach over the standard FEM approach. Also Fig. 15(b) shows a close similarity in the results of the GFEM and FEMr solutions. The absolute differences in the maximum temperatures are vey minimal between the GFEM and FEMr results. The FEMf results on the other hand show noticeable variation as compared to FEMr, with the FEMc results further deviated from the FEMr solution. 5. Conclusions The coarse meshed GFEM approach is used to solve timedependent heat diffusion problems in 3D domains and is assessed against the FEM solution on refined mesh grids. For GFEM, the solution space is enriched with 3D Gaussian or hyperbolic tangent functions describing the heat diffusion decay with various rates. To show the advantage of the GFEM approach over the standard FEM, three different example problems are considered. To quantify the errors in both approaches, a problem with a known analytical solution is considered as a first example. It is concluded that for a comparable accuracy, GFEM requires less DOFs than the conventional FEM. Numerical results show that with GFEM, the total DOFs are less than 5% of that used for the standard FEM. The coarse meshes used with GFEM made the computational time only a small fraction of that for the FEM solution. Also the time-independent nature of the enrichment functions resulted in remarkable reduction in the total computation time for large number of steps. To investigate the performance of GFEM further, two more example problems are considered for which the analytical solutions are not known. Example Problem 2 considered a cubic domain while Example Problem 3 assessed the method on an annular computational domain. Exponential Gaussian functions are used for Example Problem 1 and 2 while hyperbolic tangent functions are used for the third example due to the different nature of the problem. Both enrichment approaches are shown to be successful in significantly reducing the computational efforts and lead to better quality results with fewer DOFs as compared to the standard FEM. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1] P. O’Hara, Finite element analysis of three-dimensional heat transfer for problems involving sharp thermal gradients Master’s thesis, University of Illinois at Urbana-Champaign, USA, 2007. [2] M. Frank, A. Klar, E.W. Larsen, S. Yasuda, Time-dependent simplified PN approximation to the equations of radiative transfer, J. Comput. Phys. 226 (2007) 2289–2305. [3] A.D. Klose, E.W. Larsen, Light transport in biological tissue based on the simplified spherical harmonics equations, J. Comput. Phys. 220 (2006) 441–470. [4] G. Thömmes, R. Pinnau, M. Seaïd, T. Götz, A. Klar, Numerical methods and optimal control for glass cooling processes, Transp. Theory Stat. Phys. 31 (2002) 513–529. [5] M. Seaïd, An Eulerian-Lagrangian method for coupled parabolic-hyperbolic equations, Appl. Numer. Math. 59 (2009) 754–768. [6] M. Seaïd, Multigrid Newton-Krylov method for radiation in diffusive semitransparent media, J. Comput. Appl. Math. 203 (2007) 498–515. [7] R. Tezaur, L. Zhang, C. Farhat, A discontinuous enrichment method for capturing evanescent waves in multiscale fluid and fluid/solid problems, Comput. Methods Appl. Mech. Eng. 197 (2008) 1680–1698. [8] J.M. Melenk, I. Babuška, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Eng. 139 (1996) 289– 314. [9] J.M. Melenk, I. Babuška, The partition of unity method, Int. J. Numer. Meth. Eng. 40 (1997) 727–758. [10] T. Huttunen, P. Monk, J.P. Kaipio, Computational aspects of the ultra-weak variational formulation, J. Comput. Phys. 182 (2002) 27–46. [11] C. Farhat, I. Harari, L.P. Franca, The discontinuous enrichment method, Comput. Methods Appl. Mech. Eng. 190 (2001) 6455–6479. [12] T. Strouboulis, I. Babuška, K. Copps, The design and analysis of the generalized finite element method, Comput. Methods Appl. Mech. Eng. 181 (2000) 43–69. [13] T. Strouboulis, I. Babuška, R. Hidajat, The generalized finite element method for Helmholtz equation: theory, computation, and open problems, Comput. Methods Appl. Mech. Eng. 195 (2006) 4711–4731. [14] T. Strouboulis, K. Copps, I. Babuska, The generalized finite element method: an example of its implementation and illustration of its performance, Int. J. Numer. Meth. Eng. 47 (2000) 1401–1417. [15] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Meth. Eng. 45 (1999) 601–620. [16] S. Wu, G. Liu, H. Zhang, G. Zhang, A node-based smoothed point interpolation method (ns-pim) for thermoelastic problems with solution bounds, Int. J. Heat Mass Transf. 52 (2009) 1464–1471. [17] X. Cui, Z. Li, H. Feng, S. Feng, Steady and transient heat transfer analysis using a stable node-based smoothed finite element method, Int. J. Therm. Sci. 110 (2016) 12–25. [18] S. Wu, G. Liu, X. Cui, T. Nguyen, G. Zhang, An edge-based smoothed point interpolation method (es-pim) for heat transfer analysis of rapid manufacturing system, Int. J. Heat Mass Transf. 53 (2010) 1938–1950. [19] E. Li, Z. He, X. Xu, An edge-based smoothed tetrahedron finite element method (es-t-fem) for thermomechanical problems, Int. J. Heat Mass Transf. 66 (2013) 723–732. [20] S. Feng, X. Cui, A. Li, G. Xie, A face-based smoothed point interpolation method (fs-pim) for analysis of nonlinear heat conduction in multi-material bodies, Int. J. Therm. Sci. 100 (2016) 430–437. [21] H. Zhang, S. Han, L. Fan, D. Huang, The numerical manifold method for 2d transient heat conduction problems in functionally graded materials, Eng. Anal. Boundary Elem. 88 (2018) 145–155. [22] I. Babuška, G. Caloz, J.E. Osborn, Special finite element methods for a class of second order elliptic problems with rough coefficients, SIAM J. Numer. Anal. 31 (1994) 945–981. [23] I. Kergren, Kenan Babuška, U. Banerjee, Stable generalized finite element method and associated iterative schemes; application to interface problems, Comput. Methods Appl. Mech. Eng. 305 (2016) 1–36. [24] T. Strouboulis, K. Copps, I. Babuška, The generalized finite element method, Comput. Methods Appl. Mech. Eng. 190 (2001) 4081–4193. [25] C. Duarte, I. Babuška, J.T. Oden, Generalized finite element methods for threedimensional structural mechanics problems, Comput. Struct. 77 (2000) 215– 232. [26] Y. Abdelaziz, A. Hamouine, A survey of the extended finite element, Comput. Struct. 86 (2008) 1141–1151. [27] E.A. Munts, S.J. Hulshoff, R. De Borst, The partition-of-unity method for linear diffusion and convection problems: Accuracy, stabilization and multiscale interpretation, Int. J. Numer. Meth. Fluids 43 (2003) 199–213. [28] R.L. Taylor, O.C. Zienkiewicz, E. Oñate, A hierarchical finite element method based on the partition of unity, Comput. Methods Appl. Mech. Eng. 152 (1998) 73–84. [29] P. O’Hara, C. Duarte, T. Eason, Generalized finite element analysis of threedimensional heat transfer problems exhibiting sharp thermal gradients, Comput. Methods Appl. Mech. Eng. 198 (2009) 1857–1871. [30] P. O’Hara, C. Duarte, T. Eason, Transient analysis of sharp thermal gradients using coarse finite element meshes, Comput. Methods Appl. Mech. Eng. 200 (2011) 812–829. [31] P. O’Hara, C. Duarte, T. Eason, J. Garzon, Efficient analysis of transient heat transfer problems exhibiting sharp thermal gradients, Comput. Mech. 51 (2013) 743–764.

M. Iqbal et al. / International Journal of Heat and Mass Transfer 149 (2020) 118969 [32] S. Soghrati, A.M. Aragón, C. Armando Duarte, P.H. Geubelle, An interfaceenriched generalized FEM for problems with discontinuous gradient fields, Int. J. Numer. Meth. Eng. 89 (2012) 991–1008. [33] S. Soghrati, H. Ahmadian, 3D hierarchical interface-enriched finite element method: implementation and applications, J. Comput. Phys. 299 (2015) 45–55. [34] F. Van der Meer, R. Al-Khoury, L. Sluys, Time-dependent shape functions for modeling highly transient geothermal systems, Int. J. Numer. Meth. Eng. 77 (2009) 240. [35] A. Cosimo, V. Fachinotti, A. Cardona, An enrichment scheme for solidification problems, Comput. Mech. 52 (2013) 17–35. [36] A. Coppola-Owen, R. Codina, Improving eulerian two-phase flow finite element approximation with discontinuous gradient pressure shape functions, Int. J. Numer. Meth. Fluids 49 (2005) 1287–1304. [37] C. Zeller, B. Surendran, M.F. Zaeh, Parameterized extended finite element method for high thermal gradients, J. Comput. Des. Eng. 5 (2018) 329–336.

13

[38] M. Mohamed, M. Seaid, J. Trevelyan, O. Laghrouche, A partition of unity FEM for time-dependent diffusion problems using multiple enrichment functions, Int. J. Numer. Meth. Eng. 93 (2013) 245–265. [39] M.S. Mohamed, M. Seaïd, J. Trevelyan, O. Laghrouche, Time-independent hybrid enrichment for finite element solution of transient conduction– radiation in diffusive grey media, J. Comput. Phys. 251 (2013) 81–101. [40] M. Malek, N. Izem, M.S. Mohamed, M. Seaid, O. Laghrouche, A partition of unity finite element method for three-dimensional transient diffusion problems with sharp gradients, J. Comput. Phys. 396 (2019) 702–717. [41] M. Iqbal, H. Gimperlein, M.S. Mohamed, O. Laghrouche, An a posteriori error estimate for the generalized finite element method for transient heat diffusion problems, Int. J. Numer. Meth. Eng. 110 (2017) 1103–1118. [42] M.S. Mohamed, Numerical aspects of the PUFEM for efficient solution of Helmholtz problems Ph.D. thesis, Heriot-Watt University, 2010.