Closed Form Integration of Singular and Hypersingular Integrals in 3 D BEM Formulations for Heat Conduction

The evaluation of the singular and hypersingular integrals that appear in three-dimensional boundary element formulations for heat diffusion, in the frequency domain, is presented in analytical form. This improves computational efficiency and accuracy. Numerical integrations using existing techniques based on standard Gaussian integration schemes that incorporate an enormous amount of sampling points are used to verify the solutions of singular integrals. For the hypersingular integrals the comparison is evaluated by making use of an analytical solution that is valid for circular domains, combined with a standard Gaussian integration scheme for the remaining boundary element domain. Closed form solutions for cylindrical inclusions with null temperatures and null heat fluxes prescribed on the boundary are then derived and used to validate the three-dimensional boundary element formulations.


Introduction
A number of models have been proposed to study heat diffusion.They range from the analytical methods 1 to purely numerical methods based on domain discretization, such as, finite elements FEMs 2 and finite differences FDMs 3, 4 or on boundary discretization, such as, the boundary element method BEM 5, 6 .More recently, some researchers are focused on the development of meshless methods which require neither domain nor boundary discretization, such as, the method of fundamental solutions MFSs 7 .
Of the numerical methods, the BEM is possibly one of the most suitable tools for modelling homogeneous unbounded or semi-infinite systems because it automatically satisfies the far field conditions, so that only the boundaries of the inclusions require discretization.The BEM thus allows a compact description of the regions, leading to fully diffusion in the vicinity of cylindrical inclusions embedded in a homogeneous medium, when heated by dynamic point sources.Analytical solutions have been derived for this problem.

3D Geometry
Consider an unbounded spatially uniform solid medium of density ρ, thermal conductivity λ, and specific heat c, with a three-dimensional embedded inclusion, bounded by a surface S see Figure 1 .
This system is subjected to a pressure point heat source placed at x s , y s , z s , f x, y, z, t δ x − x s δ y − y s δ z − z s e iωt , 2.1 where δ x − x s , δ y − y s , and δ z − z s are Dirac-delta functions, and ω is the frequency of the source.The response of this heat source can be expressed by in which K is the thermal diffusivity defined by λ/ρc, ω is the oscillating frequency, i √ −1,

Boundary Integral Formulations
This section describes the formulation of the BEM and the TBEM used to obtain the heat field generated when the inclusion is heated by the incident three-dimensional source.
The solutions are subsequently verified against a numerical standard Gaussian quadrature scheme.

Boundary Element Formulation (BEM)
The transient heat transfer by conduction to calculate the heat T at any point of the spatial 3D homogeneous solid domain is governed by the diffusion equation: x, y, z and k c −iω/K.The boundary integral equation, formulated in the frequency domain, can be constructed by applying the reciprocity theorem, leading to where G and H are, respectively, the fundamental solutions Green's functions for the temperature T and heat flux q , at a point x x, y, z on the boundary S, due to a virtual point heat source at x 0 x 0 , y 0 , z 0 ; n n1 represents the unit outward normal along the boundary S, at x x, y, z ; b is a constant defined by the shape of the boundary, taking the value 1/2 if x 0 x 0 , y 0 , z 0 ∈ S and 1 otherwise.The required Green's functions for temperature and heat flux in an unbounded medium, in Cartesian coordinates, are given by G x, x 0 , ω e −ik c r 4λπr ,

Null Heat Fluxes along Its Boundary
When null heat fluxes are prescribed along the boundary, 3.

Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the following equation:

The Normal-Derivative Integral Equation (TBEM)
The normal-derivative integral equation can be derived by applying the gradient operator to the boundary integral equation 3.2 , which can be seen as assuming the existence of dipole heat sources.When the boundary of the inclusion is loaded with dipoles dynamic doublets the required integral equations can be expressed as:

3.6
The Green's functions G and H are defined by applying the gradient operator to G and H, which can be seen as the derivatives of these former Green's functions, to obtain heat fluxes and temperatures.In these equations, n n2 is the unit outward normal to the boundary S at the collocation points x 0 x 0 , y 0 , z 0 , defined by the vector n n2 .In this equation, the factor a is null for piecewise planar boundary elements.
The required three-dimensional Green's functions for an unbounded space are now defined as: In 3.6 the incident heat field is computed as

Null Heat Fluxes along Its Boundary
Null heat fluxes are prescribed along the boundary, which leads to the following equation:

Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the equation: bq x 0 , n n1 , ω S q x, n n1 , ω G x, n n2 , x 0 , ω ds T inc x 0 , n n2 , x s , ω . 3.11

