Coupling Heat Conduction and Radiation by an Isogeometric Boundary Element Method in 2-D Structures

We propose an ecient isogeometric boundary element method to address the coupling of heat conduction and radiation in homogeneous or inhomogeneous materials.e isogeometric boundary element method is used to construct irregular 2Dmodels, which eliminate errors in model construction. e physical unknowns in the governing equations for heat conduction and radiation are discretized using an interpolation approximation, and the integral equations are nally solved by Newton–Raphson iteration; it is noteworthy that we use the radial integration method to convert the domain integrals to boundary integrals, and we combine the numerical schemes for heat conduction and radiation.e results of the three numerical cases show that the adopted algorithm can improve the computational accuracy and eciency.


Introduction
Since Hughes et al. [1] proposed the isogeometric boundary element method, which has been developed like a mushroom, and its principle is using the B-spline interpolation function instead of the traditional Lagrangian interpolation basis function. e traditional nite element and boundary element techniques use the Lagrangian interpolation basis function; therefore, it is necessary to convert the geometric model created by computer aided design (CAD) into a mesh rst, and then we calculate the mesh; such an operation increases the error of mesh reconstruction and waste of time; the isogeometric boundary element method inherits the advantages of dimensionality reduction and deals with in nite domains in the traditional eld. Simultaneously, the CAD model can be directly imported into CAE for calculation, simplifying the process of mesh reconstruction and reducing the error of the geometric model caused by tedious preprocessing. In addition, the IGABEM can provide highorder continuity and exible re nement schemes, and it can be extended to many 2D modeling elds. erefore, the IGABEM has been widely studied by many scholars and has been successfully applied in many elds such as elasticity [2,3], uid-structure coupling [4][5][6], shape optimization [7,8], and acoustic [9][10][11][12][13][14][15][16]. e heat conduction and radiation require coupling the two physical elds, which can be signi cantly deviated if only heat conduction or radiation is considered. Many scholars have studied the coupled heat conduction and radiation at home and abroad. Furmanski and Banaszek [17] proposed a method based on the combination of nite element space discretization and iterative skills and used it to study the coupled heat conduction and radiation in a rectangular region.
Kong and Viskanta [18] learned the associated heat transfer and radiation in two-dimensional rectangular glass media by using the discrete coordinate method. Lacroi et al. [19] studied the radiation coupling heat transfer of a rectangular translucent medium irradiated in a speci c direction. Zhou et al. [20] combined the nite volume method with the spectral band model to solve the heat transfer in the absorbing and scattering nongray medium. Luo et al. [21] solved the two-dimensional coupled heat transfer and radiation in the isotropic rectangular medium by the ray trace node analysis method. Mondal and Mishra [22] combined the Boltzmann method with the finite volume method to deal with the coupled heat transfer and radiation with heat flow and temperature boundary conditions. Gu et al. [23][24][25][26][27][28][29][30] studied the heat transfer of coating structure based on the isogeometric boundary element method, and they also used the isogeometric boundary element method by calculating the effective property of steady-state thermal conduction in 2D heterogeneities with a homogeneous interphase. Fu et al. [31][32][33] studied the boundary collocation method for anomalous heat conduction analysis in functionally graded materials. Nie et al. [34][35][36][37][38] studied the noniterative inversion of heat flow boundary conditions and thermal stress estimation of gradient materials based on the refined integration finite element method. Chen et al. [39][40][41][42][43] studied the heat conduction analysis of two-dimensional and threedimensional geometric boundary element methods. e study in this paper is the first to apply the isogeometric boundary element method to coupled heat conduction and heat radiation. Compared with the traditional boundary element method, this method not only has unparalleled advantages in model construction but also is closer to the software solution in terms of computational accuracy. e remainder of this paper is organized as follows: e theory of the isogeometric analysis is given in Section 2. e coupled heat conduction and radiation are put in Section 3. en, Section 4 gives the discretization of the integral equation, and Section 5 gives how to discrete equations and form the system equations. Section 6 gives some numerical examples, and Section 7 gives details of the summary.

