Partition of Unity Finite Element Analysis of Nonlinear Transient Diffusion
Problems Using p-Version Refinement

We propose a high-order enriched partition of unity finite element method for linear and nonlinear time-dependent diffusion problems. The solution of this class of problems often exhibits non-smooth features such as steep gradients and boundary layers which can be very challenging to recover using the conventional low-order finite element methods. A class of steady-state exponential functions has been widely used for enrichment and its performance to numerically solve these challenges has been demonstrated. However, these enrichment functions have been used only in context of the standard h-version refinement or the so-called q-version refinement. In this paper we demonstrate that the p-version refinement can also be a very attractive option in terms of the efficiency and the accuracy in the enriched partition of unity finite element method. First, the transient diffusion problem is integrated in time using a semi-implicit scheme and the semi-discrete problem is then integrated in space using the p-version enriched finite elements. Numerical results are presented for three test examples of timedependent diffusion problems in both homogeneous and heterogeneous media. The computed results show the significant improvement when using the p-version refined enriched approximations in the finite element analysis. In addition, these results support our expectations for a robust and high-order accurate enriched partition of unity finite element method.


Introduction
Transient diffusion equations have been used in many applications in physics and engineering, for instance, to describe cooling down of molten glass or heat transfer in enclosures. In glass manufacturing, a hot melt of glass is cooled down to room temperature. This annealing must be monitored carefully to avoid excessive temperature differences, which may affect the quality of the product or even lead to cracking [1,2]. To control the annealing process, the transient diffusion equations may be used to predict accurately the temperature evolution in the glass. In addition, transient diffusion equations have also been used to model several problems in thermal radiation heat transfer [3] and optical tomography [4] among others. In general, thermal radiation has to be modeled by equations that involve the direction and frequency dependent thermal radiation field due to the energy transport by photons. However, using asymptotic expansions, the full radiative transfer equation can be replaced by a class of time-dependent diffusion equations equipped with Robin boundary conditions that depend on space and time but not direction. Such practical transient diffusion problems are not trivial to simulate since the geometry can be complex and internal source/sink terms may produce steep gradient propagating along the computational domain. It is well known that unstructured grids can be highly advantageous based on their ability to provide local mesh refinement near important diffusion features and structures. As a consequence, the ability to provide local mesh refinement where it is needed leads to improved accuracy for a given computational cost as compared to methods that use structured meshes.
Over the past two decades the Partition of Unity Finite Element (PUFE) method was developed within the general framework of the generalized finite element method, see for example [5,6] and further references are therein. In these references, it has been shown that the use of interpolation functions with the conventional partition of unity property allows the approximation space in a numerical approximation, such as the Finite Element (FE) method, to be enriched by the inclusion of other functions suitable for the problem under study. This approach has also been classified as a variety of the well-established Trefftz methods. The PUFE method is considered to be a powerful technique when solving engineering problems such as wave propagation in both steady-state regimes [7] and time-dependent domains [8,9]. The PUFE method was first introduced in [10] for solving boundary-value problems using enrichment techniques. The method allows functions describing the local behavior of the numerical solution to be incorporated in the approximation without violating the boundary conditions. These techniques have demonstrated major improvements in the convergence rate as well as a significant reduction in the total number of degrees of freedom for the problem under study, see for instance [11]. The robustness and efficiency of the PUFE method for twoand three-dimensional time-dependent problems were discussed in [12,13] among others.
The standard high-order finite element methods have the potential to exhibit a great level of accuracy and robustness for solving partial differential equations in complex geometries. However, for problems with sharp gradients and/or steep boundary layers, the FE method would require high computational effort especially for time-dependent problems, see for example [14]. On the other hand, the main idea of the PUFE method is to combine the classical FE method with special functions (known by enrichment functions) while retaining the best properties of the FE method, compare [10]. Several research studies have been carried out in the literature to develop a varied choice of special enrichment functions for different applications in science and engineering. For example, the harmonic functions have been used in [5] for solving the Laplace and Helmholtz equations using the PUFE method. For the wave equation, the FE method has been enriched in space while using implicit or explicit time integration schemes [8,9]. Multiple enrichment functions have been proposed in [12] for solving the two-dimensional heat equation using linear finite elements and the obtained results have been compared to those obtained using the hversion FE method. The extension of this PUFE method to three-dimensional problems has been presented in [13] using unstructured tetrahedral meshes. Other recent works on the PUFE method include application to transient conduction-radiation problems [15,16], to nonlinear boundary-value problems [17] and uncertainty quantification in acoustic waves [18]. The PUFE method presented in these studies can be applied as an effective tool to solve many practical engineering problems with respect to complex geometries and complicated solution features such internal and boundary layers. more advanced for conventional FE methods than for PUFE methods. In the present work, we implement the enrichment functions for the P 3 finite elements to reconstruct a high-order PUFE method for linear and nonlinear transient diffusion problems. A notable feature of the approach is that the same (enriched) approximation space is used at each time step. The more rapidly decaying Gaussian enrichment functions are useful for early time steps around localized diffusion sources, while other slowly decaying enrichment functions become important both in the far field and more generally as a steady-state is approached. Therefore, there is no need to time-dependent shape functions to handle the transient nature of the diffusion problem under study. The main focus of our work is to investigate the impact of using highorder Lagrangian polynomials for solving the diffusion problem in two space dimensions using both linear and nonlinear problems. To the best of our knowledge, solving boundary-value problems using high-order partition of unity finite element methods is presented for the first time. Several numerical examples are presented to demonstrate the performance of the proposed PUFE method to accurately solve the nonlinear transient diffusion equations. Numerical results presented in the current study demonstrate high resolution of the proposed PUFE and confirm its capability to provide highly accurate solutions for linear and nonlinear diffusion problems. We demonstrate through several numerical examples that the proposed p-version PUFE method can recover the sharp solution gradients on coarse meshes and with much fewer degrees of freedom compared to the standard FE method. Thus, a significant reduction in the computational requirements is achieved without compromising on the solution accuracy of the diffusion problem at hand. This paper is organized as follows. In Section 2 we present the governing equations for transient diffusion problems, the initial and boundary conditions and the weak formulation. The p-version finite element partition of unity method is formulated in Section 3. The time-independent Gaussian enrichment functions are also discussed in this section. Numerical test examples are considered in Section 4 to investigate the robustness and the efficiency of the proposed approach. We present computational results for a diffusion problem with known analytical solution to quantify the errors for the proposed method and for a diffusion problem in cracked domains to examine the ability of the method to handle complex domains. We also solve a nonlinear heat transfer in a functionally graded material formed from Zirconium dioxide (ZrO 2 ) and Titanium alloy (Ti 6Al-4V). Finally, conclusions are summarized in Section 5.

