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

Authors: Wei Jiang; Jiancheng Xu; Li Kunpeng; Ye Ouyang; Jinzhou Yan
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

We propose an efficient isogeometric boundary element method to address the coupling of heat conduction and radiation in homogeneous or inhomogeneous materials. The isogeometric boundary element method is used to construct irregular 2D models, which eliminate errors in model construction. The physical unknowns in the governing equations for heat conduction and radiation are discretized using an interpolation approximation, and the integral equations are finally 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. The results of the three numerical cases show that the adopted algorithm can improve the computational accuracy and efficiency.


---

## Body

## 1. 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. The traditional finite 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 first, 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 infinite domains in the traditional field. 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 high-order continuity and flexible refinement schemes, and it can be extended to many 2D modeling fields. Therefore, the IGABEM has been widely studied by many scholars and has been successfully applied in many fields such as elasticity [2, 3], fluid-structure coupling [4–6], shape optimization [7, 8], and acoustic [9–16].The heat conduction and radiation require coupling the two physical fields, which can be significantly 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 finite 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 specific direction. Zhou et al. [20] combined the finite 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–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–33] studied the boundary collocation method for anomalous heat conduction analysis in functionally graded materials. Nie et al. [34–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–43] studied the heat conduction analysis of two-dimensional and three-dimensional geometric boundary element methods.The 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. The remainder of this paper is organized as follows: The theory of the isogeometric analysis is given in Section2. The coupled heat conduction and radiation are put in Section 3. Then, 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.


## 2. 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. Therefore, 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Ξ=ξ1,ξ2,…,ξn+p+1,ξa∈R.Where a is the knot index,p is the curve degree, and n is the number of basis functions or control points Pa. Each ξa∈Ξ is called a knot. The B-spline basis functions of degree p are defined recursively; starting with p = 0, we define:(1)Na,0=1,ζa≤ζ<ζa+10,otherwise,(2)Na,pζ=ζ−ζaζa+p−ζaNa,p−1ζ+ζa+p+1−ζζa+p+1−ζa+1Ni+1,p−1ζThus, 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 Pa. 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. Therefore, in this work, the descriptions of boundary geometry or physical quantities are given in a knot span withp + 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:(3)xξ=∑a=1p+1Na,pξPa,where Pa 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. The NURBS geometry is a weighted form of the B-spline definition, i.e.,(4)xξ=∑a=1p+1Ra,pξPa,where Ra,p are the NURBS basis functions, which are defined by(5)Ra,pξ=Na,pξWaWξWξ=∑a=1p+1Na,pξwa,where a denotes the control point index, Na,p is the B-spline basis function from equation (2), and wa is a weight associated with control point Pa when all the weights are of equal value.The 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.


## 3. Basic Governing Equations of Heat Conduction and Radiation

The 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. This paper expounds on using the isogeometric boundary element method to solve this problem.


### 3.1. Heat Conduction Integral Equation

For the steady-state linear heat conduction and radiation problems, it is assumed that the thermal conductivityk is constant for homogeneous material in the calculation region Ω. The heat conduction control equation with a radiative heat source term can be expressed as follows:(6)k∇2Tx+qVrx=0,where Tx,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; qVrx 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:(7)∫ΩG˜k∇2TxdΩ+∫ΩG˜qVrxdΩ=0.Then, we use the partial integration method and Gauss divergence theorem to solve Equation (8). The integral equation can be obtained by partial integration operation as follows, and the specific derivation process can be seen in [39]:(8)cxkTx+G˜qcx=−kG˜x+G˜xqVrxx,where the coefficient cx=1 when the source point x is located in the domain, and cx=0.5 when the source point is on the boundary of the structure; the integrals in the above formula are expressed as follows:(9)G˜qcx≔∫ΓG˜x,yqcxdΓykG˜x,yx≔−∫Γ∂G˜x,y∂xikTxdΓyG˜x,yqVrx≔∫ΩG˜x,yqVrydΩy,among them, qc=−k∂T/∂n is the heat flux caused by heat conduction; for two-dimensional problems, G˜=1/2πln1/r, for three-dimensional problems, G˜=1/4πr, and to solve integral Equation (9), we must know the qVr in the equation, so we introduce the radiation integral equation.


