Resolving Domain Integral Issues in Isogeometric Boundary Element Methods
via Radial Integration: A Study of Thermoelastic Analysis

The paper applied the isogeometric boundary element method (IGABEM) to thermoelastic problems. The Non-Uniform Rational B-splines (NURBS) used to construct geometric models are employed to discretize the boundary integral formulation of the governing equation. Due to the existence of thermal stress, the domain integral term appears in the boundary integral equation. We resolve this problem by incorporating radial integration method into IGABEM which converts the domain integral to the boundary integral. In this way, IGABEM can maintain its advantages in dimensionality reduction and more importantly, seamless integration of CAD and numerical analysis based on boundary representation. The algorithm is verified by numerical examples.


Introduction
Isogeometric analysis (IGA) has received considerable attention in engineering science ever since the seminal works of Hughes et al. [1]. IGA adopts spline functions used to construct the geometric model in CAD, for example Non-Uniform Rational B-splines (NURBS), as the basis function to discretize partial differential equations. The main advantage of IGA lies in its capability of performing numerical analysis from CAD models directly without needing meshing, and therefore avoids cumbersome preprocessing procedure and geometric errors as encountered by traditional numerical methods such as finite element methods (FEM) and boundary element methods (BEM). Other benefits include high order continuity, flexible refinement scheme, etc. IGA has found wide applications in a number of areas, such as elasticity mechanics [2,3], fluid dynamics [4], fluid-structure interaction [5][6][7], shape optimization [8][9][10], structural vibration [11], shell analysis [12], contact problems [13], composite layering [14], etc. In addition, advancements were made in IGA by utilizing T-splines [4], PHT-splines [15,16], and subdivision surfaces [17] that allow for the analysis of complicated geometries. The introduction of Bézier extraction technique further eases the implementation of IGA with the existing FEM codes [18].
In elasticity mechanics, thermal stress plays a critical role in many engineering applications. However, it is not a trivial task to perform thermoelasticity analysis with (IGA)BEM because an additional domain integral term appears which must be treated carefully [55,56]. Although one may discretize the analysis domain as in FEM, the advantages of (IGA)BEM in dimensionality reduction is sacrificed. Hence, to maintain the boundary representation properties of (IGA)BEM, the domain integral should be converted to boundary integral with the aid of some special techniques. Fang et al. [57] utilized Galerkin vector in IGABEM for thermoelasticity but the availability of this technique only holds for steady heat conduction, that is, the variation of temperature should satisfy the Laplace equation. The radial integration method proposed by [58] might be a feasible alternative and is investigated in the present work. To the authors' best knowledge, this is the first time that the radial integration method is integrated into IGABEM, which is not only significant for the effective use of IGABEM in thermoelasticity analysis, but more importantly, paves the way towards a general approach of addressing domain integrals in IGABEM. The remainder of the paper is organized as follows. Section 2 outlines the NURBS. Section 3 introduces the boundary integral equation for thermoelasticity. Section 4 details implementation of the isogeometric boundary element method in the field of thermoelasticity. Section 5 provides some numerical examples to verify the present algorithm, followed by the conclusion in Section 6.

