Application of the topological derivative to post-processing infrared time-harmonic thermograms for defect detection

Abstract

This paper deals with active time-harmonic infrared thermography applied to the detection of defects inside thin plates. We propose a method to post-process raw thermograms based on the computation of topological derivatives which will produce much sharper images (namely, where contrast is highly enhanced) than the original thermograms. The reconstruction algorithm does not need information about the number of defects, nor the size or position. A collection of numerical experiments illustrates that the algorithm is highly robust against measurement errors in the thermograms, giving a good approximation of the shape, position and number of defects without the need of an iterative process.

1 Introduction

Nowadays, structural health monitoring is of paramount interest in a wide range of fields. Structures may be monitored to evaluate the damage suffered after an earthquake, where information is needed almost in real time. Airplanes are periodically inspected looking for critical defects due to fatigue, halting them for long times and the list goes on. This means that faster, more reliable, and cheaper inspection methods are eagerly sought.

The standard methods used in the airspace industry are based on ultrasonic tests, which are extremely accurate but time-consuming as the inspected area per unit of time is quite small. In this paper we study the capabilities of the infrared thermography as an alternative inspection method.

An infrared thermogram is an image of a surface where each pixel has been assigned a temperature. This temperature is estimated by means of an infrared thermal camera, which senses the amount of infrared radiation coming from every direction in view (see [21] for further details). Thermography has been successfully used in many fields [43], including not only aircraft structural health monitoring [12, 16, 60], but also building insulation inspection [36], medical diagnosis [15, 33, 44], nano-composites monitoring [35], drying problems [57], et cetera.

Motivated by some experiments [61, 62] that were carried out at the Spanish Aeronautical Technology Center (CTA), we consider in this paper thin rectangular plates, where the illuminated side coincides with the one where the temperature is measured, i.e., where the thermogram is taken, as sketched on Fig. 1. This could represent part of the fuselage of an airplane which can only be inspected from the outside.

Figure 1

Experiments setup. Both the camera and the lamp are situated at the same side of the inspected plate

Compared with the ultrasonic test, the main advantage of the infrared thermography is that it can take measurements over a whole surface at the same time. On top of that it is a contactless (or nearly contactless) method, which can be very important when the unknown state of damage of the structure can be dangerous for the person that performs the inspection [20, 62].

The main disadvantage is that being heat transfer a short range phenomenon, the amount of information that can be extracted from the thermogram is less than with the same amount of data in an ultrasonic case. Another problem is that the temperature on the surface of the structure is more easily influenced by the surrounding temperature and humidity, and consequently, it should be carefully used in controlled environments. This paper is devoted to active thermography, which consists in heating the medium actively, in contrast to passive thermography, where the medium is not excited. There will be a heat source (in our case an infrared lamp generating time-harmonic excitations) controlled by the person performing the inspection. This way the amount of user-controlled conditions is increased compared with the passive infrared thermography.

The usage of time-harmonic thermal excitations, i.e. of thermal waves, goes back to A.J. Ångström [3] who was the pioneer in the theoretical and experimental study of this kind of waves. However, the main advances in their use are rather recent, as reported in [1, 40] and references therein. As indicated in [5], the use of thermal waves is preferable to other thermal imaging methods, in particular, to steady-state infrared thermography.

From the thermal point of view, defects are manifested as sharp changes in the thermal properties of the structure. These changes affect the way temperature is distributed on the whole structure. From the measured temperature in one of the outer surfaces (i.e., from a thermogram) one tries to devise the presence of such internal defects. However, simple inspection or simple noise-filtering methods are not enough in most of the cases to identify the defects, and post-processing thermograms via very efficient tools is required. This will be illustrated in the forthcoming Figs. 5 and 6.

From the mathematical point of view [32, 45], following J. Hadamard’s definition, a problem is said to be well-posed if satisfies the following three conditions: (i) there is at least one solution (existence), (ii) there is at most one solution (uniqueness) and (iii) the solution depends continuously on the data (stability). Otherwise, the problem is called ill-posed. The problem of obtaining the expected temperature measurements on a surface of the inspected medium, given the boundary conditions, initial conditions, and thermal coefficients is called the direct problem. It is a well-posed problem in the Hadamard’s sense: it has a unique solution that depends continuously on the data. In particular, small perturbations in the inputs (thermal coefficients, boundary conditions, etc.) give rise to small perturbation in the outputs (temperature on the surfaces). Furthermore, the heat equation acts as a smoothing operator, meaning that in general the size of the perturbation in the output will be smaller than the perturbation in the inputs. The greater the time difference between excitation and measurement or distance between the discontinuity of the thermal properties and the measuring location, the smaller the perturbation in the outputs. This means that the inverse problem, i.e. finding the thermal coefficients distribution (which will give away the presence of defects) given the boundary and initial conditions, and the temperature distribution on part of its sides is an ill-posed problem. Small differences in the input (the temperature on part of its sides) can be the consequence of really different thermal coefficient distributions, and in the extreme case, more than one solution can be possible, that is, more than one distribution of thermal coefficients can give rise to the same temperature distribution on one of the sides of the inspected medium.

On top of that, if we take into account that we do not know with infinite precision the values of the thermal coefficients of the material, nor the rest of the conditions of the experiment, and that the thermogram is not the temperature on the side of the structure but an estimation which contains measurement errors, we can almost assure that the thermogram will not be a solution of the mathematical problem that models the direct problem for any given defect distribution. That is, there is in general no solution to the inverse problem as it is formulated, because, in general we cannot find a set of defects for which the solution to the direct problem on the measurement side exactly agrees with the given thermogram.

Because of this, in this paper, the inverse problem will be reformulated in a weakened form, as a minimization problem of a shape functional, where we will try to find the set of defects that gives rise to a temperature distribution on the inspected side which minimizes its $L^{2}$ distance to the given thermogram. This kind of reformulation is not new in the literature and there are many different methods to solve this minimization problem. Most of these reconstruction methods are iterative, namely, they begin with an “a priori” distribution of defects and proceed by deforming them step by step trying to minimize the aforementioned distance. Early methods were based on classical shape deformations, using perturbations of the identity operator, as the ones proposed in [27, 34, 55], and more sophisticated ones deal with the computation of the so-called shape derivative [13, 59], which provides the direction of deformation that minimizes this distance the most. However, the main drawback of these methods is that they need to start with an initial guess having the true number of defects, which in practise is also an unknown for the inverse problem. This difficulty was solved by using a new type of deformations that allow for topological changes: the level-set methods [39, 49, 58]. In any case, these methods, being iterative, are generally very time-consuming and require a rather large number of iterations unless a good initial guess of the defects is used.