### 3.2. Heat Radiation Integral Equation

To solve the boundary equation of the heat conduction problem, the radiant heat sourceqVr 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:(10)qrPεP+EbP=∫ΓEbQ+1−εQEbQqrQKQPQ,PdΓQ+∫ΩaqEbqK0q,pdΩq.The integral equation of radiant heat source is as follows:(11)qVrPaP+c′EbP=∫ΓEbQ+1−εQEbQqrQKQQ,pdΓQ+∫ΩaqEbqK0q,pdΩq.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].


## 3.1. Heat Conduction Integral Equation

For the steady-state linear heat conduction and radiation problems, it is assumed that the thermal conductivityk is constant for homogeneous material in the calculation region Ω. The heat conduction control equation with a radiative heat source term can be expressed as follows:(6)k∇2Tx+qVrx=0,where Tx,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; qVrx 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:(7)∫ΩG˜k∇2TxdΩ+∫ΩG˜qVrxdΩ=0.Then, we use the partial integration method and Gauss divergence theorem to solve Equation (8). The integral equation can be obtained by partial integration operation as follows, and the specific derivation process can be seen in [39]:(8)cxkTx+G˜qcx=−kG˜x+G˜xqVrxx,where the coefficient cx=1 when the source point x is located in the domain, and cx=0.5 when the source point is on the boundary of the structure; the integrals in the above formula are expressed as follows:(9)G˜qcx≔∫ΓG˜x,yqcxdΓykG˜x,yx≔−∫Γ∂G˜x,y∂xikTxdΓyG˜x,yqVrx≔∫ΩG˜x,yqVrydΩy,among them, qc=−k∂T/∂n is the heat flux caused by heat conduction; for two-dimensional problems, G˜=1/2πln1/r, for three-dimensional problems, G˜=1/4πr, and to solve integral Equation (9), we must know the qVr in the equation, so we introduce the radiation integral equation.


## 3.2. Heat Radiation Integral Equation

To solve the boundary equation of the heat conduction problem, the radiant heat sourceqVr 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:(10)qrPεP+EbP=∫ΓEbQ+1−εQEbQqrQKQPQ,PdΓQ+∫ΩaqEbqK0q,pdΩq.The integral equation of radiant heat source is as follows:(11)qVrPaP+c′EbP=∫ΓEbQ+1−εQEbQqrQKQQ,pdΓQ+∫ΩaqEbqK0q,pdΩq.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].


## 4. 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:(12)∫ΩfxdΩ=∫Γ1rαQ,p∂r∂n∫0rQ,pfxrαdr.Equation (13) is the calculation formula of the radial integration method, α=1 for two-dimensional problems and α=2 for three-dimensional problems.


### 4.1. Discretization of the Heat Conduction Equation

Due to the heat source term,qVrx was contained in domain integral Equation (8). Therefore, we use the interpolation function to express it as follows:(13)qVrx=MqTqVr=MIqqVrH,where qVrH is the value of the radiation source term on point H,and MIq is the global interpolation function given as follows:(14)MIq=ϕJqϕJI−1,where the repetition index represents periodic summation, and the range of I and j is 1∼MS. Bringing Equation (13) into Equation (9) and then using the radial integral method Equation (12), we can get(15)∫ΩG˜x,yqVrxdΩx=qVrH∫Γ1rQ,p∂r∂n∫0rQ,pG˜x,yMIqrx,ydrdΓQ.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 Ne th boundary elements and intradomain layout Ni th points. The total number of nodes is Ne+Ni. We take all nodes as source points in turns, and we can use the radial integration method and can get the following equation:(16)∫ΩG∗qVrdΩ=MNAqVr,where MNA is the matrix of NA∗NA , 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:(17)HU−Gq=MNAqVr,where U is the normalized temperature kT; H is a matrix of NA∗NA; 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 NA∗NB.


