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

Authors: A. Tadeu; J. Prata; N. Simões
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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 populated systems of equations, contrary to the sparse systems given by the FDM and FEM techniques.  One disadvantage of the BEM is that it can only be applied to more general geometries and media when the relevant fundamental solution is known. But this may not always be possible. In addition, it is well known that the boundary integrals may become singular or nearly singular, depending on the distance between the source point and the node being integrated. It is also known that when the heterogeneity to be modelled tends towards zero thickness the direct BEM degenerates and becomes inaccurate. Various techniques have been proposed, such as, the use of dual boundary element techniques (DBEM), which lead to hypersingular integrals.Knowing that the simulation of three-dimensional phenomena involves high computational costs, different authors have proposed a number of BEM formulations. Ma et al. [8] applied a BEM formulation to study transient heat conduction in 3D solids with fiber inclusions. Jabłoński [9] solved 3D Laplace and Poisson equations by proposing the analytical evaluation of the surface integrals appearing in BEMs. Qin et al. [10] implemented changes to the conventional distance transformation technique to evaluate nearly singular integrands on 3D boundary elements, including planar and curved surface elements and very irregular elements of slender shape.The correct integration of the singular and hypersingular integrals is one of the big challenges of the BEM and DBEM.  Various methods have been proposed by the BEM research community to overcome some difficulties posed by singularities, as described by Zhou et al. [11]. As well as analytical and semianalytical methods, the alternatives include nonlinear transformation techniques [12–14] or distance transformation techniques [15–17].Different approaches have been proposed to deal with hypersingular integrals that arise in DBEM [18]. Solutions for specific 2D problems may be found in Cruse [19], V. Sládek and J. Sládek [20], Prosper [21], Prosper and Kausel [22], and Mendes and Tadeu [23]. Contrary to 2D problems, for which some closed form solutions for singular and hypersinguar integrals may be found, the authors are not aware that analytical solutions are known for singular and hypersingular integrals for 3D for heat diffusion problems. In fact, for 3D problems, singular integrals are mostly solved using numerical schemes based on Gaussian integration schemes. However, as the accuracy of the BEM is highly dependent on the precision of those integrals, some researchers are looking for semianalytical solutions or sophisticated approaches to solve specific problems [24]. Niu et al. [25] proposed a semianalytical algorithm for 3D elastic problems that require the evaluation of nearly strongly singular and hypersingular integrals on the triangular and quadrilateral elements. Applying a scheme of integration by parts, the nearly singular surface integrals are transformed to a set of line integrals along the boundary for which standard numerical quadrature can be used. Tadeu and António [26] have recently presented a set of analytical solutions for singular and hypersingular integrals for 3D acoustic problems.This paper first presents an analytical evaluation of the singular and hypersingular integrals that appear in three-dimensional boundary element formulations in the frequency domain for heat diffusion, when the element being integrated is the loaded one. The proposed analytical expressions are then compared with the numerical integration results by means of matured techniques based on standard Gaussian quadrature, which uses an enormous amount of sampling points, to validate the singular integrals presented. The proposed hypersingular integral solutions are verified against solutions obtained by combining an analytical solution defined for a circular domain and a standard Gaussian quadrature integration scheme applied to the remaining boundary element domain.   These functions are subsequently incorporated in a BEM and used to simulate the three-dimensional heat 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.


## 2. 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).Figure 1
Three-dimensional geometry of the problem.This system is subjected to a pressure point heat source placed at(xs,ys,zs),
(2.1)f^(x,y,z,t)=δ(x-xs)δ(y-ys)δ(z-zs)eiωt,
where δ(x-xs),δ(y-ys), and δ(z-zs) are Dirac-delta functions, and ω is  the frequency of the source. The response of this heat source can be expressed by
(2.2)Tinc(x,y,z,xs,ys,zs,ω)=e-i-(iω/K)r02λr0,
in which K is the thermal diffusivity defined by λ/ρc, ω is the oscillating frequency, i=-1, and r0=(x-xs)2+(y-ys)2+(z-zs)2.


## 3. 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.


