A polygonal scaled boundary finite element method for solving heat conduction problems

This paper presents a steady-state and transient heat conduction analysis framework using the polygonal scaled boundary finite element method (PSBFEM) with polygon/quadtree meshes. The PSBFEM is implemented with commercial finite element code Abaqus by the User Element Sub-routine (UEL) feature. The detailed implementation of the framework, defining the UEL element, and solving the stiffness/mass matrix by the eigenvalue decomposition are presented. Several benchmark problems from heat conduction are solved to validate the proposed implementation. Results show that the PSBFEM is reliable and accurate for solving heat conduction problems. Not only can the proposed implementation help engineering practitioners analyze the heat conduction problem using polygonal mesh in Abaqus, but it also provides a reference for developing the UEL to solve other problems using the scaled boundary finite element method.


Introduction
The analysis of heat conduction has crucial importance in many fields, such as aerospace engineering, mechanical, electrical, and civil engineering.Over the last decades, analytical and numerical methods have been utilized to investigate heat conduction problems.The analytical approaches are limited by the complex geometry and boundary conditions.Numerical methods have become a powerful tool to simulate heat conduction problems.
Due to the equation of transient heat conduction exist space and time variables, we need to discretize the space and the time by the numerical methods.The domain of complex geometry is partitioned into a finite number of non-overlapping subdomains of simplex shapes by introducing the concept of discretization [1].The most widely used computational techniques in discretized space are the finite element method (FEM) [2].
These methods must discretize the whole computational domain.It is often unavoidable that the generation of highly distorted meshes for complex geometries, which leads to difficulty in solving many complex problems.
The polygonal scaled boundary finite element method (PSBFEM) is a novel approach integrating the standard SBFEM and the polygonal mesh technique.This method is flexible in meshing complex geometries, and the use of polygons to discretize the computational domain naturally complements the SBFEM.Polygonal element with more than four sides involves more nodes in their interpolation compared with a conventional FEM.
Simultaneously, they generally exhibit superior solution accuracy [1].Moreover, it is more flexible in the discretization of complex geometry.Therefore, these advantages have further motivated polygonal elements as an alternative to conventional FEM using triangles or quadrilaterals.
Recently, the quadtree mesh generation technique [22][23] is applied to discrete domains.Mesh generation and adaptive refinement of quadtree meshes are straightforward [22].However, the quadtree discretization will lead to additional nodes, usually called 'handing nodes'.Due to hanging nodes between two adjacent elements of different sizes, it is problematic that quadtree meshes are directly used to simulate within the finite element method's framework.However, the SBFEM only discretizes in the boundary of geometry.
Hence, each cell in a quadtree mesh, regardless of hanging nodes, is treated as a generic polygon.The enables the structure of the quadtree to be exploited for efficient computations.The ability to assume any number of sides also enables the SBFEM to discretize curved boundaries better.
In this work, we present a framework for a heat conduction analysis using the SBFEM with polygon/quadtree meshes.We implement steady-state and transient heat conduction analysis of SBFEM by developing a User Element subroutine (UEL).This paper is organized as follows.The governing equations for the heat conduction problem are described in Section 2. The SBFEM of heat conduction problem concept is stated in Section 3. Section 4 outlines the solution procedures.Section 5 describes the implementation of PSBFEM in the heat conduction problems by the Abaqus UEL subroutine.Then, several benchmarks examples are presented in Section 6.Finally, the main concluding remarks of this paper are presented in Section 7.

Governing equations for the heat conduction problem
In this section, the governing equations for two-dimensional transient heat conduction problems are considered.The governing equations without heat sources are written as with the initial conditions 0 ( , , 0) ( , ) and the boundary conditions where  is the density, c is the heat capacity, T is the temperature, T is the temperature change rate,  is the gradient operator, k is the heat conductivity,  is the computational domain, n is outward normal vector to  .
S is the boundary.Moreover, T and 2 q are the prescribed boundary temperature and heat flux, respectively.
h is the convection heat transfer coefficient and T  is the ambient temperature.
By applying the Fourier transform on Equation (1), the governing equation is transformed into the frequency domain as where T is the Fourier transform of T . is the frequency.When =0  , the problems are transformed into the steady-state heat conduction problem.

The SBFEM for the transient heat conduction problem
As illustrated in Figure 1, the SBFEM presents a local coordinate system ( , )  .The coordinates of a point ( , ) xy along the radial line and inside the domain can be expressed as follows [24]: where , are the scaled boundary coordinates in two-dimensions. is a radial coordinate.
 is the circumferential coordinate.The temperature ( , ) T  at any point in SBFEM coordinates is written as: where () T  is radial temperature functions along a line connecting the scaling center  and a node at the boundary.  By using Equation ( 8) and ( 12), the heat flux ( , ) q  can be expressed as where, Using the weighted residual method and Green theorem and introducing the boundary conditions, the following SBFEM equation can be written as follows [25]  where， The coefficient matrices   that depend only on the geometry and material properties.