### 4.2. 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. The black body radiation force Eb in the domain integral is unknown. We use the interpolation basis function to represent the following equation:(18)Ebq=NIqEbI,where NIq is the global interpolation function; EbI 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(19)∫ΩaqEbqK0q,pdΩq=EbI∫Γ1rq,p∂r∂nFQ,pdΓq,(20)FQ,p=∫0rQ,paqK0q,pNIqrQ,pdrq.The domain integral in equation (11) can be obtained using the same method, and the only difference is K0q,p instead of Kq,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 intoNe th boundary elements and arrange Ni internal points in the domain. The total number of points is Ne+Ni. 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:(21)∫ΩaqEbKq,PdΩq=ZEb,(22)∫ΩaqEbK0q,PdΩq=Z′Eb.In Equations (21) and (22), the matrix Z and Z′ are NA∗NA order matrix; the algebraic equations set each boundary node i as the source point P in turns. The algebraic equations can be obtained by calculating the boundary integral in the formula in each boundary elements as follows:(23)HEb+Gqr+ZEbQm=0,where Eb and qr are NA th order sequence vector composed of the boundary blackbody radiation force and boundary radiation heat flux; EbQm is NA, composed of the medium blackbody radiation force rank vector. For the coefficient matrix of the element expression of square matrix of orderNA, we can refer to the work of Peng et al [45]; then we use the same method in equation (12). The following system of algebraic equations can be obtained:(24)qVr=ZEb+Nqr+WEbQm.Equation (23) is the boundary node equation, and Equation (24) is applied for all nodes. In this way, the heat radiation equation is discretized.


## 4.1. Discretization of the Heat Conduction Equation

Due to the heat source term,qVrx was contained in domain integral Equation (8). Therefore, we use the interpolation function to express it as follows:(13)qVrx=MqTqVr=MIqqVrH,where qVrH is the value of the radiation source term on point H,and MIq is the global interpolation function given as follows:(14)MIq=ϕJqϕJI−1,where the repetition index represents periodic summation, and the range of I and j is 1∼MS. Bringing Equation (13) into Equation (9) and then using the radial integral method Equation (12), we can get(15)∫ΩG˜x,yqVrxdΩx=qVrH∫Γ1rQ,p∂r∂n∫0rQ,pG˜x,yMIqrx,ydrdΓQ.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 Ne th boundary elements and intradomain layout Ni th points. The total number of nodes is Ne+Ni. We take all nodes as source points in turns, and we can use the radial integration method and can get the following equation:(16)∫ΩG∗qVrdΩ=MNAqVr,where MNA is the matrix of NA∗NA , 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:(17)HU−Gq=MNAqVr,where U is the normalized temperature kT; H is a matrix of NA∗NA; 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 NA∗NB.


## 4.2. 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. The black body radiation force Eb in the domain integral is unknown. We use the interpolation basis function to represent the following equation:(18)Ebq=NIqEbI,where NIq is the global interpolation function; EbI 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(19)∫ΩaqEbqK0q,pdΩq=EbI∫Γ1rq,p∂r∂nFQ,pdΓq,(20)FQ,p=∫0rQ,paqK0q,pNIqrQ,pdrq.The domain integral in equation (11) can be obtained using the same method, and the only difference is K0q,p instead of Kq,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 intoNe th boundary elements and arrange Ni internal points in the domain. The total number of points is Ne+Ni. 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:(21)∫ΩaqEbKq,PdΩq=ZEb,(22)∫ΩaqEbK0q,PdΩq=Z′Eb.In Equations (21) and (22), the matrix Z and Z′ are NA∗NA order matrix; the algebraic equations set each boundary node i as the source point P in turns. The algebraic equations can be obtained by calculating the boundary integral in the formula in each boundary elements as follows:(23)HEb+Gqr+ZEbQm=0,where Eb and qr are NA th order sequence vector composed of the boundary blackbody radiation force and boundary radiation heat flux; EbQm is NA, composed of the medium blackbody radiation force rank vector. For the coefficient matrix of the element expression of square matrix of orderNA, we can refer to the work of Peng et al [45]; then we use the same method in equation (12). The following system of algebraic equations can be obtained:(24)qVr=ZEb+Nqr+WEbQm.Equation (23) is the boundary node equation, and Equation (24) is applied for all nodes. In this way, the heat radiation equation is discretized.