In our paper we study a fairly recent reconstruction method based on the computation of the topological derivative of the shape functional [18, 47, 59]. The topological derivative at a given point is the sensitivity of the shape functional to introducing an infinitesimal defect at that point of the domain. It is a scalar function that can be interpreted (as we will do) as an indicator function, that classifies each of these points as belonging to a defect or to the surrounding medium. The main advantages of this method are: (i) it is a one-step method, (ii) it does not require any a priori information about the number, size, shape or location of the defects, (iii) it has a very small computational cost. Topological derivatives have been successfully used for solving a wide variety of problems related with defect identification in many fields, like in acoustics [6, 9, 19], electromagnetism [37, 42], elasticity [4, 25, 41], and electrical impedance tomography problems [2, 8, 11], to mention a few. Related work dealing with thermal problems was published in [7, 10], where a two dimensional unsteady thermal propagation problem in an unbounded media was studied. In this paper we deal with the identification of defects inside thin three dimensional plates by post-processing time-harmonic thermograms via the topological derivative method. The simplification to the steady-state case was very recently considered in [52, 53] for two-dimensional plates and in [28] for three-dimensional ones. The two-dimensional time-harmonic case was already studied in [53]. The current paper extends our short paper published in the Proceedings of the 20th European Conference on Mathematics and Industry [54]. The study of the general unsteady problem (with no time-harmonic dependence) will be the subject for our future work.

Throughout a gallery of numerical experiments, it will be illustrated that reconstructions highly depend on both, the frequency of excitation and on the location of the lamp. Defects that are located close to the lamp are easier to identify, almost independently of their size. The frequency of excitement makes more or less visible the presence of defects depending on their size. To overcome these difficulties, we propose to simultaneously post-process several thermograms corresponding to different frequencies and different locations of the lamp by considering a weighted sum of the individual $L^{2}$ distances, in the spirit of the previous work [19] where measurements at different frequencies in an acoustic problem were processed. We will see that this promotes a clear improvement on the reconstructions.

In comparison with simple steady excitations (considered in [28]), time-harmonic excitations have two main advantages. First, when considering steady excitations, only variations in the thermal conductivity are relevant, while changes in the density or in the specific heat capacity do not play a role. In case of time-harmonic excitations, the three thermal properties are taken into account. Secondly, the usage of different frequencies improves the identification of small defects, that are in general not detected by steady thermograms when in the configuration some defects bigger in size are also present.

For the numerical experiments presented in this paper, only syntectic data (namely, numerically generated) will be considered. Processing experimental data is left as future work.

The remaining of the paper is structured as follows. First we mathematically formulate in detail the direct and inverse problems. Next we introduce the concept of topological derivative, which is central to our inversion method, and present then an efficient way (a closed form formula) to implement it by solving two related heat conduction problems that take place in a healthy plate, i.e., in a plate without defects. The reconstructed defects will be defined by clustering all the points where the topological derivative attains the more pronounced negative values. As a way to assess the method we present numerical results obtained on a representative configuration. Finally, we extract some conclusions and propose possible improvements to the method.

2 Mathematical formulation of the direct and inverse problems

In this section we mathematically describe both, the heat conduction process in a thin plate that models the thermographic inspection basis, and the problem to be solved, namely, the identification of damage inside the plate from measurements of the temperature distribution at one side of such plate.

For simplicity, we assume the thin plate is rectangular and occupies the region

$$\mathscr{R}=(-\ell _{x}/2,\ell _{x}/2)\times (-\ell _{y}/2,\ell _{y}/2) \times (-\ell _{z}/2, \ell _{z}/2),\quad \text{with }\ell _{x}\ll\ell _{y}, \ell _{y}\sim \ell _{z}. $$

The assumption $\ell _{x}\ll\ell _{y}$, $\ell _{y}\sim \ell _{z}$ is motivated by some laboratory experiments performed at the Aeronautical Technology Center (CTA) at Miñano, Spain (see the description of the experiments in [61, 62]).

We will denote by $\varGamma _{\mathrm{front}}$ the side of $\mathscr{R}$ where the thermogram is taken and the radiation from the lamp is received (see Fig. 1), described by the set $\varGamma _{\mathrm{front}}=\{x=-\ell _{x}/2\}\cap \overline{\mathscr{R}}$. The opposite side will be denoted as $\varGamma _{\mathrm{back}}$, namely, $\varGamma _{\mathrm{back}}=\{x=\ell _{x}/2\}\cap \overline{\mathscr{R}}$, and the remaining of the sides will be denoted collectively by $\varGamma _{\mathrm{sides}}=\partial \mathscr{R}\setminus ( \varGamma _{\mathrm{front}}\cup \varGamma _{\mathrm{back}})$. Notice that $\varGamma _{\mathrm{sides}}$ stands for the sides of the plate with the smallest areas (see Fig. 2). The internal defects conform the subset $\mathscr{D}\subset \mathscr{R}$ which is not necessarily connected (it can consist of more than one disjoint defect).

Figure 2

Domain. Mathematical description of the domain

In this paper, we study homogeneous metallic plates where the presence of defects is defined by a sharp change in the value of the thermal properties. The extension to heterogeneous plates only requires small adjustments in the results. The thermal conductivity κ, density ρ, and specific heat capacity c will therefore be described by constant values inside the defects (denoted by the subindex “i”, to recall that are interior parameters) and different constant ones in the remaining of the plate (denoted by the subindex “e”, meaning exterior), that is:

$$\begin{aligned}& \kappa (\textbf{x} )= \textstyle\begin{cases} \kappa _{\mathrm{e}}, & \textbf{x}\in \mathscr{R}\setminus \overline{ \mathscr{D}}, \\ \kappa _{\mathrm{i}}, & \textbf{x}\in \mathscr{D}, \end{cases}\displaystyle \qquad \rho (\textbf{x} )= \textstyle\begin{cases} \rho _{\mathrm{e}}, & \textbf{x}\in \mathscr{R}\setminus \overline{ \mathscr{D}}, \\ \rho _{\mathrm{i}}, & \textbf{x}\in \mathscr{D}, \end{cases}\displaystyle \\& c (\textbf{x} )= \textstyle\begin{cases} c_{\mathrm{e}}, & \textbf{x}\in \mathscr{R}\setminus \overline{ \mathscr{D}}, \\ c_{\mathrm{i}}, & \textbf{x}\in \mathscr{D}. \end{cases}\displaystyle \end{aligned}$$

The temperature distribution $T(\textbf{x},t)$ inside the plate satisfies then the homogeneous linear heat equation

$$ \rho (\textbf{x} )c (\textbf{x} )\partial _{t}T (\textbf{x},t )+ \nabla \cdot \bigl(-\kappa (\textbf{x} ) \nabla T (\textbf{x},t ) \bigr)=0, \quad \textbf{x}\in \mathscr{R}\setminus \partial \mathscr{D}, t>0. $$

On the boundary of the defects $\partial \mathscr{D}$ we impose continuity of temperature and normal thermal fluxes. By using the superscripts “i” and “e” to refer respectively to limit values from the interior and the exterior of the defects, and denoting by n the exterior unit vector perpendicular to the boundary, this is modelled by the boundary conditions

$$ \textstyle\begin{cases} T^{\mathrm{{e}}} (\textbf{x},t ) - T^{\mathrm{{i}}} (\textbf{x},t ) =0,& \textbf{x} \in \partial \mathscr{D}, t>0, \\ \kappa _{\mathrm{e}}\partial _{\textbf{{n}}}T^{\mathrm{{e}}} (\textbf{x},t )- \kappa _{\mathrm{i}}\partial _{\textbf{{n}}}T^{\mathrm{{i}}} (\textbf{x},t )=0,& \textbf{x}\in \partial \mathscr{D}, t>0. \end{cases} $$