Isogeometric Analysis
In this section, we focus our attention on methods of modeling construction by using the isogeometric boundary element method, which is a variant of the BEM based on the use of NURBS basis functions. Due to the precise representation of geometry, the isogeometric boundary element method can be applied to the field of heat conduction and radiation. erefore, a brief introduction of B-splines and NURBS is given in this part. For the interested reader, a complete description can be found in the study conducted by Piegl and Tiller [44].
A B-spline is defined using a group of piecewise polynomials, which are determined by the following "knot vector" of nondecreasing values Where a is the knot index, p is the curve degree, and n is the number of basis functions or control points P a . Each ξ a ∈ Ξ is called a knot. e B-spline basis functions of degree p are defined recursively; starting with p � 0, we define: us, as is given in [44], a whole B-spline curve can be defined by n basis functions (in Equations (1) and (2)) and control points P a . However, in the implementation of BEM, a discretized form of the integral boundary equation is commonly used, in which computations are focused on a.
A single knot span corresponds to a single element in conventional BEM implementations. erefore, in this work, the descriptions of boundary geometry or physical quantities are given in a knot span with p + 1 nonzero basis functions, and the values of which can be obtained by Equations (1) and (2). A pth degree B-spline curve in a knot span (with p + 1 non-zero basis functions) is constructed by mapping from the parameter space to physical space, as shown in the following equation: where P a denotes the set of control point coordinates and x � (x, y, z) is the location of the physical curve corresponding to the spatial coordinate ξ in parametric space. NURBS, a dominant tool used to describe curves and surfaces in CAD systems, is developed from B-splines and can offer significant advantages due to their ability to describe circular arcs and other conic sections exactly. e NURBS geometry is a weighted form of the B-spline definition, i.e., where R a,p are the NURBS basis functions, which are defined by where a denotes the control point index, N a,p is the B-spline basis function from equation (2), and w a is a weight associated with control point P a when all the weights are of equal value. e NURBS curve will degenerate to a B-spline curve. At this point, it is useful to emphasize that the control points do not all lie on the spline curve.

Basic Governing Equations of Heat
Conduction and Radiation e solution to the heat conduction and radiation problem needs to consider two basic equations simultaneously. Due to the complexity of radiation heat transfer, the main difficulty of solving this problem is transferred to the thermal radiation source term. We first consider the steady-state linear heat conduction and radiation problem. is paper expounds on using the isogeometric boundary element method to solve this problem. 2 Mathematical Problems in Engineering

Heat Conduction Integral Equation.
For the steady-state linear heat conduction and radiation problems, it is assumed that the thermal conductivity k is constant for homogeneous material in the calculation region Ω. e heat conduction control equation with a radiative heat source term can be expressed as follows: where T(x, t) represents the temperature of the point x, and because it is a steady-state problem, no time terms exists; k is the thermal conductivity; q r V (x) is the heat source value of the point x. Green's function is introduced, and the weighted integral is carried out in the calculation area by introducing the weight function G. Equation (7) can be obtained by taking the weighted integral in the calculation area, which is as follows: en, we use the partial integration method and Gauss divergence theorem to solve Equation (8).
e integral equation can be obtained by partial integration operation as follows, and the specific derivation process can be seen in [39]: where the coefficient c(x) � 1 when the source point x is located in the domain, and c(x) � 0.5 when the source point is on the boundary of the structure; the integrals in the above formula are expressed as follows: among them, q c � −k(zT/zn) is the heat flux caused by heat conduction; for two-dimensional problems, G � (1/2π)ln(1/r), for three-dimensional problems, G � (1/4πr), and to solve integral Equation (9), we must know the q r V in the equation, so we introduce the radiation integral equation.

Heat Radiation Integral Equation.
To solve the boundary equation of the heat conduction problem, the radiant heat source q r V should be solved first, so the radiation heat transfer of the participatory medium would be considered. We give the integral equation of radiation heat transfer on the boundary, and its boundary medium has the ability to absorb and emit radiation as follows: e integral equation of radiant heat source is as follows: Equation (11) can be applied when the source point is located at the boundary and inside, while Equation (10) is only applicable to the edge. Definition of each parameter can be seen in [45].