High-Order p-Version Finite Element Method
The aim in this study is to develop a highly accurate p-version finite element method for solving nonlinear diffusion problems of this form @u @t À r Á j u ð Þru ð Þ¼f ðt; xÞ; ðt; xÞ 2 ½0; T Â ; au þ j u ð Þru Á n ¼ gðt; xÞ; ðt; xÞ 2 ½0; T Â À; uð0; xÞ ¼ u 0 ðxÞ; x 2 ; (1) where & R 2 is the computational domain with piecewise smooth boundary Γ, [0,T] is the time interval, x = (x,y) ┬ is the space coordinates, t is the time variable, and n is the outward unit normal on Γ. In (1), u(t, x) is the unknown solution, κ(u) is the nonlinear diffusion coefficient, α is the boundary coefficient, u 0 (x) is the initial condition, f(t,x) and g(t,x) are respectively, the internal and the boundary source terms which may depend on the solution u as well. The boundary-value problem has been widely used in the literature to model many applications in science and engineering including heat transfer, image processing and many reactiondiffusion processes in biology and chemistry. In addition, numerical solutions of the nonlinear diffusion problem described by (1) can be computationally demanding especially with the presence of steep gradients, boundary layers and sharp moving fronts, see for example [19].
To integrate the Eq. (1) in time, we divide the time interval [0,T] into equi-distributed subintervals [t n ,t n + 1 ] with length Δt = t n + 1 − t n for n ¼ 0; 1; . . . . We also use the notation w n to denote the value of a generic function w at time t n . Hence, applied to the system (1) a semi-implicit time stepping scheme yields u nþ1 Àu n Dt À r Á j u n ð Þru nþ1 ð Þ¼f nþ1 ðxÞ; x 2 ; au nþ1 þ j u n ð Þru nþ1 Á n ¼ g nþ1 ðxÞ; x 2 À; u 0 ðxÞ ¼ u 0 ðxÞ; Note that in the considered semi-implicit scheme (2), the nonlinear terms are calculated at the current time t n and therefore only linear systems of algebraic equations are required to update the solution at the next time t n+1 . Other semi-implicit time integration schemes can also be used in the proposed p-version finite element method with major conceptual modifications.
The starting point in the finite element method is the weak formulation of the corresponding semidiscrete Eq. (2). We proceed as in the conventional finite element formulations by multiplying the first equation in (2) by a weighting function v and then integrating over Ω. This results in Z Using the divergence theorem, Eq. (3) can be reformulated as Z Substituting the boundary condition from (2) into (4), we obtain Z Thus, the statement of weak formulation of the diffusion problem (1) states: find u ∈ H 1 (Ω) such that Z where H 1 (Ω) denotes the set of square integrable functions whose first derivatives are also square integrable. To solve the above weak form using the finite element method, first the domain Ω is discretized. To perform this step, we generate a triangulation h & of a total number N e of triangular elements T i . Here, the index i referring to the ith element T i in the computational domain h ¼ [ N e i¼1 T i . The conforming finite element space for the solution is given by where P k (T i ) is the set of polynomials of degree ≤ k defined on the element T i . Next, we formulate the finite element solution to U as u n i f i ðxÞ; (8) where N d is the number of solution mesh points in the partition Ω h and u n i are the corresponding nodal values of u n h ðxÞ defined by u n i ¼ u n h ðx i Þ where fx i g Nd i¼1 are the set of solution mesh points in the partition Ω h . In (8), are the set of global nodal basis functions of V h characterized by the property f i (x j ) = δ ij with δ ij denoting the Kronecker symbol. We introduce fx 1 ; . . . ; x N p g as the set of N p nodal points in the element T i . Hereafter, unless otherwise stated, the subscripts h and i are used to refer to coefficients associated with the whole mesh Ω h and a mesh element T i , respectively. The approximation space is then defined as Hence, assembling the global finite element matrices, the fully-discrete formulation of (6) can be written in a compact form as where u nþ1 È É is the N p -valued vector of unknown solutions with entries u nþ1 i , u n f g is the N p -valued vector of known solutions at time t n with entries u n i , M ½ is the N p × N p -valued mass matrix with entries m ij , S ½ is the N p × N p -valued stiffness matrix with entries s ij , and R ½ is the N p × N p -valued remainder matrix with entries r ij defined as In (9), f n+1 is the N p -valued vector of the source term and r n+1 is the N p -valued remainder vector the entries of which are given by Z f nþ1 f i d; and It is evident that, to update the solution u nþ1 È É from (9) one has to solve at each time step the linear system of algebraic equation where In the present study we consider Lagrangian polynomials of order up to k = 3. Using the barycentric coordinates η 1 , η 2 and η 3 = 1 − η 1 − η 2 of the reference element b T , the basis functions associated with P 1 , P 2 and P 3 are defined as follows: • Linear shape functions (k = 1): • Quadratic shape functions (k = 2): • Cubic shape functions (k = 3): Note that, in all the cases considered in this paper, the quadratic elements are formed from linear elements but adding a vertex/node at the mid edge of all three sides while cubic elements by adding two vertices 1 3 and 2 3 at each edge and a further one to the center of the linear element, see Fig. 1 for an illustration. Needless to mention that the partition of unity method proposed in this study can handle curved elements without major conceptual modifications in its implementation.

High-Order p-Version Partition of Unity Method
As stated in [12,13,17], solving diffusion problems with sharp moving fronts and steep boundary layers, the conventional FE method needs large number of nodes in order to be able to capture these solution features. Alternatively, to improve the polynomial-based approximations, we may enrich the P k discretization using some well-defined functions with capabilities to capture sharp moving fronts and steep boundary layers in the computational domain. Examples of these enrichment functions include exponential Gaussian functions [12] and hyperbolic tangent functions [16]. The basic idea of these techniques is to build basis functions that are the product of the polynomial functions and the considered enrichment functions. In this case no extra nodes are added to the element in the computational mesh. However, the number of degrees of freedom at each node is increased by the number of enrichment functions used. Hence, given G q,k (x) a class of well defined functions to be used to enrich the FE approximation (8), the partition of unity method consists of expanding the nodal values u n i in (8) as a product of polynomial functions P N p i¼1 f i ðxÞ and a sum of enrichment functions P Q q¼1 G q;k ðxÞ as where Q is the total number of enrichment functions, k the order of the Lagrangian polynomials in the P k method, and U n i;q;k the degree of freedom associated with the product f i (x)G q,k (x). In the current study, we propose exponential enrichment functions defined by T for the P k finite elements with k = 1 (left), k = 2 (middle) and k = 3 (right) used in the proposed PUFE method where x c is a control point used to define the location of the enrichment function G q,k (x) in the computational domain, γ and σ are parameters used to define the amplitude and variance of the function G q,k (x) accordingly. Note that the control points x c and the parameters γ and σ need to be selected a priori depending of sharp gradients and steep boundary layers in the diffusion problem under study. It is also pointed out that the proposed approach is flexible such that the functions G q,k (x) in the expression (13) can easily be replaced with other functions without major modifications in the implementation. Note that similar enrichment functions have been considered in [20] using a posteriori error estimate for the partition unity finite element method. In our simulations, a sensitivity analysis of the results with different enrichment parameters is performed a priori and results are presented only for the most suitable selections.
Hence, the finite element approximation space can be defined as For convenience, we denote the multiplication of the enrichment functions by the nodal Lagrangian polynomial shape functions by Φ i which we treat as a new shape function in the high-order p-version partition of unity method. Hence, the new basis functions Φ i satisfy the relation It is also clear that the total number of degrees of freedom in the proposed method is M = N p × Q. Hence, using the new basis functions in the weak form (6) and assembling the global partition of unity finite element matrices, the fully-discrete approximation can be formulated in a compact form as where U nþ1 È É is the M-valued vector of unknown solutions with entries u nþ1 i;q , U n f g is the M-valued vector of known solutions at time t n with entries u n i;q , M ½ is the M × M-valued mass matrix with entries m ij,q , S ½ is the M × M-valued stiffness matrix with entries s ij,q , and R ½ is the M × M-valued remainder matrix with entries r ij,q defined as In (16), F nþ1 is the M-valued vector of the source term and R nþ1 is the M-valued remainder vector the entries of which are given by Z f nþ1 È i d; and As in the conventional FE method, to update the solution U nþ1 È É from (16) one has to solve at each time step the linear system of algebraic equation where Similar to the conventional finite element approximations (10) and (11), the entries (17) and (18) are evaluated using standard Gauss quadratures but with a relatively high number of integration points. In this paper the convergence of the integration scheme is ensured by increasing the number of integration points and checking the convergence in the linear solver. It should also be stressed that the enrichment functions are applied in a global sense at all nodes in the computational domain. Hence, the continuity is guaranteed and the implementation can performed in an already existing FE code with minimum changes.

Numerical Results and Applications
In this section we examine the accuracy and performance of the proposed high-order p-version partition of unity method using three test examples. The first example solves a linear transient diffusion problem with an available exact solution that can be used to quantify errors in the PUFE method. The second and third examples consider the problem of transient diffusion heat transfer in a cracked plate and a nonlinear diffusion in a plate with several holes whereas, the third example studies a nonlinear heat transfer in a functionally graded material formed from Zirconium dioxide (ZrO 2 ) and Titanium alloy (Ti 6Al-4V). This last example is considered to illustrate the ability of the developed PUFE method to deal with more complicated transient diffusion problems. The results obtained using the PUFE method are compared to those obtained with the conventional FE method using linear P 1 , quadratic P 2 and cubic P 3 elements. To evaluate elementary matrix entries for PUFE and FE methods, all the integrals are evaluated numerically using standard Gaussian quadrature. The number of integration points is chosen to be 20 integration points per spatial direction in the PUFE method which small enough so that the results are not affected by the integration errors. Notice that in practice, the conventional FE method does not require high numbers of integration points. A direct solver using LU decomposition is used to solve the resulting linear systems of algebraic equations. All the computations are performed on an Intel ® Core(TM) i5-4200U @ 1.60 GHz with 8 GB of RAM. The codes only take the default optimization of the machine, i.e., they are not parallel codes.

Accuracy Test
To evaluate the accuracy of the proposed PUFE method in terms of accuracy and efficiency we consider a test example of linear transient diffusion problems with known analytical solution introduced in [12]. Here, we solve the Eq. (1) in the squared domain Ω = [0, 2] × [0, 2] with constant coefficient κ(u). The source term f(t, x), the boundary function g(t, x) and the initial condition u 0 (x) are explicitly calculated such that the exact solution of the problem (1) is given by uðt; x; yÞ ¼ x 20 y 20 ð2 À xÞ 20 ð2 À yÞ 20 ð1 À expðÀjtÞÞ; (20) To quantify the errors in this problem, we consider the relative L 2 -error defined as where ‖·‖ L 2(Ω) is the L 2 norm, u h and u are respectively, the computed and exact solutions. In all presented computations the parameter α = 1, the time step Δt = 0.1 and results are presented at two different instants t = 0.5 and 2. To study the combination of an increased enrichment number with a refined polynomial order, different meshes are used. We consider a coarse mesh of 16 linear elements and 13 nodes where we refer to the element size in this mesh as h 1 . We then preform a h-refinement to create other two meshes namely, h 1 /2 and h 1 /4. Thus, the number of elements and nodes for h 1 /2 is increased to 40 and 29, respectively. The respective numbers for h 1 /4 are 184 elements and 109 nodes. All the considered linear meshes are displayed in Fig. 2.
Next, a further p-refinement is also performed on the meshes h 1 , h 1 /2 and h 1 /4 so that k is increased to 2 and then to 3 while the number of elements is unchanged. As can been seen in Tab. 2, performing h-or p-refinements always lead to improved accuracy. This is consistently observed for the FE and PUFE methods equally. The convergence rate in the FE method is consistent with the standard convergence rates for k = 1, 2 and 3 so that the rate is higher for an increased k. It is also clear in Tab. 2 that adding enrichment functions always improves the computed error. In general, the improvement seems to be larger when enriching higher order P k elements. For example, at t = 0.5 with linear elements and for the coarse Mesh 1 (h 1 ) the error is reduced from 6.87E-01 to 2.33E-01 when adding two enrichment functions whereas for cubic elements the error is reduced from 1.03E-01 to 1.14E-02. This is consistently observed across all the cases and at all the considered time instants. In fact the error with cubic elements is consistently improved by an order of magnitude when adding one enrichment function while the error is roughly halved with linear elements. Furthermore, for the same mesh, when no enrichment is considered the p-refinement reduces the error at a smaller rate than the enriched cases. Adding three enrichment functions to Mesh 3 (h 1 /4), the error is   1  13  16  29  40  109  184  2  41  16  97  40  401  184  3  85  16  205  40  877  184 reduced by two orders of magnitude when refining k. For example, at t = 0.5 the error is reduced from 1.28E-02 to 2.86E-04 when increasing k from 1 to 2 and it is reduced again to 5.07E-05 when k is increased to 3. The results suggest that adding enrichment functions can be more efficient in reducing the error compared to h-and p-refinements. For example, refining Mesh 1 (h 1 ) to Mesh 3 (h 1 /4) increases the number of nodes from 13 to 109 while the error is reduced from 6.87E-01 to 9.52E-02. Similarly, the p-refinement increases the number of nodes to 85 while the error is reduced to 1.03E-01. However, adding 3 enrichment functions, it increases the number of degrees of freedom only to 39 which is the corresponding number to the number of nodes in the previous cases, while dropping the error to 1.33E-01. In general a comparable reduction in the error is achieved in all cases but the increase in the number of degrees of freedom is significantly smaller with the enrichment.
The computational cost referred to by CPU time in Tab. 3 is split into the time needed to build the linear system and that needed to solve the linear system. Obviously, in both h and p-refinements, adding enrichment functions leads to an increase in the CPU time needed to solve the problem which is clearly reflected in Tab. 3. For smaller systems the solution total time is dominated by the time needed to build the linear system. As the size of the linear system becomes larger the solve time becomes more dominant. It should be noted that it is possible to significantly reduce the time needed to build the linear system from the second time step as the system matrix remains unchanged and only the right-hand side must be updated from the second time step onward. Similarly, it is possible to decompose the linear system at the first time step and then repeatedly reuse the same decomposition at later time steps after updating the right-hand side. However, in this work we repeatedly build and solve the linear system as in the following examples we also consider nonlinear problems where it becomes necessary to rebuild the system at every time step.
In Tab. 3 we also show the condition numbers for all the spatial discretizations considered in this test example and summarized in Tab. 2. Both types of the considered mesh refinements as well as the addition of enrichment functions leads to higher condition numbers of the associated matrices in the linear systems. However, the condition number associated with an enriched system is significantly higher than the standard FE method system. For the Mesh 3 (h 1 /4) and with k = 3, adding three enrichment functions increases the condition number from 16.09 to 9.52E+5. Indeed, adding more enrichment functions will further deteriorate the condition number and it may lead to ill-conditioned systems. This aspect of the enrichment is well studied in the literature, see for example [12,21]. Table 2: Relative errors obtained using the PUFE method with k = 1, 2 and 3 on the considered meshes for the accuracy test problem using different numbers of enrichments Q at two different instants t = 0.5 and 2 using the time step Δt = 0.1

Diffusion in a Cracked Domain
In the second numerical test, we evaluate the linear heat diffusion (1) in a squared plate with a crack. The considered thermal conductivity is κ = 0.01 kg m/K s 2 while the time step is Δt = 0.1 s. The configuration of the computational domain along with the location and size of the crack is shown in Fig. 3. Initially, the plate is at the ambient temperature u 0 = 300 K and an instantaneous heat source is defined by Here, the source is rectangular located below the crack and the considered simulation time span is [0, 10 s]. On the boundary α = 1 kg/K s 2 and g(t, x) = 300 K. Since there is no analytical solution for this example, we solve the problem on a highly refined mesh grid for which the solution convergence is ensured. The reference finite element mesh is composed of 16948 linear elements and 8674 nodes. To solve this problem we consider the PUFE method with two enrichment functions using the coarse mesh of 139 linear elements and 89 nodes illustrated in Fig. 3 alongside the problem configuration. To check the PUFE method convergence we use p-refinement where the polynomial order is increased to 2 and then 3 so that the number of nodes is increased to 316 and then 682. The results are then compared to those obtained using the FE method on the reference mesh. Fig. 4 shows the comparisons between the numerical results obtained for the temperature using the PUFE method with k = 1, 2 and 3 and the reference solution. The results are compared at the time instants t = 1 s, 5 s and 10 s and at two different cross sections namely, line AB and line CD as indicated on the domain configuration in Fig. 3. It is clear that, as the polynomial order increases, the solution obtained using the PUFE method converges toward the reference solution. Once k = 3 is reached, solution obtained using the PUFE method and reference solution are well matched. For demonstration reasons, Fig. 5 depicts the temperature patterns obtained using the FE method with k = 3, the PUFE method with k = 3 and the reference solution at the time instants t = 1 s, 5 s and 10 s. It should also be noted that the CPU time to construct and solve the reference solution is on average 329 s compared to only 11.61 s for the PUFE method with k = 3.

Nonlinear Diffusion Problems
In the final class of numerical test examples, we study the nonlinear problem (1) using the p-refinement approach with the PUFE method. This problem studies the heat transfer in a functionally graded material formed from Zirconium dioxide (ZrO 2 ) and Titanium alloy (Ti 6Al-4V) which is common in biological applications, see for example [17]. The considered domain is a square with a side length of 20 cm and a central circular hole with a 4 cm radius and four smaller circular holes with 1 cm radius arranged uniformly around the domain center with one hole located at each corner. The problem configuration is illustrated in Fig. 6 with the dimensions included. Because of the symmetry, only a quarter of the domain is used in our simulations. Here, the problem statement consists on solving the following set of equations where κ(u) is the conductivity, c(u) the heat capacity and ρ(u) is the density of the medium. As can be seen from above equations, all these parameters are now temperature-dependent and hence, appears the nonlinearity in the system. The initial temperature and ambient temperatures are u 0 = 300 K and g = 300. The domain is heated with a steady source defined by f x; y ð Þ ¼ 5 Â 10 6 1 þ sinð10pxÞ sinð10pyÞ ð Þ : Figure 6: Configuration of the domain used for nonlinear diffusion problem along with finite element meshes used for the problem of nonlinear heat conduction in functionally graded materials. Numbers of elements and nodes are displayed under each mesh The heterogeneous medium properties are given for the two considered materials by Material 1 (ZrO 2 ): jðuÞ ¼ 1:71 þ 2:1 Â 10 À4 u þ 1:16 Â 10 À7 u 2 ; cðuÞ ¼ 2:74 Â 10 2 þ 7:95 Â 10 À1 u À 6:19 Â 10 À4 u 2 þ 1:71 Â 10 À7 u 3 ; qðuÞ ¼ 3657 ð1 þ aðu À 300ÞÞ 3 ; where aðuÞ ¼ 1:331 Â 10 À5 À 1:89 Â 10 À8 u þ 1:27 Â 10 À11 u 2 : Material 2 (Ti 6Al-4V): kðuÞ ¼ 1:1 þ 1:7 Â 10 À2 u; cðuÞ ¼ 3:5 Â 10 2 þ 8:78 Â 10 À1 u À 9:74 Â 10 À4 u 2 þ 4:43 Â 10 À7 u 3 ; where aðuÞ ¼ 7:43 Â 10 À6 þ 5:56 Â 10 À9 u À 2:69 Â 10 À12 u 2 : The volume fraction of each constituent varies linearly with the radial distance from ZrO 2 to Ti 6Al-4V. Similar material properties were also investigated in [22] using an enriched finite element method. It should be stressed that this problem is more challenging than the previous test examples not only because of the nonlinearity of the governing equations but also due to the presence of moving thermal fronts in the computational domain.
Since no analytical solutions are available for this class of equations, the problem is first solved using the FE method on a highly refined mesh to obtain a reference solution. The problem is then solved again on a coarse mesh using the PUFE method using Q = 2 and k = 1. The two meshes are shown in Fig. 6. As in the previous test example, a p-refinement procedure is used to ensure the convergence of the solution obtained using the PUFE method to the reference solution, where the value of k is increased to 2 and then 3 so that the number of nodes is increased to 190 then to 414. Fig. 7 compares the PUFE method with k = 1, 2 and 3 to the reference solution at two different cross-sections and at three different time instants t = 60 s, 90 s and 120 s. Again it is clear that the PUFE method with k = 3 has converged to the same solution as the one obtained using the FE method on a fine mesh. For comparison purposes, we display in Fig. 8 the two-dimensional snapshots of the temperature obtained using the PUFE method with k = 3 compared to the FE method k = 3 and the reference solutions. As expected the temperature builds up at the domain parts further away from the hole and drops rapidly in the area near the hole. Thermal fronts forms in the domain due to the difference between the heat isolation boundary conditions and the high heat release rates near the source. These moving fronts and their change in time are efficiently captured using the proposed PUFE method despite using a coarse mesh grid. From the computed results we can observe that the complicated heat structures in the computational domain captured by the enriched finite element method. It is worth remarking that all these features have been achieved using triangular meshes far coarser than those required for conventional finite element methods.

Concluding Remarks
A high-order enriched partition of unity finite element method is presented for linear and nonlinear timedependent diffusion problems in unstructured triangular meshes. Unlike, the standard h-version refinement or the so-called q-version refinement, the p-version refinement can be a very attractive alternative in terms of the efficiency and the accuracy in the enriched partition of unity finite element method. In this study, a class of approximate Gaussian functions describing the diffusion decay is embedded in the P k finite element shape functions with k varies up to 3. Increasing the number of the enrichment functions in the q-version refinement would ameliorate the computational results. However, the enrichment functions can be increased only to certain limit after which the system becomes strongly ill-conditioned. In contrast, the p-version refinement can significantly reduce the computational costs when compared to other available enrichment techniques where the full linear system must be updated and resolved repeatedly as those discussed in [23,24] among others. In addition, the advantage of this method compared to h-version refinement is its ability to locally refine the solution by adapting the enrichment without need to generate new meshes. The favorable performance of the developed high-order partition of unity finite element method has been demonstrated using a series of numerical examples, including a nonlinear heat transfer in a functionally graded material formed from Zirconium dioxide (ZrO 2 ) and Titanium alloy (Ti 6Al-4V). One obvious next step is to test the accuracy of the proposed method for high order P k finite elements with k > 3. Future work will also concentrate on the extension of the high-order enriched partition of unity finite element method to computational heat transfer problems in three spatial dimensions implemented on unstructured tetrahedral meshes.