NURBS
Non-Uniform Rational B-splines (NURBS) are used in the present work since they are ubiquitous in the CAD industry. For the sake of completeness, the fundamentals of NURBS are briefed in this section. Before the introduction of NURBS, it is necessary to start with the concept of B-splines. A B-spline curve is expressed by where x is the Cartesian coordinate of a point located on the B-spline curve, n the number of control points, p the polynomial order, P i the control point, and N i,p B-spline basis functions. Given a knot vector Ξ = [ξ 0 ,ξ 1 , …,ξ m ], the B-spline basis function is defined as follows: N i;p ðnÞ¼ n À n i n iþp À n i N i;pÀ1 ðnÞ þ n iþpþ1 À n n iþpþ1 À n iþ1 N iþ1;pÀ1 ðnÞ As an extension of B-splines by adding the weight coefficient w i , NURBS curve is parameterized as follows 586 CMES, 2020, vol.124, no.2 where R i,p (ξ) is the NURBS basis function which is written as The advantage of NURBS over B-splines is that NURBS can exactly represent all conic sections. Take a quarter circle as an example. As shown in Fig. 1, the weight coefficients w 0 = w 2 = 1, and the value of w 1 varies. As can be seen from the figure, the larger the weight coefficient w 1 , the closer the NURBS curve is to the control point. When a quarter arc is recovered [59]. The h-refinement (knot insertion) refers to inserting a new knot " n on the basis of the original knot vector. This process adds control points without changing the original geometry, and finally achieves the purpose of  increasing the model degree of freedom. Assuming the newly inserted knot " n is in the interval [ξ k , ξ k + 1 ), the weight coefficient corresponding to the new control point can be defined as where The corresponding control point " P i can be expressed as It is noteworthy that although the continuity of the basis is reduced by one for each repetition of a give knot value, the continuity of the curve is preserved.

Boundary Integral Equation for Thermoelasticity
The constitutive relationship in thermoelasticity is expressed by where T represents the variation of temperature and k is the thermal expansion coefficient. The second term on the right hand side of the above equation denotes the strain induced by thermal stress. The above equation can be reformulated as The boundary integral equation for evaluating the displacement of a pointp is expressed as [56,60,61]: in whichp is the source point, Q the field point on the boundary, and A = 1 − 2υ. c ij = δ ij for inner point and c ij = δ ij /2 if the point is located on the boundary.
The boundary integral equation for evaluating stress at pointp is given by [56,60,61], where n is the normal vector, and and We can find the domain integral term exists in Eqs. (15) and (19). It will be treated with radial integration method [58].
Using the radial basis operation, any integrand in a domain can be transformed into a integral over its boundary where f(x) is a integrand function, rðQ;pÞ denotes the distance from the source pointp to the boundary point Q, and F(x) is a radial integral function as follows As shown in Fig. 3, the Cartesian coordinate x needs to be expressed as a function of the integral variable r: where Upon substitution of the domain integral term ∫ Ω Ψ i T dΩ into the radial integral formula, the domain integral is converted into a boundary integral: Z where By substituting Eq. (28) into the Eq. (15), we arrive at the displacement boundary integral equation, In a similar manner, the domain integral for stress in Eq. (19) can be expressed as: Substituting Eq. (31) into Eq. (19) can eliminate the domain integral term in the boundary integral formulation for stress evaluation as follows

IGABEM Formulation for Thermoelasticity Analysis
In this work, the NURBS basis function is used for both geometry generation and physical field approximation. The boundary of the domain is first discretized into N e non-overlapping NURBS elements. The global coordinates of any point in every NURBS element can be obtained through linear combination of the NURBS basis functions and the coordinates of the control points x l xðnÞ ¼ where l is the number of control points, and ξ is the local parametric coordinate varying in the range [−1, 1].
In order to perform Gaussian numerical integration, variables need to be mapped from physical space to local coordinate space. Therefore, the Jacobian transformation J ðnÞ consists of two parts: the mapping from physical space to parametric space and that from parametric space to local coordinate space. The corresponding mathematical expression is where ξ c and ξ d represent the papametric coordinates of the first and last point in a NURBS element, respectively. Thus, the Jacobian can be expressed as The approximation of the physical field with NURBS basis functions can be expressed by where u l j and t l j are nodal coefficients related to displacement and traction. Different from the Lagrangian interpolation used in traditional boundary element methods (TBEM), NURBS do not possess Kroneckerdelta property, so the control points may not lie on the boundary and u l j and t l j have no physical interpretations. We use collocation scheme to generate a series of equations, and adopt Greville abscissae [62,63] to obtain the coordinates of these collocation points, where n 0 a is the parameteric coordinate of the a-th collocation point. Substituting Eqs. (38) and (39) to boundary integral equation, the first term of Eq. (30) is expressed as follows where ē represents the element containing the collocation pointp 0 , andn 0 represents the local coordinate of the collocation point.
Similarly, the first term of the right hand side of Eq. (30) is discretized as where u le j and t le j represent nodal coefficients in element e related to displacement and traction. For the domain integral in the displacement boundary integral equation, the physical coordinates are approximated by NURBS. Given a known temperature field T, the radial integral can be evaluated using Eq. (26). For the simple temperature field distribution the radial integral can be calculated analytically, but in the general cases it needs to be calculated using the Gaussian numerical integration. The conversion relationship of the integral variable is Since the domain integral has been approximated by the radial integration method, only the geometric field approximation is required. Using Eqs. (26) and (43) or collected as a matrix form: where the matrix H contains all integrals kernel functions u Ã ij and the jump terms c ij , the matrix G contains all integration kernel functions t Ã ij , and the vectors u and t contain the nodal coefficient of displacements and traction field. y p is a column vector formed by the domain integral caused by thermal stress (Eq. (44)). After imposing boundary conditions, Eq. (47) transforms to a set of linear equations as with x being a vector of unknown degrees of freedom.

