Modeling of active thermography through uncertainty quantification of parameters of the heat transfer equation

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.

Introduction

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 [4]. 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 [5].

Problem possing and objectives

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) [6]. 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

$\frac{\partial Q}{\partial t}=\Phi -q\text{ (1)}$

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 [6], which can be expressed as

$\frac{\partial \theta }{\partial t}+\frac{h}{{\rho }_m\cdot C_p\cdot th}\theta -\frac{\Phi }{{\rho }_m\cdot C_p\cdot V}=0\text{ (2)}$

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, C_p 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 [6].

${\theta }_{\uparrow }=\frac{\Phi }{Ah}\left(1-\text{exp }\left(-\frac{t}{\tau }\right)\right)\text{ (3)}$

where A is the total area of the sample and

$\tau =\frac{th\cdot {\rho }_m\cdot C_p}{2h}\text{ (4)}$

is the relaxation (transient) time of the system. Eq.3 has been used to calculate the specific heat capacity ·C_p for thin slabs of known thickness [6], 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 [6].

Successful as it is the approximation mentioned in [6] (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))[1], 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.,) [2]. 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.

Methodology

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.

Theoretical fundaments

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 [6].

Conceptual model

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.

Mathematical derivation of the model

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,

${\rho }_m\frac{\partial E}{\partial t}=\nabla \cdot \text{q+}q\text{o, (5)}$

the relationship between heat and temperature of Eq.35, and Eq.31, and assuming no heat sources, it is obtained

${\rho }_m\cdot C_p\frac{\partial T}{\partial t}=-\nabla \cdot (-{\text{λ}}_m\nabla T)\text{ (6)}$

or, to show explicitly the parabolic form of this partial differential equation,

${\rho }_m\cdot C_p\frac{\partial T}{\partial t}=-\nabla \cdot ({\text{λ}}_m\nabla T)\text{ (7)}$

Boundary conditions

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 q₀, which in general is

n · (λm∇T) = q0 + hconv · (Tenv − T(x,y,z,t)) + σεm(Tenv4− T4) (8)

For the boundaries with q₀=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. q₀ is represented by the laser beam through the extremes of the slab, as shown in Figure 1, namely Φ1 Φ2.

Figure 1: (A) Schematics showing how heat is injected in the laboratory experiments.

Initial conditions

The initial condition, as mentioned, is T(0) = temperature of air.

Summary of notation of mathematical model

T(x,y,z,t) Temperature [ K ]

Tenv Air temperature [ K ]

Ρm Density of material [kg/m3]

Φ Incident Power[W · m⁻²]

hconv λm σ εm Convective heat transfer coefficient [W/(m²·K)]

Thermal conductivity [W/(m·K)]

Stefan-Boltzmann constant [W·m⁻²· K⁻⁴]

Total emissivity of the material [1]

ε_mb Total emissivity of the material darked side [1]

C_p heat capacity of the material (at constant pressure) [J/(kg · K)]

Uncertainty in the coefficients

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)[11], of the Partial Differential Equation, namely for λ_m,h_conv,ε_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 [12].

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) [12].

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 [12].

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)[12].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.

Discretization of mathematical model (numerical model)

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 [13].

The FEM naturally incorporates a broader spatial extent, thus FEM can use a coarser mesh, compared with FD [14].

The basic mixed finite element (MFE) method [15], 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 [16].

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

Computational implementation-comsol multiphysics r

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 [17]. In particular, using the PDE mode for time dependent analysis in the coefficient form. An alternative can be using the FEniCS implementation [18]. For concreteness, it is exposed here how has been done the implementation in Comsol v 3.5 a.

Parameters

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

Figure 2: Constants of model, as defined in Comsol 3.5a. In particular, notice the collection of emissivities.

T(x,y,z,t) = u

T_env = T env

ρ_m = rhom

Φ = Phi

h_conv = h conv

λ_m = lambda m

σ = sigma

ε_mi = Epsilonm

ε_mb = Epsilon-mb

C_p = Cp

Notice that i ∈ [2,11] is an index for the emissivities taken from the Uniform distribution.

Geometry

The selected geometrical object represents an experimental metal sheet (Figures 3,4)

Figure 3: Definition of hexahedron for simulation.

Figure 4: Sheet as a result of hexahedron definition.

Scalar expressions

The definition of the coefficients of our equation are contained in the section scalar expressions, as can be seen in Figure 5.

Figure 5: Scalar expressions.

Comsol equation

Eq.7 is transformed into the comsol form, using the scalar expressions to define the coefficients, as can be shown in Figure 6.

Figure 6: Subdomain settings is where the coefficients of the Comsol equation are written, using the definitions of the scalar expressions.

Initial conditions, and lenght of simulation

The initial temperature of the material is considered the same as the surrouding air, that is T(0) = T_env. The length of the simulation is for 5.3 seconds, which was seen enough to obtain a stationary solution.

Boundary conditions

There is only one type of conditions: Robin. However, in most boundaries T(t = 0) = T_env and only in boundaries 2 and 5 T = Φ, for the duration of the numerical experiment, as can be seen in Figures 7,8.

Figure 7: Robin boundary conditions, Φ =6 0.

