Isogeometric Boundary Element Analysis for 2D Transient Heat Conduction Problem with Radial Integration Method

Leilei Chen,Kunpeng Li,Xuan Peng,Haojie Lian,Xiao Lin and Zhuojia Fu

1College of Architecture and Civil Engineering,Xinyang Normal University,Xinyang,464000,China

2Department of Mechanical Engineering,Suzhou University of Science and Technology,Suzhou,215009,China

3Key Laboratory of In-Situ Property-Improving Mining of Ministry of Education,Taiyuan University of Technology,Taiyuan,030024,China

4Institute for Computational Engineering,Faculty of Science,Technology and Communication,University of Luxembourg,Luxembourg

5The York Management School,University of York,York,YO10 5GD,UK

6College of Mechanics and Materials,Hohai University,Nanjing,211100,China

ABSTRACT This paper presents an isogeometric boundary element method(IGABEM)for transient heat conduction analysis.The Non-Uniform Rational B-spline(NURBS)basis functions,which are used to construct the geometry of the structures,are employed to discretize the physical unknowns in the boundary integral formulations of the governing equations.B´ezier extraction technique is employed to accelerate the evaluation of NURBS basis functions.We adopt a radial integration method to address the additional domain integrals.The numerical examples demonstrate the advantage of IGABEM in dimension reduction and the seamless connection between CAD and numerical analysis.

KEYWORDS Isogeometric analysis; NURBS; boundary element method; heat conduction; radial integration method

1 Introduction

Isogeometric analysis (IGA) has drawn extensive attentions in engineering science community since the seminal works of Hughes et al.[1].As an alternative to the traditional analysis based on Lagrange polynomial, IGA uses the same spline functions that are used in the geometric expression in CAD, for example Non-Uniform Rational B-splines (NURBS), as the basis functions to approximate the unknown field.Hence, IGA is able to perform numerical analysis directly from CAD without cumbersome preprocessing procedure and geometric errors.Additionally, IGA offers high-order continuity and flexible refinement schemes which are particularly attractive in numerical simulation.Due to the aforementioned salient features, IGA has been successfully applied in many fields [2-6].To facilitate the implementation of IGA, Bézier extraction [7] technique is proposed that enables one to develop IGA codes in existing finite element codes, without changing the program framework.By using the extraction operator, the complicated iterative process for evaluating B-spline basis functions is eliminated, which improves the efficiency of simulation significantly.

Although the concept of IGA was originally proposed in the context of finite element methods, its application in engineering practice is restricted because FEM relies on volume parameterization which conflicts with the boundary representation in CAD field.Mierzwiczak et al.[8]and Fu et al.[9] applied the isogeometric boundary element method to the heat conduction problem, and proposed the origin intensity factor to eliminate the singularity of the integral.The order reduction of POD model for transient heat transfer was studied by Li et al.[10].Wang et al.[11] solved the three-dimensional elasto-plastic problem by boundary element method.On the contrary, boundary element method (BEM) [12-19] only needs boundary parameterization and thus is compatible with CAD models naturally.The isogeometric analysis with boundary element method (IGABEM) was first applied to potential problems [20] and elasticity analysis [21,22].Since then, IGABEM has been applied successfully to a wide range of fields, such as heat conduction [23], linear elasticity [24,25], fracture mechanics [26-28], electromagnetics [29],acoustics [30-35], structural optimization [34,36-40], etc.

In this paper, we are aimed to extend the application of IGABEM to heat conduction problems which are an increasingly important topic in many engineering fields.The radial integration method proposed by Gao [41] is incorporated to solve the additional domain integrals.The problems of constant and variable coefficients, transient and steady state of temperature field are studied.The rest of the paper is arranged as follows:Section 2 gives an overview of the B-spline and NURBS basis.Bézier extraction is detailed in Section 3.In Section 4, the boundary integral equation for heat conduction is introduced and the discrete formula derivation of the boundary integral equation, as well as the time marching method are introduced in detail.Section 5 provides two sets of numerical examples to verify the algorithm, followed by the conclusion in Section 6.

2 NURBS

For the sake of completeness, the fundamentals of NURBS are briefed in this section.NURBS are generalized from B-splines, whose basis functions are defined over a knot vector which is a monotonic increasing real number sequence, denoted byU={u1,u2,...,um} whereuiis thei-th knot in the vectorU.B-spline basis function is evaluated iteratively as:

