Numerical Computation of 2D Domain Integrals in Boundary Element Method by ( α , β ) Distance Transformation for Transient Heat Conduction Problems

: When the time-dependent boundary element method, also termed the pseudo-initial condition method, is employed for solving transient heat conduction problems, the numerical evaluation of domain integrals is necessitated. Consequently, the accurate calculation of the domain integrals is of crucial importance for analyzing transient heat conduction. However, as the time step decreases progressively and approaches zero, the integrand of the domain integrals is close to singular, resulting in large errors when employing standard Gaussian quadrature directly. To solve the problem and further improve the calculation accuracy of the domain integrals, an ( α , β ) distance transformation is presented. Distance transformation is a simple and efficient method for eliminating near-singularity, typically applied to nearly singular integrals. Firstly, the ( α , β ) coordinate transformation is introduced. Then, a new distance transformation for the domain integrals is constructed by replacing the shortest distance with the time step. With the new method, the integrand of the domain integrals is substantially smoothed, and the singularity arising from small time steps in the domain integrals is effectively eliminated. Thus, more accurate results can be obtained by the ( α , β ) distance transformation. Different sizes of time steps, positions of source point, and shapes of integration elements are considered in numerical examples. Comparative studies of the numerical results for the domain integrals using various methods demonstrate that higher accuracy and efficiency are achieved by the proposed method.


Introduction
Numerical methods are extensively employed across various engineering applications [1][2][3].In many practical engineering scenarios, the problems are transient, which means that they are space-and time-dependent.Since it is usually very difficult to find an analytical solution, numerical methods are widely employed for these problems, such as the finite element method (FEM) [4], the finite difference method (FDM) [5][6][7], and the boundary element method (BEM) [8][9][10].Compared to other methods, the BEM is identified as being a significantly simpler and more efficient tool for solving transient heat conduction problems.Yu et al. [11][12][13] introduced the isogeometric dual reciprocity boundary element method to solve the problems of transient heat conduction.Erhart et al. [14] adopted a parallel domain decomposition boundary element method approach to solve large-scale transient heat conduction problems.Feng et al. [15,16] proposed an analytically integrated radial integration BEM and a meshless interface integral BEM for solving transient heat conduction.Zhang et al. [17,18] developed the boundary face method for analyzing transient heat conduction.
In the analysis of transient heat conduction, the pseudo-initial method, based on the time-dependent fundamental solution, is prevalently employed within the BEM.The method employs the temperature calculated in the previous step as the initial condition for the current step, necessitating careful consideration of the domain integral associated with this pseudo-initial condition.Although the domain integrals typically exhibit regularity, the integrand (time-dependent fundamental solution) of these integrals may vary dramatically near the source point as the time step approaches zero, with the domain integrals exhibiting properties similar to singular integrals.The numerical integration of functions possessing the above characteristics leads to numerical instability in the BEM solution, as evidenced in references [19][20][21].Such conditions render the Gaussian quadrature method inadequate for achieving accurate calculations.Therefore, eliminating the singularity caused by the time step is essential for accurately computing domain integrals, which is crucial for solving transient heat conduction problems using the pseudo-initial condition method.
Various strategies have been developed to counteract the adverse effects arising from the use of small time steps.Dong et al. [22,23] introduced a method combining a coordinate transformation with an element subdivision technique.However, this technique is timeconsuming and requires more complex coding due to the need to divide the integration element into enough sub-elements.Gao et al. [24][25][26] proposed the radial integration method, which converted the domain integrals into equivalent boundary integrals, but its efficiency is compromised in large-scale problem applications.In recent work, Dong et al. expanded upon the sinh transformation [27] and the distance transformation [28], originally applied to the nearly singular integrals [29][30][31][32][33][34], to calculate the domain integrals.These transformations are implemented based on the polar coordinate transformation.
To improve calculation accuracy and simplify the implementation process, an (α, β) distance transformation is proposed to compute the 2D domain integrals in the BEM for transient heat conduction problems.The (α, β) coordinate transformation is initially introduced to smooth the integrand to some degree.The shortest distance in the nearly singular integral is replaced by the time step, and the new distance transformation is constructed in the α direction.In this paper, a new distance transformation is represented and combined with the (α, β) coordinate transformation.The singularity arising from small time steps in the domain integral is effectively eliminated by the Jacobian of the (α, β) distance transformation.Comparative analyses of the numerical examples demonstrated that the proposed method is simpler and more effective than the other four existing methods.Thus, the new method can be effectively applied to accurately compute the 2D domain integrals.
The structure of the paper is organized as follows.The boundary integral equation, the domain integral, and the element subdivision technique are discussed in Section 2. The implementation of the proposed method is detailed in Section 3. Several numerical examples are given to demonstrate the validity of the method in Section 4. The conclusion is provided in Section 5.