Analytical Integration of Singular and Hypersingular Integrals
The global solution is found by solving 3.2 or 3.6 , which requires the discretization of the interface S into N planar boundary elements, with one nodal point in the middle of each element.
The integrations in 3.2 and 3.6 are evaluated using a Gaussian quadrature scheme when the element to be integrated is not the loaded element.For the loaded element the singular element , however, the integrands exhibit a singularity and the integration can be carried out in closed form, as will be demonstrated.
Consider the singular rectangular element of width W in the x direction and length L in the z direction shown in Figure 2.

Singular Integrals in 3.2
Since in this case r is perpendicular to the normal e.g., rn n1 0 , the singular term −W/2 H x, n n1 , x 0 , ω dx dz disappears.On the other hand, the integration of the Green's function This integration is evaluated by first expressing G x, x 0 , ω as the sum of twodimensional Green's functions with varying spatial wavenumbers.This is accomplished by first applying a spatial Fourier transformation along the z direction to the three-dimensional Green's function G x, x 0 , ω .The application of a Fourier-transformation to 3.3 in that direction leads to a line heat field, whose amplitude varies sinusoidally in the third dimension z , where k z is the wavenumber in the z direction.
Assuming the presence of an infinite set of equally spaced sources in the z direction, the former Green's function can be recast as: where L vs is the spatial source interval, and k zm √ 2π/L vs m.This equation converges very fast and can be approximated by a finite sum of terms M .The distance L vs needs to be large enough to avoid spatial contamination.Given that the frequency increment defines the time window for the analysis, 1/Δf, the minimum distance needs to be at least 2.0 √ K 1/Δf.The number of terms depends on this distance L vs .However, a larger value of L vs corresponds to a higher number of terms M .
Note that the use of complex frequencies further reduces the influence of the neighboring fictitious sources.The 3D Green's field can therefore be computed as the pressure irradiated by a sum of harmonic steady-state line loads, whose amplitude varies sinusoidally in the z dimension.
This procedure allows the integration of −W/2 G x, x 0 , ω dx dz to be obtained in a similar manner, thus

3.16
where S ns • • • are Struve functions of order ns.

Hypersingular Integrals in 3.6
Since in this case r is perpendicular to the normal e.g., rn n2 0 , the singular term However, the integration of the Green's function , ω dx dz with n n2 n n1 leads to a hypersingular term.Using the procedure above, the evaluation of this integration can performed by writing H x, n n1 , n n2 , x 0 , ω as the sum of two-dimensional Green's functions with varying spatial wavenumbers.This can be accomplished by applying the traction operator to the Green's function: H 1 kr 0 r 0 .

3.18
This procedure allows the integration of where −W/2 ik/8πλ −kH 2 kr 0 ∂r 0 /∂x ∂x/∂n nl ∂r 0 /∂y ∂y/∂n nl 2 H 1 kr 0 /r 0 dx.The integration of I 2 kr 0 is performed indirectly by isolating a semicylinder just above the boundary element and by considering its thermal equilibrium see Figure 3 .The integration of I 2 kr 0 is performed indirectly by isolating a semicylinder just above the boundary

3.21
Both integrands to the right of the equal sign in 3.21 are well behaved and can be evaluated directly,

Verification of the Analytical Integrations
The above analytical expressions are next implemented and compared with the numerical integration by means of a standard Gaussian quadrature, using a large amount of sampling Mathematical Problems in Engineering points in the case of singular integrals and making use of analytical expressions deduced for circular domains in the case of hypersingular integrals.

Verification of the Singular Integral Integration
The integrations To illustrate the convergence of 3.14 , Figure 5 presents the imaginary and real parts of the result obtained as the number of terms is successively taken into account in the integration k zm √ 2π/L vs m, m 1, . . ., M , that is, as k z increases, when the integration is performed for a frequency of 0.5 × 10 −7 Hz.It can be observed that the integration converges very fast as k z increases.Similar behavior was found throughout the frequency range.
For the purpose of verification, the analytical results for the frequency range 0.0, 1 × 10 −5 Hz are illustrated in Figure 6 by solid lines.The numerical integration results, given  by a standard Gaussian quadrature with 2000 × 2000 points, are plotted by marked lines in this figure.The two solutions exhibit similar results.Even so, the accuracy of the result of the numerical integration depends on the number of sampling points.The variation of the error, as the number of Gauss points increases from 500 × 500 to 2000 × 2000 sampling points, is given in Figure 7.As expected, the accuracy of the numerical integration improves as the number of sampling points increases.

