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

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.

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 𝑐, with a three-dimensional embedded inclusion, bounded by a surface 𝑆 (see Figure 1).

Figure 1 

Three-dimensional geometry of the problem.

This system is subjected to a pressure point heat source placed at (𝑥_𝑠,𝑦_𝑠,𝑧_𝑠), 𝑓𝑥,𝑦,𝑧,𝑡=𝛿𝑥−𝑥_𝑠𝛿𝑦−𝑦_𝑠𝛿𝑧−𝑧_𝑠𝑒^𝑖𝜔𝑡, (2.1) where 𝛿(𝑥−𝑥_𝑠),𝛿(𝑦−𝑦_𝑠), and 𝛿(𝑧−𝑧_𝑠) are Dirac-delta functions, and 𝜔 is the frequency of the source. The response of this heat source can be expressed by 𝑇_inc𝑥,𝑦,𝑧,𝑥_𝑠,𝑦_𝑠,𝑧_𝑠=𝑒,𝜔^{−𝑖−(𝑖𝜔/𝐾)𝑟0}2𝜆𝑟₀, (2.2) in which 𝐾 is the thermal diffusivity defined by 𝜆/𝜌𝑐, 𝜔 is the oscillating frequency, 𝑖=√−1, and 𝑟₀=(𝑥−𝑥_𝑠)²+(𝑦−𝑦_𝑠)²+(𝑧−𝑧_𝑠)².

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 (𝑇) at any point of the spatial 3D homogeneous solid domain is governed by the diffusion equation: 𝜕²𝜕𝑥²+𝜕²𝜕𝑦²+𝜕²𝜕𝑧²𝑇(𝐱,𝜔)+𝑘_𝑐²𝑇(𝐱,𝜔)=0, (3.1) in which 𝐱=(𝑥,𝑦,𝑧) and 𝑘_𝑐=−𝑖𝜔/𝐾.

The boundary integral equation, formulated in the frequency domain, can be constructed by applying the reciprocity theorem, leading to 𝑏𝑇𝐱_𝟎,𝜔=_𝑆𝑞𝐱,𝐧_𝑛1,𝜔𝐺𝐱,𝐱₀,𝜔d𝑠−_𝑆𝐻𝐱,𝐧_𝑛1,𝐱₀,𝜔𝑇(𝐱,𝜔)d𝑠+𝑇_inc𝐱₀,𝐱_𝑠,𝜔, (3.2) where 𝐺 and 𝐻 are, respectively, the fundamental solutions (Green’s functions) for the temperature (𝑇) and heat flux (𝑞), at a point 𝐱=(𝑥,𝑦,𝑧) on the boundary 𝑆, due to a virtual point heat source at 𝐱₀=(𝑥₀,𝑦₀,𝑧₀); 𝐧_𝑛1 represents the unit outward normal along the boundary 𝑆, at 𝐱=(𝑥,𝑦,𝑧); 𝑏 is a constant defined by the shape of the boundary, taking the value 1/2 if 𝐱₀=(𝑥₀,𝑦₀,𝑧₀)∈𝑆 and 1 otherwise.

The required Green’s functions for temperature and heat flux in an unbounded medium, in Cartesian coordinates, are given by 𝐺𝐱,𝐱₀=𝑒,𝜔^{−𝑖𝑘𝑐𝑟},𝐻4𝜆𝜋𝑟𝐱,𝐧_𝑛1,𝐱₀=𝑒,𝜔^{−𝑖𝑘𝑐𝑟}−𝑖𝑘_𝑐𝑟−14𝜆𝜋𝑟²𝜕𝑟𝜕𝐧_𝑛1, (3.3) with 𝑟=(𝑥−𝑥₀)²+(𝑦−𝑦₀)²+(𝑧−𝑧₀)².

Null Heat Fluxes along Its Boundary
When null heat fluxes are prescribed along the boundary, (3.2) is simplified resulting in 𝑏𝑇𝐱₀,𝜔=−_𝑆𝐻𝐱,𝐧_𝑛1,𝐱₀,𝜔𝑇(𝐱,𝜔)𝑑𝑠+𝑇_inc𝐱₀,𝐱_𝑠,𝜔. (3.4)

Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the following equation: 0=_𝑆𝑞𝐱,𝐧_𝑛1,𝜔𝐺𝐱,𝐱₀,𝜔𝑑𝑠+𝑇_inc𝐱₀,𝐱_𝑠,𝜔. (3.5)

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: 𝑎𝑇𝐱₀,𝜔+𝑏𝑞𝐱₀,𝐧_𝑛1,𝜔=_𝑆𝑞𝐱,𝐧_𝑛1,𝜔𝐺𝐱,𝐧_𝑛2,𝐱₀−,𝜔𝑑𝑠_𝑆𝐻𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔𝑇(𝐱,𝜔)𝑑𝑠+𝑇_inc𝐱₀,𝐧_𝑛2,𝐱_𝑠,𝜔. (3.6)

