Efficient numerical analysis of transient heat transfer by Consecutive-Interpolation and Proper Orthogonal Decomposition

 Abstract—The consecutive-interpolation technique has been introduced as a tool enhanced into traditional finite element procedure to provide higher accurate solution. Furthermore, the gradient fields obtained by the proposed approach, namely consecutive-interpolation finite element method (CFEM), are smooth, instead of being discontinuous across nodes as in FEM. In this paper, the technique is applied to analyze transient heat transfer problems. In order increase time efficiency, a modelreduction technique, namely the proper orthogonal decomposition (POD), is employed. The idea is that a given large-size problem is projected into a small-size one which can be solved faster but still maintain the required accuracy. The optimal POD basis for projection is determined by mathematical operations. With the combination of the two novel techniques, i.e. consecutive-interpolation and proper orthogonal decomposition, the advantages of numerical solution obtained by CFEM are expected to be maintained, while computational time can be significantly saved.


Efficient numerical analysis of transient heat transfer by Consecutive-Interpolation and Proper Orthogonal Decomposition
Nguyen Ngoc Minh, Nguyen Thanh Nha, Truong Tich Thien * , Bui Quoc Tinh  Abstract-The consecutive-interpolation technique has been introduced as a tool enhanced into traditional finite element procedure to provide higher accurate solution.Furthermore, the gradient fields obtained by the proposed approach, namely consecutive-interpolation finite element method (CFEM), are smooth, instead of being discontinuous across nodes as in FEM.In this paper, the technique is applied to analyze transient heat transfer problems.In order increase time efficiency, a modelreduction technique, namely the proper orthogonal decomposition (POD), is employed.The idea is that a given large-size problem is projected into a small-size one which can be solved faster but still maintain the required accuracy.The optimal POD basis for projection is determined by mathematical operations.With the combination of the two novel techniques, i.e. consecutive-interpolation and proper orthogonal decomposition, the advantages of numerical solution obtained by CFEM are expected to be maintained, while computational time can be significantly saved.

INTRODUCTION
NE may encounter heat transfer problems in many human activities.For example, all three types of heat transfer can be found in cooking, i.e conduction, convection and radiation.Design of air-conditioning system is usually based on knowledge of heat convection.Day by day, the Earth is receiving heat from the Sun by thermal radiation.In industry, heat transfer analysis is required in many fields of engineering, such as mechanical engineering, electrical engineering, aeronautical engineering, etc.However, analytical solutions are only available for some specific problems, most of which are described with relatively simple geometry and boundary conditions.When it comes to deal with complicated geometries and/or boundary conditions, which are usually the cases of engineering applications, numerical analysis seems to be a more practical approach.Currently, the standard finite element (FEM) [1] has been widely used for heat transfer problems due to its simplicity and reasonable accuracy.However, several shortcomings of the method have been pointed out, see [2].The FEM shape function is C0 continuous, resulting in non-physical discontinuity of gradient fields, e.g.temperature gradient in case of heat transfer problems.
As alternatives to FEM, various other methods have been proposed for heat transfer analysis, such as the Boundary Element Method (BEM) [3] and the class of meshfree method [4,5].On the other hand, amendments that can be integrated into FEM was also suggested to overcome the weakness while keeping the familiar FEM framework.In recent years, the consecutive-interpolation procedure (CIP) has been introduced as an enhancement for traditional FEM, to develop the O so-called Consecutive-interpolation Finite Element Method (CFEM).In CFEM, the continuity of gradient fields is improved by taking the averaged nodal gradients into interpolation.Interestingly, the number of degrees of freedom are equal to that of FEM, given the same mesh.The CFEM was first investigated for two-dimensional linear elastic problems [6,7] and later was further developed for heat transfer analyis in two-dimensional space [8] and three-dimensional space [9].A general formulation which allows application of CFEM in a wide range of finite elements was also proposed in [9].
As a model reduction method, the Proper Orthogonal Decomposition (POD) was introduced to reduce the computational time by projecting the problem of interest to another one which is much smaller in size.Hence, computer memory and elapsed time can be greatly saved.POD has been applied to structural vibration analysis based on experiments [10].Investigation on combination of POD with finite element analysis of heat transfer problems is discussed by [11].
In this study, POD is combined with CFEM to effectively save computational time in the context of three-dimensional transient thermal analysis, such that the applicability of CFEM is further expanded.The proposed procedure is named by CFEM-POD for brevity while the CFEM without POD is mentioned by CFEM.
The paper is organized as follows.After the introduction, a brief review on application of CIP to three-dimensional element is presented in Section 2. The integration of POD into analysis is discussed in Section 3. In Section 4, the efficiency and accuracy of the proposed formulation are investigated by several numerical examples.Concluding remarks are given in Section 5.

CONSECUTIVE-INTERPOLATION FOR
HEAT TRANSFER PROBLEMS

