SciELO - Scientific Electronic Library Online

vol.65 issue1Growth of a TiO2 nanotubular layer without presence of nanograss in a short timeThe Riemann-Silberstein vector in the Dirac algebra author indexsubject indexsearch form
Home Pagealphabetic serial listing  

Services on Demand




Related links

  • Have no similar articlesSimilars in SciELO


Revista mexicana de física

Print version ISSN 0035-001X

Rev. mex. fis. vol.65 n.1 México Jan./Feb. 2019  Epub Nov 09, 2019



Axisymmetric modelling of transient thermal response in solids for application to infrared photothermal radiometry technique

J. Hernández-Wonga 

J.B. Rojas-Trigosb 

U. Nogala 

V. Suárezc  b 

E. Marínb 

A. Calderónb 

aConsejo Nacional de Ciencia y Tecnología-Instituto Politécnico Nacional, Centro de Investigación en Ciencia Aplicada y Tecnología Avanzada, Legaria No. 694, 11500 Ciudad de México.

bInstituto Politecnico Nacional, Centro de Investigación en Ciencia Aplicada y Tecnología Avanzada, Legaria No. 694, 11500 Ciudad de México. e-mail:

cConsejo Nacional de Ciencia y Tecnología-Universidad Autonoma Metropolitana-Iztapalapa, Departamento de Química, ECOCATAL. Av. San Rafael Atlixco No. 186, 09340 Ciudad de México.


To induce temperature changes on the sample surface by the incidence of a monochromatic modulated light beam and detect the changes produced in the thermal radiation emission, is the basic principle of the infrared photothermal radiometry technique. Until now, in order to analyze the thermal response mathematical models based in an one-dimensional model were used considering a sample with a finite thickness and an infinite incidence surface, as well as the linear approximation of the Stefan-Boltzmann Law in the calculus of the heat losses due to thermal radiation. In this work, analytical and numerical models for the 2D heat diffusion in homogenous finite solid samples, are presented. These models were obtained by solving the heat diffusion equation, under cylindrical symmetry, considering mixed boundary conditions to include radiation and convection heat losses through the surfaces of the sample, and a monochromatic Gaussian excitation beam impinging on the front of the sample. The analytical models were obtained by solving the governing equations, considering the well-known linear approximation of the Stefan-Boltzmann law in the calculus of the heat losses due to thermal radiation. To analyse the effects of the nonlinearity of the heat losses by thermal radiation on the thermal transient response, in the numerical model it was taken into account the full expression of the Stefan-Boltzmann law, and the transport equation was solved numerically by means of the Finite Element Method (FEM). The analytical solution for the oscillatory thermal response reveals the close dependence of the thermal response on the ratio of thickness to the radius of the sample, represented by the form factor sf. Both, the analytical and the numerical solutions were employed to simulate the thermal response of homogenous materials, and compared with experimental results reported elsewhere by part of our same research group.

Finally, the difference between the thermal response predictions, from the analytical and numerical models, were analyzed.

Keywords: Heat transfer; infrared photothermal radiometry; thermal properties; Robin boundary conditions; finite element method

PACS: 44.05.+e; 66.70.-f; 44.25.+f; 44.40.+9; 02.70.Dh

1. Introduction

There are several techniques for the thermal characterization of condensed matter samples, which can be divided in two classes: Steady state and dynamical methods [1]. Among the latter one can mention the Photothermal (PT) techniques [12], which are based in the generation and detection of thermal waves. A particular type of PT technique is the one proposed in 1979 by Nordal and Kanstad for the spectroscopic analysis of solid and semi-solid materials [3], which they called Photothermal Radiometry (PTR), although a preliminary presentation was published a year earlier by these authors, for the analysis of powders [4]. The basic principle of this technique is to induce temperature changes on the sample surface by the incidence of a monochromatic modulated light beam and detect the changes produced in the thermal radiation emission. Since then, several variants of PTR technique have been reported, but its principle remains the same. Nowadays, this technique is named Modulated Photothermal Radiometry (MPTR) and, if a pulse of light replaces the exciting light beam the technique is named Pulsed Photothermal Radiometry (PPTR) [5,6]. Within the reports on the quantitative description of the modulated PTR technique, it is found that, in the eighties Santos and Miranda [7] presented a quantitative derivation of the modulated PTR production in both, frequency and time domain, using the one-dimensional model for the heat flow. From the model developed, Santos and Miranda successfully found the surface temperature fluctuation, neglecting the heat lost through the lateral walls of the sample and convection heat lost. A few years ago, Fuente et al. analyzed the ability of modulated PTR to retrieve simultaneously and accurately the optical absorption coefficient and the thermal diffusivity in homogeneous slabs [8]. Their analysis was based on considering the oscillating component of the temperature in the absence of heat losses obtained from the solution of one dimensional (1D) heat diffusion equation. Salazar et al. [9] reported the solution of 1D heat diffusion equation with non-adiabatic boundary conditions and concluded that the front and rear surface temperatures are affected by heat losses at low frequencies. Recently, Mart´ınez et al. [10] obtained a theoretical model for the PT signal in frequency domain in a 1D configuration considering convection and radiation heat losses (CRHL) in order to determine the thermal diffusivity of some low thermal conductivity materials, showing that CRHL should be taken into account for poor heat conductors at low frequencies.