These conditions model conduction as the heat transfer method acting inside the defects $\mathscr{D}$, while convection and radiation effects are neglected due to the expected small dimensions of the defects (consisting of air) in comparison with the size of the plate.

On the sides of the plate, there will be heat transfer by convection and radiation between the surface of the plate and the surrounding air at temperature $T_{\mathrm{air}}$. On the one hand, convection will be modelled by Newton’s law of cooling, being the heat transferred per unit area equal to $h (T (\textbf{x},t )-T_{ \mathrm{air}} )$, where $h>0$ is a constant. On the other hand, the radiation term would have two terms, the emitted heat $\varepsilon \sigma T^{4} (\textbf{x},t )$, and the absorbed heat $\alpha \sigma T_{\mathrm{air}}^{4}$, where $\sigma >0$ is the Stefan–Boltzmann constant and $\varepsilon >0$ and $\alpha >0$ are the hemispherical emissivity and absorptance of the surface, which for simplicity are assumed to be constant [30]. Due to the small thickness of the plate and the small difference between the temperature on the surfaces of the plate we can perform some simplifications in the model. First, we consider thermal insulation in the smallest sides $\varGamma _{\mathrm{sides}}$, which is modelled by imposing a zero normal flux (a homogeneous Neumann boundary condition):

$$-\kappa _{e}\partial _{\textbf{n}}T(\textbf{x},t)=0,\quad \textbf{x} \in \varGamma _{\mathrm{sides}}, t>0. $$

In the biggest sides, $\varGamma _{\mathrm{front}}$ and $\varGamma _{\mathrm{back}}$, we simplify the emitting term by linearising it around $T_{\mathrm{air}}$, resulting in a term of the form

$$4\varepsilon \sigma T_{\mathrm{air}}^{3}T (\textbf{x},t )- { (3 \varepsilon +\alpha )\sigma T_{ \mathrm{air}}^{4}}, \quad \textbf{x} \in \varGamma _{\mathrm{front}}\cup \varGamma _{\mathrm{back}}, t>0, $$

and taking into account both radiation and convection we have

$$\bigl(4\varepsilon \sigma T_{\mathrm{air}}^{3}+h \bigr)T ( \textbf{x},t )- { \bigl( (3\varepsilon +\alpha )\sigma T_{ \mathrm{air}}^{3}+h \bigr)T_{\mathrm{air}}}, \quad \textbf{x}\in \varGamma _{\mathrm{front}}\cup \varGamma _{\mathrm{back}}, t>0, $$

which we will rewrite as

$$a T (\textbf{x},t )-b, $$

where ${a = 4\varepsilon \sigma T_{\mathrm{air}}^{3}+h}$ and ${b= ( (3\varepsilon +\alpha )\sigma T _{\mathrm{air}}^{3}+h )T_{\mathrm{air}}}$ are two positive constants. In addition, on $\varGamma _{\mathrm{front}}$ we have a term accounting for the radiation coming from the lamp. This will be modeled as an isotropic radiator situated at a point s (with first component $s_{x}<-\ell _{x}/2$), and with intensity $I (t )$. If $\theta _{\mathrm{inc}} (\textbf{x} )$ is the incidence angle for a point $\textbf{x}\in \varGamma _{\mathrm{front}}$ with respect to the source s, then we have that the incoming heat can be described by the function

$$ q_{\textbf{s}} (\textbf{x},t )=\frac{I (t )}{4 \pi } \frac{\cos (\theta _{\mathrm{inc}} (\textbf{x} ) )}{ \vert \textbf{x}-\textbf{s} \vert ^{2}},\quad \textbf{x}\in \varGamma _{\mathrm{front}}, t>0. $$

(1)

With all of this in mind the full non-stationary model would be:

$$\textstyle\begin{cases} \rho _{\mathrm{e}}c_{\mathrm{e}}\partial _{t}T (\textbf{x},t )- \kappa _{\mathrm{e}}\Delta T (\textbf{x},t )=0, & \textbf{x}\in \mathscr{R}\setminus \overline{\mathscr{D}}, t>0, \\ \rho _{\mathrm{i}}c_{\mathrm{i}}\partial _{t}T (\textbf{x},t )- \kappa _{\mathrm{i}}\Delta T (\textbf{x},t )=0, & \textbf{x}\in \mathscr{D}, t>0, \\ T^{\mathrm{e}} (\textbf{x},t ) = T^{\mathrm{i}} (\textbf{x},t ), & \textbf{x}\in \partial \mathscr{D}, t>0, \\ -\kappa _{\mathrm{e}}\partial _{\textbf{{n}}}T^{\mathrm{e}} (\textbf{x},t )=- \kappa _{\mathrm{i}}\partial _{\textbf{{n}}}T^{\mathrm{i}} (\textbf{x},t ),& \textbf{x}\in \partial \mathscr{D}, t>0, \\ -\kappa _{\mathrm{e}}\partial _{\textbf{n}}T (\textbf{x},t )=0, & \textbf{x}\in \varGamma _{\mathrm{sides}}, t>0, \\ -\kappa _{\mathrm{e}}\partial _{\textbf{n}}T (\textbf{x},t )=a T (\textbf{x},t )-b, & \textbf{x}\in \varGamma _{ \mathrm{back}}, t>0, \\ -\kappa _{\mathrm{e}}\partial _{\textbf{n}}T (\textbf{x},t )=a T (\textbf{x},t )-b-\alpha q_{\textbf{s}} (\textbf{x},t ), & \textbf{x}\in \varGamma _{\mathrm{front}}, t>0, \\ T (\textbf{x},0 )=T_{0} (\textbf{x} ), & \textbf{x}\in \mathscr{R}, \end{cases} $$

where $T_{0}(\textbf{x})$ is the initial temperature distribution in the plate.

As already mentioned in the introduction, we are interested in the case where the intensity of the lamp is a time-harmonic function. That would require an experiment where the intensity of the lamp has an average value and oscillates around it, that is:

$$I (t ) = \hat{I}+\mathfrak{Re} \bigl(\mathcal{I}e^{-i \omega t} \bigr), $$

where $\mathcal{I}$ is a complex number carrying information about both amplitude and phase of the oscillation, Î stands for the mean value and $\omega >0$ is the frequency of the oscillation. The corresponding heat coming from the lamp can also be expressed as:

$$q (\textbf{x},t ) = \hat{q}_{\textbf{s}} (\textbf{x} ) + \mathfrak{Re} \bigl(\mathcal{Q} _{\textbf{s}} (\textbf{x} )e^{-i\omega t} \bigr), $$

where the function $\mathcal{Q}_{\textbf{s}}$ is defined as its counterpart time-dependent one (see (1)), namely,

$$\mathcal{Q}_{\textbf{s}} (\textbf{x} )=\frac{\mathcal{I}}{4 \pi } \frac{\cos (\theta _{\mathrm{inc}} (\textbf{x} ) )}{ \vert \textbf{x}-\textbf{s} \vert ^{2}}. $$