## 5. 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. The column vectorEb is known so that it can directly obtain the unknown boundary quantity qr and then bring it to get the radiative heat source qVr 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. The unknown quantities in Equation (24) on the boundary are obtained by solving the linear equations. We inverse to solve the unknown qr, which is as follows:(25)qr=−ZEbQm−WEbG.Finally, the matrix equation can be obtained by merging. The 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 temperatureT (also known as the absolute temperature).(26)Ax+BT4=y,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. Then, 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:(27)Tx=T¯xx∈Γ1qx=q¯xx∈Γ2qx=hTx−Tfx+εσT4x−Tf4xx∈Γ3..For the matricesA 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. The elements in matrices A and B and vector y are denoted as Aij, Bij, and yi, and when i and j lie on the boundary, we can get(28)Aij=−GijBij=0yi=yi−HijkTj−σQijTj4..Equation (28) is the point j in the first type boundary conditions.(29)Aij=kjHijBij=σQijyi=yi+Gijqj..Equation (29) is the point j in the second type boundary conditions.(30)Aij=−GijBij=0yi=yi−HijkTj−σQijTj4..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].


## 6. Numerical Examples

### 6.1. Rectangular Example

A rectangular domain with a size of2L∗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. The bottom temperature of the rectangle is maintained at T0=1000K. Other wall temperatures are T1=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. The temperature distribution of the inner points is calculated under the condition of thermal coupling radiative heat transfer. The results of the IGABEM are compared with the results of the COMSOL software Figure 1. The algorithm’s accuracy is verified by comparing the IGABEM results with those of COMSOL software.Figure 1
2D rectangular domain boundary conditions and boundary cell inner points. (a) 2D rectangular domain and boundary conditions. (b) 2D rectangular boundary cell and inner points.
(a)(b)Figure1(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.Figure 2
2D rectangular domain boundary conditions and boundary cell inner points. (a) Dimensionless temperature distribution curve along black interior points in 2D rectangular. (b) Red point temperature in 2-D.
(a)(b)Table1 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. The 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.Table 1
Error table of 2-D rectangular black interior points.
x = L, y coordinate0.125L0.25L0.375L0.5L0.625L0.75L0.875LIGABEM926.45853.76784.07718.76658.30602.34549.96COMSOL926.56854.98786.72722.56662.56606.17552.40Error1.1E − 41.4E − 33.4E − 35.3E − 36.5E − 36.4E − 34.4E − 3Table 2
Error table of 2-D rectangular red interior points.
y = 0.5L, x coordinate0.1L0.4L0.7L1.0L1.3L1.6L1.9LIGABEM548.73659.82709.55722.56709.55659.82548.73COMSOL547.75657.07706.02706.02706.02657.07547.75Error1.8E − 34.2E − 35.0E − 32.3E − 35.0E − 34.2E − 31.8E − 3We arrange the red points and black points in the rectangular area as shown in Figures1(b), 2(a), and 2(b),which are the temperature distribution curve of red points and black points. The error results only reach 2% when x = 1.0L and y = 0.5L, which is lower than 1% at the other points. The main reason is that the number of discrete parts of the boundary is too small. In the next example, we can improve the accuracy by increasing the number of discrete parts.Figures3(a) and 3(b) show the temperature distribution diagram and isotherm in COMSOL. The figure shows that the closer to the bottom boundary, the higher the temperature and the closer to the upper boundary, the lower the temperature. The 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.Figure 3
COMSOL model temperature distribution diagram and its isotherm. (a) COMSOL model temperature distribution diagram. (b) Isotherm of the COMSOL software.
(a)(b)


### 6.2. NURBS Curve Model for Homogeneous and Heterogeneous Material

The isogeometric boundary element method can eliminate model errors, so its advantages can be applied to the field of thermal radiation, which we use for model construction. A blackbody wall is composed of a square domain with a size of1∗1, and a semicircular domain is composed of a quadratic NURBS curve. The interior is filled 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 T0=800K, other wall temperatures are T1=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). The temperature distribution of each internal point is calculated under the condition of thermal conductivity coupled with radiation heat transfer. The algorithm’s accuracy is verified by comparing the IGABEM results with COMSOL results.Figure 4
Two-dimensional square semicircle region and boundary conditions and vertical inner points of boundary elements. (a) 2-D square semicircle region and boundary condition element. (b) Black points of boundary.
(a)(b)Comparing Figures5(a) and 5(b), the NURBS curve can be seen with different control points but the same weight coefficients; 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 different weight coefficients 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 difference is relatively small. In practical engineering, the standard model components do not exist, and the local modification is more conducive to modeling and establishing an accurate two-dimensional structural model. Based on this, we build a two-dimensional model with a square bottom, as shown in Figure 4(a), and some internal points can be seen in Figure 4(b).Figure 5
NURBS curves constructed with different control points and weight coefficients. (a) Control points (0, 0), (0.2, 0.7), (0.8, 0.7), and (1, 0), with the weight coefficient (1, 1, 1, 1); (b) control points (0, 0), (0, 1), (1, −0.5), and (1, 0.5), with the weight coefficient (1, 1, 1, 1); (c) control points (0, 0), (0, 1), (1, −0.5), and (1, 0.5), with the weight coefficient (1, 2, 1, 1).
(a)(b)(c)Figure6(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. Therefore, 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 refinement 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 refined 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 find 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. 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 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 coefficient is (1, 1, 1, 1).Figure 6
(a) Initial. (b) Refine1. (c) Refine2. (d) Refine3. (e) Refine4. (f) Refine5.
(a)(b)(c)(d)(e)(f)We give the temperature error table of the IGABEM and COMSOL by the black points, as shown in Table3. 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. The error is less than 1%, and the convergence of the result is better than that of the example with few discrete boundaries.Table 3
Black point error table of the square and semicircle model.
x = 0.5, y coordinate0.10.30.50.70.91.11.3IGABEM682.95484.77356.04282.61242.39220.24207.57COMSOL682.68485.59358.16285.09244.57221.78208.38Error3.9E − 41.7E − 36.0E − 38.8E − 39.0E − 37.0E − 33.9E − 3Table 4
Red point error table of the square and semicircle model.
y = 0.7, x coordinate0.050.20.350.50.650.80.95IGABEM213.10249.29274.30283.15274.30249.293213.10COMSOL213.55250.57276.09285.09276.09250.57213.55Error2.1E − 35.1E − 36.5E − 36.9E − 36.5E − 35.1E − 32.1E − 3Figure7 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 first derivative of the temperature function is less than zero. The reason is that the smaller the value of y, the higher the temperature and the faster the diffusion. When the points is close to the top, the temperature is less diffused and slower, which is corresponding to the laws of thermodynamics.Figure 7
Error graph of black and red points in example 2. (a) Dimensionless temperature distribution of black points of the square and semicircle model. (b) Error between the IGABEM and COMSOL of red points of square.
(a)(b)Figure8 shows the COMSOL model temperature distribution diagram and isotherm figure. From Figures 8(a) and 8(b), we can see that the temperature distribution is in line with the laws of thermodynamics. The 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.Figure 8
COMSOL model temperature distribution diagram and isotherm. (a) COMSOL model temperature distribution diagram. (b) Isotherm.
(a)(b)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 ofy, that is, ky=400−0.2Ty. We discuss the concerned temperature error of the inner point distribution . The boundary conditions and internal points are shown in Figure 9.Figure 9
Two-dimensional square semicircle region and boundary conditions (a); vertical inner points of boundary elements (b). (a) Boundary conditions. (b) Inner point of boundary elements.
(a)(b)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 verified.Figure10 and Table 5 list the temperature distribution between the circular domain and square domain constructed by a two-dimensional Bézier curve. The 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.Figure 10
Compared the COMSOL and the IGABEM for nonhomogeneous materials.Table 5
Comparison of nonhomogeneous error tables.
x = 0.5, y coordinate0.10.30.50.70.91.11.3IGABEM928.95786.53644.80505.57372.69254.87163.03COMSOL933.29802.25665.51521.70379.38255.82163.51Error4.7E − 32.0E − 23.2E − 23.2E − 21.8E − 33.7E − 32.9E − 3Table6 gives the efficiency comparison of two calculation methods; from the table, we can see that the computational efficiency of the IGABEM is better than that of COMSOL under the same discrete conditions.Table 6
Comparison time of nonhomogeneous error tables.
Time10 interior points (S)20 interior points (S)30 interior points (S)40 interior points (S)IGABEM8.369.9512.0214.22COMSOL11.0913.9517.2019.52


### 6.3. 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. The wall is a blackbody wall with the absorption ratioa = 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. The 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 coefficients are (1, 1, 1, 1).Figure 11
Boundary conditions and different inner points of boundary elements. (a) Boundary conditions. (b) Inner point of boundary elements.
(a)(b)Figure12(a) shows the temperature of the red points in Figure 11(b). Since the 2D model is axisymmetric along the y-axis, 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 first 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 radiation distribution, and Figure 13(d) shows the isoradiance diagram.Figure 12
Inner point temperature diagram. (a) Temperature diagram of internal points in horizontal chart distribution. (b) Vertical distribution inner point temperature.
(a)(b)Figure 13
Inner point temperature and radiation diagram. (a) Temperature pattern of the example circle subtract square. (b) Isotherm of the example circle subtract square. (c) Internal radiation of the example circle subtract square. (d) Isoradiance curve of the example circle subtract square.
(a)(b)(c)(d)Tables7 and 8 show the tabulation of Figures 12(a) and 12(b). 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.Table 7
Comparison of temperature of red points in example3.
y = −0.6, x coordinates−0.7−0.5−0.3−0.10.10.30.50.7IGABEM236.84301.64347.48370.57370.50347.28301.34236.66COMSOL237.98304.22350.80374.20374.20350.79304.22237.99Error4.8E − 38.5E − 39.4E − 39.7E − 39.9E − 31.0E − 29.5E − 35.6E − 3Table 8
Comparison of temperature of black points in example3.
x = −0.7, y coordinates−0.6−0.4−0.200.20.40.6IGABEM236.66310.56396.33479.49509.39473.61352.37COMSOL237.99313.74401.06484.51512.13475.24352.49Error5.6E − 31.0E − 21.2E − 21.0E − 25.4E − 33.4E − 23.2E − 4


## 6.1. Rectangular Example

A rectangular domain with a size of2L∗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. The bottom temperature of the rectangle is maintained at T0=1000K. Other wall temperatures are T1=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. The temperature distribution of the inner points is calculated under the condition of thermal coupling radiative heat transfer. The results of the IGABEM are compared with the results of the COMSOL software Figure 1. The algorithm’s accuracy is verified by comparing the IGABEM results with those of COMSOL software.Figure 1
2D rectangular domain boundary conditions and boundary cell inner points. (a) 2D rectangular domain and boundary conditions. (b) 2D rectangular boundary cell and inner points.
(a)(b)Figure1(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.Figure 2
2D rectangular domain boundary conditions and boundary cell inner points. (a) Dimensionless temperature distribution curve along black interior points in 2D rectangular. (b) Red point temperature in 2-D.
(a)(b)Table1 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. The 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.Table 1
Error table of 2-D rectangular black interior points.
x = L, y coordinate0.125L0.25L0.375L0.5L0.625L0.75L0.875LIGABEM926.45853.76784.07718.76658.30602.34549.96COMSOL926.56854.98786.72722.56662.56606.17552.40Error1.1E − 41.4E − 33.4E − 35.3E − 36.5E − 36.4E − 34.4E − 3Table 2
Error table of 2-D rectangular red interior points.
y = 0.5L, x coordinate0.1L0.4L0.7L1.0L1.3L1.6L1.9LIGABEM548.73659.82709.55722.56709.55659.82548.73COMSOL547.75657.07706.02706.02706.02657.07547.75Error1.8E − 34.2E − 35.0E − 32.3E − 35.0E − 34.2E − 31.8E − 3We arrange the red points and black points in the rectangular area as shown in Figures1(b), 2(a), and 2(b),which are the temperature distribution curve of red points and black points. The error results only reach 2% when x = 1.0L and y = 0.5L, which is lower than 1% at the other points. The main reason is that the number of discrete parts of the boundary is too small. In the next example, we can improve the accuracy by increasing the number of discrete parts.Figures3(a) and 3(b) show the temperature distribution diagram and isotherm in COMSOL. The figure shows that the closer to the bottom boundary, the higher the temperature and the closer to the upper boundary, the lower the temperature. The 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.Figure 3
COMSOL model temperature distribution diagram and its isotherm. (a) COMSOL model temperature distribution diagram. (b) Isotherm of the COMSOL software.
(a)(b)


## 6.2. NURBS Curve Model for Homogeneous and Heterogeneous Material

The isogeometric boundary element method can eliminate model errors, so its advantages can be applied to the field of thermal radiation, which we use for model construction. A blackbody wall is composed of a square domain with a size of1∗1, and a semicircular domain is composed of a quadratic NURBS curve. The interior is filled 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 T0=800K, other wall temperatures are T1=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). The temperature distribution of each internal point is calculated under the condition of thermal conductivity coupled with radiation heat transfer. The algorithm’s accuracy is verified by comparing the IGABEM results with COMSOL results.Figure 4
Two-dimensional square semicircle region and boundary conditions and vertical inner points of boundary elements. (a) 2-D square semicircle region and boundary condition element. (b) Black points of boundary.
(a)(b)Comparing Figures5(a) and 5(b), the NURBS curve can be seen with different control points but the same weight coefficients; 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 different weight coefficients 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 difference is relatively small. In practical engineering, the standard model components do not exist, and the local modification is more conducive to modeling and establishing an accurate two-dimensional structural model. Based on this, we build a two-dimensional model with a square bottom, as shown in Figure 4(a), and some internal points can be seen in Figure 4(b).Figure 5
NURBS curves constructed with different control points and weight coefficients. (a) Control points (0, 0), (0.2, 0.7), (0.8, 0.7), and (1, 0), with the weight coefficient (1, 1, 1, 1); (b) control points (0, 0), (0, 1), (1, −0.5), and (1, 0.5), with the weight coefficient (1, 1, 1, 1); (c) control points (0, 0), (0, 1), (1, −0.5), and (1, 0.5), with the weight coefficient (1, 2, 1, 1).
(a)(b)(c)Figure6(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. Therefore, 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 refinement 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 refined 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 find 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. 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 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 coefficient is (1, 1, 1, 1).Figure 6
(a) Initial. (b) Refine1. (c) Refine2. (d) Refine3. (e) Refine4. (f) Refine5.
(a)(b)(c)(d)(e)(f)We give the temperature error table of the IGABEM and COMSOL by the black points, as shown in Table3. 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. The error is less than 1%, and the convergence of the result is better than that of the example with few discrete boundaries.Table 3
Black point error table of the square and semicircle model.
x = 0.5, y coordinate0.10.30.50.70.91.11.3IGABEM682.95484.77356.04282.61242.39220.24207.57COMSOL682.68485.59358.16285.09244.57221.78208.38Error3.9E − 41.7E − 36.0E − 38.8E − 39.0E − 37.0E − 33.9E − 3Table 4
Red point error table of the square and semicircle model.
y = 0.7, x coordinate0.050.20.350.50.650.80.95IGABEM213.10249.29274.30283.15274.30249.293213.10COMSOL213.55250.57276.09285.09276.09250.57213.55Error2.1E − 35.1E − 36.5E − 36.9E − 36.5E − 35.1E − 32.1E − 3Figure7 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 first derivative of the temperature function is less than zero. The reason is that the smaller the value of y, the higher the temperature and the faster the diffusion. When the points is close to the top, the temperature is less diffused and slower, which is corresponding to the laws of thermodynamics.Figure 7
Error graph of black and red points in example 2. (a) Dimensionless temperature distribution of black points of the square and semicircle model. (b) Error between the IGABEM and COMSOL of red points of square.
(a)(b)Figure8 shows the COMSOL model temperature distribution diagram and isotherm figure. From Figures 8(a) and 8(b), we can see that the temperature distribution is in line with the laws of thermodynamics. The 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.Figure 8
COMSOL model temperature distribution diagram and isotherm. (a) COMSOL model temperature distribution diagram. (b) Isotherm.
(a)(b)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 ofy, that is, ky=400−0.2Ty. We discuss the concerned temperature error of the inner point distribution . The boundary conditions and internal points are shown in Figure 9.Figure 9
Two-dimensional square semicircle region and boundary conditions (a); vertical inner points of boundary elements (b). (a) Boundary conditions. (b) Inner point of boundary elements.
(a)(b)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 verified.Figure10 and Table 5 list the temperature distribution between the circular domain and square domain constructed by a two-dimensional Bézier curve. The 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.Figure 10
Compared the COMSOL and the IGABEM for nonhomogeneous materials.Table 5
Comparison of nonhomogeneous error tables.
x = 0.5, y coordinate0.10.30.50.70.91.11.3IGABEM928.95786.53644.80505.57372.69254.87163.03COMSOL933.29802.25665.51521.70379.38255.82163.51Error4.7E − 32.0E − 23.2E − 23.2E − 21.8E − 33.7E − 32.9E − 3Table6 gives the efficiency comparison of two calculation methods; from the table, we can see that the computational efficiency of the IGABEM is better than that of COMSOL under the same discrete conditions.Table 6
Comparison time of nonhomogeneous error tables.
Time10 interior points (S)20 interior points (S)30 interior points (S)40 interior points (S)IGABEM8.369.9512.0214.22COMSOL11.0913.9517.2019.52


## 6.3. 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. The wall is a blackbody wall with the absorption ratioa = 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. The 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 coefficients are (1, 1, 1, 1).Figure 11
Boundary conditions and different inner points of boundary elements. (a) Boundary conditions. (b) Inner point of boundary elements.
(a)(b)Figure12(a) shows the temperature of the red points in Figure 11(b). Since the 2D model is axisymmetric along the y-axis, 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 first 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 radiation distribution, and Figure 13(d) shows the isoradiance diagram.Figure 12
Inner point temperature diagram. (a) Temperature diagram of internal points in horizontal chart distribution. (b) Vertical distribution inner point temperature.
(a)(b)Figure 13
Inner point temperature and radiation diagram. (a) Temperature pattern of the example circle subtract square. (b) Isotherm of the example circle subtract square. (c) Internal radiation of the example circle subtract square. (d) Isoradiance curve of the example circle subtract square.
(a)(b)(c)(d)Tables7 and 8 show the tabulation of Figures 12(a) and 12(b). 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.Table 7
Comparison of temperature of red points in example3.
y = −0.6, x coordinates−0.7−0.5−0.3−0.10.10.30.50.7IGABEM236.84301.64347.48370.57370.50347.28301.34236.66COMSOL237.98304.22350.80374.20374.20350.79304.22237.99Error4.8E − 38.5E − 39.4E − 39.7E − 39.9E − 31.0E − 29.5E − 35.6E − 3Table 8
Comparison of temperature of black points in example3.
x = −0.7, y coordinates−0.6−0.4−0.200.20.40.6IGABEM236.66310.56396.33479.49509.39473.61352.37COMSOL237.99313.74401.06484.51512.13475.24352.49Error5.6E − 31.0E − 21.2E − 21.0E − 25.4E − 33.4E − 23.2E − 4


## 7. Conclusions

This work applies for the IGABEM to solve two-dimensional steady-state heat conduction and radiation problems. The NURBS method is used to construct a smooth geometric model, which eliminates geometric errors and improves the computational accuracy. The radial integration method converts the domain integrals due to variable coefficients into boundary integrals. Numerical results show that the developed algorithm effectively solves the nonhomogeneous steady-state coupled heat conduction and heat radiation problems with variable coefficients. 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.



---

*Source: 4209720-2022-09-07.xml*