### 3.1. 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:
(3.1)(∂2∂x2+∂2∂y2+∂2∂z2)T(x,ω)+(kc)2T(x,ω)=0,
in which x=(x,y,z) and kc=-iω/K.The boundary integral equation, formulated in the frequency domain, can be constructed by applying the reciprocity theorem, leading to(3.2)bT(x0,ω)=∫Sq(x,nn1,ω)G(x,x0,ω)ds-∫SH(x,nn1,x0,ω)T(x,ω)ds+Tinc(x0,xs,ω),
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 x0=(x0,y0,z0); nn1 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 x0=(x0,y0,z0)∈S and 1 otherwise.The required Green’s functions for temperature and heat flux in an unbounded medium, in Cartesian coordinates, are given by(3.3)G(x,x0,ω)=e-ikcr4λπr,H(x,nn1,x0,ω)=e-ikcr(-ikcr-1)4λπr2∂r∂nn1,
with r=(x-x0)2+(y-y0)2+(z-z0)2.Null Heat Fluxes along Its Boundary
When null heat fluxes are prescribed along the boundary, (3.2) is simplified resulting in
(3.4)bT(x0,ω)=-∫SH(x,nn1,x0,ω)T(x,ω)ds+Tinc(x0,xs,ω).Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the following equation:(3.5)0=∫Sq(x,nn1,ω)G(x,x0,ω)ds+Tinc(x0,xs,ω).


### 3.2. 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)aT(x0,ω)+bq(x0,nn1,ω)=∫Sq(x,nn1,ω)G¯(x,nn2,x0,ω)ds-∫SH¯(x,nn1,nn2,x0,ω)T(x,ω)ds+T¯inc(x0,nn2,xs,ω).The Green's functionsG¯ 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, nn2  is the unit outward normal to the boundary S at the collocation points x0=(x0,y0,z0), defined by the vector nn2. 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:(3.7)G¯(x,nn2,x0,ω)=e-ikcr(-ikcr-1)4λπr2∂r∂nn2,H¯(x,nn1,nn2,x0,ω)=∂H∂x∂x∂nn2+∂H∂y∂y∂nn2+∂H∂z∂z∂nn2,
with
(3.8)∂H∂x(x,nn1,nn2,x0,ω)=14πλ{A[(∂r∂x)2∂x∂nn1+∂r∂x∂r∂y∂y∂nn1+∂r∂x∂r∂z∂z∂nn1]+B[∂x∂nn1]},∂H∂y(x,nn1,nn2,x0,ω)=14πλ{A[∂r∂x∂r∂y∂x∂nn1+(∂r∂y)2∂y∂nn1+∂r∂y∂r∂z∂z∂nn1]+B[∂y∂nn1]},∂H∂z(x,nn1,nn2,x0,ω)=14πλ{A[∂r∂x∂r∂z∂x∂nn1+∂r∂y∂r∂z∂y∂nn1+(∂r∂z)2∂z∂nn1]+B[∂z∂nn1]},A=-kc2e-ikcrr+3ikce-ikcrr2+3e-ikcrr3,B=-ikce-ikcrr2-e-ikcrr3.
In (3.6) the incident heat field is computed as
(3.9)T¯inc(x,nn2,xs,ω)=e-ikcr0(-ikcr0-1)2λr02(∂r0∂x∂x∂nn2+∂r0∂y∂y∂nn2+∂r0∂z∂z∂nn2).Null Heat Fluxes along Its Boundary
Null heat fluxes are prescribed along the boundary, which leads to the following equation:(3.10)aT(x0,ω)=-∫SH¯(x,nn1,nn2,x0,ω)T(x,ω)ds+T¯inc(x0,nn2,xs,ω).Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the equation:(3.11)bq(x0,nn1,ω)=∫Sq(x,nn1,ω)G¯(x,nn2,x0,ω)ds+T¯inc(x0,nn2,xs,ω).


### 3.3. 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 widthW (in the x direction) and length L (in the z direction) shown in Figure 2.Figure 2
Scheme of a planar boundary element.


#### 3.3.1. Singular Integrals in (3.2)