Due to the linearity of the problem, the temperature can then be decomposed in the sum of three terms,

$$T (\textbf{x},t )=\tilde{T} (\textbf{x},t )+ \hat{T} (\textbf{x} )+\mathfrak{Re} \bigl(\mathcal{T} (\textbf{x} )e ^{-i\omega t} \bigr), $$

where T̃ is a transient term that accommodates the initial condition and vanishes in a short time, T̂ is a time-steady temperature distribution and $\mathcal{T}:\mathscr{R}\to \mathbb{C}$ is the term corresponding to the complex amplitude of oscillation of the temperature at every point. In this paper, we focus on the effects of the latter term $\mathcal{T}$. Taking into account that the initial conditions are fulfilled by the transient term T̃, and that the steady terms are absorbed by T̂, the time-harmonic one ${\mathcal{T}}$ solves the following (stationary) problem:

$$ \textstyle\begin{cases} \kappa _{\mathrm{e}}\Delta \mathcal{T} (\textbf{x} )+i \omega \rho _{\mathrm{e}}c_{\mathrm{e}}\mathcal{T} (\textbf{x} )=0, & \textbf{x}\in \mathscr{R}\setminus \overline{\mathscr{D}}, \\ \kappa _{\mathrm{i}}\Delta \mathcal{T} (\textbf{x} )+i \omega \rho _{\mathrm{i}}c_{\mathrm{i}}\mathcal{T} (\textbf{x} )=0, & \textbf{x}\in \mathscr{D}, \\ \mathcal{T}^{\mathrm{e}} (\textbf{x} ) - \mathcal{T}^{ \mathrm{i}} (\textbf{x} )=0, &\textbf{x}\in \partial \mathscr{D}, \\ \kappa _{\mathrm{e}}\partial _{\textbf{{n}}}\mathcal{T}^{\mathrm{e}} (\textbf{x} )-\kappa _{\mathrm{i}}\partial _{\textbf{{n}}} \mathcal{T}^{\mathrm{i}} (\textbf{x} )=0, &\textbf{x} \in \partial \mathscr{D}, \\ \kappa _{\mathrm{e}}\partial _{\textbf{n}}\mathcal{T} (\textbf{x} )=0, & \textbf{x}\in \varGamma _{\mathrm{sides}}, \\ \kappa _{\mathrm{e}}\partial _{\textbf{n}}\mathcal{T} (\textbf{x} )+a \mathcal{T} (\textbf{x} )=0, & \textbf{x}\in \varGamma _{\mathrm{back}}, \\ \kappa _{\mathrm{e}}\partial _{\textbf{n}}\mathcal{T} (\textbf{x} )+a \mathcal{T} (\textbf{x} )=\alpha \mathcal{Q}_{\textbf{s}} (\textbf{x} ), & \textbf{x}\in \varGamma _{\mathrm{front}}, \end{cases} $$

(2)

where no dependence in time is present.

The direct problem, which models the generation of a thermogram is then the following: knowing the dimensions of the plate $\mathscr{R}$, the damaged region $\mathscr{D}$, the involved physical parameters ($\kappa _{e}$, $\kappa _{i}$, $\rho _{e}$, $\rho _{i}$, $c_{e}$, $c_{i}$, ε, α, σ, h), the surrounding air temperature $T_{\mathrm{air}}$, and the location s, frequency ω and intensity $\mathcal{I}$ of the lamp, find the complex thermal amplitude of the temperature $\mathcal{T}_{\textbf{s},\boldsymbol{\omega}} ^{\mathrm{front}}(\textbf{x})$ for $\textbf{x}\in \varGamma _{ \mathrm{front}}$ which is the solution $\mathcal{T}$ of problem (2) on the boundary $\varGamma _{\mathrm{front}}$. It is important to remark here that infrared cameras rather than the temperature (or the complex thermal amplitude of the temperature) they measure thermal radiation, and from it, they calculate the temperature distribution. In this paper, we assume that the thermal radiation obtained by the camera has already been processed and therefore the thermogram contains directly the values of the amplitude of the temperature oscillations $\mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{ \mathrm{front}}$. The associated inverse problem is then the following: knowing the thermogram $\mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{ \mathrm{front}}$, the dimensions of the plate $\mathscr{R}$, the involved physical parameters ($\kappa _{e}$, $\kappa _{i}$, $\rho _{e}$, $\rho _{i}$, $c_{e}$, $c_{i}$, ε, α, σ, h), the surrounding air temperature $T_{\mathrm{air}}$, and the location s, frequency ω and intensity $\mathcal{I}$ of the lamp, find the damaged region $\mathscr{D}$ such that the solution $\mathcal{T}$ of problem (2) for that region $\mathscr{D}$ agrees on $\varGamma _{\mathrm{front}}$ with the thermogram $\mathcal{T} _{\textbf{s},\boldsymbol{\omega}}^{\mathrm{front}}$, i.e., such that

$$ \mathcal{T}(\textbf{x})=\mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{ \mathrm{front}}( \textbf{x}),\quad \textbf{x}\in \varGamma _{ \mathrm{front}}. $$

(3)

3 Inversion method

This section is devoted to the solution of the inverse thermal problem by means of a topological derivative-based algorithm. The starting point is to rewrite our inverse problem in a less demanding way as follows. For a given temperature distribution $\mathcal{T}_{\textbf{s},\boldsymbol{\omega}} ^{\mathrm{front}}$ (that is, a thermogram), instead of imposing the identity (3) to find the damaged region $\mathscr{D}$, we consider a weaker formulation, where the aim is to find the set of defects $\mathscr{D}$ that minimize the quadratic error cost function

$$ J_{\textbf{s},\omega }(\mathscr{R}\setminus \overline{\mathscr{D}})= \frac{1}{2} \int _{\varGamma _{\mathrm{front}}} \bigl\vert { \mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{\mathscr{D}}} {( \textbf{x})}-\mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{\mathrm{front}}{( \textbf{x})} \bigr\vert ^{2}\,{d\gamma _{\textbf{x}}}, $$

(4)

where ${\mathcal{T}_{\textbf{s},\boldsymbol{\omega}}^{\mathscr{D}}}$ is the solution of (2) when the defects are occupying the region $\mathscr{D}$, and $d\gamma _{\textbf{x}}$ stands for surface integration. Notice that in this formulation the defects $\mathscr{D}$ are the variable of the cost function, while the time-harmonic heat conduction problem plays the role of a constraint.

The reconstruction algorithm to be applied consists now of computing the so-called topological derivative or topological gradient of the cost function $J_{\textbf{s},\omega }$, which is the scalar function $D_{\textbf{s},\omega }$ defined at each point $\textbf{x} \in \mathscr{R}$ by the asymptotic expansion [47, 59]

$$ J_{\textbf{s},\omega } \bigl(\mathscr{R}\setminus \overline{B_{\epsilon }(\textbf{x})} \bigr)= J_{\textbf{s},\omega } (\mathscr{R} )+ f(\epsilon ) D_{ \textbf{s},\omega }(\textbf{x}) + o \bigl(f(\epsilon ) \bigr), \quad \mbox{as }\epsilon \to 0, $$

