3D heat diffusion simulation using 3D and 1D heat sources – Temperature and phase contrast results for defect detection using IRT

Applied Mathematical Modelling 40 (2016) 1576–1587

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

3D heat diffusion simulation using 3D and 1D heat sources – Temperature and phase contrast results for defect detection using IRT C. Serra a,∗, A. Tadeu b, N. Simões b a ITeCons – Instituto de Investigação e Desenvolvimento Tecnológico em Ciências da Construção, Rua Pedro Hispano s/n, 3030-289 Coimbra, Portugal b Department of Civil Engineering, University of Coimbra, Pólo II, Rua Luís Reis Santos, 3030-788 Coimbra, Portugal

a r t i c l e

i n f o

Article history: Received 30 May 2013 Revised 17 April 2015 Accepted 26 August 2015 Available online 26 September 2015 Keywords: Transient heat conduction Phase contrast Infrared thermography 3D heat source 1D heat source Normal-derivative integral equations (TBEM)

a b s t r a c t Infrared thermography (IRT) has proven to be a powerful non-destructive technique and it has been successfully used for detecting defects in a wide range of applications and sectors. Numerical modeling of heat transfer and diffusion in the presence of defects combined with experimental IRT results may be used both to detect and characterize existing defects with specific geometries and depths. However, the three-dimensional (3D) nature of defects combined with the need to simulate heat transfer and diffusion phenomena in transient regime often presents many challenges for researchers. The study presented in this paper is motivated by such difficulties and is intended to contribute to the interpretation of quantitative data results collected in experimental defect detection IRT tests and help with the definition of IRT experimental set-up parameters. In this study, 3D heat diffusion by conduction in the proximity of a 3D defect is modeled using a boundary element method (BEM) formulated in the frequency domain. The defect is a crack lodged in an unbounded solid medium with null thickness. In order to overcome difficulties that occur in the presence of null thickness elements, the BEM formulation was written in terms of normal-derivative integral equations (TBEM) and known analytical solutions were used to solve the resulting hypersingular integrals. The focus of this paper is to study the influence using either a 3D (point) energy source or a 1D (planar) energy source in heat diffusion simulations performed for IRT defect detection studies. Heat field and thermal wave phase results were computed in the presence of a defect and for when there is no defect present and a comparative analysis of the results obtained for a point and a planar source was carried out. In order to contribute to defects characterization studies, the influence of the crack’s characteristics such as its size, shape and placement (depth and position) was analyzed using a phase contrast approach. Other features that may be relevant in IRT experiments, such as the nature of the stimulus provided and its distance from the surface were also studied. The major findings achieved from varying those parameters are presented in the conclusions. © 2015 Elsevier Inc. All rights reserved.

∗

Corresponding author. Tel.: +351 910304560. E-mail addresses: [email protected] (C. Serra), [email protected] (A. Tadeu), [email protected] (N. Simões).

http://dx.doi.org/10.1016/j.apm.2015.08.007 S0307-904X(15)00489-8/© 2015 Elsevier Inc. All rights reserved.

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

1577

1. Introduction Analyzing thermal pattern images obtained using infrared thermography (IRT) has shown to be an effective non-destructive testing (NDT) technique which is used in many sectors for many applications, including the detection of hidden defects. Defects appear as anomalies in the thermal patterns of thermographic images since their presence affects the heat and moisture diffusion phenomena. Passive IRT studies are performed with materials in their natural thermal state and are mostly done using a static approach. However if a temperature gradient that allows the identification of anomalies in the thermal patterns is not naturally present, it may be necessary the use an additional heat source in order to produce it. This technique is known as active IRT [1]. Additionally, if results are obtained in a transient regime, and the nature of the heat source generated in active IRT is well defined, a quantitative characterization of defects may be performed by solving heat transfer problems using thermographic data. For this reason, active IRT has been used in NDT in a great number of areas, including in civil engineering [2]. In particular, the technique has been successfully used to test the integrity of composites and structures, having detected surface moisture and located defects in concrete composites up to a depth of about 10 cm, as reported by Wiggenhauser and Maierhofer [3]. Fig. 1 is a general representation of the experimental apparatus for an active IRT test. A known thermal wave generated by a heat source (stimulus) is applied to a test specimen, leading to heat flow inside it. The reflected thermal wave is altered by defects inside the specimen, which produces disturbances in the temperature patterns on the surface. These are recorded by the IRT camera as a function of time. Data is stored in a computer and can then be subject to a range of processing techniques (data analysis). A number of active IRT techniques have been categorized, depending on the thermal stimulation scheme used to produce a relevant temperature gradient on the surface of a test specimen, and on the data processing technique employed. In [4] Maldague considers several methods, such as pulse thermography, step-heating and lock-in thermography. Pulse thermography (PT) is one of the most popular thermal stimulation methods. It consists of briefly heating a test specimen and recording its temperature decay curve. Since a subsurface defect reduces the diffusion rate of the thermal front propagating by diffusion under the surface, the presence of such defect appears as an area of higher surface temperature. PT data results are thermal images corresponding to the mapping of the emitted thermal infrared power. In step heating (SH), a long heat pulse is generated and the increase in surface temperature is monitored. In lock-in thermography (LT), modulated heating at frequency ω generates a temperature field inside the material (near the surface), which is recorded by the equipment. This allows the observation of both the amplitude and phase of the resulting thermal waves. In LT experimental results, both amplitude and phase images are available. Phase images have the added advantage of being relatively independent of local optical and thermal surface features. Each of the various schemes has its own advantages and limitations and choosing the most suitable approach depends on the intended application. In any case, several heat distributions are possible: point inspections using laser or focused light beams; line inspection using line lamps, heated wires, linear air jets or a scanning laser, and surface inspections in reflection or transmission [5]. Ibarra-Castanedo et al. [6] explored various data analysis methods for processing raw thermal image results in subsurface defect detection and characterization studies. Maldague and Marinetti [7] first introduced pulsed phase thermography (PPT) as a processing technique which combines the advantages of both PT and LT. Test specimens are pulse-heated but thermal waves are unscrambled by performing the Fourier transform of the temperature decay and enabling the computation of phase images [7]. More recently, Arndt [8] introduced the term square-pulse thermography (SPT) in the frequency domain in an adaptation of PPT for civil engineering applications. Arndt concluded that SPT is an ideal approach for IRT use in civil engineering because phase