Brief on consecutive-interpolation
Let us consider a 3D body in the domain Ω bounded by Г = Гu + Гt và Гu ∩ Гt = {ø}.In finite element analysis, the domain Ω is discretized into non-overlapping sub-domains Ωe called elements.The points interconnected by the elements called nodes.Each node is associated with a shape function.Any function u(x) defined in Ω can be approximated by a linear combination as (1) Here n is the number of nodes, is the vector containting nodal values and N is the vector of shape functions.By assigning the approximated value at node i as , and the vector of shape functions evaluated at node i as , the average nodal derivatives (similarly for and ) can then be determined by [6,7,9] (2) in which the vector of averaged derivative is calculated by In Eq. ( 3), denotes the derivative of computed in element e. Si is the set containing all the elements connected to node i, and we is a weight function defined by [9] , with being the volume of element .
One well-known shortcoming of the standard FEM is the non-physical discontinuity of gradient fields, e.g.temperature gradient in case of heat transfer analysis.Such drawback can be overcome by taking both the averaged nodal derivatives (and and ) and the nodal values u [i] and into interpolations, following the consecutiveinterpolation procedure (CIP) [9].By means of CIP scheme, the approximation in Eq. ( 1) can be rewritten as In Eq. ( 5), the CIP shape functions is given by , in which , , and iz are the auxiliary functions dependent on the element type.Determination of auxiliary functions used to be bottleneck in application of CIP into finite element analysis, i.e.CFEM.However, a general formulation recently suggested by [9] can be used to determine auxiliary functions for a wide range of standard finite element types.For the sake of completeness, the formulation will be briefly presented here.Let us denote the following terms and ( 6) where n is the number of nodes within the element of interest and Li is the Lagrange shape function associated with the i th node of the element.The functions and can be written by In Eq. ( 8), xi and xj denote the x-coordinate of node i and node j, respectively.Functions and are obtained analogously by replacing xcoordinate in Eq. ( 8) with y-coordinate and zcoordinate, respectively.Figure 1 illustrates the application of CIP approach into the four-node quadrilateral (Q4) element, which results in the namely CQ4 element.Without loss of generality, the scheme is described particularly in an irregular finite element mesh.As shown in Figure 1, the supporting nodes for the point of interest x include all the nodes in the four sets Si, Sj, Sk, Sm, which contain all the adjacent elements that share the nodes i, j, k, m, respectively.Thus, the support domain for a point x in CQ4 element is larger than that in the standard Q4, since it includes not only the nodes of the element in interest but also the nodes of the neighboring elements.Similar observation is reported by [9] for the case of tetrahedral element, as shown in Figure 2. In Eqs. ( 9) to ( 13), k = diag(kxx, kyy, kzz) is the tensor of thermal conductivities; T is the temperature field; Q is the body heat flux; ρ is the density; c is the specific heat capacity; h is the convective coefficient; and Ta is the ambient temperature.Multiplying both sides of Eq. ( 9) with an arbitrary test function δT, then applying Green's theorem and integration by parts, the variational form is obtained as follows (17) in which M is the matrix related to the specific heat capacitance, K is the matrix related to conductivity and convection terms, while F is the vector accounts for heat source and thermal interaction with external environment. (20)

PROPER ORTHOGONAL DECOMPOSITION (POD).
The Proper Orthogonal Decomposition (POD) was initially developed to statistically analyze experimental data.Firstly, a series of snapshots are generated.Each snapshot is actually a vector containing data of system response at a specified time period.An orthogonal basis is then obtained from the snapshots.The orthogonal basis is constructed such that it reduces the size of the problem to be solved, but the required accuracy is still kept.Due to the reduction of problem size, computational cost can be greatly saved.Finally, the full system can then be reproduced from the reduced system without much loss of accuracy.
Denoting the column vector Ti, i = 1, 2, …, d, as the response at the i th time step and d is the total number of time steps, the set of snapshots can be expressed by an n x d matrix, with n being the total number of degrees of freedom Denoting , the snapshot matrix is rewritten by (23) Since is orthogonal, V can be calculated by The snapshot matrix can be approximated by a truncated basis where is the first k columns of (25) The truncation error of approximation is determined by (26) Due to orthogonality, is an identity matrix.The key point in POD procedure is to determine k such that truncation error less than a given tolerance.Similar to [12], the cumulative "energy" coefficient is defined by Here, e(k) represents the ratio of "energy" in the total first k modes with respect to the total "energy".As k increases, the truncation error reduces.Once the POD basis is selected, the following reduced problem can be obtained 4 NUMERICAL EXAMPLES.
In this section, three numerical examples are investigated to demonstrate the effectiveness of the proposed procedure.We denote Q4 for the standard four-node quadrilateral element and TH4 for the tetrahedral element, while CQ4 and CTH4 are the CIP-version of Q4 and TH4, respectively.