(5)

where $B_{\epsilon }(\textbf{x})$ is the ball of radius ϵ and center x, and $f(\epsilon )$ is a positive increasing function selected to guarantee that the expansion (5) holds. In our case, adapting the results in [6, 24] we find that it can be taken as the volume of the ball, namely $f(\epsilon )=4\pi \epsilon ^{3}/3$.

In view of the asymptotic expansion (5), the topological derivative is a measure of the sensitivity of the cost function to considering an infinitesimal ball-like defect at each point x of the plate. Consequently, we can interpret it as a damage indicator function since the cost function is expected to decrease at the points where this derivative attains large negative values. Then, our reconstruction strategy just consists of identifying the points where the topological derivative attains pronounced negative values as belonging to a defect. That is, we approximate the true set of defects $\mathscr{D}$ by the set of points

$$ \mathscr{D}_{\mathrm{approx}}= \Bigl\{ \textbf{x}\in \mathscr{R}; D_{\textbf{s},\omega }(\textbf{x})< \lambda \min_{\textbf{y}\in \mathscr{R}}D_{\textbf{s},\omega }( \textbf{y}) \Bigr\} , $$

(6)

where $0<\lambda <1$ is a parameter that can be calibrated.

It is clear now that the asymptotic expansion (5) motivates our reconstruction strategy. However, from the numerical point of view it does not give a practical way to compute the topological derivative. Using the strategy proposed in [46] one can obtain closed-form formulae for the topological derivative of cost functionals of the form (4). In particular, this was used by one of us to obtain in [6] a closed-form formula for a time-harmonic problem in two dimensions analogous to the one in the present paper, and in [37, 38] for analogous time-harmonic electromagnetic problems in three dimensions. Performing suitable adjustments in the derivations in [6, 37, 38] to deal with the boundary conditions in the direct problem (2) it can be proven the following result.

Theorem

For the cost functiondefined in (4), the topological derivative at each pointis given by the formula

(7)

whereis the solution to the direct problem corresponding to a healthy plate:

(8)

andis the solution to a related adjoint problem

(9)

The previous theorem provides a very simple formula to implement the reconstruction method, we just need to solve problems (8) and (9), being (8) the first one to be solved since the source term at the last equation in the adjoint problem depends on the solution to (8), and compute then the gradients of such solutions to implement the formula (7). Then, the approximate domain is obtained easily by identifying the points where the topological derivative attains the largest negative values, as explained in formula (6). Notice also that both problems (8) and (9) take place in the healthy plate $\mathscr{R}$ and consequently, no a priori information about the number, size or shape of the searched defects is needed. The information about the defects (the measured thermogram $\mathcal{T}^{ \mathrm{front}}_{\textbf{s},\boldsymbol{\omega}}$) is incorporated in the topological derivative by means of the right hand side in the last equation in the adjoint problem (9).

In practice, and as we will illustrate in the forthcoming section devoted to the numerical experiments, processing just one thermogram corresponding to one lamp position and one frequency could be inadequate to properly identify the damaged region. To overcome this problem, one can either consider different lamp locations and process all the associated thermograms, or consider different thermograms obtained at different frequencies, or combine both strategies. Then, when having a set of thermograms $\mathcal{T}_{\textbf{s}_{i},\omega _{j}}^{ \mathrm{front}}$ corresponding to locating the lamp at $N_{ \mathrm{lamps}}$ different points $s_{i}$, $i=1,\dots , N_{ \mathrm{lamps}}$ ($N_{\mathrm{lamps}}=1$ in case just one position is considered), and use $N_{\mathrm{freq}}$ different oscillation frequencies $\omega _{j}$, for $j=1,\dots , N_{ \mathrm{freq}}$ ($N_{\mathrm{freq}}=1$ if just one frequency is used), one can consider a new cost functional, corresponding to a linear combination of each individual cost functional (4), namely,

$$ J(\mathscr{R}\setminus \mathscr{D})=\sum _{i=1}^{N_{\mathrm{lamps}}} \sum_{j=1}^{N_{\mathrm{freq}}} p_{ij} J_{\textbf{s}_{i},\omega _{j}}( \mathscr{R}\setminus \mathscr{D}), $$

(10)

where $p_{ij}>0$ are weights to be selected. Then, by linearity (see the definition (5)), we find that for each $\textbf{x}\in \mathscr{R}$ the topological derivative of the functional (10) is

$$ D(\textbf{x})=\sum_{i=1}^{N_{\mathrm{lamps}}} \sum_{j=1}^{N_{ \mathrm{freq}}} p_{ij} D_{\textbf{s}_{i},\omega _{j}}(\textbf{x}), \quad \textbf{x}\in \mathscr{R}. $$

(11)

Then, the true domain can be approximated by the set

$$ \mathscr{D}_{\mathrm{approx}}= \Bigl\{ \textbf{x}\in \mathscr{R}; D(\textbf{x})< \lambda \min_{\textbf{y}\in \mathscr{R}}D(\textbf{y}) \Bigr\} , $$

(12)

where, as in (6), the parameter $0<\lambda <1$ has to be selected.

The simplest choice for the weights in (10) is $p_{ij}=1$ for all i, j. However, by doing this the contribution of some lamps and/or frequencies could be disregarded while the contribution of some other ones could be too stressed (recall that we are interested in the points where the topological derivative attains the largest negative values). To somehow make uniformly important all the contributions, we propose, as already done in [19, 38, 50] for some related multifrequency data, to use the weights

$$ p_{ij}= \Bigl\vert \min_{\textbf{y}\in \mathscr{R}} D_{\textbf{s}_{i}, \omega _{j}}(\textbf{y}) \Bigr\vert ^{-1} $$

(13)

for which it follows that the contribution of each term in (11) satisfies

$$\min_{\textbf{y}\in \mathscr{R}} \bigl(p_{ij} D_{\textbf{s}_{i}, \omega _{j}}( \textbf{y}) \bigr)=-1. $$

Selecting the weights $p_{ij}$ as in (13) for the definition (12) of the approximate domains $\mathscr{D}_{\mathrm{approx}}$ provides good approximations of the true defects, as illustrated in the next section. Nevertheless, some other choices of $p_{ij}$ are possible and could improve the results. For related selections of weights, either a priori or a posteriori (as done in this paper) can be found in [23, 29, 31, 51, 56].

4 Numerical results

In this section we present some numerical results to illustrate the performance of the method.

Both the computation of the topological derivative, i.e., the solution of the state (8) and adjoint (9) problems; as well as the generation of synthetic thermograms via solving the forward problem (2) was done by means of the finite element solver Freefem++ which allows for variational formulation of generic problems and solving them with different types of elements [26]. In our case P2 elements (piecewise polynomials of degree two) defined on a tetrahedral mesh were used.

For all our experiments we consider a plate of dimensions $\ell _{x}=0.01 \mbox{ m}$, $\ell _{y}=0.50\mbox{ m}$ and $\ell _{z}=1.00\mbox{ m}$ containing three different defects:

• A box-shaped defect of height $\frac{6\ell _{x}}{4}=15\times 10^{-3}\mbox{ m}$, width $\frac{3\ell _{x}}{4}=7.5\times 10^{-3}\mbox{ m}$ and thickness $\frac{\ell _{x}}{2}=5\times 10^{-3}\mbox{ m}$ with centroid at $(-\frac{\ell _{x}}{6},- \frac{\ell _{y}}{4},\frac{\ell _{z}}{4})= (-1.67\times 10^{-3},0.125,0.25 ) \mbox{ m}$.

• A spherical defect of radius $\frac{\ell _{x}}{3}=3.33\times 10^{-3}\mbox{ m}$ and center at $(0,\frac{\ell _{y}}{4},0 ) = (0,0.125,0 )\mbox{ m}$.

• A groove-shaped defect of height $\frac{30\ell _{x}}{4}=75\times 10^{-3} \mbox{ m}$, width $\frac{\ell _{x}}{4}=2.5\times 10^{-3}\mbox{ m}$ and thickness $\frac{\ell _{x}}{2}=5\times 10^{-3}\mbox{ m}$, with centroid at $(-\frac{\ell _{x}}{6},- \frac{\ell _{y}}{4},-\frac{\ell _{z}}{4})= (-1.67\times 10^{-3},0.125,-0.25 ) \mbox{ m}$.

The faulty plate is represented in Fig. 3.

Figure 3

Faulty plate. Plate and three zoomed regions that contain the defects

The faulty plate was meshed in gmsh [22] using tetrahedrons obtaining 70,287 nodes and 405,089 elements. For solving the state and adjoint problems, the associated healthy plate (i.e, the plate without defects) was meshed obtaining 85,261 nodes and 438,944 elements. As can be seen in Fig. 4 these meshes not only differ in the number of nodes and elements, but also in the structure: the mesh for the faulty plate concentrates more elements near the defects, while in the mesh for the healthy plate the size of all the elements is uniform. Furthermore, this second mesh is finer in order to be able to capture objects of small size.

Figure 4

Meshes. Left: mesh corresponding to the faulty plate. Right: mesh corresponding to the healthy plate. Both meshes have approximately the same number of nodes, however for the faulty plate, these are concentrated near the three defects

The exterior thermal coefficients in the model were chosen as $\kappa _{\mathrm{e}} = 200 \mbox{ W/(m }\cdot \mbox{ K)}$, $\rho _{ \mathrm{e}} = 2700 \mbox{ kg/m}^{3}$ and $c_{\mathrm{e}} = 900 \mbox{ J/(kg }\cdot \mbox{ K)}$ which correspond to an aluminium alloy, while the defects are supposed to be air cavities with associated interior thermal coefficients $\kappa _{\mathrm{i}} = 0.025 \mbox{ W/(m }\cdot \mbox{ K)}$, $\rho _{\mathrm{i}} = 1 \mbox{ kg/m}^{3}$ and $c_{\mathrm{i}} = 1000 \mbox{ J/(kg }\cdot \mbox{ K)}$.

For the boundary conditions, the air is assumed to be at room temperature $T_{\mathrm{air}} = 290 \mbox{ K}$ and $h = 15 \mbox{ W/(m} ^{2}\cdot \mbox{ K)}$ which accounts for the possible natural convection. The hemispherical emissivity and absorptance of the aluminium surface at the expected temperatures are $\varepsilon =0.08$ and $\alpha =0.4$ respectively. Finally $\mathcal{I}=6000\mbox{ W}$, that is, the intensity of the lamp oscillates with amplitude equal to 6000 watts and zero phase.

With this configuration and the lamp situated 0.15 m apart from the surface of the plate, we get maximum temperature variations of $|\mathcal{T}|\approx 0.5\mbox{ K}$ for the lower frequencies (0.8 Hz) to $|\mathcal{T}|\approx 0.25\mbox{ K}$ for the higher ones (2 Hz).

Inverse crimes may occur when the same model is used to generate and process synthetic data [14, 45]. We avoid them both by solving the forward and inverse problems in different meshes (i.e., by considering different meshes to generate the synthetic data in the faulty plate, and to compute the topological derivative which requires solving two problems in the healthy plate), as well as by adding some random noise to simulate measurement errors of the form:

$$ N (\textbf{x} )=\eta \max \bigl( \vert \mathcal{T} \vert \bigr) \bigl(U_{ (0,1 )} (\textbf{x} )-0.5+ \bigl(U _{ (0,1 )} ( \textbf{x} )-0.5 \bigr)i \bigr), \quad \text{for each node }\textbf{x}, $$

(14)

where $U_{ (0,1 )}$ stands for the uniform distribution between 0 and 1 and the dependency on x means that for each node of the mesh two realizations of the distribution are made, one for the real part and one for the imaginary one. The uniform distribution was generated using the randreal1 function from FreeFem++ [26]. That way, if $\mathcal{T}$ is the simulated temperature obtained by solving (2) for the true defects $\mathscr{D}$, the corresponding noisy thermogram will be

$$\mathcal{T}^{\mathrm{front}} (\textbf{x} ) = \mathcal{T} (\textbf{x} ) + N ( \textbf{x} ),\quad \textbf{x}\in \varGamma _{\mathrm{front}}. $$

The constant η in (14) allows us to control the level of measurement noise, namely, the maximum magnitude of the difference $|\mathcal{T}^{\mathrm{front}}- \mathcal{T}|$. For example, setting $\eta =0.1$ means having a level of noise of 0.05 degrees for the lowest frequencies and 0.025 for the higher ones, which agrees with the precision of a high-end infrared camera (0.02 degrees [33]). If we compare the ratio noise to information, namely, the quotient between the noise we are adding and the useful information we have (where the useful information is defined as the difference between the thermogram of a faulty plate and the one we would get on a healthy one), we find that for $\eta =0.1$ the noise is comparable to the 70% of the information for the higher frequencies and up to the 90% for the lower ones. See Table 1 for a more in depth explanation.

Table 1 Comparative between the amount of information and added noise. For six linearly-spaced frequencies between 0.8 and 2 Hz we present in the first row the maximum amplitude of a thermogram, on the second row the maximum amplitude of the difference between the thermogram of a faulty plate and the one corresponding to a healthy one, that is the amount of useful information, and on the last row the quotient between the added noise and that amount of information. For the considered frequency range, this ratio decreases as we increase the frequency

In Figs. 5 and 6 we show the thermogram obtained when a lamp is situated at $\textbf{s}= (-0.15,0,0 ) \mbox{ m}$ and its intensity oscillates at frequencies of 0.8 and 2 Hz respectively. We can see that with that amount of noise is virtually impossible to use the raw thermograms as an indicator of the position, size or shape of the defects.

Figure 5

Thermograms with $\omega =0.8\mbox{ Hz}$. Left: amplitude of the temperature $|\mathcal{T}^{\mathrm{front}}|$. Right: difference between thermograms corresponding to a healthy plate and a faulty one. At low frequencies defect detection is impossible at plain sight

Figure 6