He Sou at rce GENERATED THERMAL WAVE

STIMULUS PROVIDED

DATA ANALYSIS

Control Unit

IRT Camera

REFLECTED THERMAL WAVE

Computer System

t Hea ce r u o S Fig. 1. Schematic representation of an experimental setup for an active IRT test.

Test Specimen

1578

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

images allow deeper probing than is possible with thermal images, which enhances defect detectability and improves the resolution of defect geometry. Furthermore, phase images are less influenced by surface characteristics, less sensitive to non-uniform heating, and do not need previous knowledge of the position of non-defect areas, which is necessary in simple thermal contrast processing. SPT in the frequency domain involves longer heating and observation times to fit the long response of construction materials and the depth of the inhomogeneities in the test specimens. A long pulse in the time domain (from a few seconds to several minutes) will yield data with clear amplitude in the frequency domain, which can be used effectively to estimate defect depth. However, this requires a sufficient recording time and frequency resolution. The best combination of heat source conditions and heating time can be determined numerically to avoid performing a large number of experiments. Experimental conditions may differ from theoretical conditions due to non-uniform emissivity, uneven heating of the test specimen surface and variations in the test material. However, numerical simulations may help define intervals for setup parameters and set the limit of the depth at which a defect is detectable [9]. In general terms, numerical modeling may be done using methods that require the discretization of the domain, such as the finite elements (FEM) [10] and finite differences methods (FDM) [11,12], on methods that require the discretization of the boundary, such as the boundary element method (BEM) [13,14], or even on meshless methods that do not need either domain discretization or boundary discretization, such as the method of fundamental solutions (MFS) [15]. In the present study, the BEM is used to model 3D heat diffusion in a homogeneous unbounded solid medium. The benefits and disadvantages of using the BEM in this type of problem are explored in the authors’ previous studies [16]. One of the major difficulties stems from the null thickness of the defect being modeled. To overcome it, the dual boundary element method (DBEM) [17] or the traction boundary element method (TBEM) written in terms of normal-derivative integral equations, may be used. However, these methods lead to hypersingular integrals [18]. For this problem, an analytical method proposed and verified by Tadeu et al. [19] for 3D heat diffusion using the BEM formulated in the frequency domain is used. The present study was motivated by an interest in assessing the potential of using IRT as an NDT to detect defects with specific geometries and depths, located within constructive elements. It has been established that the numerical modeling of heat transfer in the presence of defects, in conjunction with experimental IRT, can be employed as to detect the presence and assess existing defects in a nondestructive way. In SPT, in particular, numerical simulations are useful for establishing the optimum experimental conditions and reducing the number of tests performed, for a specific defect depth. But this modeling is complex. More specifically, the contribution of this study is that it examines the effect of the nature of the heat source on IRT defect depth studies. It does so by comparing numerical phase contrast calculation results obtained from using a point heat source with results from applying a planar heat source. In the sections which follow, we define the problem and briefly outline the boundary element formulation used (the normalderivative integral equations – TBEM), formulated in the frequency domain. This model is used to simulate heat diffusion generated by planar or point energy sources placed in the vicinity of a 3D crack. The resulting hypersingular integrals are solved using analytical solutions. The methodology used to obtain results in the time domain from the calculations performed in the frequency domain is briefly described, as well as the calculation required for generating thermal wave phase results. Then, the proposed approach is used to study the influence of applying either a point (3D) or planar (1D) heat source when modeling phase contrast results obtained for varying experimental parameters and defect properties which are relevant for the quantitative analysis of defects using IRT. These include the defect location or depth and position (inclination), along with its size and shape, as well as the distance at which the artificial heat source is placed. 2. Methods 2.1. Problem definition This work simulates the 3D heat diffusion that occurs around a 3D crack embedded in a unbounded spatially uniform solid medium of density ρ , thermal conductivity λ and specific heat c. Computations are made for heat diffusion generated by either a 1D or 3D heat source, as shown in Fig. 2. The crack is assumed to be of null thickness and null heat fluxes are prescribed along its surface S. 2.1.1. Incident field generated by a point heat source The heat field response of 3D heat source placed at xs = (xs , ys zs ) can be expressed by [16]:

√ iω Pe−i − K r0 Tinc (x, xs , ω) = , 2λr0

(1)

in which x = (x, y, z), ω is the frequency of the source, P is the amplitude of the heat source, K is the thermal diffusivity given by λ , i = √−1 and r = (x − x )2 + (y − y )2 + (z − z )2 . s s s 0 ρc 2.1.2. Incident field generated by a planar heat source The incident heat field generated by a 1D heat source can be computed by applying a double Fourier transformation in the y and z directions to Eq. (1). This allows the 3D incident field to be written as a summation of plane waves, and thus only the plane

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

1579

perpendicular to the x axis is taken into account, leading to:

Tinc (x, xs , ω) =

√ iω −iPe−i − K |x−x0 |



2λ − iKω

.

(2)

2.2. The normal-derivative integral equation (TBEM) Heat transfer by conduction in transient regime can be expressed by the diffusion equation:



∂2 ∂2 ∂2 + + ∂ x2 ∂ y 2 ∂ z 2 

in which kc =



T (x,ω) + (kc ) T (x, ω) = 0, 2

(3)

−iω K .

Given the null thickness of the defect, the following normal-derivative integral equation, derived from the classical boundary integral equation formulated in the frequency domain, can be used [19]:

a T (x0 , ω) = −



S

H (x, nn1 , nn2 , x0 , ω) T (x, ω) ds + T inc (x0 , nn2 , xs , ω),

(4)

in which:

H (x, nn1 , nn2 , x0 , ω) =

∂H ∂x ∂H ∂y ∂H ∂z + + , ∂ x ∂ nn2 ∂ y ∂ nn2 ∂ z ∂ nn2

(5)

  2

 1 ∂H ∂x ∂r ∂r ∂y ∂r ∂r ∂z ∂x ∂r with (x, nn1 , nn2 , x0 , ω) = A + + +B , ∂x 4π λ ∂ x ∂ nn1 ∂ x ∂ y ∂ nn1 ∂ x ∂ z ∂ nn1 ∂ nn1

   2  ∂r ∂r ∂x 1 ∂H ∂y ∂r ∂r ∂z ∂r ∂y (x, nn1 , nn2 , x0 , ω) = A + + +B , ∂y 4π λ ∂ x ∂ y ∂ nn1 ∂ y ∂ nn1 ∂ y ∂ z ∂ nn1 ∂ nn1  

 2  ∂r ∂r ∂x ∂H ∂r ∂r ∂y ∂z 1 ∂r ∂z (x, nn1 , nn2 , x0 , ω) = A + + +B , ∂z 4π λ ∂ x ∂ z ∂ nn1 ∂ y ∂ z ∂ nn1 ∂ z ∂ nn1 ∂ nn1 2

A=−

kc e−ikc r 3ikc e−ikc r 3e−ikc r + + 2 r r r3

and B = −

ikc e−ikc r e−ikc r − 3 . 2 r r

In Eqs. (4) and (5), nn1 represents the unit outward normal to boundary S at x = (x, y, z), nn2 is the unit outward normal to boundary S at the collocation points x0 = (x0 , y0 , z0 ), the factor a is null for piecewise planar boundary elements and 

r=

(x − x0 )2 + (y − y0 )2 + (z − z0 )2 .

In the case of a point heat source, the incident heat field in Eq. (4) is computed as

Pe−ikc r0 (−ikc r0 − 1) T inc (x, nn2 , xs , ω) = 2λr0 2



 ∂ r0 ∂ y ∂ r0 ∂ z ∂ r0 ∂ x + + . ∂ x ∂ nn2 ∂ y ∂ nn2 ∂ z ∂ nn2

(6)

In the case of a planar heat source, the incident heat field in Eq. (4) is given by:

−Pe−ikc |x−x0 | T inc (x, nn2 , xs , ω) = 2λ

a



 ∂x . ∂ nn2

(7)

b

Fig. 2. Three-dimensional geometry of the problem: (a) point heat source and (b) planar heat source.

1580

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

Fig. 3. 3D view of a 3D null thickness crack hosted in an unbounded solid medium.