Verification of the Hypersingular Integral Integration
To illustrate the convergence of 3.19 , Figure 8 shows the integration result as k z increases, that is, as the number is successively taken into account, when the integration is performed for a frequency × 10 −7 Hz.The were computed for frequencies c ω − iη with η 0.7Δω, Δω 2 × π × Δf, and Δf 0.5 × 10 −7 Hz .
Figures 8 a real part and 8 b imaginary part give the analytical results for the integration performed for the full quadrangular region, 0.2 m wide.As k z increases the imaginary part converges very fast, while the real part oscillates between two fixed values.It can be seen that the mean value between these two values is the correct integration result indicated by the horizontal line in the plot a .Similar behavior was found for all the cases studied.
A reference solution is required for the verification of the proposed analytical solution.It happens that the integrations −W/2 H x, n n1 , n n2 , x 0 , ω dx dz cannot be computed numerically using a standard Gaussian quadrature scheme because the solutions do not converge, even for a large amount of sampling points.This problem is overcome by comparing the proposed analytical integration with a reference integration solution constructed by splitting the former domain into two subdomains.The first assumes a circular part circumscribed by the full domain, centered on the position of the virtual load, and the second is the rest of the boundary element domain see Figure 9 .The integration over the outer part uses a standard Gaussian quadrature scheme, while the integration over the circular part is obtained analytically, following the procedure described above.
Thus, the integration of H x, n n1 , n n2 , x 0 , ω over the circular domain is performed indirectly by isolating a semisphere just above the circular domain and by considering its thermal equilibrium see Figure 10 .
As before, both integrands to the right of the equals sign in 4.1 are well behaved, and can be evaluated directly, cos θ sin θ dθ dφ  This equation resembles the one derived for acoustics domains by Terai 29 to compute the exterior sound fields around scattering objects and presents techniques for treatment of singularities.
The analytical integrations were verified by computing the solution of this integral over the planar quadrangular boundary element defined above.
The analytical integration uses 4.3 , which is derived from Terai's work, and it assumes that the domain is delimited by the circular domain, with a radius of 0.1 m.The numerical integration along the outer part of the boundary element domain was performed with a zero function value assigned sampling points inside the central circular region 2000 Gaussian quadrature sampling points were used .
The computations are obtained in the frequency range 0.0, 1 10 −5 Hz with a frequency increment of 0.5 × 10 −7 Hz.The thermal diffusivity and L vs are assumed to be, respectively, 4.96278 × 10 −7 m 2 s −1 and m.
The results obtained with the proposed analytical expressions as continuous lines and with the constructed reference solution marked lines are illustrated in Figure 11.Analysis of this figure confirms that the results are very similar.

Verification of the BEM Formulations
In this section, the proposed algorithms are validated using a circular cylindrical inclusion, with radius a, aligned along the z axis see Figure 12 , for which analytical solutions can be derived.A point heat source placed at x s , y s , z s is assumed to excite the medium.Two different conditions are prescribed on the boundary, that are, null heat fluxes and null temperatures.To enable comparison with the 3D BEM model, the length of the inclusion is limited by imposing null heat fluxes on sections z 0.0 m and z L t .
The analytical solution for this problem is obtained by first applying a spatial Fourier transformation in the z direction, which allows the solution to be obtained as the sum of two-dimensional solutions with a varying spatial wavenumber in that direction.The null normal heat fluxes at sections z 0.0 m and z L t are accomplished by adding the temperature field generated by the real source to that produced by virtual sources image sources , which are placed in the z direction and act as mirrors of the real source, so as to ensure the required boundary conditions.
Consider the following: with k zm √ 2π/L vs m m 1, . . ., NS and where The number of virtual sources NS z to be used in the calculations is defined so that the signal responses can be correctly computed in the time frame, which is determined by the frequency increment 1/Δf.This procedure does not introduce any type of error into the computed time impulse response within the time window defined.Notice that L vs should be at least twice the distance between the real source and the farthest virtual source.Each two-dimensional problem is solved using the separation of variables procedure with the Helmholtz equation and enforcing the boundary conditions throughout the boundary surface, using the Bessel series form.The following equations can be derived if the origin of the coordinate system coincides with the center of the circle, cross section of the cylinder, and the source that lies on the x axis.