Boundary Integral Equation
Assuming the material is isotropic, the governing equation for transient heat conduction problems without heat sources is initially given as follows: where u (X, t) denotes the temperature in the coordinate X at the time t; and k is the diffusion coefficient.The boundary integral equation [2,3] for transient heat conduction in an isotropic, homogeneous medium Ω bounded by Γ can be expressed as follows: Axioms 2024, 13, 490 where Y and X denote the source point and the field point, respectively.c(Y) is an internal angle function that varies according to the geometrical shape of the boundary at Y. u 0 (X, t 0 ) is the initial temperature.t 0 and t F represent the initial time and the end time of the time step, respectively.u* and q* denote the time-dependent fundamental solution and the normal derivative of u*, respectively.
The time-dependent fundamental solution u* [2,3] can be expressed as follows: where τ = t F − t 0 , which is the time step.The distance between the source Y and the field point X is denoted by r.

Domain Integral
The domain integral [22,23,28,29] of the boundary integral equation is obtained as follows: where u 0 (X, t 0 ) is a regular function that can be calculated exactly by standard Gaussian quadrature.The influence of the time step size on the integral is evident from (4).
As the time step τ approaches zero, the integrand tends to infinity, which consequently undermines the accuracy of calculating the integral I using standard Gaussian quadrature [22,23,29].The variation of u* with respect to r for different time steps is depicted in Figure 1.
Axioms 2024, 13, x FOR PEER REVIEW 3 of 14 The boundary integral equation [2,3] for transient heat conduction in an isotropic, homogeneous medium Ω bounded by Γ can be expressed as follows: where Y and X denote the source point and the field point, respectively.c(Y) is an internal angle function that varies according to the geometrical shape of the boundary at Y. u0(X, t0) is the initial temperature.t0 and tF represent the initial time and the end time of the time step, respectively.u* and q* denote the time-dependent fundamental solution and the normal derivative of u*, respectively.The time-dependent fundamental solution u* [2,3] can be expressed as follows: where τ = tF − t0, which is the time step.The distance between the source Y and the field point X is denoted by r.

Domain Integral
The domain integral [22,23,28,29] of the boundary integral equation is obtained as follows: where u0(X, t0) is a regular function that can be calculated exactly by standard Gaussian quadrature.The influence of the time step size on the integral is evident from (4).As the time step τ approaches zero, the integrand tends to infinity, which consequently undermines the accuracy of calculating the integral I using standard Gaussian quadrature [22,23,29].The variation of u* with respect to r for different time steps is depicted in Figure 1.

The Subdivision of the Integration Element
To improve the calculation accuracy of the domain integral, the subdivision of the quadrilateral integration element is employed prior to applying the proposed method, as depicted in Figure 2. Based on the position of the source point (x 0 , y 0 ) on the quadrilateral, the integration element is subdivided into different numbers of sub-triangles by connecting the source to the vertices of the element, as illustrated in Figure 2a-c.