The aim of this paper is to obtain analytical models closer to the physical reality for the transient thermal response in solids for application in infrared photothermal radiometry technique that takes into account the dimensions of the sample and heat losses by convection and radiation. Also, analyze the validity range of the linear approximation of the Stefan-Boltzmann’s Law, for the heat losses by thermal radiation, used in the deduction of the analytical model. This latter, by comparing results with the numerical model, obtained by means the finite element method considering the full expression of the Stefan-Boltzmann Law.

2. Theoretical Model

Consider a disc shaped sample (Fig. 1) of radius a and thickness ls that is heated by the absorption of a modulated monochromatic Gaussian beam. The generated heat flux density can be expressed as follows:

Φin(ρ,t)=(1-RS)I0exp(-2|ρw0|2)Ψmod(t). (1)

FIGURE 1 Schematic representation of the physical system used for the mathematical model. 

Where, I 0 and w 0 are the irradiance and the waist size of the excitation beam at the sample’s surface respectively, R s is the sample’s reflectivity, ρ is the radial cylindrical coordinate and Ψmod is a function describing the modulation. An electromagnetic-into-heat energy conversion efficiency equal to unity was assumed, as well as optical opacity of the sample. If the sample is homogenous and isotropic, then the homogenous parabolic heat diffusion equation (under cylindrical symmetry) models the heat transport, i.e.

2ΔT-1αstΔT=0. (2)

Here, ∆T = TT 0 is the variation of sample temperature, T, from room temperature, T 0. The solution of Eq. (2) is constrained by the following initial and boundary conditions:



Φcond+Φconv+Φrad|Si=Fi (3)

In Eq. (3), F i , represents the heat flux through the ith sample surface Si(i = 1,2,3 as schematically showed in Fig. 1). Also in Eq. (3), Φcond represents the conductive heat flux, given by the Fourier heat conduction law, Φconv is the convective heat flux given by the Newton’s law of cooling, and Φrad stands for the thermal radiation flux, described by the Stefan Boltzmann law.

2.1. Analytical model

The Stefan-Boltzmann law, a non-linear expression that can be written as follows, gives the thermal radiation flux.

Φrad=σbεs(T4-T04) (4)

Where σ b is the Stefan-Boltzmann constant, and ε s is the optical emissivity of the surface, at absolute temperature T. If ∆T ¿ T 0, by expanding Φrad as a Taylor series around T 0 a linear approximation of the Stefan-Boltzmann law is obtained [2]

Φrad=4σbεsT03ΔT (5)

In this way, the boundary conditions can be rewritten as:


limρ0ΔTexists (6)

