Analysis of Unsteady Heat Transfer Problems with Complex Geometries Using Isogeometric Boundary Element Method

: Numerical analysis of unsteady heat transfer problems with complex geometries by the isogeometric boundary element method (IGABEM) is presented. The IGABEM possesses many desirable merits and features, for instance, (a) exactly represented arbitrarily complex geometries, and higher-order continuity due to non-uniform rational B-splines (NURBS) shape functions; (b) using NURBS for both field approximation and geometric description; (c) directly utilizing geometry data from computer-aided design (CAD); and (d) only boundary discretization. The formulation of IGABEM for unsteady heat transfer is derived. The domain discretization in terms of IGABEM for unsteady heat transfer is required as that in traditional BEM. The internal values however are obtained with the analytical formula according to the values on the boundaries, and its computations are therefore mainly dependent on the discretization of the boundaries. The coordinates of internal control points are obtained with the coordinates of control points on the boundaries using Coons body interpolation method. The developed approach is tested with several numerical examples from simple to complicated geometries. Good agreement is gained with reference solutions derived from either analytical or finite element methods.

with more complex geometry and boundary conditions. In that circumstance, advanced numerical methods are preferable. Isogeometric analysis (IGA) uses the NURBS (non-uniform rational B-spline) basis functions in CAD for modeling geometry as shape function in finite element method (FEM) [Hughes, Cottrell and Bazilevs (2005); Shojaee, Valizadeh, Izadpanah et al. (2012)] and references therein. Compared with traditional FEM, IGA owns many desirable features such as exact representation of arbitrarily complex geometries, higherorder continuity and simple mesh refinement. Due to those aforementioned desirable properties, IGA has received extensive attention in recent years, and it gradually becomes a widely used means to solve many engineering problems [Yu, Yin, Bui et al. (2017); Bazilevs, Calo, Zhang et al. (2006); Liu, Yu, Bui et al. (2017); Lai, Yu, Bui et al. (2017); Cottrell, Reali, Bazilevs et al. (2006) ;Nguyen, Bui, Yu et al. (2014); Lorenzis, Temizer, Wriggers et al. (2011);Cottrell, Reali, Bazilevs et al. (2006); Wall, Frenzel and Cyron (2008); Seo, Kim and Youn (2010)]. Garcia et al. [Garcia, Bartoň and Pardo (2017); Garcia, Pardo, Dalcin et al. (2017)] proposed a refined isogeometric analysis (rIGA) with the use of highly continuous finite element spaces interconnected with low continuity hyperplanes to maximize the performance of direct solvers, both the solution time and best approximation errors are simultaneously improved. At present, the isogeometric concept is mainly applied to finite element implementation. However, it can also be applied to other methods, and some other hybrid approaches have been recently developed including the extended isogeometric analysis ; Yu, Bui, Yin et al. (2016)], isogeometric collocation method [Auricchio, Beirão DA Veiga, Hughes et al. (2017); Manni, Reali and Speleers (2015)], isogeometric meshless method [González, Cueto and Doblaré (2009) ;Valizadeh, Bazilevs, Chen et al. (2015)], isogeometric boundary element method [Zhou, Liu, Wang et al. (2017)], scaled boundary isogeometric analysis [Natarajan, Wang, Song et al. (2015); Li, Liu and Lin (2017)]. The geometry built by CAD is its boundary curves/surfaces, and the computation of boundary element method (BEM) is mainly performed on the boundary, which means the combination of IGA and BEM is a natural fit, and actually many researches have paid lots of attention on it. Simpson et al. [Simpson, Bordas, Lian et al. (2013); Simpson, Bordas, Trevelyan et al. (2012)] proposed IGABEM for 2-D elastostatic analysis. Subsequently IGABEM has been applied to many other fields such as 3-D potential [Gu, Zhang and Li (2012)], elastostatic problem [Gu, Zhang, Sheng et al. (2011)], shape optimization [Lian, Kerfriden and Bordas (2016) ;Sun, Yu, Nguyen et al. (2018)], elastoplastic inclusion problems [Beer, Marussig, Zechner et al. (2016); Bee, Mallardo, Ruocco et al. (2017)], gradient elasticity ], crack problems [Nguyen, Tran, Anitescu et al. (2016)], and Helmholtz problems [Peake, Trevelyan and Coates (2013) ;Coox, Atak, Vandepitte et al. (2014)]. Scott et al. [Scott, Simpson, Evans et al. (2013)] coupled IGABEM and T-spline to reduce the number of NURBS patches and improve their smoothness. In order to reduce the computational complexity and computational time, Takahasha and Matsumoto introduced fast multipole method into IGABEM for two-dimnesional Laplace equation [Takahashi and Matsumoto (2012)]. Marussig et al. applied fast IGABEM to elasticity [Marussig, Zechner, Beeret al. (2015)], and Campos et al. [Campos, Albuquerque and Wrobel (2017)] applied the boundary conditions on control points, further improving the applicability of the method. Simpson et al. [Simpson and Liu (2016)] employed a black-box FMM and T-spline to accelerate computation. Beer et al. [Beer, Marussig, Zechner et al. (2014); Wang, Benson and Nagy (2015)] introduced trimmed NURBS technique into IGABEM. Gong et al. [Gong, Dong and Bai (2017)] solved nearly singular integral problems using exponential transformation. The main objective of the present contribution is to study unsteady heat transfer problems by an effective computational method. In general, the numerical implementation of transient heat transfer is more challenging than that of the steady heat problems. Different from the steady heat transfer ], the transient heat transfer, which is studied here, has the following features: (1) since the transient temperature problem is time-dependent, the time domain has to be discretized while considering the spatial domain; (2) the differential equation possesses a time derivative term, dealing with the temperature in domain interior is required, i.e., the domain interior is also discretized, while only the boundaries are discretized for the steady heat problem; (3) an exponential integral function is included in the fundamental solution, and special method is required to deal with the computation of the exponential integral function. Due to the inherent features of NURBS, the boundary and domain of geometry can be represented exactly by data in CAD. This method retains the advantage of mesh refinement on parameter space in IGA. It is important to point out that, for the steady heat transfer using the IGABEM, the discretization is performed only on the boundary, whereas the discretization of domain is still required for unsteady heat transfer as that in traditional BEM. In the IGABEM, the internal values are obtained with the analytical formula according to the values on the boundaries, and the computations mainly focus on the discretized calculations on the boundaries. The coordinates of internal control points are obtained with the coordinates of control points on the boundary using Coons body interpolation method. The body of this paper is structured as follows. Section 2 briefly presents transient heat transfer problem. In Section 3, the IGABEM formulation for solving transient heat conduction problem is derived in detail. Numerical examples are presented and discussed in Section 4, in which the obtained results are compared with analytical solutions or other methods such as FEM. Some conclusions are drawn in Section 5.