Solving the normal-derivative integral equation requires the discretization of S into N planar boundary elements, with a nodal point in each element. A Gaussian quadrature scheme is used to integrate elements when they are not the loaded element. However, for the loaded element the integrands exhibit a singularity (hypersingular element). The hypersingular integrals were solved using the expressions described and validated by Tadeu et al. [19] to solve singular and hypersingular integrals. 3. Results and discussion 3.1. Cases studied The proposed model is used to simulate heat diffusion around a 3D crack with null thickness which is represented by plane S in Fig. 3. The crack is a heterogeneity lodged in a uniform unbounded solid medium of density ρ = 2300 kg m−3 , thermal conductivity λ = 1.40 W m−1 ◦ C−1 and specific heat c = 880 J kg−1 ◦ C−1 . The transient temperature response and phase contrast images corresponding to the resulting thermal waves are calculated along three grids of receivers (G1, G2 and G3 shown in Fig. 3). The grid G1 emulates the surface of an object from which thermographic data is obtained using IRT. Various cases with changing defect shape, size, inclination and position relative to the energy source and to the surface are studied in order to investigate the influence that defect characteristics and experimental parameters have on heat diffusion and phase contrast image results. All computations are performed using both a point heat source and a planar heat source and the final results are compared. The initial case study refers to a vertical defect with width l1 = 0.2 m and height l2 = 0.2 m, placed parallel to the yz plane at x = 0.6 m. A point heat source is located at the origin (xs = 0.0 m, ys = 0.0 m, zs = 0.0 m), or in the case of a planar source, it is placed parallel to the yz plane, at xs = 0.0 m. Grids G1 and G2 are placed parallel to the yz plane, G1 over a plane at x = 0.5325 m (0.0675 m distanced from the crack) and grid G2 over the plane at x = 0.6725 m (see Fig. 4b). G3 is placed parallel to the xy plane at z = 0.0 m. In all three grids, receivers are spaced at equal intervals of x = 0.005 m, y = 0.00625 m and z = 0.006 m. In order to assess the influence of the size of a crack, three different geometries were considered (0.1 × 0.1 m2 , 0.2 × 0.2 m2 and 0.3 × 0.3 m2 ). A 0.05 × 0.2 m2 geometry was also modeled in order to study the influence of changing the shape of a crack. By tilting the vertical crack around an axis parallel to the z axis at an angle of α = 15° and α = 40°, the influence of the position of the crack (in terms of inclination) was also analyzed. To further comprehend how the outline of each defect is translated into the temperature distribution and the phase images, results were also analyzed for specific receivers placed at the center (REC A) and at the extremities of the different defects (REC B to REC F). Fig. 4a shows a front view of the case study systems with the different defect sizes that were considered and the specific placement of receivers REC A to REC F. Fig. 4b illustrates a vertical cross section of the systems modeled and shows the tilting of the defect by an angle α . The changes produced by moving the crack closer to the simulated object’s surface, meaning the grid G1, were studied by moving G1 0.04 m closer to the crack, to x = 0.5725 m. By moving the point source to (xs = 0.4525 m, ys = 0.0 m, zs = 0.0 m), or the planar source to xs = 0.4525 m, the effects of changing the distance between the heat source and the defect were studied. In addition to the geometries defined earlier, results were also computed for the case of an infinitely sized crack (named in this study as the reference case). This was done by means of an image source technique, in which responses are computed by placing an additional virtual source as the mirror image of the real source in relation to the surface of the crack. 3.2. Simulation results This section presents the results in terms of heat field in the time domain and thermal wave phase.

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

a

1581

b

Fig. 4. Geometry of the systems modeled: (a) front view with placement of receivers REC A to REC F and (b) vertical cross section.

Fig. 5. Snapshots of temperature distribution, in °C, at t = 20 h: (a) in the presence of a 0.2 × 0.2 m2 vertical crack; (b) medium without cracks and (c) temperature differences between results in (a) and (b). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

3.2.1. Temperature in time The time-domain temperature response was obtained by applying a numerical inverse fast Fourier transform in the frequency domain. Aliasing was dealt by introducing complex frequencies with a small imaginary part ωc = ω − iη (where η = 0.7ω and ω is the frequency step). This shift was subsequently taken into account in the time domain by means of an exponential window (eη t ) applied to the response. The source can have any time variation and the frequency domain solution can be determined by applying a time Fourier transformation to a frequency range from 0.0 Hz to quite high frequencies. However, there is no need to compute the highest frequencies in the range since the heat response falls rapidly for higher frequencies. The 0.0 Hz frequency is the static response. The use of complex frequencies allows this response to be computed since the argument of the Hankel functions in the integral equations is −iη, i.e., other than zero. Computations were performed in the frequency range [0.0, 1.024 × 10−3 ] Hz with a frequency increment of  f = 0.5 × 10−6 Hz, resulting in a time window of 1/(0.5 × 10−6 ) s. The imaginary part of the frequency is η = 0.7 × 2 × π ×  f . The medium was assumed to be at an initial temperature of 20.0 °C. The source emitted energy from instant t = 0.5 h to 1.5 h following a rectangular heating function with amplitude P. This amplitude was defined so that a maximum temperature increase of 15.0 °C was registered at the center of grid G1, at REC A (see Fig. 4) for the case when the crack is not present in the medium. An amplitude of P = 1843.2 was employed in the case of the point heat source, and an amplitude of P = 205.9 in the case of a planar source. In order to study the effects of the presence of a defect on the temperature diffusion pattern recorded by IRT, calculations were made for a medium with no cracks and a medium with an embedded crack. The difference between both emulated the temperature differences that appear when IRT is used to detect defects. This procedure is illustrated in Figs. 5 and 6. Fig. 5 shows a snapshot of the time domain temperature distribution recorded on the grid of receivers placed at xG1 = 0.5325 m, xG2 = 0.6725 m and zG3 = 0.0 m, taken at the instant t = 20 h for the case of a 0.2 × 0.2 m2 crack placed vertically at x = 0.6 m, subject to heat generated by a 3D point source located at (xs = 0.0 m, ys = 0.0 m, zs = 0.0 m). The results on the