The Green's functions 𝐺 and 𝐻 are defined by applying the gradient operator to 𝐺 and 𝐻, which can be seen as the derivatives of these former Green's functions, to obtain heat fluxes and temperatures. In these equations, 𝐧_𝑛2 is the unit outward normal to the boundary 𝑆 at the collocation points 𝐱₀=(𝑥₀,𝑦₀,𝑧₀), defined by the vector 𝐧_𝑛2. In this equation, the factor 𝑎 is null for piecewise planar boundary elements.

The required three-dimensional Green’s functions for an unbounded space are now defined as: 𝐺𝐱,𝐧_𝑛2,𝐱₀=𝑒,𝜔^{−𝑖𝑘𝑐𝑟}−𝑖𝑘_𝑐𝑟−14𝜆𝜋𝑟²𝜕𝑟𝜕𝐧_𝑛2,𝐻𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔=𝜕𝐻𝜕𝑥𝜕𝑥𝜕𝐧_𝑛2+𝜕𝐻𝜕𝑦𝜕𝑦𝜕𝐧_𝑛2+𝜕𝐻𝜕𝑧𝜕𝑧𝜕𝐧_𝑛2, (3.7) with 𝜕𝐻𝜕𝑥𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀=1,𝜔𝐴4𝜋𝜆𝜕𝑟𝜕𝑥²𝜕𝑥𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑦𝜕𝑦𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑧𝜕𝑧𝜕𝐧_𝑛1+𝐵𝜕𝑥𝜕𝐧_𝑛1,𝜕𝐻𝜕𝑦𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀=1,𝜔𝐴4𝜋𝜆𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑦𝜕𝑥𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑦²𝜕𝑦𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑦𝜕𝑟𝜕𝑧𝜕𝑧𝜕𝐧_𝑛1+𝐵𝜕𝑦𝜕𝐧_𝑛1,𝜕𝐻𝜕𝑧𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀=1,𝜔𝐴4𝜋𝜆𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑧𝜕𝑥𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑦𝜕𝑟𝜕𝑧𝜕𝑦𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑧²𝜕𝑧𝜕𝐧_𝑛1+𝐵𝜕𝑧𝜕𝐧_𝑛1,𝑘𝐴=−^2𝑐𝑒^{−𝑖𝑘𝑐𝑟}𝑟+3𝑖𝑘_𝑐𝑒^{−𝑖𝑘𝑐𝑟}𝑟²+3𝑒^{−𝑖𝑘𝑐𝑟}𝑟³,𝐵=−𝑖𝑘_𝑐𝑒^{−𝑖𝑘𝑐𝑟}𝑟²−𝑒^{−𝑖𝑘𝑐𝑟}𝑟³. (3.8) In (3.6) the incident heat field is computed as 𝑇_inc𝐱,𝐧_𝑛2,𝐱_𝑠=𝑒,𝜔^{−𝑖𝑘𝑐𝑟0}−𝑖𝑘_𝑐𝑟₀−12𝜆𝑟²₀𝜕𝑟₀𝜕𝑥𝜕𝑥𝜕𝐧_𝑛2+𝜕𝑟₀𝜕𝑦𝜕𝑦𝜕𝐧_𝑛2+𝜕𝑟₀𝜕𝑧𝜕𝑧𝜕𝐧_𝑛2. (3.9)

Null Heat Fluxes along Its Boundary
Null heat fluxes are prescribed along the boundary, which leads to the following equation: 𝑎𝑇𝐱₀,𝜔=−_𝑆𝐻𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔𝑇(𝐱,𝜔)𝑑𝑠+𝑇_inc𝐱₀,𝐧_𝑛2,𝐱_𝑠,𝜔. (3.10)