wherepis the order of the polynomial in B-spline basis functions, andnis the number of the basis functions.

The derivative expression of the basis function can be derived as follows:

The B-spline curve is described by the linear combination of NURBS basis functions and the corresponding control points,

where Piis the Cartesian coordinates of control points.

Non-uniform rational B-spline (NURBS) basis functions are obtained by introducing the weight coefficientωto rationalize the B-spline basis,

where

Similar to B-splines, the NURBS curve is expressed by:

Since the NURBS basis function adopts the rational form composed of the B-spline basis function and the weight factorω, it is able to exactly express conic curves, such as ellipses and circles.

3 Bézier Extraction Operation

The computation of IGABEM can be accelerated using Bézier extraction.The idea is to extract linear operators and map the Bernstein polynomial basis functions of Bézier elements to global B-spline basis functions.The linear transformation is defined by a matrix called an extraction operator.The transpose of the extraction operator maps the control points of the global B-spline to that of Bézier curves.The process of Bézier extraction operation is proposed in this work, also see [7].

To perform Bézier decomposition, all internal knots of a knot vector are repeated until their multiplicity equalsp.For example, we insertuafter the k-th knot of the original knot vectorΞ={u1,u2,...,un+p+1}.The new knotξ={u1,u2,...,uk,u,uk+1,un+p+1}withk ＞p.The introduction of new knots will bring about changes in the basis functions and the generation of new control points.The relationship between the original control point P and the new control point Q can be deduced:

wherem=n+1.The expression for calculating factorαiin Eq.(8) is as follows:

After the knot is inserted, it brings about a change in the control point.Let {u1,u2,...,um}be the set of insertion knots required for the Bézier decomposition of the B-spline.Then for each new knot,,j=1,2,...,m, using Eq.(9), we define,i=1,2,...,n+j, to be theithα.So the expression for the conversion matrix C is as follows:

Then, the corresponding control points can be calculated by:

We rewrite Eq.(11) with matrix form:

To decompose a B-spline curve defined in the knot vector {0,0,0,0.25,0.5,0.75,1,1,1}with Bézier extraction, we will insert the knots {0.25,0.25,0.5,0.5,0.75,0.75} into the existing knot vector.Without changing the geometry, the B-spline curve is expressed using the Bézier basis function as follows:

Using Eqs.(12), (13) can be expressed as follows:

Since P is arbitrary, we can get the following equation

It should be noted that only the knot vector determines the C matrix.In other words, the extraction operator is a product of parameterization and does not depend on control points or basis functions.Therefore, we can apply the extraction operator directly to NURBS.Define weights as diagonal matrices

Converting Eq.(5) to matrix form, we have

Using Eq.(15), the NURBS basis function can be expressed in terms of the Bernstein basis as

The NURBS curve is expressed as:

In Eq.(19), B(u)can be used to express the NURBS curve.Letωb=CTω, and define Wbto be the diagonal matrix

Eq.(19) can be written as follows:

The expression of the Bézier control points can refer to the method of homogeneous coordinates:

Both sides are multiplied by Wbat the same time:

After solving the weight factorωband the control pointPbafter the Bézier operation, we arrive at the final Bézier representation of the NURBS:

Thus, we can express NURBS equivalently with a series ofC0Bézier elements.

4 Governing Equation

The governing partial differential equation of transient heat conduction in isotropic medium with internal heat source is expressed by:

whereT(x,t) represents the temperature at point x at timet,kis the thermal conductivity,b(x,t) shows the heat source value at point x at timet,ρis the density of the material, andcxis specific heat capacity.

The following three kinds of boundary conditions are considered:

whereγTrepresents the dirichlet boundary,γqthe neumann boundary andγfconvective boundary condition.n(x) is the unit external normal direction vector at the boundary point.

The transient problem is time-dependent, so the initial conditions need to be prescribed:

whereT0indicates initial temperature on the domain.

4.1 Boundary Element Method of Transient Temperature Field

Through weighted residual operation and partial integration operation, Eq.(25) can be expressed as the following boundary-domain integral equation:

where the coefficientc(x)=1 when the source point x is located in the domain,c(x)=0.5 when the source point x is on the boundary of structure, and the integrals are