1582

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

Fig. 6. Snapshots of temperature distribution, in °C, at t = 58.1 h: (a) in the presence of a 0.2 × 0.2 m2 vertical crack; (b) medium without cracks and (c) temperature differences between results in (a) and (b). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

0.1

fb

(rad)

0 -7 10

10

-6

10

-5

10

-4

10

-3

10

-2

-0.1

fch

-0.2

REC A (3D) REF (3D) REC A (1D) REF (1D)

fch -0.3

Frequency (Hz) Fig. 7. Phase contrast curves obtained in REC A for a 0.2 × 0.2 m2 vertical crack placed at x = 0.6 m and for an infinitely sized crack (curve REF) subjected to a point (3D) or to a planar (1D) heat source placed at the origin.

left are for the case when there is a defect within the medium, while the central plot shows the results obtained when the medium is sound (no defect). The difference between the two results is given in the plot on the right. Fig. 6 presents a snapshot of the time domain temperature distribution taken at the instant t = 58.1 h for the same defect, only it was submitted to a planar energy source placed at xs = 0.0 m. Both figures were taken at the instant when the thermal response peaked (a temperature of 35° is recorded at REC A when the medium has no defect). In Figs. 5 and 6 we can see that both heat sources are able to produce significant temperature differences when there is a defect, however the thermal pattern observed in the two cases diverges. It can be said that the point source produces a more concentrated thermal response, while the planar source gives a more disperse response. Furthermore, the temperature difference generated by the point source is more significant and makes the defect more visible.

3.2.2. Phase contrast results Thermal wave phase is calculated directly in the frequency domain. Computations are done considering harmonic sources with unitary amplitude, corresponding to an ideal Dirac pulse in the time domain which has an infinite flat spectrum in the frequency domain. The thermal wave phase is the arctangent of the quotient of the imaginary part of the frequency domain temperature response divided by the real part. Phase contrast is the difference between the thermal wave phase responses obtained for a sound medium and for a defected medium (with an embedded crack). Fig. 7 shows the phase contrast results obtained for the receiver located at the center of the grid G1 (REC A) for the initial case study and the reference case (infinite crack), subject to a 3D and a 1D heat source. All curves show a defined peak in the frequency spectrum. The frequency at which the maximum phase contrast |φ max | occurs is known in defect detection studies as the characteristic frequency fch . The blind frequency fb occurs when phase contrast is null, indicating the defect detection threshold. Fig. 7 shows that, for the initial case study, values for fch and fb are similar for both point and planar energy sources. However, it may be observed that the maximum phase contrast is slightly higher when subject to a planar heat source. In the case of a crack of infinite size, the REF phase contrast curves lead to lower fch values and the difference between point and planar

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

0 -7 10

(rad)

b

0.1

10

-6

10

-5

10

-4

10

-3

10

-0.1 REF 2 0.3x0.3 m 2 0.2x0.2 m 2 0.1x0.1 m 2 0.05x0.2 m

-0.2

0.1

0 -7 10

-2

(rad)

a

1583

10

-6

10

-5

10

-4

10

-3

10

-2

-0.1 REF 2 0.3x0.3 m 2 0.2x0.2 m 2 0.1x0.1 m 2 0.05x0.2 m

-0.2

-0.3

-0.3

Frequency (Hz)

Frequency (Hz)

Fig. 8. Phase contrast curve results obtained at REC A for varying crack sizes, for: (a) 3D point heat source and (b) 1D planar heat source.

a 0.1

10

-6

10

-5

10

-4

10

-3

10

0 -7 10

-2

(rad)

(rad)

0 -7 10

b 0.1

-0.1 REF 2 0.3x0.3 m 2 0.2x0.2 m 2 0.1x0.1 m 2 0.05x0.2 m

-0.2

10

-6

10

-5

10

-4

10

-3

10

-2

-0.1 REF 2 0.3x0.3 m 2 0.2x0.2 m 2 0.1x0.1 m 2 0.05x0.2 m

-0.2

-0.3

-0.3

Frequency (Hz)

Frequency (Hz)

Fig. 9. Phase contrast curve results obtained at REC D for varying crack sizes, for: (a) 3D point heat source and (b) 1D planar heat source.

sources increases. In none of the curves did the blind frequency fb change significantly, therefore neither did the defect detection threshold.