Since in this caser is perpendicular to the normal (e.g., rnn1=0), the singular term ∫-L/2L/2∫-W/2W/2H(x,nn1,x0,ω)dxdz disappears. On the other hand, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz leads to a singular term.This integration is evaluated by first expressingG(x,x0,ω) as the sum of two-dimensional 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,x0,ω). 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),
(3.12)G^(x,y,x0,y0,kz,ω)=-i8πλH0(kr-0),
in which Hn(⋯) are second kind Hankel functions of the order n, k=-iω/K-kz2, with Im(k)<0,r-0=(x-x0)2+(y-y0)2, where kz is the wavenumber in the z direction.Assuming the presence of an infinite set of equally spaced sources in thez direction, the former Green's function can be recast as:
(3.13)G(x,x0,ω)=2πLvs∑m=-∞∞G^(x,y,x0,y0,kz,ω)e-ikzmz,
where Lvs is the spatial source interval, and kzm=(2π/Lvs)m.This equation converges very fast and can be approximated by a finite sum of terms(M). The distance Lvs 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.0K1/Δf. The number of terms depends on this distance Lvs. However, a larger value of Lvs 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 thez dimension.This procedure allows the integration of∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz to be obtained in a similar manner, thus
(3.14)∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz=-i4λLvs∑m=-MM∫-L/2L/2I1(kr-0)e-ikzmzdz,
with
(3.15)∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)Lifm=0,∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I1(kr-0)=∫-W/2W/2H0(kr-0)dx.I
1
(
k
r
-
0
) is calculated analytically, following the expressions in Tadeu et al. [27, 28],
(3.16)∫0W/2H0(kr-0)dx=W2H0(kW2)+πL4[H1(kW2)S0(kW2)-H0(kW2)S1(kW2)],
where Sns(⋯) are Struve functions of order ns.


#### 3.3.2. Hypersingular Integrals in (3.6)

Since in this caser is perpendicular to the normal (e.g., rnn2=0), the singular term ∫-L/2L/2∫-W/2W/2G¯(x,nn2,x0,ω)dxdz disappears. However, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz (with nn2=nn1) leads to a hypersingular term.Using the procedure above, the evaluation of this integration can performed by writingH¯(x,nn1,nn2,x0,ω) 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:
(3.17)H^(x,y,nn1,x0,y0,kz,ω)=ik8πλH1(kr-0)(∂r-0∂x∂x∂nn1+∂r-0∂y∂y∂nn1),
which leads to
(3.18)H-(x,y,nn1,nn1,x0,y0,kz,ω)=ik8πλ[-kH2(kr-0)(∂r-0∂x∂x∂nnl+∂r-0∂y∂y∂nnl)2+H1(kr-0)r-0].
This procedure allows the integration of ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω) to be obtained as
(3.19)∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz=2πLvs∑m=-MM∫-L/2L/2I2(kr-0)e-ikzmzdz,
with
(3.20)∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)Lifm=0,∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I2(kr-0)=∫-W/2W/2(ik/8πλ)[-kH2(kr-0)((∂r-0/∂x)(∂x/∂nnl)+(∂r-0/∂y)(∂y/∂nnl))2+(H1(kr-0)/r-0)]dx.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary element and by considering its thermal equilibrium (see Figure 3).Figure 3
Dynamic equilibrium on the singular element.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary
(3.21)∫SBEH-dS=∫SW/2H-dS-∫V1K∂H^∂tdV.
Both integrands to the right of the equal sign in (3.21) are well behaved and can be evaluated directly,
(3.22)∫SW/2H-dS=ik8πλ∫0πW2∂(H1(kr-0))∂r-0|W/2sin⁡⁡θdθ=ik24πλ[W2H0(kW2)-1kH1(kW2)],∫V1K∂H^∂tdV=-ik38πλ∫0π∫0W/2r-0H1(kr-0)sin⁡⁡θdr-0dθ=ik24πλ[W2H0(kW2)-∫0W/2H0(kr-0)dr-0].
Substituting (3.22) into (3.21), one obtains
(3.23)I2(kr-0)=ik24πλ[∫0W/2H0(kr-0)dr-0-1kH1(kW2)].