Discretization of the Integral Equation
Equations (9)-and (11) are the basic integral equation for solving the steady-state heat conduction-coupled radiation problems by the isogeometric boundary element. Peng et al. [45] solved the domain integral and converted them into the boundary integral by the radial basis function method which is similar to the traditional internal grid division. In this work, the radial basis function method is used to deal with all the domain integrals, and the Gaussian integral formula is used to calculate the converted integral: Equation (13) is the calculation formula of the radial integration method, α � 1 for two-dimensional problems and α � 2 for three-dimensional problems.

Discretization of the Heat Conduction Equation.
Due to the heat source term, q r V (x) was contained in domain integral Equation (8).
erefore, we use the interpolation function to express it as follows: where q rH V is the value of the radiation source term on point H,and M I (q) is the global interpolation function given as follows: where the repetition index represents periodic summation, and the range of I and j is 1 ∼ M S . Bringing Equation (13) into Equation (9) and then using the radial integral method Equation (12), we can get Mathematical Problems in Engineering

(15)
By bringing Equation (11) into Equation (9), the domain integral can be transformed into the boundary integral and only needs to discretize the elements in the calculation domain. For the linear element of a two-dimensional problem, the number of boundary elements is the same as that of boundary nodes. It is assumed that the boundary is discretized into N e th boundary elements and intradomain layout N i th points. e total number of nodes is N e + N i . We take all nodes as source points in turns, and we can use the radial integration method and can get the following equation: where M N A is the matrix of N A * N A , and similar to the calculation of the boundary integral, the following equations can be obtained by Gaussian numerical integration and system set of all boundary elements in the unified integral equation: where U is the normalized temperature kT; H is a matrix of N A * N A ; for the matrix elements corresponding to the internal point temperature, only the diagonal has a value of 1, and the other elements are 0; G is a matrix of N A * N B .

Discretization of the Heat Conduction Equation.
Since the fundamental equation of thermal radiation we gave in Equations (11) and (12), the method of converting the domain integral to the boundary integral is the radial integral method. e black body radiation force E b in the domain integral is unknown. We use the interpolation basis function to represent the following equation: where N I (q) is the global interpolation function; E I b is the black body radiation force of the node; by putting Equation (18) into Equation (12) and using the radial integration method, we can get (20) e domain integral in equation (11) can be obtained using the same method, and the only difference is K 0 (q, p) instead of K(q, p); a is the absorption ratio, which can be solved directly by the radial integration method if it is a constant or an already given known function. if a is just the value of the distribution given at some nodes, we can solve equation 20 and represent it by a global interpolation function.. Since all the parts or constants are known within the integral, it can be solved analytically using Gaussian integration for the calculation.
Without domain integration, the integral boundary equation can be fully obtained by bringing the transformed integral equation into the unified integral equation. Simply by discretizing the boundary unit into linear or quadratic boundary units and using the Gaussian basic formula on each boundary integral, it can be calculated on each boundary cell. Since the radial basis function is used to approximate the unknown quantity in conversion, we often arrange some points inside to get the exact result.
We discretize boundaries into N e th boundary elements and arrange N i internal points in the domain. e total number of points is N e + N i . Taking each node and internal point as the source point in turn, the boundary integral in the unified equation can be calculated to obtain the following matrix equation: In Equations (21) and (22), the matrix Z and Z ′ are N A * N A order matrix; the algebraic equations set each boundary node i as the source point P in turns. e algebraic equations can be obtained by calculating the boundary integral in the formula in each boundary elements as follows: where E b and q r are N A th order sequence vector composed of the boundary blackbody radiation force and boundary radiation heat flux; E b (Q m ) is N A, composed of the medium blackbody radiation force rank vector. For the coefficient matrix of the element expression of square matrix of orderN A , we can refer to the work of Peng et al [45]; then we use the same method in equation (12). e following system of algebraic equations can be obtained: Equation (23) is the boundary node equation, and Equation (24) is applied for all nodes. In this way, the heat radiation equation is discretized.