Null Temperatures along Its Boundary
Null temperatures on the surface of Inclusion are now prescribed, which leads to the equation: 𝑏𝑞𝐱₀,𝐧_𝑛1,𝜔=_𝑆𝑞𝐱,𝐧_𝑛1,𝜔𝐺𝐱,𝐧_𝑛2,𝐱₀,𝜔𝑑𝑠+𝑇_inc𝐱₀,𝐧_𝑛2,𝐱_𝑠,𝜔. (3.11)

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 𝑆 into 𝑁 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 𝑊 (in the 𝑥 direction) and length 𝐿 (in the 𝑧 direction) shown in Figure 2.

Figure 2 

Scheme of a planar boundary element.

3.3.1. Singular Integrals in (3.2)

Since in this case 𝑟 is perpendicular to the normal (e.g., 𝐫𝐧_𝑛1=0), the singular term ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐻(𝐱,𝐧_𝑛1,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 disappears. On the other hand, the integration of the Green’s function ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐺(𝐱,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 leads to a singular term.

This integration is evaluated by first expressing 𝐺(𝐱,𝐱₀,𝜔) 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 𝑧 direction to the three-dimensional Green's function 𝐺(𝐱,𝐱₀,𝜔). 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 (𝑧), 𝐺𝑥,𝑦,𝑥₀,𝑦₀,𝑘_𝑧,𝜔=−𝑖𝐻8𝜋𝜆⁰𝑘𝑟₀, (3.12) in which 𝐻_𝑛(⋯) are second kind Hankel functions of the order 𝑛, 𝑘=−𝑖𝜔/𝐾−𝑘²_𝑧, with Im(𝑘)<0,𝑟₀=(𝑥−𝑥₀)²+(𝑦−𝑦₀)², where 𝑘_𝑧 is the wavenumber in the 𝑧 direction.

Assuming the presence of an infinite set of equally spaced sources in the 𝑧 direction, the former Green's function can be recast as: 𝐺𝐱,𝐱₀,𝜔=2𝜋𝐿^∞_{𝑣𝑠𝑚=−∞}𝐺𝑥,𝑦,𝑥₀,𝑦₀,𝑘_𝑧,𝜔𝑒^{−𝑖𝑘𝑧𝑚𝑧}, (3.13) where 𝐿_𝑣𝑠 is the spatial source interval, and 𝑘_𝑧𝑚=(√2𝜋/𝐿_𝑣𝑠)𝑚.

This equation converges very fast and can be approximated by a finite sum of terms (𝑀). The distance 𝐿_𝑣𝑠 needs to be large enough to avoid spatial contamination. Given that the frequency increment defines the time window for the analysis, 1/Δ𝑓, the minimum distance needs to be at least 2.0√𝐾1/Δ𝑓. The number of terms depends on this distance 𝐿_𝑣𝑠. However, a larger value of 𝐿_𝑣𝑠 corresponds to a higher number of terms (𝑀).

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 𝑧 dimension.

This procedure allows the integration of ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐺(𝐱,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 to be obtained in a similar manner, thus ^𝐿/2_−𝐿/2^𝑊/2_−𝑊/2𝐺𝐱,𝐱₀,𝜔𝑑𝑥𝑑𝑧=−𝑖4𝜆𝐿^𝑀_{𝑣𝑠𝑚=−𝑀}^𝐿/2_−𝐿/2𝐼₁𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧, (3.14) with ^𝐿/2_−𝐿/2𝐼₁𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧=𝐼₁𝑘𝑟₀𝐿if𝑚=0,_{𝐿/2−𝐿/2}𝐼₁𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧=𝐼₁𝑘𝑟₀2sin𝑘_𝑧𝑚(𝐿/2)𝑘_𝑧𝑚if𝑚≠0, (3.15) where 𝐼₁(𝑘𝑟₀)=∫^𝑊/2_−𝑊/2𝐻₀(𝑘𝑟₀)𝑑𝑥.

𝐼₁(𝑘𝑟₀) is calculated analytically, following the expressions in Tadeu et al. [27, 28], ₀^𝑊/2𝐻₀𝑘𝑟₀𝑊𝑑𝑥=2𝐻₀𝑘𝑊2𝐿+𝜋4𝐻₁𝑘𝑊2𝑆₀𝑘𝑊2−𝐻₀𝑘𝑊2𝑆₁𝑘𝑊2, (3.16) where 𝑆_𝑛𝑠(⋯) are Struve functions of order 𝑛𝑠.

3.3.2. Hypersingular Integrals in (3.6)

Since in this case 𝑟 is perpendicular to the normal (e.g., 𝐫𝐧_𝑛2=0), the singular term ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐺(𝐱,𝐧_𝑛2,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 disappears. However, the integration of the Green’s function ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐻(𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 (with 𝐧_𝑛2=𝐧_𝑛1) leads to a hypersingular term.

Using the procedure above, the evaluation of this integration can performed by writing 𝐻(𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔) 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: 𝐻𝑥,𝑦,𝐧_𝑛1,𝑥₀,𝑦₀,𝑘_𝑧,𝜔=𝑖𝑘𝐻8𝜋𝜆¹𝑘𝑟₀𝜕𝑟₀𝜕𝑥𝜕𝑥𝜕𝐧_𝑛1+𝜕𝑟₀𝜕𝑦𝜕𝑦𝜕𝐧_𝑛1, (3.17) which leads to 𝐻𝑥,𝑦,𝐧_𝑛1,𝐧_𝑛1,𝑥₀,𝑦₀,𝑘_𝑧,𝜔=𝑖𝑘8𝜋𝜆−𝑘𝐻₂𝑘𝑟₀𝜕𝑟₀𝜕𝑥𝜕𝑥𝜕𝐧_𝑛𝑙+𝜕𝑟₀𝜕𝑦𝜕𝑦𝜕𝐧_𝑛𝑙²+𝐻₁𝑘𝑟₀𝑟₀. (3.18) This procedure allows the integration of ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐻(𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔) to be obtained as ^𝐿/2_−𝐿/2^𝑊/2_−𝑊/2𝐻𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔𝑑𝑥𝑑𝑧=2𝜋𝐿^𝑀_{𝑣𝑠𝑚=−𝑀}^𝐿/2_−𝐿/2𝐼₂𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧, (3.19) with ^𝐿/2_−𝐿/2𝐼₂𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧=𝐼₂𝑘𝑟₀𝐿if𝑚=0,_{𝐿/2−𝐿/2}𝐼₂𝑘𝑟₀𝑒^{−𝑖𝑘𝑧𝑚𝑧}𝑑𝑧=𝐼₂𝑘𝑟₀2sin𝑘_𝑧𝑚(𝐿/2)𝑘_𝑧𝑚if𝑚≠0, (3.20) where 𝐼₂(𝑘𝑟₀)=∫^𝑊/2_−𝑊/2(𝑖𝑘/8𝜋𝜆)[−𝑘𝐻₂(𝑘𝑟₀)((𝜕𝑟₀/𝜕𝑥)(𝜕𝑥/𝜕𝐧_𝑛𝑙)+(𝜕𝑟₀/𝜕𝑦)(𝜕𝑦/𝜕𝐧_𝑛𝑙))²+(𝐻₁(𝑘𝑟₀)/𝑟₀)]𝑑𝑥.

The integration of 𝐼₂(𝑘𝑟₀) 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 of 𝐼₂(𝑘𝑟₀) is performed indirectly by isolating a semicylinder just above the boundary _𝑆BE𝐻𝑑𝑆=_𝑆𝑊/2𝐻𝑑𝑆−_𝑉1𝐾𝜕𝐻𝜕𝑡𝑑𝑉. (3.21) Both integrands to the right of the equal sign in (3.21) are well behaved and can be evaluated directly, _𝑆𝑊/2𝐻𝑑𝑆=𝑖𝑘8𝜋𝜆^𝜋₀𝑊2𝜕𝐻₁𝑘𝑟₀𝜕𝑟_0𝑊/2=sin𝜃𝑑𝜃𝑖𝑘²𝑊4𝜋𝜆2𝐻₀𝑘𝑊2−1𝑘𝐻₁𝑘𝑊2_𝑉,1𝐾𝜕𝐻𝜕𝑡𝑑𝑉=−𝑖𝑘³8𝜋𝜆^𝜋₀₀^𝑊/2𝑟₀𝐻₁𝑘𝑟₀sin𝜃𝑑𝑟₀=𝑑𝜃𝑖𝑘²𝑊4𝜋𝜆2𝐻₀𝑘𝑊2−₀^𝑊/2𝐻₀𝑘𝑟₀𝑑𝑟₀. (3.22) Substituting (3.22) into (3.21), one obtains 𝐼₂𝑘𝑟₀=𝑖𝑘²4𝜋𝜆₀^𝑊/2𝐻₀𝑘𝑟₀𝑑𝑟₀−1𝑘𝐻₁𝑘𝑊2. (3.23)

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 ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐺(𝐱,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 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⁻⁷ Hz in the frequency range (0.0,1×10⁻⁵ Hz). The thermal properties assigned to the medium are those of the mortar: specific heat 𝑐=780 J·kg⁻¹·°C⁻¹, density 1860 kg·m⁻³ and thermal conductivity 0.72 W·m⁻¹·°C⁻¹ (thermal diffusivity of 4.96278×10⁻⁷ m² s⁻¹). 𝐿_𝑣𝑠 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 (𝑘_𝑧𝑚=(√2𝜋/𝐿_𝑣𝑠)𝑚, 𝑚=1,…,𝑀), that is, as 𝑘_𝑧 increases, when the integration is performed for a frequency of 0.5×10⁻⁷ Hz. It can be observed that the integration converges very fast as 𝑘_𝑧 increases. Similar behavior was found throughout the frequency range.

(a)

(b)

Figure 5 

Convergence of the proposed analytical integration of the singular integrals: (a) real part; (b) imaginary part.

For the purpose of verification, the analytical results for the frequency range (0.0,1×10⁻⁵ 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 𝑘_𝑧 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⁻⁷ Hz. The responses were computed for complex frequencies 𝜔_𝑐=𝜔−𝑖𝜂 (with 𝜂=0.7Δ𝜔, Δ𝜔=2×𝜋×Δ𝑓, and Δ𝑓=0.5×10⁻⁷ Hz).

(a)

(b)

Figure 8 

Convergence of the proposed analytical integration solution: (a) real part; (b) imaginary part.

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 𝑘_𝑧 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 ∫^𝐿/2_−𝐿/2∫^𝑊/2_−𝑊/2𝐻(𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔)𝑑𝑥𝑑𝑧 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 of 𝐻(𝐱,𝐧_𝑛1,𝐧_𝑛2,𝐱₀,𝜔) 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). _𝑆BE𝐻𝑑𝑆=_𝑆𝑊/2𝐻𝑑𝑆−_𝑉1𝐾𝜕𝐻𝜕𝑡𝑑𝑉. (4.1) As before, both integrands to the right of the equals sign in (4.1) are well behaved, and can be evaluated directly, _𝑆𝑊/21𝐻𝑑𝑆=4𝜋𝜆₀^2𝜋₀^𝜋/2𝑊2²𝜕𝜕𝑟₀𝑒^{−𝑖𝑘𝑐𝑟}−𝑖𝑘_𝑐𝑟−1𝑟²_𝑊/2=1cos𝜃sin𝜃𝑑𝜃𝑑𝜙4𝜆−𝑘²_𝑐𝑊2𝑒^{−𝑖𝑘𝑐(𝑊/2)}+2𝑖𝑘_𝑐𝑒^{−𝑖𝑘𝑐(𝑊/2)}+4𝑊𝑒^{−𝑖𝑘𝑐(𝑊/2)},_𝑉1𝐾𝜕𝐻𝜕𝑡𝑑𝑉=−𝑘²4𝜋𝜆₀^2𝜋₀^𝜋/2₀^𝑊/2𝑒^{−𝑖𝑘𝑐𝑟}−𝑖𝑘_𝑐=1𝑟−1cos𝜃sin𝜃𝑑𝑟𝑑𝜃𝑑𝜙4𝜆−𝑘²_𝑐𝑊2𝑒^{−𝑖𝑘𝑐(𝑊/2)}+2𝑖𝑘_𝑐𝑒^{−𝑖𝑘𝑐(𝑊/2)}−2𝑖𝑘_𝑐. (4.2) Substituting (4.2) into (4.1), one obtains _𝑆BE1𝐻𝑑𝑆=𝜆1𝑊𝑒^{−𝑖𝑘𝑐(𝑊/2)}+𝑖𝑘_𝑐2. (4.3) 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⁻⁵ Hz) with a frequency increment of 0.5×10⁻⁷ Hz. The thermal diffusivity and 𝐿_𝑣𝑠 are assumed to be, respectively, 4.96278×10⁻⁷ m² s⁻¹ 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 Figure 11. 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 radius 𝑎, aligned along the 𝑧 axis (see Figure 12), for which analytical solutions can be derived. A point heat source placed at (𝑥_𝑠,𝑦_𝑠,𝑧_𝑠) 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 𝑧=0.0 m and 𝑧=𝐿_𝑡.

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 the 𝑧 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 𝑧=𝐿_𝑡 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 𝑧 direction and act as mirrors of the real source, so as to ensure the required boundary conditions.

Consider the following: 𝑇𝐱,𝐱_𝑠,𝜔=2𝜋𝐿^∞_{𝑣𝑠𝑚=−∞}𝑇^′𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝑘_𝑧𝑚,𝜔𝑒^{−𝑖𝑘𝑧𝑚𝑧}, (5.1) with 𝑘_𝑧𝑚=(√2𝜋/𝐿_𝑣𝑠)𝑚(𝑚=1,…,𝑁𝑆) and 𝑇^′𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝑘_𝑧𝑚,𝜔=𝑇𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝑘_𝑧𝑚,𝜔𝑒^{−𝑖𝑘𝑧𝑚𝑧𝑚0}+^{𝑁𝑆𝑧4}_𝑚=1𝑗=1𝑒^{−𝑖𝑘𝑧𝑚𝑧𝑚𝑗}, (5.2) where 𝑧_𝑚0=𝑧−𝑧_𝑠, 𝑧_𝑚1=𝑧+𝑧₀−2𝐿_𝑡𝑚, 𝑧_𝑚2=𝑧+𝑧₀+2𝐿_𝑡(𝑚−1), 𝑧_𝑚3=𝑧−𝑧₀−2𝐿_𝑡𝑚, and 𝑧_𝑚4=𝑧−𝑧₀+2𝐿_𝑡𝑚. The number of virtual sources 𝑁𝑆_𝑧 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/Δ𝑓. This procedure does not introduce any type of error into the computed time impulse response within the time window defined. Notice that 𝐿_𝑣𝑠 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 𝑥 axis.

Null Temperatures along Its Boundary
Consider the following: 𝑇𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝑘_𝑧𝑖,𝜔=−𝐻4𝜆⁰𝑘𝑟₀+^∞_𝑛=0𝐴_𝑛𝐻_𝑛(𝑘𝑟)cos(𝑛𝜃), (5.3) with 𝐴_𝑛=(𝑖/4𝜆)(−1)^𝑛𝜀_𝑛𝐻_𝑛(𝑘𝑟₀₀)𝐽_𝑛(𝑘𝑎)/𝐻_𝑛(𝑘𝑎), 𝑟=𝑥²+𝑦², 𝐽_𝑛(⋯) are Bessel functions of order 𝑛, 𝜀_𝑛=^1if₂^𝑛=0_if𝑛≠0, 𝑟₀₀=|𝑥₀| and 𝜃=arctan(𝑦/𝑥).

Null Heat Fluxes along Its Boundary
𝑇𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝑘_𝑧𝑖,𝜔=−𝐻4𝜆⁰𝑘𝑟₀+^∞_𝑛=0𝐵_𝑛𝐻_𝑛(𝑘𝑟)cos(𝑛𝜃), (5.4) with 𝐵_𝑛=(𝑖/4𝜆)(−1)^𝑛𝜀_𝑛𝐻_𝑛(𝑘𝑟₀₀)[(𝑛/𝑎)𝐽_𝑛(𝑘𝑎)−𝑘𝐽_𝑛+1(𝑘𝑎)]/[(𝑛/𝑎)𝐻_𝑛(𝑘𝑎)−𝑘𝐻_𝑛+1(𝑘𝑎)].

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 𝑛=𝑁+20 and for 𝑛=𝑁, 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 (𝑦=0.0m), grid 2 (𝑧=1.0m)), in the vicinity of a cylinder inclusion with a radius of 𝑎=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⁻⁷ Hz. The responses were computed for complex frequencies 𝜔_𝑐=𝜔−𝑖𝜂 (with 𝜂=0.7Δ𝜔, Δ𝜔=2×𝜋×Δ𝑓 and Δ𝑓=0.5×10⁻⁷ Hz).

Figure 13 

Scheme of the receiver grids placed in the vicinity of the cylindrical inclusion subjected to a heat source located at 𝑂. Grid 1 (𝑦=0.0m) and grid 2 (𝑧=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 𝑧 direction) and 50×32 (32 in the 𝑧 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.

(a)

(b)

(c)

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).

(a)

(b)

(c)

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).

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.

(a)

(b)

(c)

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).

(a)

(b)

(c)

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).

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. 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).

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.

Acknowledgments

The research work presented herein was supported by FEDER funds through the Operational Programme for Competitiveness Factors (COMPETE) and by national funds through the FCT (Portuguese Foundation for Science and Technology), under research Project PTDC/ECM/114189/2009.