Fi={0,i1Φin,i=1. (7)

Here, n^i are the unitary outward vectors normal to the surfaces, S i , k s is the thermal conductivity of the sample, and h s = h conv + h rad is the total thermal exchange coefficient, being hrad=4ε8σbT03 [2] the radiative heat transfer coefficient. The value for the convective heat transfer coefficient, h conv, was set on 4 × 10−4 W·cm−2·K−1 for all samples, consistent with the value reported by Martinez et al. [10]. To solve the model uprising from Eq. (2) and (3), considering Eq. 4-7, the dimensionless variables r = ρ · a −1 and ς = z · l s −1 will be employed, along with the parameter σ 0 = w 0 · 2 −1/2 α −1. Therefore, using the Fourier-Bessel orthonormal basis, ∆T can be spanned as:


Rn(r)=2j0(vnr)j02(vnj12(vn) (8)

In Eq. (8), Θ n are the expansion coefficients of ∆T in the Fourier-Bessel orthonormal basis, defined through the set of the Bessel functions of first kind, J 0(v n r) and the associated eigenvalues, v n , are given by the roots of the next transcendental equation:


Biρ=ahs2ks (9)

Here, Bi ρ stands for the radial Biot number. Substituting Eq. (8) into the previous Eqs. (2-6) leads to the subsequent reduced problem:

2ζ2Θn-1πfctΘn-λn2Θn=0limt0Θn =0 (10)


ζΘn+BizΘn|ζ=1=0. (11)

In Eq. (11), f c α s · ( πl82)−1 is the characteristic frequency that represents the modulation frequency at which the thermal diffusion length µ s = (α s /πf) 1/2 matches the sample thickness l s [11]; Bi z = l s h s · k s −1 defines the axial Biot number and λ n = s f v n , being s f l s · α −1 the so-called shape factor of the cylindrical sample. Additionally, the coefficients Hn0 are determined as follows:

Hn0-(1-Rs)lsI0ks×[01Rn(r)exp(-|rσ0|2)rdr] (12)

2.1.1. Oscillatory regime

If the modulation function Ψ mod is harmonic, then the reduced 1D-problem (Eq. (10)) can be re-written in frequency domain, using the Unitary Fourier Transform (UFT), F, leading to the following partial differential equation (PDE) problem:

Θ~n=Θ~nosc=mZΘ~nm(ωm,ζ)δ(ω-ωm) (13)

2ζ2Θ-nm-γnm2Θ-nm=0 (14)

ζΘ-nm-BizΘ-nm|ζ=0=Hn0FΨmod;           ζΘ-nm+BizΘ-nm|ζ=1=0      (15)

The accent marks refer to the dependant variables in frequency domain. The damping coefficients γnm2=λn2+2imffc-1 are closely related to the complex diffusion coefficients σm=(1+i)μm-1, being µ m = µ s · m −1/2 the thermal diffusion length corresponding to the mth harmonic contribution. After a straightforward calculation, the following results were obtained:




Dnm=CmHn0 (16)

ΔTosc=nNmZDnmUnm(r,ζ)exp(iωmt) (17)

In Eq. (16) and (17), {Cm} are the expansion coefficients of Ψmod in the Fourier basis (ω m = 2πmf). By definition, the Biot number - for a cylinder-shape sample - is related to the radial and axial Biot numbers through Eq. (18):

Bi=[1Biρ+2Biz]-1 (18)

When the ratio BizBip-11, or equivalently s f ≪ 1, the damping (and therefore the thermal regime), is solely characterized by ffc-1 as usual, consistent with a semiinfinite sample. Otherwise, the effects of the radial heat diffusion cannot be neglected. To exemplify the effect of different values of the shape factor s f , in Fig. 2, the magnitude of γ nm coefficients are calculated for a balsa wood sample (α s = 23 × 10−6 m2·s−1; k s = 0.11 W·m−1·K−1).

FIGURE 2 Calculated values of |γnm|for a balsa wood sample, for: (a) sf = 1, (b) sf = 0.4, (c) sf = 0.2 and (d) sf = 0.133. Here the modulation frequency was set on 100 mHz. 

As can be seen from Fig. 2, by increasing the shape factor s f , the maximal magnitudes of the damping coefficients increase -from 40 to 300, for a given value of f c . The previous implies that when the shape factor increases, the lateral surface heat transfer (heat losses) increases too, and the radial heat flow through the sample tends to acquire a greater importance.

2.1.2. Transient regime

Employing the Duhamel’s Theorem [12], and considering an arbitrary -but integrable- modulation function, the solution of the 1D reduced problem, Eq. (10), is given as follows:

Θn=Θntran=Hn0m=1EnmZm(ζ)0tΨmod(τ)texp(-πfcϖnm2t)|t-τdτ (19)

Where: ϖnm2=βm2+λn2, and:

Zm=βmcosβmζ+Bizsinβmζ (20)

Enm=01Zmχndζ (21)

The Xn function stands for the stationary solution of the reduced problem, specified by the following expression:


η(±)=ηn±Biz (22)

The eigenvalues {β m }, corresponding to the eigenfunctions {Z m }, are determined by the zeros of the following transcendental equation:

tanβm=2βmBizβm2-Biz2;          βmR+/{Biz} (23)

This equation most be solved numerically, however, it can be demonstrated that β m as long as Bi z ≪ 1. If Ψ mod is continuous and can be spanned in the Fourier basis, then the convolution product in Eq. (19) is calculated to be:

0tΨmod(τ)texp(-πfcϖnm2t)|t-τdτ=-πfcϖnm2×kZCk[exp(iωkt)-exp(-πfcϖnm2t)ϖnm2+2ikffc-1] (24)

Again, {C k } are the expansion coefficients of Ψ mod in the Fourier basis. Otherwise, if Ψ mod has discontinuity points, the integral in Eq. (19) must be solved by integration by parts over the continuity subdomains.

2.2. Numerical Model

In this work, COMSOL Multiphysics (CMP) [13] has been used to solve numerically the problem described earlier in §2, by the Finite Element Method (FEM). To solve numerically Eq. (2) it is necessary to define, not only, the physical model (equations, and boundary and initial conditions), but also global parameters, the geometry (coordinate system, symmetries, physical boundaries, etc.), the material properties, and auxiliary functions (i.e. the power density distribution and the modulation function). To simulate the power density distribution in CMP, a modulated Gaussian function was defined.

Φmodsq=(1-Rs)I0e-2((x2+y2)/w02)Ψmodsq(t) (25)

In Eq. (25) Rs is the reflectivity, I 0 is the irradiance, w 0 is the laser spot radius, x and y are the independent variables, and ψmodsq is the unitary rectangular wave train, multiplied by the left-continuous unitary step function. Notice that Φmodsq in Eq. (25), is a particular case of Φ in , defined in Eq. (1). For solving the presented physical problem using CMP, the heat transfer in solids module was employed, in which the dynamical equation is given by the convection-diffusion equation [14].

ρCptΔT+ρCpu(T)=(k(T))+Q (26)

In Eq. (26), ρ, C p and k are the mass density, the specific heat (at constant pressure), and the thermal conductivity of the material, respectively, u is the translational motion of the surrounding fluid (here, u is assumed to be null); and Q the heat source term. Here, both u and Q are assumed to be null, so that for a homogeneous and isotropic sample with constant thermal conductivity, the above equation becomes the wellknown homogeneous heat diffusion equation. The following boundary conditions have been considered for the numerical solution by FEM:

n^i(kT)+hconv(T-T0)+εsσb(T4-T04)|Si=Fi (27)

In Eq. (27), the full expression of the Stefan-Boltzmann law was considered, instead of the linear approximation of the previous section, so that the boundary conditions become non-linear. For the FEM processing, a tetrahedral mesh with two different element sizes was constructed over the domain: One element size was calibrated for “general physics” with a “normal” distribution for the cylinder’s area, and the second one using a finer element distribution for both surfaces (z = 0,l s ) as showed in Fig. 3.

FIGURE 3 Meshing of the elements which defines the spatial domain for the numerical solution (units are not shown). 

Once the model and the meshing were defined, the numerical solver was selected to generate the simulation. In this case, a direct multifrontal massively parallel sparse direct (MUMPS) solver was chosen, due to its high calculus resolution capabilities and efficient memory management in parallel computing [15]. The MUMPS solver was implemented in combination with a time stepping feature using a backward differentiation formula (BDF) [16] method in strict mode, considering second-order and fourth-order BDF to ensure a fast and accurate numerical convergence in the numerical solution of the model; obtaining a solution at the edges of the temporal subintervals, during the stipulated simulation time range. Finally, a parametric sweep was used to obtain the solutions of the model, described in the present section, at different frequencies with the same boundary and initial conditions.

3. Results and Discussion

In the present study, four different samples were considered: Two of them were low thermal conductivity balsa wood (BW) and high-density polyethylene (HDPE) materials, and the other two were high thermal conductivity cooper (Cu) and aluminium (Al) samples. The geometrical characteristics and reported physical properties of the samples are displayed in Table I [17,18].

TABLE I Geometrical and physical parameters considered for calculations. 

Sample Material α
×10−3 m
×10−3 m
εs Rs αs
M1 BW 5 1 0.91 0.20 0.22 0.11 4.0 0.070
M2 HDPE 5 1 0.93 0.20 0.21 0.52 4.0 0.067
M3 Al 5 1 0.03 0.96 93 238 4.0 29.61
M4 Cu 5 1 0.05 0.60 116 400 4.0 36.92

In all cases, and for the analytical and numerical models, the modulation function Ψmod considered was the unitary rectangular wave train, with oscillation period τ mod = f −1 and wide τ mod /2, multiplied by the left-continuous unitary step function.

3.1. Analytical calculations: Oscillatory regime

The temperature distributions (surface thermograms) at ζ = 1 were calculated - using the analytical model, Eq. (10b) - for the samples defined in Table I considering, a modulation frequency, f, of 100 mHz, a value of 0.2 for the shape factor, s f , and the following values for the other input parameters: I 0 = 12.7 × 104 W·m−2; σ 0 = 4.24 × 10−2; and T 0 = 300 K. In Fig. 4, the obtained results in temperature difference are shown at four different simulation times (0,(1/4)f −1 ,(1/2)f −1 ,(3/4)f −1).

FIGURE 4 Simulated thermograms for: a) M1, b) M2, c) M3, and d) M4 samples, at four different simulation times. In each figure, 2.5 s equals f −1 /4. The units in colour bars are K. 