## 3.1. 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:
(3.1)(∂2∂x2+∂2∂y2+∂2∂z2)T(x,ω)+(kc)2T(x,ω)=0,
in which x=(x,y,z) and kc=-iω/K.The boundary integral equation, formulated in the frequency domain, can be constructed by applying the reciprocity theorem, leading to(3.2)bT(x0,ω)=∫Sq(x,nn1,ω)G(x,x0,ω)ds-∫SH(x,nn1,x0,ω)T(x,ω)ds+Tinc(x0,xs,ω),
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 x0=(x0,y0,z0); nn1 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 x0=(x0,y0,z0)∈S and 1 otherwise.The required Green’s functions for temperature and heat flux in an unbounded medium, in Cartesian coordinates, are given by(3.3)G(x,x0,ω)=e-ikcr4λπr,H(x,nn1,x0,ω)=e-ikcr(-ikcr-1)4λπr2∂r∂nn1,
with r=(x-x0)2+(y-y0)2+(z-z0)2.Null Heat Fluxes along Its Boundary
When null heat fluxes are prescribed along the boundary, (3.2) is simplified resulting in
(3.4)bT(x0,ω)=-∫SH(x,nn1,x0,ω)T(x,ω)ds+Tinc(x0,xs,ω).Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the following equation:(3.5)0=∫Sq(x,nn1,ω)G(x,x0,ω)ds+Tinc(x0,xs,ω).


## 3.2. 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)aT(x0,ω)+bq(x0,nn1,ω)=∫Sq(x,nn1,ω)G¯(x,nn2,x0,ω)ds-∫SH¯(x,nn1,nn2,x0,ω)T(x,ω)ds+T¯inc(x0,nn2,xs,ω).The Green's functionsG¯ 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, nn2  is the unit outward normal to the boundary S at the collocation points x0=(x0,y0,z0), defined by the vector nn2. 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:(3.7)G¯(x,nn2,x0,ω)=e-ikcr(-ikcr-1)4λπr2∂r∂nn2,H¯(x,nn1,nn2,x0,ω)=∂H∂x∂x∂nn2+∂H∂y∂y∂nn2+∂H∂z∂z∂nn2,
with
(3.8)∂H∂x(x,nn1,nn2,x0,ω)=14πλ{A[(∂r∂x)2∂x∂nn1+∂r∂x∂r∂y∂y∂nn1+∂r∂x∂r∂z∂z∂nn1]+B[∂x∂nn1]},∂H∂y(x,nn1,nn2,x0,ω)=14πλ{A[∂r∂x∂r∂y∂x∂nn1+(∂r∂y)2∂y∂nn1+∂r∂y∂r∂z∂z∂nn1]+B[∂y∂nn1]},∂H∂z(x,nn1,nn2,x0,ω)=14πλ{A[∂r∂x∂r∂z∂x∂nn1+∂r∂y∂r∂z∂y∂nn1+(∂r∂z)2∂z∂nn1]+B[∂z∂nn1]},A=-kc2e-ikcrr+3ikce-ikcrr2+3e-ikcrr3,B=-ikce-ikcrr2-e-ikcrr3.
In (3.6) the incident heat field is computed as
(3.9)T¯inc(x,nn2,xs,ω)=e-ikcr0(-ikcr0-1)2λr02(∂r0∂x∂x∂nn2+∂r0∂y∂y∂nn2+∂r0∂z∂z∂nn2).Null Heat Fluxes along Its Boundary
Null heat fluxes are prescribed along the boundary, which leads to the following equation:(3.10)aT(x0,ω)=-∫SH¯(x,nn1,nn2,x0,ω)T(x,ω)ds+T¯inc(x0,nn2,xs,ω).Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the equation:(3.11)bq(x0,nn1,ω)=∫Sq(x,nn1,ω)G¯(x,nn2,x0,ω)ds+T¯inc(x0,nn2,xs,ω).


## 3.3. 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 widthW (in the x direction) and length L (in the z direction) shown in Figure 2.Figure 2
Scheme of a planar boundary element.


### 3.3.1. Singular Integrals in (3.2)