Thick-Walled Cylinder Subjected to Internal and External Temperature Difference
A cylinder structure model with an inner radius of 1 m and an outer radius of 2 m is heated from the inner boundary surface, as shown in Fig. 4. This stress is induced by a large temperature gradient between the inner and outer surfaces. The temperature of the inner boundary rises by 1°C, and that of the outter boundary rises by 0°C. Linear expansion coefficient k = 1.2 × 10 −5 K −1 , modulus of elasticity E = 210 GPa, and Poisson's ratio υ = 0.3. This problem can be simplified to be a plane strain problem. Because the structure and temperature loads are axisymmetric, only a quarter of the structural model is needed for numerical analysis. This model is constructed by NURBS curve, as shown in Fig. 5.
From Fig. 5, we can find that only a few control points are needed to describe the geometric model exactly when using NURBS curves. However, for approximation of the physical field, they are not enough to capture the unknown variables. It may cause a large error for the calculation results. One of the most favourable properties of IGA is that the mesh can be refined while retaining accurate geometry at all times. Fig. 6 depicts the NURBS curve and control points with different levels of subdivision.
It can be seen that through h-refiment operation, more control points are obtained and the consistency of the structural model is still maintained.
The heat transfer equation is applied to obtain the temperature field distribution on the model satisfying initial conditions and boundary conditions, as follows where " r represents the radial distance from a point to the center of the circle.
The analytical solution of the circumferential thermal stress distribution on the structure is expressed as: A model with 200 degrees of freedom is used for numerical analysis. The distribution function of temperature field determines the values of the final displacement and stress. For a simple temperature  field distribution function, the corresponding radial integral function in Eqs. (29) and (32) can be obtained analytically, but for complex situations, Gaussian quadrature should be used.
Tab. 1 lists the numerical and analytical solutions of the hoop stresses at several internal points with different radial distances. It can be seen that the calculation results for the algorithm proposed in this work are consistent with the analytical solution.
Among them, the minimum error at the point with r = 1.5 m is −9.0766E-04, and the maximum error at the point with r = 1.96 m is −2.3843E-02.
The stress values at more points are calculated by IGABEM and compared with the analytical solution to further demonstrate the correctness of the algorithm, as shown in Fig. 7. Because the thermal stress boundary integral equation contains singular integrals, the accuracy of the boundary element method for calculating internal thermal stress decreases as the calculation point approaches the boundary. Thus, It can be seen that, the result has a small jump around r = 2 m. In this example, for the same degrees of freedom, IGABEM achieves higher accuracy than that of traditional BEM, and the accuracy of IGABEM becomes even more significant when the computing point is close to the structural boundary. For example, the errors of traditional BEM at r = 1.96 m and 1.98 m are 3.5205E-02 and 3.0584E-02, respectively. As can be seen in the Tab. 1, the errors of IGABEM are −2.3843E-02 and 8.8156E-03, respectively. The high accuracy of IGABEM can be attributed to the high order basis functions and the geometric accuracy.
To further explore the accuracy and stability of the algorithm, convergence analysis is performed, and the model is further refined, including 400, 800 and 1000 degrees of freedom. Fig. 8 illustrates the convergence curves of hoop stress values at two interior points with different radial distances that are evaluated by IGABEM and traditional BEM using quadratic elements, respectively.
It can be observed that the numerical solution of IGABEM converges rapidly and has smaller errors than the quadratic BEM.
It is noteworthy that although IGABEM improves the calculation accuracy, the NURBS employed in IGABEM need to be evaluated recursively, which is far more time consuming than the polynomial basis functions used in conventional algorithms. However, this problem is alleviated by adopting Bézier  Figure 7: Hoop stress curve 596 extraction operation that transforms NURBS to Bernstein polynomials and thus avoids the recursive calculation and considerably improves the calculation efficiency.