To have a better insight on the time evolution of the radial temperature distribution, the full-width at half-maximum (FWHM) were calculated. They are shown in Fig. 5 as a function of time, and for values of 0.2 and 0.4 for the shape factor s f .

FIGURE 5 FWHM values, at different times, of: M1 (empty circles), M2 (empty squares), M3 (empty-up triangles) and M4 (empty-down triangles) samples. (a) s f =0.2 and (b) s f =0.4. In both cases f =100 mHz. 

An important feature in the surface thermograms in all samples is that the radial distribution of temperature has a Gaussian-like profile, inherit from the power distribution. However, at t = 7.5 s (3/4 fold τ mod), the broadening of the radial temperature distribution covers the hole surface in M3 and M4 samples, for the particular choice of the shape factor s f = 0.2. This effect cannot be attributed solely to the radial heat flux, but also to a combination of the magnitudes of the radial and axial heat fluxes, and finally, to the values of Bi p , and ffc-1 since the shape factor is the same for all samples. Contrary, for M1 and M2 samples, the maximum FWHM values occur just at one-half of τ mod, Fig. 5a). However, the minimal FWHM value in sample M1 does not occur at integer multiples of τ mod, but coincidently at 3/4 fold τ mod; while the evolution of the FWHM values of the M2 sample is quite symmetric. Figure 5b) shows a significant change in the behaviour of FWHM, in all cases, with an important decrease in the FWHM value of more than 20% when the shape factor changes from 0.2 to 0.4, which shows the importance of the shape factor in the temperature distribution in the sample. In Figs. 6 and 7, the radial temperature distribution of samples M1 and M4 (corresponding to the highest and lowest Bi ρ values) are displayed. In all cases, s f = 0.2, and f = 100 mHz.