Since in this caser is perpendicular to the normal (e.g., rnn1=0), the singular term ∫-L/2L/2∫-W/2W/2H(x,nn1,x0,ω)dxdz disappears. On the other hand, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz leads to a singular term.This integration is evaluated by first expressingG(x,x0,ω) as the sum of two-dimensional 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,x0,ω). 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),
(3.12)G^(x,y,x0,y0,kz,ω)=-i8πλH0(kr-0),
in which Hn(⋯) are second kind Hankel functions of the order n, k=-iω/K-kz2, with Im(k)<0,r-0=(x-x0)2+(y-y0)2, where kz is the wavenumber in the z direction.Assuming the presence of an infinite set of equally spaced sources in thez direction, the former Green's function can be recast as:
(3.13)G(x,x0,ω)=2πLvs∑m=-∞∞G^(x,y,x0,y0,kz,ω)e-ikzmz,
where Lvs is the spatial source interval, and kzm=(2π/Lvs)m.This equation converges very fast and can be approximated by a finite sum of terms(M). The distance Lvs 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.0K1/Δf. The number of terms depends on this distance Lvs. However, a larger value of Lvs 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 thez dimension.This procedure allows the integration of∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz to be obtained in a similar manner, thus
(3.14)∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz=-i4λLvs∑m=-MM∫-L/2L/2I1(kr-0)e-ikzmzdz,
with
(3.15)∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)Lifm=0,∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I1(kr-0)=∫-W/2W/2H0(kr-0)dx.I
1
(
k
r
-
0
) is calculated analytically, following the expressions in Tadeu et al. [27, 28],
(3.16)∫0W/2H0(kr-0)dx=W2H0(kW2)+πL4[H1(kW2)S0(kW2)-H0(kW2)S1(kW2)],
where Sns(⋯) are Struve functions of order ns.


### 3.3.2. Hypersingular Integrals in (3.6)

Since in this caser is perpendicular to the normal (e.g., rnn2=0), the singular term ∫-L/2L/2∫-W/2W/2G¯(x,nn2,x0,ω)dxdz disappears. However, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz (with nn2=nn1) leads to a hypersingular term.Using the procedure above, the evaluation of this integration can performed by writingH¯(x,nn1,nn2,x0,ω) 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:
(3.17)H^(x,y,nn1,x0,y0,kz,ω)=ik8πλH1(kr-0)(∂r-0∂x∂x∂nn1+∂r-0∂y∂y∂nn1),
which leads to
(3.18)H-(x,y,nn1,nn1,x0,y0,kz,ω)=ik8πλ[-kH2(kr-0)(∂r-0∂x∂x∂nnl+∂r-0∂y∂y∂nnl)2+H1(kr-0)r-0].
This procedure allows the integration of ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω) to be obtained as
(3.19)∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz=2πLvs∑m=-MM∫-L/2L/2I2(kr-0)e-ikzmzdz,
with
(3.20)∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)Lifm=0,∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I2(kr-0)=∫-W/2W/2(ik/8πλ)[-kH2(kr-0)((∂r-0/∂x)(∂x/∂nnl)+(∂r-0/∂y)(∂y/∂nnl))2+(H1(kr-0)/r-0)]dx.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary element and by considering its thermal equilibrium (see Figure 3).Figure 3
Dynamic equilibrium on the singular element.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary
(3.21)∫SBEH-dS=∫SW/2H-dS-∫V1K∂H^∂tdV.
Both integrands to the right of the equal sign in (3.21) are well behaved and can be evaluated directly,
(3.22)∫SW/2H-dS=ik8πλ∫0πW2∂(H1(kr-0))∂r-0|W/2sin⁡⁡θdθ=ik24πλ[W2H0(kW2)-1kH1(kW2)],∫V1K∂H^∂tdV=-ik38πλ∫0π∫0W/2r-0H1(kr-0)sin⁡⁡θdr-0dθ=ik24πλ[W2H0(kW2)-∫0W/2H0(kr-0)dr-0].
Substituting (3.22) into (3.21), one obtains
(3.23)I2(kr-0)=ik24πλ[∫0W/2H0(kr-0)dr-0-1kH1(kW2)].


## 3.3.1. Singular Integrals in (3.2)