Figure 8: Robin boundary conditions, Φ = 0.

Mesh

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

Figure 9: Domain mesh.

Solver Parameters: The chosen solver es UMFPACK.

Testing of simulation

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.

Figure 10: Qualitative validation.

Validation of model

The model of Eq.7, with the boundary conditions shown in Eqs.8 and 9, with the initial condition T(t = 0) = T_env = 300.15K, with h_conv = 10[W/(m² · K)], and the weak condition of injection heat Φ = 1000[W/m²] was set up, according to [1]. The numerical solution (shown in Figure 11) has the expected radial behavior as described in [1].

Figure 11: Result of testing of computational model.

Study case and results

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.

Figure 12: The experimental setup and the numerical measurement.

Results and Discussion

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.

Figure 13: Compvarison of the family of simulations (continuos line) with the experiment of Prof. Rumen (blue circles). The error in all experiments was ≈ 7 %.

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

Conclusions

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⁻⁵%. 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.

Appendix, Systematic Modelling, ANNEX FOLLOWS.

1. Benzerrouk S (2011) Active and passive thermography for the detection of defects in green-state powdermetallic compacts. (Unpublished doctoral dissertation). Worcester Polytechnic Institute. Link: http://bit.ly/2KqVJDV
2. Cannas B, Carcangiu S, Concu G, Trulli N (2012) Modeling of active infrared thermography for defect detection in concrete structures. In Proceedings of the 2012 comsol conference. Link: http://bit.ly/2r1OwTQ
3. Ciampa F, Mahmoodi P, Pinto F, Meo M (2018) Recent advances in active infrared thermography for non-destructive testing of aerospace components. Sensors 18: 609. Link: http://bit.ly/2XpxtY3
4. Yao Y, Sfarra S, Laguela S, Ibarra-Castanedo C, Wu JY, et al. (2018) Active thermography testing and data analysis for the state of conservation of panel paintings. International Journal of Thermal Sciences 126: 143-151. Link: http://bit.ly/35fBqRV
5. Borchartt TB, Conci A, Lima R, Resmini R, Sanchez A (2013) Breast thermography from an image processing viewpoint: A survey. Signal Processing 93: 2785-2803. Link: http://bit.ly/2OhvmkU
6. Marin E (2009) Linear relationships in heat transfer. Lat Am J Phys Educ 3. Link: http://bit.ly/2qYaEyD
7. Valiente H, Delgado-Vasallo O, Abdelarrague R, Calderon A, Marin E (2006) Specific heat measurements by a thermal relaxation method: Influence of convection and conduction. International journal of thermophysics 27: 1859–1872. Link: http://bit.ly/2OfyugY
8. Marin E, Delgado-Vasallo O, Valiente H (2003) A temperature relaxation method for the measurement of the specific heat of solids at room temperature in student laboratories. American Journal of Physics 71: 1032–1036. Link: http://bit.ly/2CRt0nF
9. Lopez-Falcon D, Diaz-Viera M, Herrera I, Rodriguez-Jauregui E (2008) Systematic formulation of continuum systems: Theoretical modeling of combustion fronts in porous media. Numerical Modeling of Coupled Phenomena in Science and Engineering. Link: http://bit.ly/2Qtv6Sp
10. Diaz-Viera MA, Ortiz-Tapia A, Hernandez-Perez JR, Castorena-Cortes G, Roldan-Carrillo T, et al. (2019) A flow and transport model for simulation of microbial enhanced oil recovery processes at core scale and laboratory conditions. International Journal of Numerical Analysis and Modeling 16: 63-96. Link: http://bit.ly/33SWhKm
11. Sullivan TJ (2015) Introduction to uncertainty quantification. 63. Link: http://bit.ly/2OnpC97
12. ISO I, OIML B (1995) Guide to the expression of uncertainty in measurement. Geneva Switzerland.
13. Sewell G (2005) The numerical solution of ordinary and partial differential equations. John Wiley and Sons 75. Link: http://bit.ly/2CYBFEx
14. Simpson M, Clement T (2003) Comparison of finite difference and finite element solutions to the variably saturated flow equation. Journal of hydrology 270: 49-64. Link: http://bit.ly/33WCB8B
15. Raviart PA, Thomas JM (1977) A mixed finite element method for 2-nd order elliptic problems. Mathematical aspects of finite element methods 292–315. Link: http://bit.ly/2Qukh2K
16. Zienkiewicz OC, Taylor RL, Zienkiewicz OC, Taylor RL (1977) The finite element method. McGraw-hill London.
17. COMSOL Multiphysics (2008) User’s guide. version 3.5 a. COMSOLAB.
18. Alnæs MS, Blechta J, Hake J, Johansson A, Kehlet B, et al. (2015) The fenics project version 1.5. Archive of Numerical Software 3. Link: http://bit.ly/2XkN9fc
19. Dıaz-Viera MA, Ortiz-Tapia A (2018) Lecture notes in mathematical, numerical and computational modeling of flow and transport in porous media. UNAM Posgrado de Ciencias de la Tierra.
20. Bergman TL, Incropera FP, Lavine AS, Dewitt DP (2011) Introduction to heat transfer. John Wiley and Sons. Link: http://bit.ly/2r3j9Iq