While sample M1 shows the maximal broadening at onehalf of the time period, sample M4 exhibits a significantly larger broadening in the radial distribution, reaching the maximal broadening at t = 7.5 s. This feature in the radial temperature distribution in M4 sample predicts a “homogenization” of the superficial temperature of the sample, consistent to a 1D-heat diffusion. As consequence of the particular radial distribution, the thermal wave front differs from one sample to another. This is an important issue for the analysis of the modulated photothermal radiometric signal. Taken the normalized magnitude of the spatial average |〈∆T〉|, the form of the calculated MPTR signal is clearly distinguishable for each sample (Fig. 8). As the Bi ρ decreases and f c increases, the MPTR signal becomes symmetrical. Notice also a temporal shift, due to the axial heat diffusion, and therefore, to the f c value.

FIGURA 6 Calculated radial temperature distribution of sample M1, at: t =0 s (solid line), t =2.5 s (dash-asterisk line), t =5 s (dash-circle line) and t =7.5 s (dashed line). 

FIGURA 7 Calculated radial temperature distribution of sample M4, at: t =0 s (solid line), t =2.5 s (dash-asterisk line), t =5 s (dash-circle line) and t = 7.5 s (dashed line). For a better visualization, in this case the vertical axis is in logarithmic scale. 

FIGURE 8 Normalized magnitude of the spatial average of the temperature oscillations of samples: (a) M1, (b) M2, (c) M3 and (d) M4. In all cases, s f =0.2 and f =100 mHz. 