Since in this caser is perpendicular to the normal (e.g., rnn1=0), the singular term ∫-L/2L/2∫-W/2W/2H(x,nn1,x0,ω)dxdz disappears. On the other hand, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz leads to a singular term.This integration is evaluated by first expressingG(x,x0,ω) as the sum of two-dimensional 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,x0,ω). 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),
(3.12)G^(x,y,x0,y0,kz,ω)=-i8πλH0(kr-0),
in which Hn(⋯) are second kind Hankel functions of the order n, k=-iω/K-kz2, with Im(k)<0,r-0=(x-x0)2+(y-y0)2, where kz is the wavenumber in the z direction.Assuming the presence of an infinite set of equally spaced sources in thez direction, the former Green's function can be recast as:
(3.13)G(x,x0,ω)=2πLvs∑m=-∞∞G^(x,y,x0,y0,kz,ω)e-ikzmz,
where Lvs is the spatial source interval, and kzm=(2π/Lvs)m.This equation converges very fast and can be approximated by a finite sum of terms(M). The distance Lvs 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.0K1/Δf. The number of terms depends on this distance Lvs. However, a larger value of Lvs 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 thez dimension.This procedure allows the integration of∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz to be obtained in a similar manner, thus
(3.14)∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz=-i4λLvs∑m=-MM∫-L/2L/2I1(kr-0)e-ikzmzdz,
with
(3.15)∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)Lifm=0,∫-L/2L/2I1(kr-0)e-ikzmzdz=I1(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I1(kr-0)=∫-W/2W/2H0(kr-0)dx.I
1
(
k
r
-
0
) is calculated analytically, following the expressions in Tadeu et al. [27, 28],
(3.16)∫0W/2H0(kr-0)dx=W2H0(kW2)+πL4[H1(kW2)S0(kW2)-H0(kW2)S1(kW2)],
where Sns(⋯) are Struve functions of order ns.


## 3.3.2. Hypersingular Integrals in (3.6)

Since in this caser is perpendicular to the normal (e.g., rnn2=0), the singular term ∫-L/2L/2∫-W/2W/2G¯(x,nn2,x0,ω)dxdz disappears. However, the integration of the Green’s function ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz (with nn2=nn1) leads to a hypersingular term.Using the procedure above, the evaluation of this integration can performed by writingH¯(x,nn1,nn2,x0,ω) 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:
(3.17)H^(x,y,nn1,x0,y0,kz,ω)=ik8πλH1(kr-0)(∂r-0∂x∂x∂nn1+∂r-0∂y∂y∂nn1),
which leads to
(3.18)H-(x,y,nn1,nn1,x0,y0,kz,ω)=ik8πλ[-kH2(kr-0)(∂r-0∂x∂x∂nnl+∂r-0∂y∂y∂nnl)2+H1(kr-0)r-0].
This procedure allows the integration of ∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω) to be obtained as
(3.19)∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz=2πLvs∑m=-MM∫-L/2L/2I2(kr-0)e-ikzmzdz,
with
(3.20)∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)Lifm=0,∫-L/2L/2I2(kr-0)e-ikzmzdz=I2(kr-0)2sin⁡⁡(kzm(L/2))kzmifm≠0,
where I2(kr-0)=∫-W/2W/2(ik/8πλ)[-kH2(kr-0)((∂r-0/∂x)(∂x/∂nnl)+(∂r-0/∂y)(∂y/∂nnl))2+(H1(kr-0)/r-0)]dx.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary element and by considering its thermal equilibrium (see Figure 3).Figure 3
Dynamic equilibrium on the singular element.The integration ofI2(kr-0) is performed indirectly by isolating a semicylinder just above the boundary
(3.21)∫SBEH-dS=∫SW/2H-dS-∫V1K∂H^∂tdV.
Both integrands to the right of the equal sign in (3.21) are well behaved and can be evaluated directly,
(3.22)∫SW/2H-dS=ik8πλ∫0πW2∂(H1(kr-0))∂r-0|W/2sin⁡⁡θdθ=ik24πλ[W2H0(kW2)-1kH1(kW2)],∫V1K∂H^∂tdV=-ik38πλ∫0π∫0W/2r-0H1(kr-0)sin⁡⁡θdr-0dθ=ik24πλ[W2H0(kW2)-∫0W/2H0(kr-0)dr-0].
Substituting (3.22) into (3.21), one obtains
(3.23)I2(kr-0)=ik24πλ[∫0W/2H0(kr-0)dr-0-1kH1(kW2)].