System of Equations for Heat Conduction and Radiation Coupling Problem and Its Solution
When there is no internal point, all conditions are temperature boundary conditions. e column vector E b is known so that it can directly obtain the unknown boundary quantity q r and then bring it to get the radiative heat source q r V and finally bring the radiative heat source, the radiative heat flux, and the heat conduction boundary condition together into the discrete equation system, and the heat conduction boundary conditions are brought together into the discrete system of equations. We can obtain all the edges by solving the linear system of equations. e unknown quantities in Equation (24) on the boundary are obtained by solving the linear equations. We inverse to solve the unknown q r , which is as follows: Finally, the matrix equation can be obtained by merging. e total power radiated per unit area of a blackbody surface (called the radiance or energy flux density of the object) is directly proportional to the fourth power of the blackbody's thermodynamic temperature T (also known as the absolute temperature).
where x in Equation (26) is the unknown temperature and the column vector composed of the heat flux of the boundary node and the internal point; y is the column vector obtained by multiplying all the known boundary quantities by the corresponding coefficients in Equation (26), and the column element corresponding to the boundary node with available temperature in matrix B is 0. en, we give three kinds of boundary conditions, namely, temperature boundary condition, heat flow boundary conditions, and hybrid boundary conditions, as shown in the respective equations: For the matrices A and B and vector y in Equation (26), the set of groups is different because the values of the elements of their matrices and vectors are related to the positions of the source point p and the field point Q. e elements in matrices A and B and vector y are denoted as A ij , B ij , and y i , and when i and j lie on the boundary, we can get Equation (28) is the point j in the first type boundary conditions.
Equation (29) is the point j in the second type boundary conditions.
. (30) Equation (30) is the point j in the third type boundary conditions. When the point i and point j are boundary points and internal points, respectively, and when they are all internal points, we can find the solutions in [45]. Since Equation (21) is a nonlinear equation system about temperature, it can be solved by the Newton-Raphson iterative method, which can be seen in [39,45,46].

Rectangular Example.
A rectangular domain with a size of 2L * L is composed of the blackbody wall, filled with translucent medium with an absorption ratio a � 1 and thermal conductivity k � 226.76W/(m · k). e bottom temperature of the rectangle is maintained at T 0 � 1000K. Other wall temperatures are T 1 � 500K, and the emissivity of the four walls is ε � 1; to analyze the heat conduction coupled radiation problem of the two-dimensional rectangular model numerically, we discretize the boundary into 60 equally spaced linear cells (20 on the long side and 10 on the short side) and uniformly arrange seven interior points along the symmetry line of x � L. e temperature distribution of the inner points is calculated under the condition of thermal coupling radiative heat transfer. e results of the IGABEM are compared with the results of the COMSOL software Figure 1. e algorithm's accuracy is verified by comparing the IGABEM results with those of COMSOL software. Figure 1(a) shows the 2-D rectangular solution domain and its boundary conditions, and Figure 1(b) shows the boundary cells and their interior points (red points and black points). Firstly, we discuss the comparison between different algorithmic solutions of COMSOL and the IGABEM with the example of the interior points of the black points, and Figure 2(a) shows that the results of the two algorithms are similar. Table 1 shows the comparison between the calculation results of the COMSOL software and the IGABEM. Since there is no analytical solution in the example, COMSOL can be seen as an international standard solution algorithm and can be approximated as an analytical solution to a certain extent. e maximum error is 6.5E − 3 and the minimum error is 1.1E − 4, which verifies the accuracy and stability of the IGABEM algorithm. Due to the more minor degrees of freedom, the computational efficiency of the IGABEM for 2-D problems is better than that of software in terms of time and accuracy. As shown in Table 2, due to the rectangle being an axisymmetric graph, the temperature of the inner points distributed along the axis of symmetry is basically the same.
We arrange the red points and black points in the rectangular area as shown in Figures 1(b), 2(a), and 2(b),which are the temperature distribution curve of red points and black points. e error results only reach 2% when x � 1.0L and y � 0.5L, which is lower than 1% at the other points. e main reason is that the number of discrete Mathematical Problems in Engineering parts of the boundary is too small. In the next example, we can improve the accuracy by increasing the number of discrete parts. Figures 3(a) and 3(b) show the temperature distribution diagram and isotherm in COMSOL. e gure shows that the closer to the bottom boundary, the higher the temperature and the closer to the upper boundary, the lower the temperature. e distribution is not uniform, and the isotherm appears concentrated on both sides of the bottom, indicating that the temperature changes more, in line with the laws of thermodynamics, and more to verify the stability and accuracy of the algorithm.