Naturally, the influence of other parameters must be taken into account, for example, the magnitude of the absorbed energy flux averaged over the surface of incidence -considered in the coefficient Hn0, Eq. (7b)-decreases as the ratio Rsks-1 decreases. This fact, in combination with smallest values of Bi ρ and higher thermal diffusivities (as it was mentioned in previous lines), provokes that the radial temperature distribution gets broader. This result confirms that the radial dependence of the temperature distribution should not be neglected if the shape factor is not small, in particular in the thermally thin regime. Even so, for metallic samples like M3 and M4, the temperature oscillations could be difficult to be measured accurately because of the highly reflective and low IR emissivity surfaces, resulting in small thermal oscillations. A possible solution could be increasing of the ratio hsks-1, by the deposition of high-emissivity, low-reflectivity conductive thin coatings on the free surfaces of the sample (just as Martínez et al. does [10]). However, a precaution must be taken when the coating and the sample respond in the same thermal regime (in this case, the thermally thin regime), being necessary to consider a three layers’ system in the heat diffusion model [19].

3.2. Numerical calculations: Transient regime

Once that the FEM simulations were obtained to discuss the behaviour of the transient thermal response of all samples, the surface point P = (0,0,l s ) was chosen to calculate the ∆T vs. t curves. In Fig. 9, the results at P, calculated by FEM for M1 and M2 samples, are shown for f = 100 mHz. These samples, having high emissivity and low thermal conductivity, respond closely to thermally thin regime, and qualitatively alike as the thermal response described by Martínez et al. [10], analysed by using the 1D model. Next, Fig. 10 shows the ∆T vs. t curves for M3 and M4 samples, at P.

FIGURE 9 T calculated at f =100 mHz, by FEM, for: (a) M1, and (b) M2. Here, ∆t =0.1 s. 

FIGURE 10 T calculated at f =100 mHz, by FEM, for: (a) M3, and (b) M4. Here, ∆t =1 s. 

The small values of the temperature variation in samples M3 and M4 in relation to those of the samples M1 and M2 are because large values of k give small resistance to heat conduction (low values of Bi) and, therefore, small temperature gradients. This is the same reason why the predicted values of the temperature variation in sample M1 are greater than those corresponding to sample M2. Unlike this, the highest values of the temperature variation of the sample M4 with respect to those of sample M3 are due to the difference between the reflectivity values Rs of these samples (see Table I).

In Figs. 11 and 12, the amplitudes of the temperature oscillations as function of the modulation frequency, are shown for all samples, calculated by means of the analytical and numerical models.

FIGURE 11 Amplitudes of the temperature variations, as function of f, for M1 sample (full circle (a), empty circle (n)), and M2 sample (full square (a), empty square (n)). Here: (a) = analytical, (n) = numerical. 

FIGURE 12 Amplitudes of the temperature variations, as function of f, for M3 sample (full up-triangle (a), empty up-triangle (n)), and M4 sample (full down-triangle (a), empty down-triangle (n)). Here: (a) = analytical, (n) = numerical. 

The results of Fig. 12 show a good agreement between the analytical and numerical models for the sample M4, for which the condition ∆TT 0 is well fulfilled and the Stefan-Boltzmann law behaves with good approximation in a linear way. However, when the condition ∆TT 0 is not adequately fulfilled, the correspondence between both models is not so good, as shown by the results for samples M1, M2 y M3 in Figs. 11 and 12. In these cases, it is convenient to analyze the problem using the numerical model to obtain reliable results. On the other hand, it can be notice that M3 and M4 samples behaves almost as thermally thin samples, since the A vs. f (inset in Fig. 12) keep a f −1 dependence, unlike the behaviour of the samples M1 and M2 whose characteristic frequencies are around 5 times smaller than those of M3 and M4, see Table I.

A graph with the comparison among the numerical and analytical models with the experimental results of balsa wood and plasticine reported in Ref. [10], is shown in Fig. 13. A good agreement among them is observed, in particular, in the regions where the behavior is as a straight line the slopes are very similar, which is where the values of thermal diffusivity are obtained.

FIGURE 13 Predicted radiometric signal, calculated by analytical (solid black line), and FEM (black dashed line). Here, the full diamonds and fill circles correspond to the experimental available data for balsa wood and plasticine, respectively [10]. 