## 4. 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 points in the case of singular integrals and making use of analytical expressions deduced for circular domains in the case of hypersingular integrals.


### 4.1. Verification of the Singular Integral Integration

The integrations∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz 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 m2 s−1).  Lvs is assumed to be 60.0 m.Figure 4
Scheme of a planar quadrangular boundary element.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 (kzm=(2π/Lvs)m, m=1,…,M), that is, as kz 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 kz increases.  Similar behavior was found throughout the frequency range.Convergence of the proposed analytical integration of the singular integrals: (a) real part; (b) imaginary part.
(a)
(b)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.Figure 6
BEM: Analytical results (as continuous lines) and numerical results using a standard Gaussian quadrature with 2000 × 2000 sampling points.Figure 7
Variation of BEM error different numbers of Gauss sampling points.


### 4.2. Verification of the Hypersingular Integral Integration

To illustrate the convergence of (3.19), Figure 8 shows the integration result as kz increases, that is, as the number of terms is successively taken into account, when the integration is performed for a frequency of 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).Convergence of the proposed analytical integration solution: (a) real part; (b) imaginary part.
(a)
(b)Figures8(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 kz 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∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz 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).Figure 9
Scheme of the verification of the analytical integrations.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 ofH¯(x,nn1,nn2,x0,ω) 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). 								(4.1)∫SBEH-dS=∫SW/2H-dS-∫V1K∂H^∂tdV.
As before, both integrands to the right of the equals sign in (4.1) are well behaved, and can be evaluated directly,
(4.2)∫SW/2H-dS=14πλ∫02π∫0π/2(W2)2∂∂r-0(e-ikcr(-ikcr-1)r2)|W/2cos⁡⁡θsin⁡⁡θdθdϕ=14λ[-kc2W2e-ikc(W/2)+2ikce-ikc(W/2)+4We-ikc(W/2)],∫V1K∂H^∂tdV=-k24πλ∫02π∫0π/2∫0W/2e-ikcr(-ikcr-1)cos⁡⁡θsin⁡⁡θdrdθdϕ=14λ[-kc2W2e-ikc(W/2)+2ikce-ikc(W/2)-2ikc].
Substituting (4.2) into (4.1), one obtains
(4.3)∫SBEH-dS=1λ[1We-ikc(W/2)+ikc2].
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.Figure 10
Dynamic equilibrium on a circular hypersingular domain.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 to the 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-5Hz) with a frequency increment of 0.5×10-7 Hz. The thermal diffusivity and Lvs are assumed to be, respectively, 4.96278×10-7 m2 s−1 and 60.0 m.The results obtained with the proposed analytical expressions (as continuous lines) and with the constructed reference solution (marked lines) are illustrated in Figure11.  Analysis of this figure confirms that the results are very similar.Figure 11
TBEM: Analytical results versus constructed reference solution.


## 4.1. Verification of the Singular Integral Integration

The integrations∫-L/2L/2∫-W/2W/2G(x,x0,ω)dxdz 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 m2 s−1).  Lvs is assumed to be 60.0 m.Figure 4
Scheme of a planar quadrangular boundary element.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 (kzm=(2π/Lvs)m, m=1,…,M), that is, as kz 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 kz increases.  Similar behavior was found throughout the frequency range.Convergence of the proposed analytical integration of the singular integrals: (a) real part; (b) imaginary part.
(a)
(b)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.Figure 6
BEM: Analytical results (as continuous lines) and numerical results using a standard Gaussian quadrature with 2000 × 2000 sampling points.Figure 7
Variation of BEM error different numbers of Gauss sampling points.


## 4.2. Verification of the Hypersingular Integral Integration