The Subdivision of the Integration Element
To improve the calculation accuracy of the domain integral, the subdivision of the quadrilateral integration element is employed prior to applying the proposed method, as depicted in Figure 2. Based on the position of the source point (x0, y0) on the quadrilateral, the integration element is subdivided into different numbers of sub-triangles by connecting the source to the vertices of the element, as illustrated in Figure 2a-c

(α, β) Distance Transformation
In this section, the (α, β) coordinate transformation is first introduced.Then, to eliminate the adverse effect arising from the small time step in the calculation of the domain integrals, a new distance transformation is constructed.

The (α, β) Coordinate Transformation
In this section, the (α, β) transformation proposed by Zhang et al. [35], akin to the polar coordinate transformation, is initially introduced.As illustrated in Figure 3, the subtriangular element is mapped to a unit quadrilateral using the (α, β) transformation, wherein the source point (x0, y0) corresponds to the line where β = 1.To implement the mapping depicted in Figure 3, the following expressions can be employed: y y y y (5)

(α, β) Distance Transformation
In this section, the (α, β) coordinate transformation is first introduced.Then, to eliminate the adverse effect arising from the small time step in the calculation of the domain integrals, a new distance transformation is constructed.

The (α, β) Coordinate Transformation
In this section, the (α, β) transformation proposed by Zhang et al. [35], akin to the polar coordinate transformation, is initially introduced.As illustrated in Figure 3, the sub-triangular element is mapped to a unit quadrilateral using the (α, β) transformation, wherein the source point (x 0 , y 0 ) corresponds to the line where β = 1.

The Subdivision of the Integration Element
To improve the calculation accuracy of the domain integral, the subdivision of the quadrilateral integration element is employed prior to applying the proposed method, as depicted in Figure 2. Based on the position of the source point (x0, y0) on the quadrilateral, the integration element is subdivided into different numbers of sub-triangles by connecting the source to the vertices of the element, as illustrated in Figure 2a-c

(α, β) Distance Transformation
In this section, the (α, β) coordinate transformation is first introduced.Then, to eliminate the adverse effect arising from the small time step in the calculation of the domain integrals, a new distance transformation is constructed.

The (α, β) Coordinate Transformation
In this section, the (α, β) transformation proposed by Zhang et al. [35], akin to the polar coordinate transformation, is initially introduced.As illustrated in Figure 3, the subtriangular element is mapped to a unit quadrilateral using the (α, β) transformation, wherein the source point (x0, y0) corresponds to the line where β = 1.To implement the mapping depicted in Figure 3, the following expressions can be employed: To implement the mapping depicted in Figure 3, the following expressions can be employed: where α and β are the two variables representing the proportions.Combining ( 5)-( 7), the coordinate of the arbitrary point (x, y) on the integration element can be obtained.
Given that the initial temperature in (4), u 0 (X, t 0 ), is a regular function, for simplicity, it is assumed to be 1 in this paper.Consequently, under the influence of the (α, β) transformation, ( 4) is written as follows: where the term 2αS is the Jacobian in the transformation from the (x, y) coordinate system to the (α, β) coordinate system and S is the area of the sub-triangular element.According to (9), the variables α and β in the (α, β) coordinate system are constrained to the interval [0, 1].This constraint simplifies and facilitates the integral calculation.Additionally, the Jacobian of the transformation can enhance the smoothness of the integrand, thereby improving the accuracy of the domain integral.