3.3. Influence of crack size The graphs in Figs. 8 and 9 show the phase contrast curves obtained at receivers REC A (x = 0.5325 m, y = 0.0 m, z = 0.0 m) and REC D (x = 0.5325 m, y = 0.0 m, z = 0.15 m) for the four different crack sizes: 0.3 × 0.3 m2 ; 0.2 × 0.2 m2 ; 0.1 × 0.1 m2 and 0.05 × 0.2 m2 . The results show that as crack size increases there is a decrease in the characteristic frequency fch and an increase in the corresponding maximum phase contrast absolute value |φ max |. Fig. 8 also shows that the amplitude of phase contrast decreases when the shape of the crack changes from the 0.1 × 0.1 m2 geometry to the 0.05 × 0.2 m2 one. As expected, as the size of the crack increases, results become closer to those obtained for the infinite size crack curve REF. In no case does the blind frequency fb suffer significant alteration and the detection threshold changes very little for different crack sizes. By comparing Fig. 8 with Fig. 9 it may be seen that the maximum phase contrast is slightly greater when subject to a planar heat source and that, as the crack size increases, so does this difference between point and planar sources. As shown in Fig. 9, for receiver REC D (placed outside the limits of the 0.1 × 0.1 m2 , 0.05 × 0.2 m2 and 0.2 × 0.2 m2 cracks) phase contrast is significantly reduced for receivers placed outside the limits of the defect under study, indicating a potential for detecting and outlining a defect’s boundaries. Figs. 10 and 11 contain snapshots of phase contrast images computed at the grid of receivers G1 for the 0.2 × 0.2 m2 vertical crack. In the central plot the image is of frequency fch which corresponds to maximum phase contrast |φ max | at receiver REC A, located at the center of the crack and indicated in these plots by a red circular mark. For both the point and planar heat sources, the phase contrast snapshots clearly demonstrate the ability to detect the presence of the crack and assess its size. It can be clearly seen that outside the limits of the cracks the phase contrast tends to zero because the effect of the reflected heat field is diminished.

1584

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

Fig. 10. Phase contrast images of grid G1 for a 0.2 × 0.2 m2 vertical crack subjected to a 3D point heat source. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Fig. 11. Phase contrast images of grid G1 for a 0.2 × 0.2 m2 vertical crack subjected to a 1D planar heat source. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

0 -7 10

(rad)

b

0.1

10

-6

10

-5

10

-4

10

-3

10

-0.1 REC A - xs=0.4525 m REC A - xs=0 m REF - xs=0.4525 m REF - xs=0 m

-0.2

-0.3

0.1

0 -7 10

-2

(rad)

a

10

-6

10

-5

10

-4

10

-3

10

-2

-0.1 REC A - xs=0.4525 m REC A - xs=0 m REF - xs=0.4525 m REF - xs=0 m

-0.2

-0.3

Frequency (Hz)

Frequency (Hz)

Fig. 12. Phase contrast results at REC A when heat source is moved from xs = 0.0 m to xs = 0.4525 m, for: (a) point heat source and (b) planar heat source.

3.4. Influence of source position The influence of changing the source position was demonstrated by moving the source away from the origin, closer to the defect, to xs = 0.4525 m, thus reducing the initial source-defect distance from 0.6 m to 0.1475 m. Phase contrast results computed at REC A(x = 0.5325 m, y =0.0 m, z = 0.0 m) in the presence of a 0.2 × 0.2 m2 vertical crack are shown in Fig. 12 to illustrate the main conclusions. Fig. 12a indicates that, for a receiver located at the center (REC A), when the distance from a point source is shortened, the characteristic frequency fch value decreases, as well as the corresponding maximum absolute phase contrast value |φ max |. The blind frequency fb does not change significantly, therefore the detection threshold is not influenced by changing the source distance either. For a planar energy source, the values of fch and fb are not influenced by changing the source distance, as shown in Fig. 12b. In this case, the heat front does not lose power as the thermal waves propagate and the difference between the incident and reflected field does not change when its source position changes. 3.5. Influence of crack depth In order to study the influence of the depth at which the defect is located, the grid of receivers G1 is moved along the x axis. The graphs in Fig. 13 show the phase contrast results for receiver REC A placed at xG1 = 0.5325 m and for REC A placed at

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

b

0.2 REC A - x=0.5725 m REC A - x=0.5325 m REF - x=0.5725 m REF - x=0.5325 m

(rad)

0.1

0 -7 10

10

-6

10

-5

10

-4

10

0.2 REC A - x=0.5725 m REC A - x=0.5325 m REF - x=0.5725 m REF - x=0.5325 m

0.1

-3

10

-2

0 -7 10

(rad)

a

1585

-0.1

-0.1

-0.2

-0.2

10

-6

10

-5

10

-4

10

-3

10

-2

-0.3

-0.3

Frequency (Hz)

Frequency (Hz)