In this work, we do not consider the source.Thus, Eq.(27) can be rewritten as

By observing Eq.(28), we can find that the last term on the right hand of Eq.(30) is domain integral term rather than boundary integral term.In this work, the radial integral method is used to transform domain integral into boundary integral to retain the advantage of dimension reduction of IGABEM.

We firstly approximate unknown temperature gradientby interpolation calculation.The commonly used interpolation function is radial basis function (RBF).We select a number of collocation points and interpolation to get the approximate formula of variables,

whereNis the number of collocation points,βiis the undetermined coefficient at thei-th collocation point,fi(r)is the radial basis function at thei-th collocation point,r=|y-xi|, and xiis the coordinates of thei-th collocation point.

It is worth noting that collocation points can be boundary points or internal points.Generally,combining boundary points and internal points to form collocation points can effectively improve the approximation accuracy.In addition, the accuracy of approximation can be further improved by combining radial basis functions and polynomials, which can be expressed in combination with quadratic polynomials,

It is particularly important to determine these coefficientsβandγin Eq.(32).Although it is not easy to directly express these coefficients, we can get the implicit expression of these coefficients through a series of transformation operations.By setting collocation point asin Eq.(32) and assembling them into matrix-vector formulation, we can get the following expression:

where

and

When the collocation point is not repeated, the matrix f is invertible.Thus, the coefficient vector v can be expressed as:

By substituting Eq.(33) into the last term on the right of Eq.(27) and using the radial integration method, the following boundary integral formula can be obtained

where

and

It can be found that the domain integral in Eq.(39) is successfully transformed into the boundary integral.For the detailed derivation process of the above results, also see [42].

It is worth noting that these coefficients in Eq.(40) are unknown, so it is impossible to obtain the results ofI0,I1, andI2directly through the boundary integral calculation.In fact, these coefficients do not need to be calculated directly.

When the source point x is set as thei-th collocation point, the termsI0,I1, andI2can be rewritten as a formulas of vector and vector multiplication, respectively.By combining the three terms, we can get the following formulas

where the vector Aicontains the result of the boundary integrals in Eq.(40).By substituting Eq.(38) into Eq.(42), we can obtain the following formulation

where Ci=Aif-1.When the point x in Eq.(43) is chosen as collocation point, by assembling the equations at all points, we could get a formulation of matrix vector multiplication.

4.2 Discretization of Boundary Integral Equation

First, we discretize the boundary intoNenon-overlapping NURBS elements,

NURBS basis functions are used for approximation.In any NURBS elemente, variables can be approximately expressed as follows:

wheremis the number of control points,ξ∈[-1,1] represents local parametric coordinates, ˆTlandqlare normalized temperature and its flux at thel-th collocation points.In this work, the Greville abscissa is used to define the positions of collocation points in parameter space, also see [39].

By substituting Eq.(45) into Eq.(27), the two boundary integral terms can be rewritten as

whereJe=dΓ/dξdenotes Jacobian.Since the last term in Eq.(30) is calculated by radial basis function approximation, only geometric interpolation is needed in the final boundary integrals.Therefore, we only need to replacedΓin Eq.(40) withJedξ.On the other hand, the integrals in the above equations are solved with Gauss integral method.In this work, the number of Gauss quadrature points is set as eight.

After collecting the equations for all collocation points and expressing them in matrix forms,one can obtain the following system of linear algebraic equations:

In this work, the finite difference method is used for solution of the above time domain equations.

It is worth noting that weakly singular integrals exist in Eq.(30).For BEM, it is very important to deal with this kind of singular integral effectively, because it directly determines the accuracy of BEM.Qu et al.[43] and Gong et al.[44] used exponential transformation to calculate accurately the weak singular integrals.Gao et al.[45] used the subtraction of singularity technique to overcome the singularity problem.In this work, the subtraction of singularity technique is also used for calculation of weak singular integrals, see [39].

4.3 Time Marching Method for Solving Transient Heat Conduction Problems

Eq.(27) contains a set of integral equations about time.First, we assume that the total time period is fromt0totm.And then, divide it intomintervals.Thus, theith moment can be expressed as:

Herein,Tidenotes the temperature atti, andTi+1denotes the temperature atti+1.When the time marching method is used to solve the equations, we can obtain the derivative of temperature with respect to time, as follows:

where the time step Δtis equal toAfter interpolation approximation, we can get the temperatureand its derivativeqatt∈[ti,ti+1], as follows:

whereβchanges from 0 to 1.Whenβ=0, it means forward difference scheme is used, whenβ=1, it means backward difference scheme is used, otherwise it means central difference scheme is used.By substituting Eq.(50) into Eq.(30), we can obtain the following formulations

Thus, Eq.(47) can be rewritten as

After merging the similar terms, we can get the following expression:

where

When solving the unknown temperatureor its derivative qi+1atti+1, the temperatureand its derivative qiattiare all known.Moving all the unknown terms of Eq.(53) to the lefthand side and all the known terms to the right-hand side by considering the boundary conditions,one finally obtain the following system of linear equations:

where B is the coefficient matrix, xi+1is the vector of unknown variable values at the collocation points, yi+1is the known vector.All the unknown state values can be obtained after Eq.(55)is solved.

5 Numerical Examples

5.1 Square Example

Consider a two-dimensional flat plate, the initial temperature is 100 K, the upper and lower boundaries are adiabatic, the left boundary temperature is maintained at 100 K.Convective boundary condition is applied to the right boundary, and the ambient temperature is 400 K, as shown in the Fig.1.Material densityρ=271 kg/m3, specific heat capacitycx=871 J/(kg·K),thermal conductivity k=202.4 W/(m·K), convective heat transfer coefficient on the right boundaryh=80 W/(m2·K) and thermal emissivityε=1.The boundary element model is as shown in the Fig.1b.Each boundary is divided into 16 equally spaced linear elements, with a total of 64 boundary elements, 64 boundary nodes and 15 internal points.

Figure 1:Plate model diagram and boundary conditions.(a) Plate boundary conditions.(b) Plate model diagram

Fig.2 shows the distribution of temperature along the lower boundary of plate whent=100, 200, 400 and 600 s.In this figure, the calculation results obtained by using IGABEM are compared with those of the finite volume method (Fluent), and it demonstrates the correctness and validity of this algorithm proposed in this work.

5.2 An Example of Elliptic Plate

Consider an elliptical plate with a long radius of 2 m and a short radius of 1 m.The initial temperature is 100 K, and the ambient temperature of the plate boundary is 400 K.Material densityρ=271 kg/m3, specific heat capacitycx=871 J/(kg·K), thermal conductivityk=202.4 W/(m·K), and boundary convective heat transfer coefficienth=80 W/(m2·K).

The shape of ellipse and the distribution of internal points are shown in Fig.3.Fig.4a shows the initial NURBS curve and the position of control points.After normalizing the knot vector, the parameter space vector of boundary element can be obtained, which is expressed as=[0,0.25,0.5,0.75,1].Therefore, four boundary elements are formed, and the parameter space intervals are [0, 0.25);[0.25,0.5);[0.5, 0.75);[0.75,1) respectively.By splitting each NURBS elements equally and inserting new knots, the new control point sequence, knot vector and weight factor vector after refinement can be obtained.For example, insert a new control point on the middle point of each initial NURBS cell parameter space to get the refined control point sequence shown in Fig.4b, as shown by the black circle point “inserted knots.” After Bézier extraction operation,a new set of control point sequence is obtained, as shown in “BE operation knots.” For “initial knots” in Fig.4a and “inserted knots” in Fig.4b, we can find that inserting a new point will also change the position of the original point, but still describe the same curve.However, Bézier extraction operation does not change the position of the original control points, but inserts a new control point at the middle point of some elements.Then, the new control point sequence and Bézier extraction sequence when the initial NURBS cells are divided into 3, 4, 5 and 6 subunits are given in turn, as shown in Figs.4c-4f respectively.

Figure 2:Temperature distribution along the lower boundary of the plate at different times

Figure 3:The NURBS curve and internal points for elliptic

Figure 4:Geometric control points and NURBS curves of elliptical model with different refinement levels.(a) Initial.(b) Refine2.(c) Refine3.(d) Refine4.(e) Refine5.(f) Refine6