New Distance Transformation
As shown in Figure 4, the distance between the source point and the field point is obtained by applying the Taylor expansion near the source point and the (α, β) transformation (Equation ( 8)).
where x k,k=1,2 and y k , k=1,2 represent the nodal coordinates of the field point and the source point, respectively.(ξ, η) and (ξ 0 , η 0 ) denote the parametric coordinates of the field point X and the source point The coordinates (ξ 1 , η 1 ) and (ξ 2 , η 2 ) correspond to the other vertices of the sub-triangle in the parametric space.When the time step involved in the integrand is very small, the numerical calculation of the domain integral is negatively impacted.To eliminate these effects, a new distance transformation is proposed that utilizes the time step as the shortest distance in the nearly singular integral.
( ) The square of the distance between the source point and the field point can be expressed as follows: where When the time step involved in the integrand is very small, the numerical calculation of the domain integral is negatively impacted.To eliminate these effects, a new distance transformation is proposed that utilizes the time step as the shortest distance in the nearly singular integral.
From ( 12), the Jacobian of the distance transformation can be obtained as follows: Based on ( 12) and ( 13), the integral in (9) can be transformed into the following: In Equation (12), 2 is provided.From ( 14), it can be observed that the term exp(2v) includes the time step τ, which can effectively eliminate the time step in the denominator of the integrand.As a result, the adverse effects of the small time step on the integral calculation are significantly eliminated.

Numerical Examples
In this section, numerical examples are presented to validate the proposed method.The following domain integral is considered: where the diffusion coefficient k is assumed to be 1 in all examples.
The calculation accuracy and stability of all methods are compared using the following definition of relative error: where I numerical and I exact are the numerical solutions obtained by various methods of the considered integral and the exact solution, respectively.The exact solution is obtained by the element subdivision method [36] with enough sub-elements, which is considered stable and inefficient.

Example 1
With different time step sizes and Gaussian points, the impact on Gaussian quadrature is considered in this example.The integration element, as shown in Figure 5, with vertex coordinates (1, 1), (0, 1), (0, 0), and (1, 0), is used.The relative errors for the domain integral in (15) are illustrated in Table 1.τ is the size of the time step.The equation 5 × 5 denotes the Gaussian points used in Gaussian quadrature.
ble and inefficient.

Example 1
With different time step sizes and Gaussian points, the impact on Gaussian quadrature is considered in this example.The integration element, as shown in Figure 5, with vertex coordinates (1, 1), (0, 1), (0, 0), and (1, 0), is used.The relative errors for the domain integral in (15) are illustrated in Table 1.τ is the size of the time step.The equation 5 × 5 denotes the Gaussian points used in Gaussian quadrature.
As shown in Table 1, even with the smaller number of Gaussian points, employing direct Gaussian quadrature can yield accurate results when the time step is large.However, when the time step is progressively decreased to less than 0.001, despite employing a large number of Gaussian points, the Gaussian quadrature method proves inadequate for accurate results.As shown in Table 1, even with the smaller number of Gaussian points, employing direct Gaussian quadrature can yield accurate results when the time step is large.However, when the time step is progressively decreased to less than 0.001, despite employing a large number of Gaussian points, the Gaussian quadrature method proves inadequate for accurate results.

Example 2
This example examines the situation where the source point is located at the center of the quadrilateral element, illustrating the impact of the time step on various methods.The integration element in Figure 5 is used.Relative errors for the domain integral in (15), calculated using various methods across different time step sizes, are listed in Table 2. 'Direct' denotes the direct Gaussian quadrature method.'Polar' is the polar coordinate transformation.'Sinh' and 'Distance' represent the sinh transformation and the distance transformation based on the (ρ, θ) coordinate system, respectively.The (α, β) distance transformation, our suggested method, is indicated by 'αβDistance'.The 20 × 20 Gaussian points are used in all methods.Table 2 demonstrates that all methods maintain high accuracy when the values of τ are large; however, significant errors are incurred with the use of direct Gaussian quadrature as the time step decreases.Superior accuracy and stability are consistently exhibited by the proposed method when compared to other methods.