Thermograms with $\omega =2\mbox{ Hz}$. Left: amplitude of the temperature $|\mathcal{T}^{\mathrm{front}}|$. Right: difference between thermograms corresponding to a healthy plate and a faulty one. At high frequencies the box defect, which is the shallower one, is barely visible

The values of the topological derivative obtained when processing the thermogram for $\omega =0.8\mbox{ Hz}$ (shown in the left plot in Fig. 5) on the planes $x=-\ell _{x}/2$ (namely, on $\varGamma _{\mathrm{front}}$), $x=-\ell _{x}/4$, $x=0$, and $x=\ell _{x}/4$ are shown in Fig. 7. We observe that the largest negative values (blue colors) are not concentrated inside any of the three defects. Indeed, we only find that it presents both local maxima (red colors) and minima (blue colors) in the points closer to the lamp which hide the actual defects. Notice that since our thermogram is highly polluted by noise, the topological derivative on $\varGamma _{\mathrm{front}}$ is non-smooth, having sharp changes from one node to its neighbours. However, the diffusivity of the heat equation makes the topological derivative to be much smoother as we move away from $\varGamma _{\mathrm{front}}$ to $\varGamma _{\mathrm{back}}$. We also observe that the topological derivative is qualitatively rather similar at all the planes, attaining the most pronounced negative values on $\varGamma _{\mathrm{front}}$. For this reason and for the sake of brevity, in the forthcoming figures we will only represent the values of the topological derivative on $\varGamma _{\mathrm{front}}$, unless otherwise stated.

Figure 7

Values of the topological derivative on several planes of the form $x=x_{0}$. First: $x=-\ell _{x}/2$. Second: $x=-\ell _{x}/4$. Third: $x=0$. Forth: $x=\ell _{x}/4$

When processing the thermogram for $\omega =2\mbox{ Hz}$ shown in the left plot in Fig. 6, we find completely analogous results, as can be seen in the right plot in Fig. 8 where the topological derivative on $\varGamma _{\mathrm{front}}$ is shown. Its counterpart for $\omega =0.8\mbox{ Hz}$ is also represented in the left plot for ease of comparison.

Figure 8

Topological derivatives on $\varGamma _{\mathrm{front}}$ for a single experiment. Left: topological derivative corresponding to $\omega =0.8\mbox{ Hz}$. Right: topological derivative corresponding to $\omega =2\mbox{ Hz}$. Any possible detection of the defects is due to cognitive bias

If we use formula (11) to combine both frequencies we get the topological derivative shown in Fig. 9 suggesting that combining frequencies may be a good way to minimize the effect of the noisy thermograms. On the zoom we can see that the topological derivative does have a local minima where the defect is, however this is not visible compared with the spurious ones.

Figure 9

Topological derivative on $\varGamma _{\mathrm{front}}$ for two combined frequencies. Left: topological derivative combining single-frequency data at 0.8 and 2 Hz. Right: zoom of the previous image around the box-shaped defect with the colormap rescaled to the visible range. The outline of the defect is also shown in black on the image on the right

This is confirmed by combining 10 linearly spaced frequencies between 0.8 and 2 Hz as shown in Fig. 10. In this figure, one can clearly see traces of the three defects. Furthermore, the global minimum of the topological derivative now corresponds to a true defect, the spherical one.

Figure 10

Topological derivative on $\varGamma _{\mathrm{front}}$ for ten combined frequencies. Left: topological derivative combining ten linearly-spaced frequencies between 0.8 and 2 Hz. Right: zoom of the previous image around the groove-shaped defect with the colormap rescaled to the visible range. The outline of the defect is also shown in black on the image on the right

The main drawback of this method is that, being performed all the experiments with a lamp located in the same position, some of the defects may be situated in a more favourable region for detection than others. This is suggested by the fact that in the last example the global minimum points out the region where the spherical defect is located although this defect is the smallest and deepest one. To corroborate this, in Fig. 11 we show the results obtained for two single-frequency experiments at $\omega =2\mbox{ Hz}$ with lamps situated a two different locations $\textbf{s}= (-0.15,0,0.25 )\mbox{ m}$ and $\textbf{s}= (-0.15,0,-0.25 )\mbox{ m}$. We can see that by placing the lamp near any of the defects we highlight that defect at the expense of loosing most of the information about the other ones. This is probably caused because the gradients of the thermogram are attained near the lamp and the farther we get from it the flatter the temperature distribution is, meaning that the effect of having a defect is smaller.

Figure 11

Topological derivative on $\varGamma _{\mathrm{front}}$ for two different lamp locations. Left: lamp near the box-shaped defect. Right: lamp near the groove-shaped defect. As one can see, reconstructions highly depend on the position of the lamp

To overcome this problem, different locations of the lamp may be combined as a way to diminish the influence of the location of the defects in its detection. In Fig. 12 on the left we show the results for the combination of three different locations for the lamp at a frequency of 2 Hz. The lamps are linearly spaced over the z axis, that is, they are situated at $\textbf{s} _{1}= (-0.15,0,-0.25 )\mbox{ m}$, $\textbf{s}_{2}= (-0.15,0,0 ) \mbox{ m}$ and $\textbf{s}_{3}= (-0.15,0,0.25 )\mbox{ m}$.

Figure 12

Topological derivative on $\varGamma _{\mathrm{front}}$ combining three different lamp locations. Left: lamps linearly spaced along the z axis at a frequency of 2 Hz. Right: counterpart at 0.2 Hz. Combining different lamp locations is a great way to improve the reconstruction of several defects at the same time, but we still have to choose the adequate frequency

We should not overestimate this result as we should recall that the frequency of 2 Hz is specially well suited for this experiment. If we had chosen $\omega =0.2\mbox{ Hz}$, the result would have been much worse as shown in Fig. 12 on the right. Results for the suitable frequency were shown first to better illustrate the difference between different lamp locations. The suitability of a particular frequency depends on the size of the defect, which is, in principle, an unknown. This motivates the combination of several frequencies.

By combining thermograms coming from different frequencies and locations of the lamp, the corresponding topological derivative provides with great precision the shape, number and position of the three defects without having to search for the most adequate lamp positions and/or frequencies. In Fig. 13 we show the topological derivative corresponding to six linearly spaced frequencies between 0.2 and 2 Hz and situated on a linearly spaced grid on the xz-plane. We can see that combining multiple frequencies helps to smooth the background and to be able to recover the defects without a priori information about their sizes, whereas combining different lamp locations makes all the defects to stand out without a priori information about their location. We could see that with a 2-by-2 grid of lamps, as shown on the left, the grid is not fine enough to capture all three defects because the spherical one happens to be worse situated than the rest. This is solved by using a finer grid, i.e., by increasing the number of lamp positions.

Figure 13

Topological derivative on $\varGamma _{\mathrm{front}}$ combining several frequencies and lamp locations. Left: topological derivative corresponding to six linearly spaced frequencies between 0.2 and 2 Hz and a 2-by-2 grid of lamps. Right: counterpart when a 3-by-3 grid of lamps is used. The appearance of spurious minima is highly reduced in both cases due to the multifrequency character, and increasing the number of lamps highly improves the reconstructions