NURBS Curve Model for Homogeneous and Heterogeneous Material.
e isogeometric boundary element method can eliminate model errors, so its advantages can be applied to the eld of thermal radiation, which we use for    model construction. A blackbody wall is composed of a square domain with a size of 1 * 1 , and a semicircular domain is composed of a quadratic NURBS curve. e interior is lled with a translucent medium with an absorption ratio of a 1.5 and a thermal conductivity of k 400.76W/(m · K), and the bottom temperature is maintained at T 0 800K, other wall temperatures are T 1 200K, and the emissivity of the four walls is ε 1, which can be seen in Figure 4(a). To numerically analyze this conductive coupled radiative heat transfer problem, we discretize it into a linear unit and arrange 14 symmetry points uniformly along the x 0.5 symmetry line in Figure 4(b). e temperature distribution of each internal point is calculated under the condition of thermal conductivity coupled with radiation heat transfer. e algorithm's accuracy is veri ed by comparing the IGABEM results with COMSOL results.
Comparing Figures 5(a) and 5(b), the NURBS curve can be seen with di erent control points but the same weight coe cients; due to the various control points, the same control points' shape varies greatly. Comparing Figures 5(b) and 5(c), we can see that the two Bézier curves with the same control points but di erent weight coe cients have approximately the same trend. Among the three curves, the direction of the second and third curves is roughly the same. Due to the same control points, the local di erence is relatively small. In practical engineering, the standard model components do not exist, and the local modi cation is more conducive to modeling and establishing an accurate twodimensional structural model. Based on this, we build a twodimensional model with a square bottom, as shown in Figure 4(a), and some internal points can be seen in Figure 4(b). Figure 6(a) shows the initial NURBS curve and the position of control points. After normalizing the "knot vector," the parameter space vector of the boundary element can be obtained. erefore, only one boundary element is formed, and the parameter space intervals are given, respectively. By splitting each NURBS element equally and inserting new knots, the new control point sequence "knot vector" and weight factor vector after re nement can be obtained. For example, we insert a new control point at the middle point of each initial NURBS cell parameter space to get the re ned control point sequence shown in Figure 6(b). After the Bézier extraction operation, a new set of control point sequences is obtained, as shown in BE operation knots. For "initial knots" in Figure 6(a) and "inserted knots" in Figure 6(b), we can nd that inserting a new point will also change the position of the original point but still describe the same curve. However, the 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. en, 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 Figures 6(c) and 6(f ), respectively. In the model, its control points are (0, 1), (0.2, 1.7), (0.8, 1.7), and (1, 1), and the weight coe cient is (1, 1, 1, 1).
We give the temperature error table of the IGABEM and COMSOL by the black points, as shown in Table 3. As the y value increases, the temperature value gradually decreases, and it can be seen that the accuracy increases with the increase of the number of discrete parts. When the red points get distributed, we give the error analysis diagram of y 0.7 as shown in Table 4. e error is less than 1%, and the convergence of the result is better than that of the example with few discrete boundaries. Figure 7 shows the error graph of black and red points, respectively. Figure 7(a) shows as the value of y increases, the rate of change of its temperature value becomes slower and slower; that is, the rst derivative of the temperature function is less than zero. e reason is that the smaller the   value of y, the higher the temperature and the faster the di usion. When the points is close to the top, the temperature is less di used and slower, which is corresponding to the laws of thermodynamics. Figure 8 shows the COMSOL model temperature distribution diagram and isotherm gure. From Figures 8(a) and 8(b), we can see that the temperature distribution is in line with the laws of thermodynamics. e closer the isotherm is to the bottom corner of both sides, the denser it is, indicating the magnitude of temperature change; in the structure's design, we can focus on this part of the design and research. It can reduce the probability of temperature stress on the structure to achieve the purpose of structural optimization design and increase the safety and reliability of the structure.
To study the heat transfer properties of heterogeneous materials, we change the thermal conductivity in the previous example into a linear function varying with the value of y, that is, k(y) 400 − 0.2T(y). We discuss the concerned temperature error of the inner point distribution .
e boundary conditions and internal points are shown in Figure 9.
Compared with the homogeneous material, the thermal conductivity is a constant, so the error is relatively small, but for heterogeneous materials whose thermal conductivity is a function of the error is 3.2%, this error is relatively large but acceptable, and the correctness and stability of the algorithm are veri ed. Figure 10 and Table 5 list the temperature distribution between the circular domain and square domain constructed by a two-dimensional Bézier curve. e error between the IGABEMand COMSOL result is no more than 3%, indicating the IGABEM has a sure accuracy in solving nonlinear coupled heat transfer problems. Table 6 gives the e ciency comparison of two calculation methods; from the table, we can see that the computational e ciency of the IGABEM is better than that of COMSOL under the same discrete conditions.