Two-dimensional heat transfer
Let us consider a square domain (see Figure 3) of size L x L, where L = π m.On all four boundaries of the square, zero temperature, T = 0 o C, is imposed.Initially, the temperature distribution is given by the following equation: Material properties are given as follows: the mass density ρ = 1 kg/m 3 , the specific heat c = 1 J/(kg o C), and the heat conductivity k = 1 W/(m o C).Under the boundary conditions specified above, the temperature tends to drop down from the initial value to zero as given by follows [13] (31) Two levels of finite element mesh are used in numerical simulation: 20x20 CQ4 elements and 40x40 CQ4 elements (i.e.441 and 1681 degrees of freedom).Firstly, the matrix of snapshots is generated for a time span of t = 0.5s with time increment Δt = 0.02s, i.e. 25 time steps.Next, singular decomposition is calculated for .As shown in Table 1, the first singular value (the largest one) dominates.Thus, it is reasonable to select a POD basis of size k = 3 to approximate the response of the system.The reduced system is obtained using Eq.(30).Finally, solution for time span of t = 3s, i.e. 150 time steps, is computed by the reduced system.The results evaluated by CQ4-POD are compared with both CQ4 and analytical solutions (see Eq. ( 32)).The temperature along the line y=π/2 evaluated at t = 1s, t = 2s and t = 3s are depicted in Figure 4.At t = 3s, temperature is almost zero at every node, indicating a steady state is reached.Note that snapshot matrix is only calculated from t = 0 to t = 0.5s.Hence, the reduced system obtained by POD is able to predict responses taking place after snapshots have been generated.Relative errors between values computed by CQ4 only and by CQ4-POD at t = 1s with respect to analytical solutions are reported in Table 2. Results show that the accuracy of CQ4-POD is almost equivalent to CQ4, despite the fact that the reduced system has only 3 degrees of freedom, much smaller than the full system.Computational time in CQ4-POD (including the time required for generating snapshots) is greatly reduced, especially when finer mesh is used, see Table 3. 4.2 Three-dimensional heat transfer  For numerical simulation, a mesh of 3428 consecutive-interpolation four-node tetrahedral elements (CTH4), i.e. 848 nodes, is used to discretize the domain.The snapshot matrix Tsnap is generated by CTH4 solutions for a time span of 75s with 75 time steps, i.e. time increment Δt = 1s.Based on singular value decomposition of Tsnap, a set of 24 POD bases is chosen (largest singular value is of magnitude 10 3 and the 24 th singular value is of magnitude 10 -9 ).POD procedure is then used to predict temperature changing from t = 0 to t = 750s.Table 4 presents the relative errors between CTH4-POD and CTH4 solutions at t = 125s, t = 500s and t = 750s.The errors are all smaller than 1%.Elapsed time of CTH4-POD is approximately 160s, quite smaller than that of CTH4, which is approximately 176s. Figure 6 depicts the x-component of heat flux, showing that heat flux computed by CTH4 elements is smooth, while the one obtained by TH4 elements (standard FEM) is non-physically discontinous.Hence, CTH4-POD preserves the desirable property of CTH4, such that the nodal gradients are continous.

Figure 2 .
Figure 2. Schematic representation of support domain of CTH4 element (14) Using the approximation in Eq. (5) for both temperature T and test function δT, one gets , (15) (16) where matrix B calculates the derivatives of shape functions R. The discrete form is obtained by substitution of Eqs.(16-17) into Eq.(15) (21) By using singular value decompositions, the matrix Tsnap can be decomposed into three parts as follows (22) where V is an orthogonal matrix of size d x d; D is the rectangular matrix of size n x d containing the singular values; while the n x n matrix stores the orthogonal eigenvectors of .In matrix D, only values along the diagonal are non-negative and named by singular values, while the rest are all zero.In practice, matrix D is sorted such that the singular values are arranged in decreasing order, i.e. , with r = min(n, d).
) which can be solved much faster than Eq.(17) due to the smaller size.The terms in Eq. (28) are determined by , , , ,

Figure 5 .
Figure 5. Example 4.2: Geometry (upper) and one-quarter model (lower) In this example, three-dimensional transient heat transfer in a square plate with a cylindrical hole at center is investigated.The plate is subjected to both convection and Dirichlet boundary conditions, as shown in Figure 5. Due to symmetry, only one-quarter of the plate is modeled.Material properties are given as follows: homogeneous conductivity k = 15W/m o C, density ρ = 7800 kg/m 3 and specific heat capacitance c =

Table 4 .Figure 6 .
Figure 6.Example 4.2: The x-component of heat flux at t = 750s, obtained by CTH4-POD (upper) and TH4 (lower) . The conductivity for this example is set to be k = 100 W/m o C. The inward heat flux is applied by q = 20000 W/m 2 on the curved surface of the middle fin.Convection takes place on the left hand side surface (x = 0) with an ambient temperature of Ta= 300 o C and convection coefficient h = 100 W/m 2 .On the right hand side surface (x = 0.5), temperature is prescribed at T = 300 o C. The density is ρ = 3000 kg/m 3 and specific heat capacitance is c = 125 J/(kg o C).Initially, temperature of the whole domain is at T = 0 o C.

Figure 8 .
Figure 8. Example 4.3: Variation of temperature at point A (see Figure 7) with respect to time

Table 2 .
Example 4.1:Relative errors between CQ4 and CQ4-POD with analytical solutions, at t = 1s