To illustrate the convergence of (3.19), Figure 8 shows the integration result as kz increases, that is, as the number of terms is successively taken into account, when the integration is performed for a frequency of 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).Convergence of the proposed analytical integration solution: (a) real part; (b) imaginary part.
(a)
(b)Figures8(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 kz 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∫-L/2L/2∫-W/2W/2H¯(x,nn1,nn2,x0,ω)dxdz 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).Figure 9
Scheme of the verification of the analytical integrations.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 ofH¯(x,nn1,nn2,x0,ω) 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). 								(4.1)∫SBEH-dS=∫SW/2H-dS-∫V1K∂H^∂tdV.
As before, both integrands to the right of the equals sign in (4.1) are well behaved, and can be evaluated directly,
(4.2)∫SW/2H-dS=14πλ∫02π∫0π/2(W2)2∂∂r-0(e-ikcr(-ikcr-1)r2)|W/2cos⁡⁡θsin⁡⁡θdθdϕ=14λ[-kc2W2e-ikc(W/2)+2ikce-ikc(W/2)+4We-ikc(W/2)],∫V1K∂H^∂tdV=-k24πλ∫02π∫0π/2∫0W/2e-ikcr(-ikcr-1)cos⁡⁡θsin⁡⁡θdrdθdϕ=14λ[-kc2W2e-ikc(W/2)+2ikce-ikc(W/2)-2ikc].
Substituting (4.2) into (4.1), one obtains
(4.3)∫SBEH-dS=1λ[1We-ikc(W/2)+ikc2].
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.Figure 10
Dynamic equilibrium on a circular hypersingular domain.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 to the 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-5Hz) with a frequency increment of 0.5×10-7 Hz. The thermal diffusivity and Lvs are assumed to be, respectively, 4.96278×10-7 m2 s−1 and 60.0 m.The results obtained with the proposed analytical expressions (as continuous lines) and with the constructed reference solution (marked lines) are illustrated in Figure11.  Analysis of this figure confirms that the results are very similar.Figure 11
TBEM: Analytical results versus constructed reference solution.


## 5. Verification of the BEM Formulations

In this section, the proposed algorithms are validated using a circular cylindrical inclusion, with radiusa, aligned along the z axis (see Figure 12), for which analytical solutions can be derived. A point heat source placed at (xs,ys,zs) 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=Lt.Figure 12
Unbounded medium with an embedded cylindrical inclusion.The analytical solution for this problem is obtained by first applying a spatial Fourier transformation in thez 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=Lt 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:(5.1)T(x,xs,ω)=2πLvs∑m=-∞∞T^′(x,y,xs,ys,kzm,ω)e-ikzmz,
with kzm=(2π/Lvs)m(m=1,…,NS)  and
(5.2)T^′(x,y,xs,ys,kzm,ω)=T^(x,y,xs,ys,kzm,ω)[e-ikzmzm0+∑m=1NSz∑j=14e-ikzmzmj],
where zm0=z-zs, zm1=z+z0-2Ltm, zm2=z+z0+2Lt(m-1), zm3=z-z0-2Ltm, and zm4=z-z0+2Ltm.  The number of virtual sources NSz 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 Lvs 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:(5.3)T^(x,y,xs,ys,kz,ω)=-i4λH0(kr-0)+∑n=0∞AnHn(kr-)cos⁡⁡(nθ),
with  An=(i/4λ)(-1)nεnHn(kr-00)Jn(ka)/Hn(ka), r-=x2+y2, Jn(⋯) are Bessel functions of order n,  εn={1ifn=02ifn≠0, r-00=|x0| and θ=arctan(y/x).Null Heat Fluxes along Its Boundary
(5.4)
T
^
(
x
,
y
,
x
s
,
y
s
,
k
z
,
ω
)
=
-
i
4
λ
H
0
(
k
r
-
0
)
+
∑
n
=
0
∞
B
n
H
n
(
k
r
-
)
cos
(
n
θ
)
,
with Bn=(i/4λ)(-1)nεnHn(kr-00)[(n/a)Jn(ka)-kJn+1(ka)]/[(n/a)Hn(ka)-kHn+1(ka)].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 Figure13 (grid 1 (y=0.0m), grid 2 (z=1.0m)), 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).Figure 13
Scheme of the receiver grids placed in the vicinity of the cylindrical inclusion subjected to a heat source located atO. Grid 1 (y=0.0m) and grid 2 (z=1.0m).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).Figures14 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.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).
(a)
(b)
(c)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).
(a)
(b)
(c)Figures16 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.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).
(a)
(b)
(c)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 (with30×20 boundary elements); (c) error using TBEM (with 50×32 boundary elements).
(a)
(b)
(c)As intheprevious 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. When30×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).


## 6. 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.



---

*Source: 647038-2012-07-03.xml*