Null Temperatures along Its Boundary
Consider the following:  Equations 5.3 and 5.4 converge very fast 1 .The number of terms is defined using a small criterion of convergence applied to the difference between successive values.In the present formulation the convergence is verified by calculating the difference between the amplitude of the response for n N 20 and for n N, using the response at the receiver placed closest to the boundary as reference.
The temperature responses are computed at two grids of receivers placed as illustrated in Figure 13 grid 1 y 0.0 m , grid 2 z 1.0 m , in the vicinity of a cylinder inclusion with a radius of a 0.5 m and 2 m in length.The thermal properties are that of the mortar, mentioned above.The 3D point source is placed at −1.5 m, 0.0 m, 1.0 m .Temperature computations are performed at frequencies 0.0 Hz and 0.5 × 10 −7 Hz.The responses were computed for complex frequencies ω c ω − iη with η 0.7Δω, Δω 2 × π × Δf and Δf 0.5 × 10 −7 Hz .
The responses have been computed both analytically and using the BEM and TBEM formulations.The behavior of the formulations was verified by using two different numbers of boundary elements to discretize the inclusion: 30 × 20 20 in the z direction and 50 × 32 32 in the z direction .
Figures 14 and 15 illustrate the responses obtained when null temperatures are imposed on the boundary surface of the cylinder, on receiver grids 1 and 2, respectively.Each figure shows the real and imaginary parts of the analytical response and the error obtained when the system is solved using the BEM.As expected, the BEM accuracy improves as the receiver is placed further away from the inclusion.In both cases it can be seen that the magnitude of the error increases as the receivers approach the inclusion's boundary.It can further be seen that the solution improves as the number of boundary elements increases, which illustrates the good accuracy of the BEM response.
Figures 16 and 17 give the response real and imaginary parts at grid 1 of receivers when null heat fluxes are prescribed along the boundary surface of the cylinder.Again, each figure illustrates the analytical response and the numerical error.Figures 16 and 17 show the errors when the problem is solved using the BEM and TBEM, respectively.
As in the previous case, the receivers placed further away from boundary exhibit lower amplitude errors.Once more, the solution is seen to improve for higher numbers of boundary elements, which illustrates the good accuracy of the two models.The errors given by the TBEM are generally higher than those given by the BEM model.When 30 × 20 boundary elements are used, the highest value of the error in grid 1 is about 7 × 10 −3• C. The responses observed for grid 2 exhibit similar features not illustrated .

Conclusions
In this paper the authors presented the analytical integration of singular and hypersingular integrals that appear while modelling three-dimensional heat transfer by conduction using the classical BEM and the normal-derivative integral equation TBEM formulations when the element to be integrated is the loaded one singular element .
The analytical integrations were first compared with numerical solutions using a standard Gaussian quadrature scheme, in the case of singular integrals.The proposed analytical integrations for hypersingular integrals were then verified against a reference solution constructed by combining two integrals: an analytical integration along a circular domain, which handles the hypersingular part of the integral and a numerical integration using a standard Gaussian quadrature scheme applied to the rest of the boundary element domain.
Then the three-dimensional BEM and TBEM solutions, incorporating the analytical solutions, were verified against analytical solutions derived for a cylindrical circular inclusion limited by two perpendicular sections, where null heat fluxes where imposed: the responses showed very good accuracy.

Figure 2 :
Figure 2: Scheme of a planar boundary element.

Figure 3 :
Figure 3: Dynamic equilibrium on the singular element.
x 0 , ω dx dz are computed along a planar quadrangular boundary element 0.2 m long see Figure 4 .Computations are performed with a frequency increment of 0.5 × 10 −7 Hz in the frequency range 0.0, 1 × 10 −5 Hz .The thermal properties assigned to the medium are those of the mortar: specific heat c 780 J•kg −1 • • C −1 , density 1860 kg•m −3 and thermal conductivity 0.72 W•m −1 • • C −1 thermal diffusivity of 4.96278 × 10 −7 m 2 s −1 .L vs is assumed to be 60.0 m.

FrequencyFigure 6 :Frequency
Figure 6: BEM: Analytical results as continuous lines and numerical results using a standard Gaussian quadrature with 2000 × 2000 sampling points.

Figure 7 :
Figure 7: Variation of BEM error different numbers of Gauss sampling points.

Figure 8 :Figure
Figure 8: Convergence of the proposed analytical integration solution: a real part; b imaginary part.

Figure 10 :
Figure 10: Dynamic equilibrium on a circular hypersingular domain.

Figure 14 :
Figure 14: Comparative analysis of analytical and BEM responses real and imaginary parts at grid 1, when null temperatures are prescribed, considering different numbers of boundary elements: a analytical response; b error using BEM with 30 × 20 boundary elements ; c error using BEM with 50 × 32 boundary elements .

Figure 15 :
Figure 15: Comparative analysis of analytical and BEM responses real and imaginary parts at grid 2, when null temperatures are prescribed, considering different numbers of boundary elements: a analytical response; b error using BEM with 30 × 20 boundary elements ; c absolute error using BEM with 50 × 32 boundary elements .

Figure 16 :
Figure 16: Comparative analysis of analytical and BEM responses real and imaginary parts at grid 1, when null heat fluxes are prescribed, considering different numbers of boundary elements: a analytical response; b error using BEM with 30 × 20 boundary elements ; c error using BEM with 50 × 32 boundary elements .

Figure 17 :
Figure 17: Comparative analysis of analytical and TBEM responses real and imaginary parts at grid 1, when null heat fluxes are prescribed, considering different numbers of boundary elements: a analytical response; b error using TBEM with 30 × 20 boundary elements ; c error using TBEM with 50 × 32 boundary elements .
00 |x 0 | and θ arctan y/x .n ε n H n kr 00 n/a J n ka − kJ n 1 ka / n/a H n ka − kH n 1 ka .