Example 3
In this example, the (α, β) distance transformation is compared with the element subdivision technique incorporating the (α, β) coordinate transformation as described in reference [22].The integration element, depicted in Figure 5, is utilized with the coordinate of the source point set to (0.5, 0.5).The relative errors of the two methods used for the domain integral calculation in Equation ( 15), as well as the computation time for 10000 iterations, are presented in Tables 3 and 4, respectively.'(α, β)' denotes the technique from reference [22] that is utilized.In this method, 20 × 20 Gaussian integration points are used for the subdivided sub-triangles, and 10 × 10 Gaussian integration points are employed for the subdivided sub-quadrilaterals.'αβDistance' indicates the (α, β) distance transformation, our suggested method, where 20 × 20 Gaussian integration points are utilized for each of the four sub-triangles.Consequently, a smaller number of Gaussian points are employed by the proposed method.15) using the method in reference [22] and the (α, β) distance transformation.15) using the technique in reference [22] and the (α, β) distance transformation running 10,000 times.From the comparison in Tables 3 and 4, the following conclusions can be drawn: The proposed method achieves similar accuracy to that of the method in reference [22] while requiring fewer Gaussian integration points.Additionally, significantly less time is consumed by the new method compared to the '(α, β)' method for calculating the domain integral.The above results derived from the element subdivision method are attributed to the requirement of dividing the integration element into enough sub-elements, a process that is time-consuming and necessitates more complex code.However, the proposed method needs only the subdivision of the integration element by connecting the source point with the vertices of the element and employing the new distance transformation to eliminate the singularity in the domain integral.Thus, it can be concluded that the (α, β) distance transformation demonstrates higher efficiency.

Example 4
The effect of different source point locations on numerical calculations is considered in this example.To demonstrate the superior results achievable with the proposed method at different source points, particular attention is given to situations where the source point coincides with a vertex, is located on an edge of the integration element, or is close to a vertex or an edge, as depicted in Figure 6a-d.The results calculated by the various methods for the domain integral are presented in Table 5.In Figure 7a-d in this example.To demonstrate the superior results achievable with the proposed method at different source points, particular attention is given to situations where the source point coincides with a vertex, is located on an edge of the integration element, or is close to a vertex or an edge, as depicted in Figure 6a-d.The results calculated by the various methods for the domain integral are presented in Table 5.In Figure 7a-d As shown in Table 5 and Figure 7, the sinh transformation and the distance transformation, which are based on the (ρ, θ) coordinate system and the (α, β) coordinate system, respectively, significantly outperform the polar coordinate transformation in terms of accuracy.Moreover, the distance transformation yields superior results compared to the sinh transformation.As the source point approaches the vertex or edge of the quadrilateral element, progressively larger calculation errors are obtained.The issue arises from the formation of irregular sub-triangles during the subdivision of integration elements.Among all four methods, the (α, β) distance transformation exhibits optimal accuracy and stability.

Example 5
In this example, the calculation efficiency of different methods for computing the domain integral in ( 15) is considered.The integration element, source point positions, and other conditions are kept consistent with those in Example 4. The computation times for (15) using various methods running 10,000 times are exhibited in Table 6.As shown in Table 5 and Figure 7, the sinh transformation and the distance transformation, which are based on the (ρ, θ) coordinate system and the (α, β) coordinate system, respectively, significantly outperform the polar coordinate transformation in terms of accuracy.Moreover, the distance transformation yields superior results compared to the sinh transformation.As the source point approaches the vertex or edge of the quadrilateral element, progressively larger calculation errors are obtained.The issue arises from the formation of irregular sub-triangles during the subdivision of integration elements.Among all four methods, the (α, β) distance transformation exhibits optimal accuracy and stability.

Example 5
In this example, the calculation efficiency of different methods for computing the domain integral in (15) is considered.The integration element, source point positions, and other conditions are kept consistent with those in Example 4. The computation times for (15) using various methods running 10,000 times are exhibited in Table 6.As revealed by the analysis in Table 7, significant calculation errors are produced by both the Gaussian quadrature method and the polar coordinate transformation when the value of τ is exceedingly small, making these methods unsuitable for such a condition.The analysis indicates that calculation accuracy is significantly improved by the sinh transformation and the new distance transformation.However, superior performance is demonstrated by the new distance transformation and the (α, β) distance transformation.Moreover, the proposed method achieves optimal results compared to all other methods.