In this work, the finite difference method is used to solve the time-domain equation.Therefore,we first investigate the influence of time step on the calculation results, as shown in Tab.1.The temperature values at nine points on thexaxis are calculated, where different time step is used, for example Δt=5, 6, 7, 8, 9, and 10.From this table, we can find that the results with different time steps have small deviation.In fact, too large step value may lead to large error of temperature gradient value, even the calculation is not correct at all.If the time step is too small, it will also lead to inaccurate calculation of temperature gradient value.Therefore, it is very important to select an appropriate time step value.It can be seen from Tabs.1 and 2 that it is appropriate to set the time step from 5 to 10.

Table 1:Temperature values at several internal points with different time step values in 200 s

In this work, we use the radial basis method to approximate the temperature gradient.Therefore, it is significant to investigate the influence of different types of radial basis functions on the numerical results.Herein, eight different types of radial basis functions are given, as shown in Tab.3.Tab.4 shows the comparison of temperature values at several calculation points when using different radial basis functions, where the time step is 5.By observing Tab.4, it can be found that the deviation of calculation results in 200 s is very small when using different types of radial basis functions, more results in 400 s are also found in Tab.5.It verifies the stability of the algorithm proposed.

Table 2:Temperature values at several internal points with different time step values in 400 s

Table 3:Types of radial basis functions

Table 4:Comparison of calculation results with different types of radial basis functions in 200 s

Table 5:Comparison of calculation results with different types of radial basis functions in 400 s

Figs.5a and 5b respectively show the comparison of temperature values at points on thexaxis andyaxis based on traditional BEM and IGABEM algorithm developed in this paper.From these two graphs, it can be found that the results of the IGABEM algorithm are consistent with the results of conventional BEM, which verifies the correctness of the algorithm developed in this work.

Figure 5:Temperature distribution on X-axis and Y-axis, respectively.(a) Temperature of X-axis symmetry point.(b) Temperature of positive half axis of Y-axis

5.3 An Example of Octagonal Leaf

In this section, we consider an example of octagonal leaf, as shown in Fig.6.The initial temperature is 100 K, and the ambient temperature of the structural boundary is 600 K.Other physical parameters are consistent with the previous example.The new control point sequence and Bézier extraction sequence when the initial NURBS cells are divided into 2, 3, 4, 5 and 6 subunits are given in turn, as shown in Figs.7b-7f respectively.

Figure 6:The NURBS curve and internal points for octagonal model

Figure 7:Generation of NURBS curve control points for octagon model.(a) Initial.(b) Refine2.(c) Refine3.(d) Refine4.(e) Refine5.(f) Refine6

Similar to the previous example, we investigate the influence of different radial basis functions on the calculation results, as shown in Tab.6 for 200 s and in Tab.7 for 400 s.It shows that there is a small deviation in the calculation results when different radial basis functions are used,which verifies the stability of the algorithm in this work.

The influence of time step parameters on the calculation accuracy can be found in Tab.8 for 200 s and in Tab.9 for 400 s.Similar to the previous results, the change of time step parameters from 5 to 10 has little effect on the calculation results, which verifies the stability of the algorithm in this work.

Table 6:Comparison of calculation results with different types of radial basis functions in 200 s for octagonal leaf

Table 7:Comparison of calculation results with different types of radial basis functions in 400 s for octagonal leaf

Table 8:Temperature at internal points with different time step values in 200 s for octagonal leaf model

Table 9:Temperature at internal points with different time step values in 400 s for octagonal leaf model

6 Conclusion

This paper applied IGABEM to the heat conduction problem of temperature field.IGABEM can accurately represent the geometric model.In addition to the advantages of dimensionality reduction, IGABEM can also realize the seamless integration of CAD modeling and numerical analysis.The use of Bézier extraction technique further simplifies the implementation of IGA.This method transformed the basis functions of NURBS into Bernstein polynomials and maintained the consistency of geometric models.The main feature of this method was that the iterative calculation of NURBS basis function was avoided in the process of element physical interpolation calculation of BEM, which can effectively improve the computational efficiency.The domain integral term was successfully converted to the boundary integral in IGABEM with radial integration method.The present method provides a powerful tool for fast and high fidelity simulation of transient heat conduction problems commonly encountered in numerous industrial sectors.

Funding Statement:This research was funded by National Natural Science Foundation of China(NSFC) under Grant Nos.11702238, 51904202, and 11902212, and Nanhu Scholars Program for Young Scholars of XYNU.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.