4. Conclusions

By solving the 2D heat diffusion equation, analytical and numerical models were obtained to describe transient and oscillatory thermal response in homogenous and finite solid samples (with cylindrical symmetry). Heat losses due to radiation and convection through front, rear and lateral surfaces were considered and a monochromatic Gaussian excitation beam that impinge on the front face of the sample. The analytical solution for the oscillatory thermal response reveals the close dependence of the thermal response on the dimensions of the sample, represented by a form factor s f , showing that when s f ≪ 1 the problem is reduced to 1D of the semi-infinite sample model. In any other case, the lateral surface heat transfer (heat losses) cannot be neglected, making its consideration necessary in a complete description of the physical situation. The results obtained for the transient thermal response are in congruence with the experimental results reported by Martínez K et al [10]. In the transient thermal response it was obtained that when the condition ∆TT 0 is well fulfilled the results obtained between the analytical and numerical models agree very well showing the utility of the Stefan-Boltzmann law linearization considered in the analytical model. However, as long as the thermal response moves away from the thermally thin regime the correspondence of the results obtained between both models becomes increasingly less, until it becomes necessary to analyze the thermal response using only the numerical model to obtain reliable results.


This work has been partially supported by Consejo Nacional de Ciencia y Tecnología (CONACyT) and Secretaría de Investigación y Posgrado (SIP) from Instituto Politécnico Nacional (IPN), both from México.


1. D.P. Almond, P.M. Patel, Photothermal Science and techniques, first ed., Chapman and Hall, (London 1996). [ Links ]

2. E. Marín, Basic Principles and Recent Developments, in: E. Marín (ed.) Thermal Wave Physics and Related Photothermal Techniques, Transworld Research, Kerala, (India 2009). [ Links ]

3. E. Nordal and S.O. Kanstad, Photothermal Radiometry, Physica Scripta 20 (1979) 659-662. [ Links ]

4. S.O. Kanstad and E. Nordal, Power Technology 22 (1978) 133-137. [ Links ]

5. A.C. Tam and B. Sullivan, Appl. Phys. Lett. 43 (1983) 333-335. [ Links ]

6. J. Martan, Rev. Sci. Instrum. 86 (2015) 014902 1-9. [ Links ]

7. R. Santos and L.C.M. Miranda, J. Appl. Phys. 52 (1981) 4194-4198. [ Links ]

8. R. Fuente, E. Apiñaniz, A. Mendioroz, and A. Salazar (2011), J. Appl. Phys. 110 033515 1-9. [ Links ]

9. A. Salazar, R. Fuente, E. Apiñaniz, A. Mendioroz, J. Appl. Phys. 110 (2011) 33516 1-8. [ Links ]

10. K. Martínez et al., Int. J. Therm. Sci. 98 (2015) 202-207. [ Links ]

11. A. Calderon, R.A. Muñoz Hernández, S.A. Tomas, A. Cruz-Orea, F. Sánchez Sinencio, J. Appl. Phys. 84 (1998) 6327-6329. [ Links ]

12. F. John, Partial Differential Equations, fourth ed., SpringerVerlag, (New York 1982). [ Links ]

13. W.R. Pryor, Multiphysics Modeling using Comsol 4: A First Principles Approach, Mercury Learning and Information, Dulles, Virginia USA (2012). [ Links ]

14. V. Suarez et al., Appl. Rad. Isot. 83 (2014) 260-263. [ Links ]

15. MUMPS: Multifrontal Massively Parallel sparse direct Solver, MUMPS: Multifrontal Massively Parallel sparse direct Solver, (2017) (accessed 27 April 2018). [ Links ]

16. Comsol Multiphysics Reference Guide V 4.3a, 1998-2012 COMSOL. [ Links ]

17. A.L. Edwards, A Compilation of Thermal Property Data for Computer Heat-Conduction Calculations, Ucrl-50589, University Of California Lawrence Radiation Laboratory USA (1969). [ Links ]

18. Y.S. Touloukian, R.W. Powell, C.Y. Ho, M.C. Nicolaou, Thermophysical Properties of Matter, Vol. 10, IFI/Plenum, (New York 1973). [ Links ]

19. K.D. Cole Int. J. Thermophys. 25 (2004) 1567-1584. [ Links ]

Received: May 04, 2018; Accepted: September 04, 2018

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License