Example 7
To underscore the broad applicability of the proposed method, the impact of interpolating the shape function on the results is evaluated in this example.The domain integral is considered as follows: )dΩ (17) where N denotes the shape function of the curved quadrilateral element.The integral is performed using the curved quadrilateral depicted in Figure 8, with the coordinate of the source point set at (0.5, 0.5).The relative errors of various methods at different time steps are compared in Table 8. 'Direct' indicates that Gaussian quadrature is applied directly.'Polar' refers to the polar coordinate transformation.'Sinh' and 'Distance' present the sinh transformation and the new distance transformation.'αβDistance' denotes the (α, β) distance transformation.The 20 × 20 Gaussian points are used in all methods.(17) on the curved element are obtained using various methods at different time steps, with the source point located at (0.5, 0.5).Errors greater than 1 are indicated as E.  As can be observed from Table 8, the interpolation of the shape function into the domain integral slightly affects the calculation results.The calculation accuracy obtained using the proposed method is superior to that achieved with other methods.The results demonstrate that the (α, β) distance transformation achieves the highest level of accuracy, confirming its wide applicability.

Conclusions
An (α, β) distance transformation is proposed in this paper to compute 2D domain integrals that arise in the boundary element method for transient heat conduction problems.By employing the proposed method, the integrand is smoothed, and computational implementation is simplified.Furthermore, the singularity resulting from small time steps in the domain integral is effectively eliminated by the Jacobian of the (α, β) distance transformation.The numerical examples, which utilize different time step sizes, locations of source points, and element shapes, effectively validate the accuracy, stability, and broad applicability of the proposed method.The relative errors of the proposed method can be maintained at 10 −5 until τ decreases to 10 −5 .Compared to the other four existing methods, the most accurate results are achieved by the new method, regardless of the time step size, the source point positions, and the shapes of the integration elements.Furthermore, the shortest calculation time for the domain integrals is required.Therefore, the (α, β) distance transformation can be efficiently applied to compute 2D domain integrals.

Figure 1 .
Figure 1.Variation of u* with respect to r for different time steps τ.Figure 1. Variation of u* with respect to r for different time steps τ.

Figure 1 .
Figure 1.Variation of u* with respect to r for different time steps τ.Figure 1. Variation of u* with respect to r for different time steps τ.

Figure 2 .
Figure 2. Subdivision of the integration element when the source point lies (a) inside the quadrilateral, (b) on an edge, or (c) on a vertex, respectively.

Figure 2 .
Figure 2. Subdivision of the integration element when the source point lies (a) inside the quadrilateral, (b) on an edge, or (c) on a vertex, respectively.

Figure 2 .
Figure 2. Subdivision of the integration element when the source point lies (a) inside the quadrilateral, (b) on an edge, or (c) on a vertex, respectively.

Figure 4 .
Figure 4.The distance between the source point and the field point in local space.

Figure 4 .
Figure 4.The distance between the source point and the field point in local space.

Figure 5 .
Figure 5. Nodal coordinates of the quadrilateral element.

Table 1 .
Relative errors of Gaussian quadrature for the domain integral in(15), with different numbers of Gaussian points.Errors greater than 1 are indicated as E.
, the results of various methods at different source point positions are shown.'Polar' denotes the polar coordinate transformation.'Sinh' and 'Distance' refer to the sinh transformation and the new distance transformation, respectively.'αβDistance' indicates the (α, β) distance method.The 20 × 20 Gaussian points are used in all methods.

Table 1 .
Relative Figure 5. Nodal coordinates of the quadrilateral element.errors of Gaussian quadrature for the domain integral in (15), with different numbers of Gaussian points.Errors greater than 1 are indicated as E.

Table 3 .
Relative errors for the domain integral in Equation (

Table 4 .
Computation time for the domain integral in Equation (