2.1. Analytical model

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

Where σ _b is the Stefan-Boltzmann constant, and ε _s is the optical emissivity of the surface, at absolute temperature T. If ∆T ¿ T ₀, by expanding Φ_rad as a Taylor series around T ₀ a linear approximation of the Stefan-Boltzmann law is obtained [²]

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

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 [²] the radiative heat transfer coefficient. The value for the convective heat transfer coefficient, h _conv, was set on 4 × 10⁻⁴ W·cm⁻²·K⁻¹ for all samples, consistent with the value reported by Martinez et al. [¹⁰]. To solve the model uprising from Eq. (2) and (3), considering Eq. 4-7, the dimensionless variables r = ρ · a ⁻¹ and ς = z · ^l _s ⁻¹ will be employed, along with the parameter σ ₀ = w ₀ · 2 ^−1/2 α ⁻¹. Therefore, using the Fourier-Bessel orthonormal basis, ∆T can be spanned as:

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 ₀(v _n r) and the associated eigenvalues, v _n , are given by the roots of the next transcendental equation:

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

In Eq. (11), f _c ≡ α _s · ( πl82)⁻¹ 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 [¹¹]; Bi _z = l _s h _s · ^k _s ⁻¹ defines the axial Biot number and λ _n = s _f v _n , being s _f ≡ l _s · α ⁻¹ the so-called shape factor of the cylindrical sample. Additionally, the coefficients Hn0 are determined as follows:

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:

The accent marks refer to the dependant variables in frequency domain. The damping coefficients γnm2=λn2+2imf∙fc-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 m − th harmonic contribution. After a straightforward calculation, the following results were obtained:

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

When the ratio Biz∙Bip-1≪1, or equivalently s _f ≪ 1, the damping (and therefore the thermal regime), is solely characterized by f∙fc-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⁻⁶ m²·s⁻¹; k _s = 0.11 W·m⁻¹·K⁻¹).

FIGURE 2 Calculated values of |γ_nm|for a balsa wood sample, for: (a) s_f = 1, (b) s_f = 0.4, (c) s_f = 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 [¹²], and considering an arbitrary -but integrable- modulation function, the solution of the 1D reduced problem, Eq. (10), is given as follows:

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

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

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

This equation most be solved numerically, however, it can be demonstrated that β _m ≈ 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:

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

In Eq. (25) Rs is the reflectivity, I ₀ is the irradiance, w ₀ 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 [¹⁴].

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:

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 [¹⁵]. The MUMPS solver was implemented in combination with a time stepping feature using a backward differentiation formula (BDF) [¹⁶] 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.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 ₀ = 12.7 × 10⁴ W·m⁻²; σ ₀ = 4.24 × 10⁻²; and T ₀ = 300 K. In Fig. 4, the obtained results in temperature difference are shown at four different simulation times (0,(1/4)f ⁻¹ ,(1/2)f ⁻¹ ,(3/4)f ⁻¹).

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 ⁻¹ /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 f∙fc-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 Rs∙ks-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 hs∙ks-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 [¹⁰]). 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 [¹⁹].

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. [¹⁰], 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 ∆T ≪ T ₀ is well fulfilled and the Stefan-Boltzmann law behaves with good approximation in a linear way. However, when the condition ∆T ≪ T ₀ 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 ⁻¹ 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. [¹⁰], 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 [¹⁰].