Eight-Leaf Plate Model
In order to demonstrate the ability of IGABEM in dealing with complex models, a two-dimensional model called "eight-leaf plate" is constructed in this work, as shown in Fig. 9. In this example, the geometric model is constructed through in-house software which allows the users to input the control point coordinates in command lines or determine their positions directly with mouse in a graphical user interface. For more complicated geometries, we can resort to the professional CAD software like Rhinoceros and import the geometric information to our codes and then perform h-refinement if needed. Linear expansion coefficient k = 1.0 ×10 -5 K −1 , Poisson's ratio υ = 0.3. The temperature field distribution is taken as a linear function (T = 60y + 60) and a quadratic function (T = 40y 2 − 60y), respectively.
Because of the simple temperature field function, the radial integral function can be calculated analytically. In this example, for the quadratic temperature field function (T = 40y 2 − 60y), the radial integral function in the displacement integral Eq. (29) is derived as " Similarly, the radial integral function in the stress integral equation can be obtained as follows: We further investigate the influence of different temperature fields on the displacement and stress. The displacements and stresses at the internal points along the x-axis are calculated subject to the linear and quadratic temperature distribution, respectively. The numerical results with IGABEM are shown in Figs. 10 and 11 with 1000 degrees of freedom.
From Figs. 10 and 11, we can observe that at the internal points furthest from the boundary, the displacement attains the smallest values, and the stress under the linear function temperature field distribution is the largest. In contrary, under the quadratic temperature field distribution the stress reaches the maximum at the boundary points. It is also noticed that the calculation results of displacement and stress show symmetry to the y-axis because the structure and temperature field distribution are symmetrical along the y-axis. x/m Figure 11: Curve of stress at internal point

Wrench Model
In this section, a wrench model as shown in Fig. 12 is analyzed with IGABEM. Fig. 13(a) shows the NURBS curve and the associated control points. This model uses the same quadratic distribution function of temperature field and material properties as the previous example.
We select several internal points (Fig. 13(b)) and evaluate their displacements and stresses with IGABEM and TBEM (Tabs. 2 and 3). The TBEM uses quadratic elements. It can be seen that the numerical results of IGABEM agree well with that of TBEM.

Conclusion
This paper applied IGABEM to thermoelastic analysis. It takes full advantage of IGABEM in the seamless integration of CAD and numerical analysis. IGABEM significantly reduces the meshing time and improves the accuracy due to exact geometry representation. To address the domain integral term arising from the thermal stress, the radial integration method was used in IGABEM to successfully convert the domain integral to the boundary integral. In the case of a simple temperature field distribution function, the radial integral is evaluated analytically, while for a complex distribution function, a numerical integration is used. This method is also applicable to 3D thermal stress problems. The next step is to extend this method to thermoelastic analysis for complicated 3D geometries with multi-patches and the future work is to solve the thermal-structure coupling problems.
Funding Statement: This study was funded by the National Natural Science Foundation of China (NSFC) (Grant Nos. 11702238, 51904202 and 11902212) and Nanhu Scholars Program for Young Scholars of XYNU.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the present study.