Formulation of transient heat transfer analysis
According to the theory of heat transfer, the differential equation of unsteady heat transfer problem in 2-D isotropic solid is given by Yu et al. [Yu, Yao and Gao (2014) where Ω is the domain of geometry, T is the transient temperature, λ the thermal conductivity coefficient, w the heat source, ρ the material density, c the specific heat capacity, t the time, 2 ∇ Laplace operator.
The initial condition and the boundary conditions are expressed as follows: where 0 T is the initial temperature in the domain; T is prescribed temperature, 2 q is prescribed heat flux, while n is the normal vector pointing outward of boundary.
Using Eqs. (1) and (2), we can obtain the following weighted residual equation: where * T is an arbitrary function, 0 t and f t are two time values.
Integration by parts about Laplace operator item twice, and time derivative item once, and applying the Green formula and boundary conditions, Eq.
(3) will be written as [Wu (2008)] In this paper, the heat source will not be considered for simplicity, i.e., 0 w = , so Eq. (4) can be rewritten as

Brief on NURBS basis functions
The knot vector ( ) { } 1 + +1 = =0,..., ,..., 1 T i n p ξ ξ ξ ξ = k is defined as a set of non-decreasing numbers that are between zero and one. Here, i is the knot index, i ξ is the i th knot, n is the number of basis functions, and p is the polynomial order. The NURBS basis function , ( ) i p R ξ is a weighted average of the B-spline basis functions, is defined as [Hughes, Cottrell and Bazilevs (2005) N ξ is the i th B-spline basis function of degree p, which is defined recursively as [Hughes, Cottrell and Bazilevs (2005) The two-dimensional NURBS basis functions can be constructed by taking the tensor product of two one-dimensional B-spline basis functions as [Hughes, Cottrell and Bazilevs (2005) The definition of ( ) η k is the same as that of ( ) ξ k . By using the NURBS basis functions, a NURBS curve of order p can be constructed as where i P are the coordinates of control point i.

IGABEM for unsteady temperature field
where ′ x and f t are regarded as the source point of space and time respectively, x and t are field point of space and time, r ′ = x -x is the distance between source point and field point. (12), we can find that we can easily get the temperature at points ′∈ Ω where the term of domain integral represents the influence of initial condition and ( ) C ′ x is the same as that in steady heat transfer problem ; Wu (2008)], which is related to the shape of boundary.
( ) The boundaries are discretized into n non-overlapping elements The boundary variables are interpolated with shape functions. Time domain is uniformly divided into f steps. In each step T ∆ , T and q are considered as be constant, in time (14) can be rewritten as Integrating * T and * q over the time domain, we can get [Wu (2008) is an exponential integral function, which can be calculated by the series method, shown as follows: where c is the Euler constant.
where m is the number of nodes in each element, l is node number, where ξ ′ is local coordinate of source point, e′ the element where source point located in, and the Jacobian of transformation ( ) e J ξ is given by where 1 ξ and 2 ξ are the coordinate of the start and end point of element in parameter space respectively.
In terms of the IGABEM, the nodal points (i.e., control points) may not be situated on the boundary, the collocation points are hence defined as [Johnson (2005) Applying the source point to the discrete nodes, Eq. (24) can be written in matrix form as follows: there is no need to match the meshes in the domain with the boundary meshes.

Numerical integration
When the source point is situated in the integral element, strongly singular integral in i H and weakly singular integral in i G will exist. In general, the singularity subtraction technique (SST) [Guiggiani and Casalini (1987); Yao and Wang (2010)] is used to evaluate strongly singular integral and the transformation approach proposed [Telles (1987)] is used to evaluate weakly singular integral. There is a domain integral term in boundary integral equation Eq. (24) It is worth noting that if the boundary condition is a constant, the boundary condition can be directly added to the control points. If the boundary condition is a function distribution, this method cannot work due to its particularities. In IGA, Lagrange multiplier method, penalty method, and Nitsche's method [Nguyen, Kerfriden, Brino et al. (2014); Gu, Yu, Lich et al. (2018)] can be adopted to solve this problem, and the L2 projection method and collocation method are used in IGABEM [Lian, Kerfriden and Bordas (2016)]. Actually, in IGABEM, for the boundary with known boundary condition, we can immediately integrate it by boundary integral equation and only discrete the unknown value of this boundary, thereby avoiding the error caused by applying boundary condition.
In the calculation of transient temperature field, BEM does not adopt finite difference approximation for time domain discretization, the integral can be obtained precisely with analytical method by considering the fundamental solution including the time effect in the equation.

Numerical implementation
Using the above theory and techniques, we summarize the entire process of calculation and show it as follows: (1) Read CAD input data including control points, knot vector of both boundary and domain.
(2) Read material parameter, and determine time step.
(3) Create boundary elements and domain elements.  It should be noticed in this example that, since the temperatures at the corners of square are known, but the values of heat conduction are unknown, special corner treatment is required, which is reported in Walker et al. [Walker and Fenner (1989)], and see the Appendix B. The knot vector, weights and collocation points, elements, control points of boundary, mesh and control points of domain, all of them are the same as those in the above example. The time step t ∆ =0.001 s is used. The distribution of temperature in the domain obtained by the developed IGABEM is compared with that derived from FEM (ANSYS) with 15×15 quadratic elements. As expected, a good agreement between both solutions is obtained.
In the present study, the time step is constant, too large time step may induce large error, while too small time step will lead to inaccurate calculation of exponential integral function. The suitable time step is determined according to the numerical test. When the boundary mesh gets finer, larger time step may be used.   { } = 0, 0, 0, 1/ 9, 1/ 9, 1/ 3, 1/ 3, 2 / 3, 2 / 3, 7 / 9, 7 / 9, 1, 1, 1 k , the weights are { } 1 1,1,1, ,1,1,1,1,1,1,1 2 = w . Two points A and B are taken into account for comparison with the results obtained by ANSYS FEM. Fig. 10 sketches the collocation points, control points and NURBS basis function of initial geometry. The mesh and control points on the boundary and in the domain are shown in Fig. 11. The mesh of 173 quadratic elements in the FEM analysis is depicted in Fig. 12. The computed results of t ∆ = 5 s are then shown in Fig. 13. In addition, the distribution of temperature at time 200s derived from both the IGABEM and FEM are plotted in Fig. 14, respectively. A good agreement is obtained between the IGABEM result and the FEM result. As shown in Fig. 13, it is interesting to see that a similar variation of the temperature in the plate at two points A and B is obtained. The developed IGABEM, once again, offers good solutions and all the computed results agree well with the reference FEM solutions. The temperature at one point in the domain varies from the initial temperature (50 o C) to the specified final temperature (200 o C) with increasing the time. Because the distance is different between one point and the hole boundary, the temperature at one point is hence different at one time. The temperature in the area closer to the face (boundary) with the highest temperature will approach faster to the maximum temperature than those areas that are far from the particular boundary, e.g., the location of Point A to the hole boundary is closer than that of Point B, so the temperature at Point A is bigger than that of Point B until the same temperature is achieved at last, as shown in Fig. 13.

Heat convection in a complex domain
In the last example, the problem of multi-boundary with complicated geometry is considered. The geometry and boundary condition of considered structure are shown in Fig. 15, in which the unmarked boundaries are regarded as adiabatic boundaries. The properties of material are as follows: thermal conductivity is 391 W m C  , density is 8940 3 kg m , and specific heat is 385.2 J kg C  . The outer boundary and the inner boundary are built respectively, and it is noted that the outer boundary is constructed in a counterclockwise way, different from the inner boundary which is built in a clockwise way. Fig. 16 represents the control points, collocation points and NURBS basis function of the initial geometric boundaries. The Gaussian integration points for domain integral are extracted by IGA, and the initial elements and control points of domain are shown in Fig. 17. In order to enhance the accuracy of the solution, or a better result, the original elements are refined into 10 new elements, which are shown in Fig. 18. In this problem, the time step is t ∆ =500 s. Also, Fig. 19 shows the mesh of FEM (1228 quadratic elements). The

Conclusion and outlook
In this paper, two-dimensional unsteady heat transfer problems have been solved with the IGABEM, in which the NURBS basis functions are used to approximate the geometry and temperature fields. Arbitrarily complex geometries can be exactly represented and higher-order continuity is obtained. In IGABEM, the discretization of domain is required for unsteady heat transfer problem. This is different from the steady heat transfer problem, in which the discretization of domain is avoided. However, the internal values are obtained with the analytical formula according to the values on the boundaries, so the computations mainly focus on the discretized calculations on the boundaries. Through numerical examples, the obtained results and validated results with reference solutions show good performance and high accuracy of the proposed method. One interesting problem that left behind this contribution would be prefer to 3-D cases, which however has been scheduled for our future works. NURBS-based IGA only works well for quadrilateral domains. For the complex domain problems, there are two ways to deal with them to some extent, i.e., subdividing the design domain of complex topology into multiple quadrilateral patches [Manh, Evgrafov, Gersborg et al. (2011); Qian and Sigmund (2011)] and using trimmed surfaces [Seo, Kim, and Youn (2010)

Appendix B
For geometry with corners, we put two overlap points with the same temperature but different heat conductions of both sides on the corner. If the temperature at corner is known, but the heat conductions on both sides of the corner are unknown, we can first calculate the heat conduction on one side, and then calculate the other side. The formula for calculating the heat conduction at corner is as follows [Walker and Fenner (1989)