Circle Subtract Square Model.
We give a circle with radius as 1 and center of the circle as the point (0, 0). We remove a square with a length of 1 from its upper part. e wall is a blackbody wall with the absorption ratio a 1 and the thermal conductivity of k 226.76; the red edge temperature of Figure 11(a) is 800K, the other edge temperature is 200K, and the thermal radiation emissivity of the wall is varepsilon 1. As shown in Figure 11, Figure 11(b) is the inner points graph of the model. e red points are the transverse distribution and the black points are the vertical distribution. We discretize it into boundary elements and solve it and compare the results with COMSOL results. We constructed the model with control points (0, 0), (0.2, 0.7), (0.8, 0.7), (1, 0), (0.8, 0.7), and (0.2, 0.7), whose weight coe cients are (1, 1, 1, 1). Figure 12(a) shows the temperature of the red points in Figure 11(b). Since the 2D model is axisymmetric along the yaxis, its temperature should be distributed symmetrically. From the two numerical calculation results, it can be seen that the calculation result of COMSOL is slightly larger than that of IGABEM, but within a reasonable error range, the maximum error is 1% to meet the requirements. Figure 12(b) shows the temperature of the black points in Figure 11(b). With the increase of y value, the temperature of the inner point will rst increase and then decrease. However, since the upper inner point is close to the boundary condition with a temperature of 800K, the temperature of the upper inner point is higher than that of the lower inner point. Figure 13(a) shows the cloud diagram of temperature distribution, Figure 13(b) shows the isotherm diagram, Figure 13(c) shows the cloud diagram of  We can see that when the red point is symmetrical along the X axis, the temperature value of its inner point is also symmetrical along the X axis. When the black point increases along the Y coordinate on one side, its temperature value decreases in accordance with the law of thermodynamics.

Conclusions
is work applies for the IGABEM to solve two-dimensional steady-state heat conduction and radiation problems. e NURBS method is used to construct a smooth geometric model, which eliminates geometric errors and improves the computational accuracy. e radial integration method converts the domain integrals due to variable coe cients into boundary integrals. Numerical results show that the developed algorithm e ectively solves the nonhomogeneous steady-state coupled heat conduction and heat radiation problems with variable coe cients. Since the theory is essentially the same, in the future, the algorithm can be extended to transient analysis also to solve three-dimensional problems. From the proposed numerical example, a stepping stone is provided for the implementation of the fully integrated CAD and CAE software.

Data Availability
Data are openly available in a public repository.