One may argue that there is no clear distinction between the spurious minima and the minima corresponding to true defects, however this is just a consequence of showing the value of the topological derivatives as colormaps. If we instead plot the values of the topological derivative on the surface $\varGamma _{\mathrm{front}}$ as the three-dimensional surface $x=D (-\frac{\ell _{x}}{2},y,z )$ we can clearly see that the true defects correspond to local minima which are much more pronounced than its surroundings. This is shown in Fig. 14 where there is no doubt that the first one (which corresponds to Fig. 13-left) is detecting two defects whereas the second one (corresponding to Fig. 13-right) is detecting the presence of the three true defects.

Figure 14

Topological derivatives shown in Figure 13 plotted as surfaces. Left: topological derivative corresponding to six linearly spaced frequencies between 0.2 and 2 Hz and a 2-by-2 grid of lamps. Right: counterpart when a 3-by-3 grid of lamps is used. True defects appear as coherent mountain-like structures opposite to the random peaks of the spurious minima

Due to the short-range character of the heat transfer equation, the topological derivative appears to always have its local minima on the inspecting boundary ($\varGamma _{\mathrm{front}}$), as can be observed in Fig. 15, where we show the values of the topological derivative at different planes of the form $x=x_{0}$ for our last example (corresponding to 9 lamp positions and 6 frequencies). Therefore, in view of this, we cannot aim to recover the correct depth of the defects by using formula (12). Nevertheless, reconstructions will be accurate in the y and z directions and rather insensitive to the selected value of λ in (12). To visualize this, we consider first the value $\lambda =0.4$ and find the reconstructed defect $\mathscr{D}_{ \mathrm{approx}}$ by using formula (12). Figures 16, 17, and 18 show the reconstructions for the box-shaped, groove-shaped and the spherical defects, respectively. The three figures share a common scale and colormap to highlight that the global minimum is attained at the box defect and that the least pronounced minimum corresponds to the sphere. On each figure, we show the values of the topological derivative along three perpendicular planes that intersect each defect ($x=0$, $y=-0.125$, and $z=0.25$ for the box, $z=-0.25$ for the groove and $z=0$ for the sphere), as well as one view in perspective of the true defect in black, and the reconstructed domain $\mathscr{D}_{\mathrm{approx}}$, in green. On the planes the outline of the sections of the true (black line) and the reconstructed defects (magenta line) are superimposed. We can see that in the three cases the reconstructions begin in the surface $\varGamma _{\mathrm{front}}$ of the plate as expected. Indeed, for the groove-shaped and the spherical defects, reconstructions do not cut the $x=0$ plane (we do not see any magenta line on the corresponding images).

Figure 15

Values of the topological derivative on several planes of the form $x=x_{0}$. First: $x=-\ell _{x}/2$. Second: $x=-\ell _{x}/4$. Third: $x=0$. Forth: $x=\ell _{x}/4$

Figure 16

Reconstruction of the box-shaped defect. Upper left: Topological derivative at $x=0$. Upper right: Topological derivative at $y=-0.125$. Lower left: Topological derivative at $z=0.25$. Lower right: Perspective showing the reconstruction for $\lambda =0.4$ and the original defect

Figure 17

Reconstruction of the groove-shaped defect. Upper left: Topological derivative at $x=0$. Upper right: Topological derivative at $y=-0.125$. Lower left: Topological derivative at $z=-0.25$. Lower right: Perspective showing the reconstruction for $\lambda =0.4$ and the original defect

Figure 18

Reconstruction of the spherical defect. Upper left: Topological derivative at $x=0$. Upper right: Topological derivative at $y=-0.125$. Lower left: Topological derivative at $z=0$. Lower right: Perspective showing the reconstruction for $\lambda =0.4$ and the original defect

However, in terms of yz location, shape and size are very accurately found, even in spite that we are using thermograms with a level of noise going from the 70% to 90% of the useful information. Furthermore, while the three defects are well recovered, spurious defects do not appear, that is, $\mathscr{D}_{\mathrm{approx}}$ only consists in the three green defects that are shown in Figs. 16, 17, and 18.

By selecting the constant $\lambda =0.3$ in formula (12), we obtain a rather similar reconstruction, as can be seen in Fig. 19. This means that the method is not very sensitive to the calibration of λ (although undudobtely, a drastical increase/decrease of it would yield to reconstructions understimating/overstimating sizes, and loosing/including some small/spurious defects).

Figure 19

Reconstruction of the three defects for $\lambda =0.3$ Left: perspective of the box-shaped defect. Center: perspective of the groove-shaped defect. Right: perspective of the spherical defect. The reconstructed defects do not depend highly on the choice of the value for λ

5 Conclusions and future work

We have shown that time-harmonic infrared thermography is a promising method for contact-less defect detection when thermograms are suitably post-processed. Our proposed method based on the computation of the topological derivative enables to easily combine thermograms obtained at different frequencies and locations, what we have shown to be crucial if we do not have a priori information on the possible size or location of the defects.

As reported by many authors for different non-destructive testing settings (see for instance [24, 37, 41, 42]), one of the main advantages of using topological derivative based methods instead of using the raw measured data is that topological derivatives can give an estimation of the three dimensional position and shape of the defects, whereas the raw data can only aspire to show the projection of the defects on the inspected two dimensional surface. In our time-harmonic heating setting, however, the topological derivative seems to give barely any information about this third dimension perpendicular to the inspected surface. This does not render the topological derivative useless as we have shown that it is an excellent post-processing tool, highly improving the detection of the defects on the inspected surface, providing their number, shape and localization without any need of iterative processes. Furthermore we have tested the topological derivative against highly polluted thermograms and it was able to perform correctly.

From a practical point of view, when inspecting a plate-like structure, the topological derivative could be used to pinpoint the yz-location of possible defects and then decide to either repair that area or use at that zone another inspection method as the ultrasonic test, greatly reducing the inspection time, as this latter method would only need to analyse a discrete set of locations (namely, a very small area in comparison with the whole excitation side $\varGamma _{\mathrm{front}}$).

From an academic point of view the topological derivative could be used as an initial guess generator for iterative methods like shape-derivative methods [13, 59], which need to know the exact number of defects, or with level-set methods [39, 49, 58] which; despite not needing a priori information on the number of defects, can converge to a solution much faster if the initial guess is closer to the solution.

Another promising option is computing the second order topological derivative [17, 48] in zones where the topological gradient has marked the presence of defects. The second order topological derivative can capture more information about the effect of placing defects of finite size, thus it may be able to discern the utility of placing defects on the surface or in the interior of the plate.

On the other hand, despite the topological derivative having always its minima on the surface, the region around the minima seems to behave differently for defects situated at different depths. This may be difficult to describe in a systematic way, but maybe some machine learning algorithm could be trained to learn to extract the depth of the objects by post-processing the topological derivatives.

Finally, in view of published work comparing time-harmonic and more general time-dependent excitations for a thermal problem in an unbounded two-dimensional media [7], we believe that by processing time-dependent thermograms via topological derivative methods it could be possible to find the correct depth of the defects. The extension of the method to the general time-dependent case is not immediate and we plan to do it in the near future.