Fig. 13. Phase contrast curves for a 0.2 × 0.2 m2 crack and reference curves for infinite size crack at REC A (0.5325 m,0.0 m, 0.0 m) and REC A (0.5725 m,0.0 m, 0.0 m) subjected to a: (a) 3D point heat source and (b) 1D planar heat source.

b

0.05

0 -7 10

10

-6

10

-5

10

-4

10

-3

10

-0.05

-2

0.05

0 -7 10

(rad)

(rad)

a

10

-6

10

-5

10

-4

10

-3

10

-2

-0.05 REC E - =40º REC E - =15º REC E - =0º

REC E - =40º REC E - =15º REC E - =0º

-0.10

-0.10

Frequency (Hz)

Frequency (Hz)

Fig. 14. Phase contrast curves modeling results obtained at REC E for a tilted 0.2 × 0.2 m crack subjected to a: (a) 3D point heat source and (b) 1D planar heat source. 2

xG1 = 0.5725 m, in the presence of a 0.2 × 0.2 m2 crack. This simulates a reduction in the depth of the defect. Fig. 13 also includes the reference curves with the receivers in the same positions. For both the 0.2 × 0.2 m2 and the infinite size crack, the results show that reducing the depth of the defect increases the characteristic frequency values fch , as well as the corresponding phase contrast peak |φ max |. Furthermore, as the buried depth is reduced, it can also be said that the blind frequency fb increases as well. Comparison between planar and point heat source results show that, as the depth of the defect is reduced, results for a receiver located at the center of the crack become more similar in terms of fch and |φ max |.

3.6. Influence of crack inclination The graphs in Fig. 14 illustrate the influence of tilting a 0.2 × 0.2 m2 crack and show the phase contrast results obtained at REC E (x = 0.5325 m, y = 0.05 m, z = 0.05 m) for three different crack inclinations: α = 0◦ , α = 15◦ and α = 40◦ . Fig. 15 shows the phase contrast results obtained at REC F (x = 0.5325 m, y = −0.05 m, z = 0.05 m). When comparing planar and point heat sources, the graphs in Figs. 14 and 15 reveal similar fch and fb results. In REC E, located in front of the upper part of the crack, the maximum phase contrast is slightly greater for the planar heat source. In REC F, the peak phase contrast occurs for higher characteristic frequency values fch , which is in accordance with the results reported in the previous section (reducing the depth of the defect increases the characteristic frequency values fch as well as the corresponding phase contrast peak |φ max |). Figs. 16 and 17 show the phase contrast images recorded at the grid of receivers G1 for the 0.2 × 0.2 m2 crack, when tilted by 40°. Snapshots are taken at three different frequencies. In the center plot, a snapshot is taken at the frequency corresponding to the maximum phase contrast recorded at REC E and REC F. Fig. 16 shows that when a crack is tilted, the amplitude of phase contrast is reduced at the lower frequencies and the maximum phase contrast occurs at higher frequencies. Furthermore, as seen for α = 40◦ , in the low frequency spectrum the response is inverted and a positive difference is recorded at REC F. The images indicate that, when the crack is tilted so that the lower part of the defect is closer to the surface and the upper part is further away, the results show a mixed response. Therefore, phase contrast images given at a specific frequency can no

1586

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

0 -7 10

(rad)

b 0.05

0.05

10

-6

10

-5

10

-4

10

-3

10

0 -7 10

-2

(rad)

a

-0.05

10

-6

10

-5

10

-4

10

-3

10

-2

-0.05

-0.10

-0.10

REC F - =40º REC F - =15º REC F - =0º

REC F =40º REC F - =15º REC F - =0º

-0.15

-0.15

Frequency (Hz)

Frequency (Hz)

Fig. 15. Phase contrast curves modeling results obtained at REC F for a tilted 0.2 × 0.2 m crack subjected to a: (a) 3D point heat source and (b) 1D planar heat source. 2

Fig. 16. Phase contrast images of grid G1 for a 0.2 × 0.2 m2 tilted by 40°, subjected to a 3D point heat source. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Fig. 17. Phase contrast images of grid G1 for a 0.2 × 0.2 m2 tilted by 40°, subjected to a 1D planar heat source. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

longer provide conclusive information regarding the size and shape of the crack. A wider range of frequencies should be swept to enable a full assessment of the defect’s size and shape.

4. Conclusions In this paper, the authors have used a 3D normal-derivative integral equation (TBEM) formulation to simulate heat diffusion by conduction that occurs around a defected unbounded medium. The used simulation model allows for the calculation of temperature and phase contrast results and is intended to make a contribution to defect detection studies using active IRT as an NDT. This formulation was used to simulate temperature and phase contrast in the proximity of a null thickness defect (crack) present in an unbounded space. Parameters which are relevant features in non-destructive testing by active IRT, such as the depth at which the crack is embedded, its size, shape, inclination, and distance to heat source were analyzed, with particular focus on the influence of the type of the energy source being modeled. Namely, 3D point heat source and 1D planar heat source results were compared.

C. Serra et al. / Applied Mathematical Modelling 40 (2016) 1576–1587

1587

