Cite this asTapia AO, Tsonchev RI, D´ıaz Viera MA, Ortiz MH (2019) Modeling of active thermography through uncertainty quantification of parameters of the heat transfer equation. Ann Math Phys 2(1): 051-057. DOI: 10.17352/amp.000009
Active thermography is an experimental technique used to analyze samples of materials or entire structures without destroying them, by means of a heat source, such as a laser beam of a given power. It is posed that such experimental procedure can be modeled mathematically through the complete equation of heat transfer. The uncertainty on the assumption of the value of the parameter emissivity of this equation is to be analyzed calculating the error between concrete experimental data and simulations where such parameter has been taken from the uniform distribution. To the extent of this research, specifically for active thermography, no previous attempt has been made for using at the same time the complete equation of heat transfer (without simplifications or linearizations), with the usage of uncertainty quantification for the specific experimental results to which the mathematical theory was applied.
Active thermography uses a heat source in order to induce absorption of such heat into a sample of a possibly unknown composition, which then will dissipate through several heat transfer mechanisms [1-4]. From the parameters of those heat mechanisms, it is possible to infer the heat properties of the sample. The need of a mathematical model, particularly one that is suitable for an inverse problem, is the central matter of this project. The inverse problem is that given experimental results, determine the possible value of a coefficient of the differential equation that models the phenomenon.
In active thermography, an external energy source is required for artificial thermal excitation. Given the possibility of controlling the intensity of the external energy source, the artificial thermal excitation can reach deeper atoms in the object, and therefore information can be obtained from more internal layers . Thus, it is used thermal data from a heated sample with a laser beam, on its original state, that is, without reducing it to dust or any other destructive process. The set of performed observations (measurements), may be thus modeled through a heath-transfer, partial differential equation (transport equation). The spatial thermal parameters of a sample may differ when this has a concrete shape and dimensions, as opposite of measuring fine grinding of the material this sample is made of. Also, there are instances where the sample cannot or must not be destroyed. This is the case of thermography for cancer research, aiming at substituting the painful mammograms with much less invading techniques .
Consider that a thin slab of a solid sample of thickness th is heated with a light beam (laser) that is uniformly focused onto one or more of its surfaces (boundaries) . On the opposite side, its temperature can be monitored as a function of time, for example with a thermocouple. The variation with time, t, of the heat generated in the sample, Q, due to the absorption of light of incident power Φ, is given by
where q is a term specified by the sum of the power losses by radiation, convection and conduction. If we want to calculate the rise of temperature, θ, of the surfaces that are not in contact with an incident laser, we must express the heat term of Eq.1 as a function of that increase , which can be expressed as
This is a linear ordinary differential equation with constant coefficients, where h is the total heat transfer coefficient (a linearized sum coming from convection, conduction and radiation), ρm is the material density, Cp is the specific heat capacity, and V is the volume of the sample; the particular units used in this work will be mentioned later on. The solution of Eq.2 is .
where A is the total area of the sample and
is the relaxation (transient) time of the system. Eq.3 has been used to calculate the specific heat capacity ·Cp for thin slabs of known thickness , because it is obtained from τ by least squares fit of experimental curves of θ versus t Eq.3, and normally semi-log plots are used to avoid uncertainties due to deviations from the theoretical model .
Successful as it is the approximation mentioned in  (and other similar to this, like in [7,8], the usage of logarithms in one coordinate implies that some information is lost, as irregularities are flatten down. Other works (like in (Benzerrouk)), describe, but not use the full heat transfer equation. Still other applications for active thermography use the full heat transfer equation, but no uncertainty quantification is presented ( like in (Cannas, et al.,) . It could be therefore useful to use a non-linear model, using a numerical approach, for which one continues to develop the mathematical model. Therefore, it is the purpose of this work to develop a full non-linear heat transfer model, and then use it to try and infer the parameters of a sample of known composition, through uncertainty quantification, thus opening the possibility for using such model for inference of temperature of a sample of unknown composition. This project attempts to model active thermography mathematically, and then measure the uncertainty of one of its parameters, namely the emissivity, by dint of comparison with some experimental data, from specific, chosen experiments. The novelty of this work comes, therefore, from the implementation of the full, non-linear heat transfer equation, coupled with using Uncertainty Quantification for the particular type of experimental data that was available for this work.
It is proposed that with experimental data from active thermography, determine both experimentally and with theoretical modeling the thermal properties of different materials and the heat loss from the material under consideration.
Modeling, as is used by Dr. Mart´ın Alberto D´ıaz Viera [9,10], consists of four main steps:
1. Conceptual Modeling: The problem is put under the context of its physical phenomenology, delimiting those observables that one wish to quantify. In this case, it is modeled the heat transfer through conduction, convection and radiation.
2. Mathematical Modeling: Phases and components are systematically analyzed, together with the constitutive laws, in order to derive the balance equations and the initial and boundary conditions. In this case, a 3D heat transfer equation is obtained, where the Fourier law is part of the equation, and convection and radiation are treated as boundary conditions. Heat injection is treated as a puntual source. The emissivity was taken from the uniform distribution (see the Appendix, for some further details).
3. Numerical Modeling: Depending on the shape of the balance equations and the mathematical restrictions of the problem, it is chosen a convenient way of discretizing those equations (for example, with the finite element method). The numerical discretization depends also on the geometry of the domain of interest. In this case, Finite elements are used for discretizing the geometry, UMPFPACK is the linear solver, and Newton-Raphson is for discretizing time.
4. Computational Modeling: Once the numerical model has been obtained, a suitable computational platform is chosen, with a computational language and a software for writing the numerical model and be able to run it for computational simulations. It is the behaviour of this last model which finally gets validated with the experimental data, and where the level of uncertainty (discrepancy between observed and simulated results) can be evaluated to decide in which way and on which of the previously mentioned stages, modifications are in order, to approach as much as possible, to some predetermined level, the optimum of the parameters for the theoretical model. The implementation in this case was on COMSOL 3.5 a.
The derivation of the needed non-linear, partial differential equations includes using balance equations, heat conservations laws, and the Stefan-Boltzman Law and Fourier’s Law for radiation properties of heat .
The conceptual model consists of all the assumptions and hypothesis made:
• We are assuming isotropicity.
• Heat is provided uniformly onto one or more surfaces.
• The exchange of energy can represented by convection, conduction and radiation only.
• There are no chemical reactions, nor loss of matter, only loss of heat, and this last is lost to the air mostly.
• Gravitational potential is negligible.
• It is assumed uncertainty over the emissivity.
In order to make exposition swifter, the detailed derivation of the mathematical model has been developed in the appendix, and the reader is encouraged to study it. Here it is used the last result: Using Equation 30,
the relationship between heat and temperature of Eq.35, and Eq.31, and assuming no heat sources, it is obtained
or, to show explicitly the parabolic form of this partial differential equation,
A parabolic partial differential equations require at least one Dirichlet boundary condition, or Robin conditions (which is the case), given by the constitutive laws and the incident heating power source q0, which in general is
n · (λm∇T) = q0 + hconv · (Tenv − T(x,y,z,t)) + σεm(Tenv4− T4) (8)
For the boundaries with q0=0, the Robin boundary conditions become
n · (λm∇T) = hconv · (Tenv − T(x,y,z,t)) + σεm(Tenv4− T4) (9)
where the emissivity εm ∈ U[0.03,0.05], except for the upper part of the sample, which is painted in black and thus has emissivity εmb. q0 is represented by the laser beam through the extremes of the slab, as shown in Figure 1, namely Φ1 Φ2.
The initial condition, as mentioned, is T(0) = temperature of air.
T(x,y,z,t) Temperature [ K ]
Tenv Air temperature [ K ]
Ρm Density of material [kg/m3]
Φ Incident Power[W · m−2]
hconv λm σ εm Convective heat transfer coefficient [W/(m2·K)]
Thermal conductivity [W/(m·K)]
Stefan-Boltzmann constant [W·m−2· K−4]
Total emissivity of the material 
εmb Total emissivity of the material darked side 
Cp heat capacity of the material (at constant pressure) [J/(kg · K)]
So far the model, say, in Eq.7 assumes a problem which can be deal with directly. In this work, for validation purposes, such approach will be implemented initially. However, the final goal is to develop a methodology for analyzing the uncertainty in the parameters (Uncertainty Quantification -UQ- of the coefficients), of the Partial Differential Equation, namely for λm,hconv,εm, that is, when the sample’s composition is unknown, but experimental measurements have been obtained. In concrete, in this work it is assumed uncertainty over εm, for most of the surface of the sample.
The uncertainty of a measured parameter can be characterized with the dispersion of the values that could reasonably be attributed to the measurand error .
Some of the error components may be evaluated from the statistical distribution of the results of series of measurements and can be characterized by experimental standard deviations (Type A). The other components of the error, which also can be characterized by standard deviations, are evaluated from assumed probability distributions based on experience or other information (Type B) .
Thus Type A standard uncertainty is obtained from a probability density function derived from an observed (experimental) frequency distribution, while a Type B standard uncertainty is obtained from an assumed probability density function based on the degree of belief that an event will occur .
Because our mathematical model may be incomplete, all relevant quantities should be varied to the fullest practicable extent so that the evaluation of uncertainty can be based as much as possible on observed data (ISO & OIML, 1995).The measurement of Type A standard uncertainty is the second part of this project, whereas modeling Type B standard uncertainty is the third part of this project.
The numerical model consists of making the appropriate choice of the numerical methods in terms of precision and the efficiency for the solution of the mathematical model.
Although finite element methods (FEM) are usually substantially more difficult to program than Finite Differences (FD), this extra effort yields approximations that are of high-order accuracy even when a partial differential equation is solved in a general (nonrectangular) multidimensional region, and even when the solution varies more rapidly in certain portions of the region so that a uniform grid is not appropriate .
The FEM naturally incorporates a broader spatial extent, thus FEM can use a coarser mesh, compared with FD .
The basic mixed finite element (MFE) method , has shown to be more than good enough for solving a transport equation.
In this case, the resulting problem is a nonlinear partial differential equation with initial and boundary conditions, or only boundary conditions (stationary case). For the numerical solution the following methods are intended to be applied:
• For the time derivative, it can be used a second order backward finite differences discretization, resulting in a totally implicit scheme in time.
• For the rest of the differential operators, concerning the spatial derivatives, it can be applied a standard Galerkin finite element discretization, where Lagrange quadratic polynomials can be used as weighting and basis functions, which in this work imply a convergence of order two .
• A regular mesh of tetrahedral elements in 3D can be used.
• For the linearization of the nonlinear system of equations, an iterative Newton-Raphson method can be applied.
• For the solution of the resulting algebraic system of linear equations, it can be used a variant of the direct LU method for sparse, unsymmetrical matrices.
In view of the scale and resolution requirements for the transport model, it could be acceptable to perform the computational implementation making use of the standard finite element framework provided in COMSOL Multiphysics . In particular, using the PDE mode for time dependent analysis in the coefficient form. An alternative can be using the FEniCS implementation . For concreteness, it is exposed here how has been done the implementation in Comsol v 3.5 a.
First of all, the parameters of the simulation are defined as shown in Figure 2. Notice that at the left hand side are shown the parameters as written in Eq.7, and in the right hand side as is written in COMSOL multiphysics
T(x,y,z,t) = u
Tenv = T env
ρm = rhom
Φ = Phi
hconv = h conv
λm = lambda m
σ = sigma
εmi = Epsilonm
εmb = Epsilon-mb
Cp = Cp
Notice that i ∈ [2,11] is an index for the emissivities taken from the Uniform distribution.
The selected geometrical object represents an experimental metal sheet (Figures 3,4)
The definition of the coefficients of our equation are contained in the section scalar expressions, as can be seen in Figure 5.
Eq.7 is transformed into the comsol form, using the scalar expressions to define the coefficients, as can be shown in Figure 6.
The initial temperature of the material is considered the same as the surrouding air, that is T(0) = Tenv. The length of the simulation is for 5.3 seconds, which was seen enough to obtain a stationary solution.
There is only one type of conditions: Robin. However, in most boundaries T(t = 0) = Tenv and only in boundaries 2 and 5 T = Φ, for the duration of the numerical experiment, as can be seen in Figures 7,8.
The meshing is Physics-controlled (as per default in Comsol, Figure 9), with a so-called normal size in the elements (as per default in Comsol).
Solver Parameters: The chosen solver es UMFPACK.
Before actually making comparisons with experimental data, the computational model was tested for performance. The computation went flawlessly, and an illustration of the distributed heat loss is shown in Figure 10.
The model of Eq.7, with the boundary conditions shown in Eqs.8 and 9, with the initial condition T(t = 0) = Tenv = 300.15K, with hconv = 10[W/(m2 · K)], and the weak condition of injection heat Φ = 1000[W/m2] was set up, according to . The numerical solution (shown in Figure 11) has the expected radial behavior as described in .
The experimental setup as it is done in the laboratory of Prof. Tsonchev (Tsonchev, n.d.) is shown in its Comsol implementation as in Figure 12.
After making numerical simulations with 10 different values of emissivity, every time a graphic was made. Such graphic compares the experimental values with the numerical simulation. Wherever the experimental values were close to those of the simulation, the distance between them was measured, and it was found that the overall error was ≈ 7%. The standard deviation of errors was on the order of ≈ ±10−5% between errors (see the annexed pages for details of the numerical comparisons between experiment and simulations). One of those graphics is shown in Figure 13.
In the annex (last pages of this work), it can be seen the numerical results of the simulation for each of the ten different values of emissivity (as mentioned before, taken from the uniform distribution). In all the simulations a graphic very similar to Figure 13, almost indistinguishable, hence only one figure is exposed, but the whole numerical results can be studied in the annex. In Figure 13 it is possible to see upshots at the boundaries, clearly a result of the numerical scheme followed. However, for the purposes of this work, those numerical errors do not affect the comparison with the experimental results, because the experimental measurements where taken within the heaten sample, not at the boundaries (which are impinged by the laser beams).
The results obtained in the final average uncertainty quantification suggest that the deviation with respect to experimental results is small, and therefore that the full heat transfer equation is an adequate model to be used for analyzing the parameters of such equation, and then use the full heat transfer equation for inverse problems, that is, for inferring unknown heat parameters of a sample, within a margin of error (uncertainty), no bigger than 10 %.
Using the systematic modeling of continuous systems, a 3D heat equation was derived, from which it was considered that the emissivity was the parameter with uncertainty, as a working hypothesis. The results show an error between simulation and experiments of the order of ≈ 7%, with a possible standard deviation of ≈ ±10−5%. This error seems to be small enough to surmise that the numerical methods used could be used to determine the probability density function of the emissivity parameter with a small uncertainty.
Subscribe to our articles alerts and stay tuned.