Steady-state solution
The SBFEM equation in temperature function () T  is solved for steady-state problems.Setting =0  in Equation ( 16), the SBFEM equations for the temperature field in steady-state heat conduction can be written as By introducing the variable, The equation can be transformed into first-order ordinary differential equations: where the coefficient matrix p Z   is a Hamiltonian matrix.The solution for a bounded domain is obtained using the positive eigenvalues of The solution of Equation ( 23) can be obtained by computing the eigenvalue and eigenvector of the matrix where the real components of eigenvalues  n and  p are negative and positive, respectively.
The general solution of equation ( 23) can be obtained as follows: where 1 {} c and 2 {} c are the integration constants.To obtain a finite solution at the scaling center ( 0 {} c must be equal to zero, the solution in bounded domain can be written as: The relationship between { ( )} where the steady-state stiffness matrix can be expressed as    For the bounded domain, the dynamic-stiffness matrix [ ( )] K  on the boundary =1  formulated in the frequency domain is written as To obtain the mass matrix   M of the bound domain, the low-frequency case in the SBFEM is address, which the dynamic-stiffness [ ( )] K  of a bounded domain is assume as The steady-state stiffness matrix [] st K is computed using Equation (31).Substituting Equation (34) into Equation (33) leads to a constant term independent of i , a term in i and higher order term in i , which are neglected.Besides, the constant term vanishes.
The coefficient matrix of i can be expressed as This is a linear equation to solver the mass matrix [] M . By using the eigenvalues and , the Equation ( 35) can be written as where, where nodal temperature () Tt is the continuous derivative of time, In general, it is not easy to solve the function solution in the time domain.In this paper, the backward difference method is adopted to solve Equation (39).The time-domain is divided into several time units, and the solution of the time node is obtained step by step from the initial conditions, and the node temperature at any time is obtained by interpolation.

Overview framework
This work proposes a framework for the heat conduction analysis using the PSBFEM with polygon/quadtree meshes, as shown in Figure 2. The framework contains three submodules: the pre-processing modules, the heat conduction analysis modules, and the postprocessing modules.In the pre-processing, we develop a python script to generate a polygonal mesh and quadtree mesh.at the first incremental step and store it in the state variable.In the next incremental step, the stiffness/mass matrix is read directly to update RHS and AMATR.

Defining the UEL polygonal elements
The PSBFEM is inherently appropriate for modeling polygons and has other promising capabilities.Because the PSBFEM is discretized only in the boundary, and an element of PSBFEM can assume more complex shapes than a finite element method, the PSBFEM element is much more flexible than the FEM element.With such a significant advantage, the PSBFEM is more applicable to complex geometries than other methods.
Figure 4 shows that the supported element type in the PSBFEM.The PSBFEM also provides the complex element type: polygonal and complex quadrilateral elements (quadtree discretization).Hence, it is an effective tool to solve complex elements in numerical analysis.Listing. 1 The input file of polygon element in the UEL (c.f. Figure 4(a))

Numerical examples
In this section, we carried out several benchmark problems to demonstrate the convergence and accuracy of the proposed framework for the heat conduction analysis.
Moreover, the results of the PSBFEM are compared with the FEM.The FEM analysis uses the commercial finite element software Abaqus.For validation, the relative errors in the temperature are investigated as follows: where T is the numerical solution, and h T is the analytical or reference solution.
6.1 Steady-state heat conduction analysis , and the one at other side is zero.The analytical solution for this problem is given by [26] 100 ( , ) sin sinh sinh(0.5 ) In this example, a is set to 10.0 m.
b is set to 5.0m.Due to the symmetry of this problem, we only select the left part of the domain to analyze.The domain is discretized with the quadrilateral and polygonal elements, as shown in Figure 5 (b) and (c).A convergence study is performed by mesh refinement.The meshes are refined successively following the sequence 1m, 0.5m, 0.25m and 0.125m.In this work, the FEM analysis uses the DC2D4 element.
Figure 6 shows the convergence of the relative error in the temperature with mesh refinement.It is observed that all techniques asymptotically converge to an exact solution with an optimal convergence rate (2:1).Moreover, the PSBFEM show slightly accurate results than FEM at the same element size.The quadtree mesh is generated by setting the same number of mesh seeds on every hole, as shown in Figure 9.It is clearly present that the mesh transition between the holes of different sizes is effectively handled.The square body with multiple holes is modeled with 1632 quadtree elements, and the total nodes are 2478.Moreover, this problem is also analyzed with a similar number of nodes (2532 nodes) using the Abaqus DC2D4 element.
Figure 10 shows that the comparison between the PSBFEM quadtree element and the Abaqus DC2D4 element in the temperature at the right edge.Simple regression is investigated to compare the results of FEM and PSBFEM.Results show that the coefficient of determination 2 R is 0.9997, and the relative error is 0.67%.In addition, the temperature distributions obtained from the PSBFEM and the FEM are shown in Figure 11.It is noted that the contour plots present a good agreement.Hence, these results demonstrate PSBFEM accuracy and reliability for the quadtree mesh.The convergence of the relative error in the temperature with mesh refinement is presented in Figure 13.It is observed that the PSBFEM converges to an analytical solution with an optimal convergence rate.The temperature-time history of point A by the FEM and PSBFEM is compared in Figure 14, where the results obtained by the two methods correspond well with the analytical solution.In addition, the history temperature distributions at different times are presented in Figure 15.It is noted that the contour plots present a good agreement.Therefore, these results demonstrate that the PSBFEM is reliable and accurate for solving the transient heat conduction problems.
To evaluate the time consumption of the PSBFEM, we compare the PSBFEM and Abaqus standard elements using the quadrilateral mesh.Since the increment size setting affects the solution time, we use the auto-increment type.The comparison is evaluated with an Intel Core i7-4710MQ CPU running at 2.50GHz and 4.0 GB of RAM.The total CPU time is normalized.Results presented in Figure 16 show that the time consumption of PSBFEM and Abaqus standard elements are comparable.The quadtree mesh is generated by setting the same number of mesh seeds on the cavity, as shown in Figure 18.It is clearly present that the mesh transition between the cavity of different sizes is effectively handled.The square body with a rabbit shape cavity is modelled with 3147 quadtree elements.Moreover, this problem is also analyzed with a similar number of elements using the Abaqus DC2D4 element, and the total elements are 3358.
It is noted that the PSBFEM obtained solutions are in a very good agreement with the Abaqus for the all points, as shown in Figure 19.Table .2 shows that the relative error of four nodes is less than 1.0%.Moreover, Figure 20

Conclusions
This paper implements a framework for the heat conduction problems using the SBFEM and the polygon/quadtree discretization within the Abaqus/Standard analysis by the UEL subroutine.This work mainly focuses on the detailed implementation of the framework, input data format, and the UEL subroutine, one of the proposed work's critical features.
The implementation of PSBFEM is validated against the FEM by solving a several benchmark problems.The polygon element of PSBFEM has a higher accuracy rate than the standard FEM element.In addition, the results also show that PSBFEM converges to an analytical solution with an optimal convergence rate.The proposed framework is robust accurate for solving the steady-state and transient heat conduction problems.In the future, we will extend the new technique to three-dimensional heat conduction problems.

Figure 1 .
Figure 1.The coordinate system of scaled boundary finite element method.The differential operator can be transformed from the Cartesian coordinate system

2
Mass matrix transient solution To determine the mass matrix   M of SBFEM, introducing the dynamic-stiffness matrix [ ( , )]

3
Transient solution and time discretizationSubstitute Equation (34) into Equation (32) followed by performing inverse Fourier transform, the nodal temperature relationship of a bounded domain is expressed as a standard time domain equation using steady-state stiffness and mass matrices as The UEL subroutine is written by FORTRAN 77.Due to the Abaqus CAE not supporting the UEL element's visualization, we extract the 'vtu' format results from the 'odb' file and visualize it in the software Paraview [38].

Figure 3 .
Figure 3. Flowchart of UEL subroutine for SBFEM To solve the stiffness matrix and mass matrix, we need to employ the eigenvalue decomposition, as shown in Equation (25).At present, many mathematical libraries to perform the eigenvalue decomposition exists.This paper uses the Intel Math Kernel Library (MKL) [36] to decompose the eigenvalue.

Figure 4 .
Figure 4. Schematic diagram of the polygon meshes; (a) random polygon meshes; (b) quadtree meshesThe input file of Abaqus contains a complete description of the numerical model, such as node coordinate, element connectivity, the freedom degree, and material properties.This information needs to be defined in the 'inp' file by the user.We show a simple polygon mesh of PSBFEM that demonstrates the definition of elements in the UEL, as shown in Figure4(a).In the input file (see Listing 1), the mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and pentagon element (U5). 1 ∼ 6 is the line number.Lines 1 ∼ 6 are used to define the pentagonal element (U5); Line 1 assigns the element type, the number of nodes, the number of element properties, and the number of freedom degrees per node; Line 2 sets the active degrees of freedom.Lines 3 ∼ 4 define the element sets' E5'; Line 6 sets the element properties (conductivity, density, specific heat) of 'E5'.

6. 1 . 1
The steady-state heat conduction of square plate A classical steady-state heat conduction problem is considered to investigate the accuracy and convergence of PSBFEM, as shown in Figure5.The temperature on the top side

Figure 5 .Figure 6 .
Figure 6.Convergence of the relative error in the temperature

Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .
Figure 8.The steady-state heat conduction problem of square panel with multiple holes and boundary conditions

Figure 13 .Figure 14 .Figure 15 .Figure 16 .
Figure 14.Comparison of temperature time history of point A presents that the temperature distribution at different time using FEM and PSBFEM.The temperature distribution are virtually the same for the FEM and PSBFEM.Therefore, the PSBFEM with quadtree meshes shows a good effect for solving complex geometric in the transient heat conduction problem.

Figure 17 .Figure 18 .Figure 19 .
Figure 17.The geometry and boundary conditions of square plate with a rabbit shape cavity

Figure 20 .
Figure 20.The temperature distribution at different time using FEM and PSBFEM.

Table . 2
The relative error for the FEM and PSBFEM