Heat diffusion was modeled using an inverse fast Fourier transform applied in the frequency domain and by assuming a rectangular heating source, as used in square pulse thermography (SPT) experiments. Phase contrast results were computed for an ideal Dirac pulse in the time domain, which has an infinite flat spectrum in the frequency domain. Temperature difference and phase contrast images were modeled by computing the heat field generated in the medium containing an embedded crack (defected medium) and in the medium without defects (sound medium), subjected to both planar and point energy sources. The numerical results obtained in this study showed that the characteristics of a crack can have a great influence on active IRT data results. Results proved the usefulness of the proposed formulation in understanding heat diffusion in defected media, both in the time and frequency domains, and in interpreting experimental IRT results. Furthermore, although the results were found to be highly influenced by the 3D nature of the crack, it may be said that they were not that much influenced by the nature of the heat source. In fact, the influence of changing the crack’s size, depth and inclination did not alter significantly when the nature of the heat source was changed. However, planar source results were found to have the advantage of not being influenced by the distance between the defect and the heat source. It is concluded that studying the heat diffusion patterns generated by the proposed models when experimental IRT variables are altered helps to interpret experimental IRT results correctly. Therefore, further applications of the proposed models could be developed to establish experimental setup parameters for active IRT defect detection studies. Acknowledgments The research work presented herein was supported by FEDER funds through the Operational Programme for Competitiveness Factors – COMPETE and QREN, under the Project FCOMP-05-0128-FEDER-022994 and by national funds through the FCT – Portuguese Foundation for Science and Technology, under research project PTDC/ECM/114189/2009. This work was supported in part by the doctoral Grant SFRH/BD/91686/2012 awarded by the FCT. This work is also within the scope of the Energy for Sustainability Initiative of the University of Coimbra. References [1] X. Maldague, Theory and Practice of Infrared Technology for Non-destructive Testing, John Wiley & Sons, 2001. [2] H. Wiggenhauser, Active IR-applications in civil engineering, Infrared Phys. Technol. 43 (2002) 233–238. [3] Ch. Maierhofer, H. Wiggenhauser, A. Brink, M. Röllig, Quantitative numerical analysis of transient IR-experiments on buildings, Infrared Phys. Technol. 46 (2004) 173–180. [4] X. Maldague, Application of infrared thermography in nondestructive evaluation, Trends in Optical Nondestructive Testing, Pramod Rastogi (Ed.), 2000, pp. 591-609, (invited chapter). [5] X. Maldague, Introduction to NDT by active infrared thermography, Mater. Eval. 6 (2002) 1060–1073. [6] C. Ibarra-Castanedo, D. González, M. Klein, M. Pilla, S. Vallerand, X. Maldague, Infrared image processing and data analysis, Infrared Phys. Technol. 46 (2004) 75–83. [7] X. Maldague, S. Marinetti, Pulse phase infrared thermography, Appl. Phys. 79 (1996) 2694–2698. [8] R.W. Arndt, Square pulse thermography in frequency domain as adaptation of pulsed phase thermography for qualitative and quantitative applications in cultural heritage and civil engineering, Infrared Phys. Technol. 53 (2010) 246–253. [9] K. Srinivas, A.O. Siddiqui, J. Lahiri, Thermographic inspection of composite materials, in: Proceedings of the National Seminar on Non-Destructive Evaluation, Hyderabad, 2006, 7–9 December. [10] K.J. Bathe, Numerical Methods in Finite Element Analysis, Prentice-Hall, New Jersey, 1976. [11] M.N. Ozisik, Finite Difference Methods in Heat Transfer, CRC Press Inc., USA, 1994. [12] Gh. Juncu, Unsteady conjugate forced convection heat/mass transfer from a finite flat plate, Int. J. Therm. Sci. 47 (2008) 972–984. [13] C.A. Brebbia, J.C. Telles, L.C. Wrobel, Boundary Elements Techniques: Theory and Applications in Engineering, Springer-Verlag, Berlin/ New York, 1984. [14] Y. Ochiai, Steady heat conduction analysis in orthotropic bodies by triple-reciprocity BEM, Comput. Model. Eng. Sci. 2 (2001) 435–446. [15] B. Šarler, Towards a mesh-free computation of transport phenomena, Eng. Anal. Bound. Elem. 26 (2002) 731–738. [16] C. Serra, A. Tadeu, J. Prata, N. Simões, Application of 3D heat diffusion to detect embedded 3D empty cracks, Appl. Therm. Eng. 61 (2013) 596–605. [17] T.J. Rudolphi, The use of simple solutions in the regularisation of hypersingular boundary integral equations, Math. Comput. Model. 15 (1991) 269–278. [18] P. Amado Mendes, A. Tadeu, Wave propagation in the presence of empty cracks in an elastic medium, Comput. Mech. 38 (2006) 183–199. [19] A. Tadeu, J. Prata, N. Simões, Closed form integration of singular and hypersingular integrals in 3D BEM formulations for heat conduction, Math. Probl. Eng. 2012 (2012), Article ID 647038, 21 pages, doi:10.1155/2012/647038.