# Correlation between hardening depth and thermal parameters based on photothermal measurements/Fototermiliste mootmiste alusel leitud korrelatsioon tahkumise sugavuse ja termiliste parameetrite vahel.

1. INTRODUCTION

The photothermal detection enables characterization of materials in near-surface layers. New capabilities in machining processes and decreasing tolerances down to the nanometer scale have lead to material treatments, affecting only the layers near the surface. These near-surface areas influence the functional behaviour, quality and lifetime of the workpiece. Methods for the assessment of these areas should be non-destructive and meet standards concerning accuracy, measuring time, robustness and costs.

During the last years the photothermal techniques (Fig. 1), including various detector options, have gained a substantial progress.

[FIGURE 1 OMITTED]

Photothermal measuring techniques exploit the time-dependent heat flow, induced by a time-varying heat source, applied to the sample. This can be appropriately realized either by periodic laser heating or pulsed illumination. Both excitation schemes initiate propagating temperature oscillations, the so-called thermal waves. Their penetration into the inspected sample and different interactions with subsurface inhomogeneities open a vast field of non-destructive photothermal testing [1]. Today, photothermal measurements are applied to various materials and composites including coatings (thickness, adhesion strength, local defects), hardness profiles, residual stresses, wear, thermal influences and surface-affecting manufacturing techniques. Photothermal imaging includes quick, reliable and non-contact near-surface inspection with pattern sizes in the mm range as well as photothermal microscopy with [micro]m resolution [2].

2. PHYSICAL PRINCIPLES

2.1. Thermal waves and their detection

Fundamental photothermal and photoacoustic effects, a suitable experimental set-up and the two most important detection principles (photothermal radiometry and optical beam deflection) have been explained in detail in [3]. The corresponding basic theory can be traced back to the late 1970s ([1] and references therein), whereas intensive research, aiming at an application in production processes, started in the early 1990s [4-7].

In all photothermal measurements, an intensity-modulated light source (e.g. a laser beam) excites the investigated sample at its surface, as shown in Fig. 1. Optically opaque materials absorb the light directly at their surface, and, subsequently, the sample heating will also be modulated, causing a propagating thermal wave. The wave spreads within the sample, leading to a temperature rise by a few K at the surface. The penetration depth of the thermal wave is determined by the thermal diffusion length [[micro].sub.th] = [square root of (2[alpha]/[omega]), which is a function of the modulation frequency [omega] and the thermal diffusivity [alpha]. For metals, the penetration depth varies from a mm at a few Hz to some tens of [micro]m at a kHz-modulation [1]. Thermal inhomogeneities scatter the thermal wave and the surface temperature will differ from the value occurring in a homogeneous material. Photothermal measuring techniques map these (modulated and time-dependent) temperature changes at the surface due to hidden structures and near-surface changes of material properties [8].

2.2. The thermal diffusion equation

The heat propagation within a sample is given by the thermal diffusion equation [1,8,9]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1) ## which contains the locally-dependent thermal parameters temperature , T heat flux h, thermal conductivity k, mass density [rho] and thermal density c = h[rho]. With a harmonic excitation Q(x, t) the thermal response T(x, t) and the photothermal signal S(x, t) will be harmonic, too. The heat source Q(x, t) and, thus, the resulting temperature distribution T(x, t) can be separated into static and dynamic components. The dynamic part includes a (locally dependent) phase shift with respect to the excitation, described by a complex temperature distribution [T.sub.dyn]. This separation splits Eq. (1) into two equations, where only the dynamic part is of interest due to the lock-in detection of photothermal signals (locked on the modulation frequency):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

3. THERMOPHYSICAL PROPERTIES AND MICROSTRUCTURE

The thermal response of a volume element can be expressed by its thermal resistance R = [Delta]d/k . and its heat energy storage capacity C = c[Delta]d, assuming a volume element A[Delta]d . (Fig. 2) and heat flux perpendicular to the cross section A.

The local thermal impedance Z(x) = 1/[square root of i[omega]kc] introduces a network model of elementary cells, each described by an appropriate thermal RC filter as shown in Fig. 2. All the rules, used for describing electric networks, can be applied to thermal networks, too.

[FIGURE 2 OMITTED]

Z(x) is connected with the energy storage capability of the investigated material, while obstructions of the heat energy carriers' propagation (electrons and phonons) affect the diffusivity.

4. TEMPERATURE OSCILLATION OF INHOMOGENEOUS SAMPLES

[FIGURE 3 OMITTED]

If the excitation focus exceeds the detection spot and the diffusion length [[micro].sub.th] significantly, an approximately plane thermal wave can be assumed (Fig. 3). This requires a detection spot, which is smaller than [[micro].sub.th]. Thus Eq. (1) is reduced to a 1D problem. Besides, it requires nearly opaque samples, where the optical penetration depth is much smaller than the thermal diffusion length [[micro].sub.th]. A 1D description may not reflect the real heat flow on cracks or localized defects which have vertical boundaries.

4.1. Horizontal, vertical and random boundaries

Horizontal interfaces are considered as infinitely extended boundaries parallel to the sample surface, e.g. layered structures with the main heat flux perpendicular to the surface. The detection of vertical interfaces requires observation of the lateral heat flow. This can be made very sensitively by photothermal means providing an appropriate lateral offset between the heating and probing spots. Some materials are characterized by statistically distributed interfaces (porous materials, ceramics, alloys, minerals, mixtures, short-fiber reinforced plastics). Here, grain size related to the thermal diffusion length is decisive for the thermal wave propagation. This paper focuses mainly on horizontal interfaces.

Today, thermal waves are widely used for the estimation of the coating thickness and for the investigation of layered structures. For a 1D step-like thermal impedance depth profile, the reflection of thermal waves at the interface between materials A (coating) and B (substrate) can be described by the thermal reflection coefficient [r.sub.AB] = ([Z.sub.B] - [Z.sub.A])/ ([Z.sub.B] + [Z.sub.A]). For an ideal thermal contact, both temperature and heat flux are continuous across these interfaces; AB r is real and independent of the modulation frequency.

If the thermal contact differs from the ideal one, an additional thermal contact resistance has to be introduced, modifying [r.sub.AB] to a complex and frequency-dependent parameter [6,7]. None of these types of inhomogeneities will be pursued any further in this paper.

For thermal waves, propagating along the x axis, the local thermal impedance Z(x) is a complex-valued and frequency-dependent quantity. Z(x) describes the attenuation as well as the phase lag of thermal waves, penetrating the solid. Any depth profile of Z(x) determines the thermal response of a periodically heated sample. In contrast to the abrupt change in material properties, also continuous transition may occur. The thermal wave reflection of a "smeared" interface (decreasing slope in the transition zone) shows significant change of phase signals [6]. Estimation of thermal inhomogeneities with "smeared" horizontal interfaces requires photothermal depth profiling.

4.2. Evaluation of gradient depth profiles

4.2.1. Analytical evaluation

To interpret photothermally measured data, the oscillating surface temperature of the inspected sample has to be calculated. For simplification, a surface absorbing sample with a polygonal depth profile of thermal conductivity k is considered, whereas the thermal density c is assumed to be constant. Furthermore, the conditions of plane thermal wave propagation are fulfilled. The region from the sample surface x = 0 to the depth , [x.sub.N], limiting the domain of the depth-dependent , k is subdivided into N virtual sublayers (Fig. 4). The value of k(x) = [k.sub.m] + [[gamma].sub.m] (x - [x.sub.m]) varies piecewise linearly with the slope [[gamma].sub.m] (m = 0, ... , N - 1). In each sublayer [x.sub.m] [less than or greater than] x [less than or greater than] [x.sub.m+1], the temperature field can be composed of a left-handed and a right-handed propagating thermal wave with amplitudes [T.sup.+] and [T.sup.-] satisfying steady temperature and heat flow conditions at all interfaces x = [x.sub.m]. The complex-valued oscillating temperature can be derived from [9]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

A complete calculation has to include the gas (air) above the sample surface and the bulk material for . N x x > 0 K and 0 I are the modified Bessel functions of zero order with a complex argument. The thermal wave propagation in such a sample can be described by a matrix formalism, where the surface temperature oscillation can be expressed by the product of N component matrices. Each of these matrices transforms thermal waves from one sublayer to the next. For this case, the surface temperature can be expressed by an analytical formula, derived in [9]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

The transforming matrices [M.sub.m,m+1] are again functions of the modified Bessel functions of the first and zero order. Since optical absorptivity and infrared emissivity changes in the modified sample, the phase data are preferred to amplitude data for k profile reconstruction [10].

[FIGURE 4 OMITTED]

4.2.2. Inversion of 1-dimensional frequency scans

The measured frequency sweep of the photothermal signal [S.sub.meas] ([florin]) = |S([florin])| [e.sup.i[Phi]([florin])], which is proportional to [T.sub.surf]([florin]), has to be inverted to find out a polygonal depth profile, corresponding to the depth profile of the microstructural change. The inversion approach "builds-up" the profile of k iteratively sublayer by sublayer (Fig. 4), starting at the surface x = 0. For example, in sublayer No. m+1, by lowering the modulation frequency from [[florin].sub.m] to [[florin].sub.m+1], the thermal waves provide additional information about the sublayer between m x and [x.sub.m+1], since their penetration depth [[micro].sub.th] is inversely proportional to [[florin].sup.1/2]. In order to estimate the unknown [k.sub.m+1] = k([x.sub.m] + dx), the measured frequency sweep is compared with the calculated one S([florin], [k.sub.0]), using the least squares method. A reasonable estimation for the "true" value of [k.sub.0] can be found by searching the minimum of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

[N.sub.[florin]] is the number of measuring frequencies. The advantage of this approach consists in the evaluation of the unknown depth profiles by performing a sequence of one-parameter fits instead of using multiparameter methods. Besides, [k.sub.0] and [k.sub.m] are the only unknown parameters in each iteration step.

4.3. Finite element and finite difference methods

Numerical evaluation methods like the finite element method (FEM) and the finite difference method (FDM) divide the sample and the air above into small finite elements (Fig. 5).

The total heat flux Q and the boundary temperatures are described by a set of equations, each describing one cell. They balance for each element the variable temperature [T.sub.i], defined at nodes, and the heat fluxes [Q.sub.i] through neighbouring cells. Variational methods and superposition lead to the total temperature distribution within the whole sample.

Extending the whole grid to ten diffusion lengths, the size of the sample and the surrounding air can be assumed as "infinite", since the propagation of a thermal wave is limited to a few thermal diffusion lengths. The step size of the grid itself must be smaller than 0.1 [[micro].sub.th] in order to describe the thermal wave with sufficient resolution.

Conventional FEM programs like ANSYS are based on the general Eq. (1). This fact limits the performance of simulation (calculation time, resolution) significantly, since the periodic modulation has to be represented by small time steps. In addition, transient effects before reaching the equilibrium state must be calculated too, even though they are not of any interest. Thus a 3D simulation must be carried out as a 4D one, leading to high calculation efforts [11]. The latter disadvantages are avoided by introducing the FDM, causing additional programming work, since ready-to-use software packages are not available as yet.

[FIGURE 5 OMITTED]

In the 1D case, the thermal quantities T(x), k(x) and Q(x) only vary with the depth coordinate x (Fig. 3). The substitution

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

transforms Eq. (2) into an ordinary differential equation

u" + [p.sub.1](x)u' + [p.sub.0](x)u = -1/k Q(x). (8)

The solution of Eq. (8) is a discrete function u([x.sub.i] = [u.sub.i] within the interval [a,b]. The points i = 0,..., N, are equally distributed. The derivatives at [x.sub.i] can be approximated by difference quotients, using neighbouring temperatures [u.sub.i-1] and 1. i u + The indices i = 0 and i = N represent the boundary conditions, which have to be treated separately. Together with two boundary conditions, sufficient linear equations of N + 1 unknown values are available [11].

The set of equations can be written as GU = g The matrix G is tridiagonal, i.e. only 3N - 1 entries on the three main diagonals are different from zero requiring a special algorithm for tridiagonal matrices. The vector u gives the depth-dependent dynamic temperature distribution of Eq. (2) within the sample and in the air above.

The computing performance depends on the geometry, the thermal parameters and the frequency. For the examples, presented later, the simulation time for one frequency took less than one second (MATLAB realization) on a conventional PC.

If point-like defects and other local material variations occur or if the lateral heat flow cannot be neglected, 3D solution of Eq. (2) is necessary. Even influence of the finite excitation beam diameter may require a spatial description of the heat propagation. If, for example, an investigated hardness layer in the mm range is considered, the modulation frequency sweep has to start at 0.1 Hz. Then the diffusion length amounts to about 6 mm for a typical steel. For such a beam diameter, the heat flux density and the signal amplitude will be too small. Therefore, a 3D calculation has to be carried out in these cases.

By designing a non-regular three-dimensional grid, almost any kind of photothermal problems can be treated. A very promising option deals with the heat propagation in objects with rotational symmetry and with axial excitation (Fig. 5). They can be described in cylindrical coordinates, reducing the TDE to a 2D problem [11].

5. RESULTS

5.1. Test of the analytical inversion

The analytical inversion technique was at first tested by simulation, assuming a polygonal k profile of arbitrary shape with typical characteristics of a casehardened steel. The surface temperature according to Eq. (4) leads to the amplitude and phase signal. A certain amount of added noise simulated real measurement conditions. The original k profiles were retrieved from the noisy amplitude and phase sweeps with reasonable accuracy, applying the suggested inversion algorithm [9].

Laser-hardened steel plates, made of C60 with two different hardness traces, were investigated. Former investigations on case-hardened steel had proved a correlation between local hardness and thermal conductivity [12]. Applying the algorithm described above, the k profiles were reconstructed. The results are shown in Fig. 6a. The non-destructively evaluated k profiles are compared with the hardness profiles in Fig. 6b, destructively obtained by a series of Vickers micro-indentation measurements on the cross-section of the cut samples. The close intercorrelation between the microstructure (hardness) and local thermal conductivity is obvious. A quantitative comparison and, thus, a complete nondestructive hardness profiling is possible using specific calibration curves (Fig. 7).

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

5.2. Test of the FDM inversion

Figure 8 illustrates a simulation for different excitation beams of 1 mm, 3 mm and infinite diameters. The first two violate the condition for the 1D case. The phase maximum appears at the same frequency, whereas the maximum value itself and the phase contrast function is smeared for higher frequencies, depending on the beam diameter. The hardness depth, i.e. the depth of the k profile, corresponds to the frequency of the phase maximum. Its decay at higher frequencies includes information about the profile shape.

For the inversion problem, the FDM can also be used; the profile is described with several degrees of freedom, using estimated starting values for the first approximation step. The theoretical photothermal signals are calculated and compared with the measured ones in terms of a convenient objective function, e.g. Eq. (6). Thus the parameters to be identified can be approximated (minimum of the objective function). As an experimental verification, a hardness measurement was carried out on a shaft of 16MnCr5. The shaft was again inspected in both ways, non-destructively by photothermics and conventionally by a series of Vickers tests along a transverse cut. The Vickers test shows a smooth arctangent-like profile (width about 0.5 mm) at a mean depth of about 0.75 mm. Then, a photothermal measurement was carried out. Its phase contrast is illustrated in Fig. 9 by circles. FDM inversion (solid line) shows the mean hardening depth to 0.8 mm with a profile width of 0.4 mm.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

6. COMPARISON OF THE FDM WITH ANALYTICAL INVERSION

The thermal wave analysis, using the photothermal measuring technique, offers a non-destructive method to evaluate subsurface material properties. For example, the correlation between hardness and thermal conductivity enables the determination of a hardness profile. Both presented inversion algorithms meet industrial purposes, since they show reasonable results on case- and laser-hardened samples, even in the presence of typical measurement noise. The analytical approach causes less implementation efforts, retrieving the unknown profile by performing a sequence of one-parameter fits instead of using multiparameter methods. The algorithm promises a higher stability, since even in the case of occasional measuring and reconstruction errors in one sublayer, the uncertainty is reduced by including all the measured data up to the current probing depth. Additionally, no significant difference between both methods in terms of spatial resolution (depth or lateral) was observed [11]. In comparison with the analytical approach, the FDM inversion scheme shows decisive advantages, if the assumption of a 1D heat flow (plane thermal wave) is no longer valid. 2D and 3D heat propagation within the sample and the surrounding air can be modelled, enabling an easy variation of the geometry as well as of several thermal parameters (particularly their spatial distribution) at the same time. Also problems with partially transparent samples can be treated. Restrictions like ideal surface absorption or transparent air layers are no longer necessary. A determination of different parameters (e.g., profile, depths) and thermal quantities is possible for a given model. Also profiles of thermal density, even simultaneously with the conductivity, can be calculated using FDM. Calculation times are reduced by two orders compared with calculations, performed by a commercial FEM program.

ACKNOWLEDGEMENT

The authors wish to thank the German Research Foundation for funding parts of this research work within the Collaborative Research Centre SFB/TR4.

Received 27 July 2007

REFERENCES

[1.] Mandelis, A. (ed.). Nondestructive Evaluation. Progress in Photothermal and Photoacoustic Science and Technology, Vol. 2. Prentice Hall, Englewood Cliff, NJ, 1994.

[2.] Kruse, D., Prekel, H. and Goch, G. Automated photothermal detection of burning and hardening depth. In Proc. 9th International Conference on Infrared Sensors & Systems. Nurnberg, 2006, 341-346.

[3.] Goch, G., Schmitz, B. and Reick, M. Photothermal sensing techniques for measuring material properties and near-surface defects. Ann. CIRP, 1994, 42, 623-626.

[4.] Goch, G., Geerkens, J., Reick, M. and Schmitz, B. Analysis of surface layer variations by photothermal means. J. Physique IV, 1994, 319-322.

[5.] Fivez, J. and Thoen, J. Thermal waves in materials with linearly inhomogeneous thermal conductivity. J. Appl. Phys., 1994, 75, 7696-7699.

[6.] Friedrich, K., Seidel, U., Walther, H. G., Karpen, W. and Busse, G. Proposal for photothermal characterization of boundaries between layer and substrate. Res. Nondestr. Eval., 1993, 5, 31-39.

[7.] Ritter, R., Reick, M., Schmitz, B. and Goch, G. Nondestructive and contactless evaluation of surface coatings and adhesion defects by photothermal radiometry. Proc. SPIE, 1996, 2782, 662-673.

[8.] Seidel, U., Lan, T. T. N., Walther, H. G., Schmitz, B., Geerkens, J. and Goch, G. Quantitative characterization of material inhomogeneities by thermal waves. Opt. Eng., 1997, 36, 376- 390.

[9.] Lan, T. T. N., Seidel, U. and Walther, H. G. Theory of microstructure depth profiling by photothermal measurements. J. Appl. Phys., 1995, 77, 4739-4741.

[10.] Lan, T. T. N. and Walther, H. G. Photothermal depth profiling using only phase data. J. Appl. Phys., 1996, 80, 5289-5291.

[11.] Reigl, M. Methoden zur Quantifizierung photothermischer Signale. PhD Thesis, Universitat Ulm, 1997.

[12.] Lan, T. T. N., Seidel, U., Walther, H. G., Goch, G. and Schmitz, B. Experimental results of photothermal microstructural depth profiling. J. Appl. Phys., 1995, 78, 4108-4111.

Dennis Kruse, Helmut Prekel, Gert Goch and Heinz G. Walther

BIMAQ--Bremen Institute for Metrology, Automation and Quality Science, 28359 Bremen, Germany; kru@biba.uni-bremen.de

The photothermal detection enables characterization of materials in near-surface layers. New capabilities in machining processes and decreasing tolerances down to the nanometer scale have lead to material treatments, affecting only the layers near the surface. These near-surface areas influence the functional behaviour, quality and lifetime of the workpiece. Methods for the assessment of these areas should be non-destructive and meet standards concerning accuracy, measuring time, robustness and costs.

During the last years the photothermal techniques (Fig. 1), including various detector options, have gained a substantial progress.

[FIGURE 1 OMITTED]

Photothermal measuring techniques exploit the time-dependent heat flow, induced by a time-varying heat source, applied to the sample. This can be appropriately realized either by periodic laser heating or pulsed illumination. Both excitation schemes initiate propagating temperature oscillations, the so-called thermal waves. Their penetration into the inspected sample and different interactions with subsurface inhomogeneities open a vast field of non-destructive photothermal testing [1]. Today, photothermal measurements are applied to various materials and composites including coatings (thickness, adhesion strength, local defects), hardness profiles, residual stresses, wear, thermal influences and surface-affecting manufacturing techniques. Photothermal imaging includes quick, reliable and non-contact near-surface inspection with pattern sizes in the mm range as well as photothermal microscopy with [micro]m resolution [2].

2. PHYSICAL PRINCIPLES

2.1. Thermal waves and their detection

Fundamental photothermal and photoacoustic effects, a suitable experimental set-up and the two most important detection principles (photothermal radiometry and optical beam deflection) have been explained in detail in [3]. The corresponding basic theory can be traced back to the late 1970s ([1] and references therein), whereas intensive research, aiming at an application in production processes, started in the early 1990s [4-7].

In all photothermal measurements, an intensity-modulated light source (e.g. a laser beam) excites the investigated sample at its surface, as shown in Fig. 1. Optically opaque materials absorb the light directly at their surface, and, subsequently, the sample heating will also be modulated, causing a propagating thermal wave. The wave spreads within the sample, leading to a temperature rise by a few K at the surface. The penetration depth of the thermal wave is determined by the thermal diffusion length [[micro].sub.th] = [square root of (2[alpha]/[omega]), which is a function of the modulation frequency [omega] and the thermal diffusivity [alpha]. For metals, the penetration depth varies from a mm at a few Hz to some tens of [micro]m at a kHz-modulation [1]. Thermal inhomogeneities scatter the thermal wave and the surface temperature will differ from the value occurring in a homogeneous material. Photothermal measuring techniques map these (modulated and time-dependent) temperature changes at the surface due to hidden structures and near-surface changes of material properties [8].

2.2. The thermal diffusion equation

The heat propagation within a sample is given by the thermal diffusion equation [1,8,9]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1) ## which contains the locally-dependent thermal parameters temperature , T heat flux h, thermal conductivity k, mass density [rho] and thermal density c = h[rho]. With a harmonic excitation Q(x, t) the thermal response T(x, t) and the photothermal signal S(x, t) will be harmonic, too. The heat source Q(x, t) and, thus, the resulting temperature distribution T(x, t) can be separated into static and dynamic components. The dynamic part includes a (locally dependent) phase shift with respect to the excitation, described by a complex temperature distribution [T.sub.dyn]. This separation splits Eq. (1) into two equations, where only the dynamic part is of interest due to the lock-in detection of photothermal signals (locked on the modulation frequency):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

3. THERMOPHYSICAL PROPERTIES AND MICROSTRUCTURE

The thermal response of a volume element can be expressed by its thermal resistance R = [Delta]d/k . and its heat energy storage capacity C = c[Delta]d, assuming a volume element A[Delta]d . (Fig. 2) and heat flux perpendicular to the cross section A.

The local thermal impedance Z(x) = 1/[square root of i[omega]kc] introduces a network model of elementary cells, each described by an appropriate thermal RC filter as shown in Fig. 2. All the rules, used for describing electric networks, can be applied to thermal networks, too.

[FIGURE 2 OMITTED]

Z(x) is connected with the energy storage capability of the investigated material, while obstructions of the heat energy carriers' propagation (electrons and phonons) affect the diffusivity.

4. TEMPERATURE OSCILLATION OF INHOMOGENEOUS SAMPLES

[FIGURE 3 OMITTED]

If the excitation focus exceeds the detection spot and the diffusion length [[micro].sub.th] significantly, an approximately plane thermal wave can be assumed (Fig. 3). This requires a detection spot, which is smaller than [[micro].sub.th]. Thus Eq. (1) is reduced to a 1D problem. Besides, it requires nearly opaque samples, where the optical penetration depth is much smaller than the thermal diffusion length [[micro].sub.th]. A 1D description may not reflect the real heat flow on cracks or localized defects which have vertical boundaries.

4.1. Horizontal, vertical and random boundaries

Horizontal interfaces are considered as infinitely extended boundaries parallel to the sample surface, e.g. layered structures with the main heat flux perpendicular to the surface. The detection of vertical interfaces requires observation of the lateral heat flow. This can be made very sensitively by photothermal means providing an appropriate lateral offset between the heating and probing spots. Some materials are characterized by statistically distributed interfaces (porous materials, ceramics, alloys, minerals, mixtures, short-fiber reinforced plastics). Here, grain size related to the thermal diffusion length is decisive for the thermal wave propagation. This paper focuses mainly on horizontal interfaces.

Today, thermal waves are widely used for the estimation of the coating thickness and for the investigation of layered structures. For a 1D step-like thermal impedance depth profile, the reflection of thermal waves at the interface between materials A (coating) and B (substrate) can be described by the thermal reflection coefficient [r.sub.AB] = ([Z.sub.B] - [Z.sub.A])/ ([Z.sub.B] + [Z.sub.A]). For an ideal thermal contact, both temperature and heat flux are continuous across these interfaces; AB r is real and independent of the modulation frequency.

If the thermal contact differs from the ideal one, an additional thermal contact resistance has to be introduced, modifying [r.sub.AB] to a complex and frequency-dependent parameter [6,7]. None of these types of inhomogeneities will be pursued any further in this paper.

For thermal waves, propagating along the x axis, the local thermal impedance Z(x) is a complex-valued and frequency-dependent quantity. Z(x) describes the attenuation as well as the phase lag of thermal waves, penetrating the solid. Any depth profile of Z(x) determines the thermal response of a periodically heated sample. In contrast to the abrupt change in material properties, also continuous transition may occur. The thermal wave reflection of a "smeared" interface (decreasing slope in the transition zone) shows significant change of phase signals [6]. Estimation of thermal inhomogeneities with "smeared" horizontal interfaces requires photothermal depth profiling.

4.2. Evaluation of gradient depth profiles

4.2.1. Analytical evaluation

To interpret photothermally measured data, the oscillating surface temperature of the inspected sample has to be calculated. For simplification, a surface absorbing sample with a polygonal depth profile of thermal conductivity k is considered, whereas the thermal density c is assumed to be constant. Furthermore, the conditions of plane thermal wave propagation are fulfilled. The region from the sample surface x = 0 to the depth , [x.sub.N], limiting the domain of the depth-dependent , k is subdivided into N virtual sublayers (Fig. 4). The value of k(x) = [k.sub.m] + [[gamma].sub.m] (x - [x.sub.m]) varies piecewise linearly with the slope [[gamma].sub.m] (m = 0, ... , N - 1). In each sublayer [x.sub.m] [less than or greater than] x [less than or greater than] [x.sub.m+1], the temperature field can be composed of a left-handed and a right-handed propagating thermal wave with amplitudes [T.sup.+] and [T.sup.-] satisfying steady temperature and heat flow conditions at all interfaces x = [x.sub.m]. The complex-valued oscillating temperature can be derived from [9]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

A complete calculation has to include the gas (air) above the sample surface and the bulk material for . N x x > 0 K and 0 I are the modified Bessel functions of zero order with a complex argument. The thermal wave propagation in such a sample can be described by a matrix formalism, where the surface temperature oscillation can be expressed by the product of N component matrices. Each of these matrices transforms thermal waves from one sublayer to the next. For this case, the surface temperature can be expressed by an analytical formula, derived in [9]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

The transforming matrices [M.sub.m,m+1] are again functions of the modified Bessel functions of the first and zero order. Since optical absorptivity and infrared emissivity changes in the modified sample, the phase data are preferred to amplitude data for k profile reconstruction [10].

[FIGURE 4 OMITTED]

4.2.2. Inversion of 1-dimensional frequency scans

The measured frequency sweep of the photothermal signal [S.sub.meas] ([florin]) = |S([florin])| [e.sup.i[Phi]([florin])], which is proportional to [T.sub.surf]([florin]), has to be inverted to find out a polygonal depth profile, corresponding to the depth profile of the microstructural change. The inversion approach "builds-up" the profile of k iteratively sublayer by sublayer (Fig. 4), starting at the surface x = 0. For example, in sublayer No. m+1, by lowering the modulation frequency from [[florin].sub.m] to [[florin].sub.m+1], the thermal waves provide additional information about the sublayer between m x and [x.sub.m+1], since their penetration depth [[micro].sub.th] is inversely proportional to [[florin].sup.1/2]. In order to estimate the unknown [k.sub.m+1] = k([x.sub.m] + dx), the measured frequency sweep is compared with the calculated one S([florin], [k.sub.0]), using the least squares method. A reasonable estimation for the "true" value of [k.sub.0] can be found by searching the minimum of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

[N.sub.[florin]] is the number of measuring frequencies. The advantage of this approach consists in the evaluation of the unknown depth profiles by performing a sequence of one-parameter fits instead of using multiparameter methods. Besides, [k.sub.0] and [k.sub.m] are the only unknown parameters in each iteration step.

4.3. Finite element and finite difference methods

Numerical evaluation methods like the finite element method (FEM) and the finite difference method (FDM) divide the sample and the air above into small finite elements (Fig. 5).

The total heat flux Q and the boundary temperatures are described by a set of equations, each describing one cell. They balance for each element the variable temperature [T.sub.i], defined at nodes, and the heat fluxes [Q.sub.i] through neighbouring cells. Variational methods and superposition lead to the total temperature distribution within the whole sample.

Extending the whole grid to ten diffusion lengths, the size of the sample and the surrounding air can be assumed as "infinite", since the propagation of a thermal wave is limited to a few thermal diffusion lengths. The step size of the grid itself must be smaller than 0.1 [[micro].sub.th] in order to describe the thermal wave with sufficient resolution.

Conventional FEM programs like ANSYS are based on the general Eq. (1). This fact limits the performance of simulation (calculation time, resolution) significantly, since the periodic modulation has to be represented by small time steps. In addition, transient effects before reaching the equilibrium state must be calculated too, even though they are not of any interest. Thus a 3D simulation must be carried out as a 4D one, leading to high calculation efforts [11]. The latter disadvantages are avoided by introducing the FDM, causing additional programming work, since ready-to-use software packages are not available as yet.

[FIGURE 5 OMITTED]

In the 1D case, the thermal quantities T(x), k(x) and Q(x) only vary with the depth coordinate x (Fig. 3). The substitution

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

transforms Eq. (2) into an ordinary differential equation

u" + [p.sub.1](x)u' + [p.sub.0](x)u = -1/k Q(x). (8)

The solution of Eq. (8) is a discrete function u([x.sub.i] = [u.sub.i] within the interval [a,b]. The points i = 0,..., N, are equally distributed. The derivatives at [x.sub.i] can be approximated by difference quotients, using neighbouring temperatures [u.sub.i-1] and 1. i u + The indices i = 0 and i = N represent the boundary conditions, which have to be treated separately. Together with two boundary conditions, sufficient linear equations of N + 1 unknown values are available [11].

The set of equations can be written as GU = g The matrix G is tridiagonal, i.e. only 3N - 1 entries on the three main diagonals are different from zero requiring a special algorithm for tridiagonal matrices. The vector u gives the depth-dependent dynamic temperature distribution of Eq. (2) within the sample and in the air above.

The computing performance depends on the geometry, the thermal parameters and the frequency. For the examples, presented later, the simulation time for one frequency took less than one second (MATLAB realization) on a conventional PC.

If point-like defects and other local material variations occur or if the lateral heat flow cannot be neglected, 3D solution of Eq. (2) is necessary. Even influence of the finite excitation beam diameter may require a spatial description of the heat propagation. If, for example, an investigated hardness layer in the mm range is considered, the modulation frequency sweep has to start at 0.1 Hz. Then the diffusion length amounts to about 6 mm for a typical steel. For such a beam diameter, the heat flux density and the signal amplitude will be too small. Therefore, a 3D calculation has to be carried out in these cases.

By designing a non-regular three-dimensional grid, almost any kind of photothermal problems can be treated. A very promising option deals with the heat propagation in objects with rotational symmetry and with axial excitation (Fig. 5). They can be described in cylindrical coordinates, reducing the TDE to a 2D problem [11].

5. RESULTS

5.1. Test of the analytical inversion

The analytical inversion technique was at first tested by simulation, assuming a polygonal k profile of arbitrary shape with typical characteristics of a casehardened steel. The surface temperature according to Eq. (4) leads to the amplitude and phase signal. A certain amount of added noise simulated real measurement conditions. The original k profiles were retrieved from the noisy amplitude and phase sweeps with reasonable accuracy, applying the suggested inversion algorithm [9].

Laser-hardened steel plates, made of C60 with two different hardness traces, were investigated. Former investigations on case-hardened steel had proved a correlation between local hardness and thermal conductivity [12]. Applying the algorithm described above, the k profiles were reconstructed. The results are shown in Fig. 6a. The non-destructively evaluated k profiles are compared with the hardness profiles in Fig. 6b, destructively obtained by a series of Vickers micro-indentation measurements on the cross-section of the cut samples. The close intercorrelation between the microstructure (hardness) and local thermal conductivity is obvious. A quantitative comparison and, thus, a complete nondestructive hardness profiling is possible using specific calibration curves (Fig. 7).

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

5.2. Test of the FDM inversion

Figure 8 illustrates a simulation for different excitation beams of 1 mm, 3 mm and infinite diameters. The first two violate the condition for the 1D case. The phase maximum appears at the same frequency, whereas the maximum value itself and the phase contrast function is smeared for higher frequencies, depending on the beam diameter. The hardness depth, i.e. the depth of the k profile, corresponds to the frequency of the phase maximum. Its decay at higher frequencies includes information about the profile shape.

For the inversion problem, the FDM can also be used; the profile is described with several degrees of freedom, using estimated starting values for the first approximation step. The theoretical photothermal signals are calculated and compared with the measured ones in terms of a convenient objective function, e.g. Eq. (6). Thus the parameters to be identified can be approximated (minimum of the objective function). As an experimental verification, a hardness measurement was carried out on a shaft of 16MnCr5. The shaft was again inspected in both ways, non-destructively by photothermics and conventionally by a series of Vickers tests along a transverse cut. The Vickers test shows a smooth arctangent-like profile (width about 0.5 mm) at a mean depth of about 0.75 mm. Then, a photothermal measurement was carried out. Its phase contrast is illustrated in Fig. 9 by circles. FDM inversion (solid line) shows the mean hardening depth to 0.8 mm with a profile width of 0.4 mm.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

6. COMPARISON OF THE FDM WITH ANALYTICAL INVERSION

The thermal wave analysis, using the photothermal measuring technique, offers a non-destructive method to evaluate subsurface material properties. For example, the correlation between hardness and thermal conductivity enables the determination of a hardness profile. Both presented inversion algorithms meet industrial purposes, since they show reasonable results on case- and laser-hardened samples, even in the presence of typical measurement noise. The analytical approach causes less implementation efforts, retrieving the unknown profile by performing a sequence of one-parameter fits instead of using multiparameter methods. The algorithm promises a higher stability, since even in the case of occasional measuring and reconstruction errors in one sublayer, the uncertainty is reduced by including all the measured data up to the current probing depth. Additionally, no significant difference between both methods in terms of spatial resolution (depth or lateral) was observed [11]. In comparison with the analytical approach, the FDM inversion scheme shows decisive advantages, if the assumption of a 1D heat flow (plane thermal wave) is no longer valid. 2D and 3D heat propagation within the sample and the surrounding air can be modelled, enabling an easy variation of the geometry as well as of several thermal parameters (particularly their spatial distribution) at the same time. Also problems with partially transparent samples can be treated. Restrictions like ideal surface absorption or transparent air layers are no longer necessary. A determination of different parameters (e.g., profile, depths) and thermal quantities is possible for a given model. Also profiles of thermal density, even simultaneously with the conductivity, can be calculated using FDM. Calculation times are reduced by two orders compared with calculations, performed by a commercial FEM program.

ACKNOWLEDGEMENT

The authors wish to thank the German Research Foundation for funding parts of this research work within the Collaborative Research Centre SFB/TR4.

Received 27 July 2007

REFERENCES

[1.] Mandelis, A. (ed.). Nondestructive Evaluation. Progress in Photothermal and Photoacoustic Science and Technology, Vol. 2. Prentice Hall, Englewood Cliff, NJ, 1994.

[2.] Kruse, D., Prekel, H. and Goch, G. Automated photothermal detection of burning and hardening depth. In Proc. 9th International Conference on Infrared Sensors & Systems. Nurnberg, 2006, 341-346.

[3.] Goch, G., Schmitz, B. and Reick, M. Photothermal sensing techniques for measuring material properties and near-surface defects. Ann. CIRP, 1994, 42, 623-626.

[4.] Goch, G., Geerkens, J., Reick, M. and Schmitz, B. Analysis of surface layer variations by photothermal means. J. Physique IV, 1994, 319-322.

[5.] Fivez, J. and Thoen, J. Thermal waves in materials with linearly inhomogeneous thermal conductivity. J. Appl. Phys., 1994, 75, 7696-7699.

[6.] Friedrich, K., Seidel, U., Walther, H. G., Karpen, W. and Busse, G. Proposal for photothermal characterization of boundaries between layer and substrate. Res. Nondestr. Eval., 1993, 5, 31-39.

[7.] Ritter, R., Reick, M., Schmitz, B. and Goch, G. Nondestructive and contactless evaluation of surface coatings and adhesion defects by photothermal radiometry. Proc. SPIE, 1996, 2782, 662-673.

[8.] Seidel, U., Lan, T. T. N., Walther, H. G., Schmitz, B., Geerkens, J. and Goch, G. Quantitative characterization of material inhomogeneities by thermal waves. Opt. Eng., 1997, 36, 376- 390.

[9.] Lan, T. T. N., Seidel, U. and Walther, H. G. Theory of microstructure depth profiling by photothermal measurements. J. Appl. Phys., 1995, 77, 4739-4741.

[10.] Lan, T. T. N. and Walther, H. G. Photothermal depth profiling using only phase data. J. Appl. Phys., 1996, 80, 5289-5291.

[11.] Reigl, M. Methoden zur Quantifizierung photothermischer Signale. PhD Thesis, Universitat Ulm, 1997.

[12.] Lan, T. T. N., Seidel, U., Walther, H. G., Goch, G. and Schmitz, B. Experimental results of photothermal microstructural depth profiling. J. Appl. Phys., 1995, 78, 4108-4111.

Dennis Kruse, Helmut Prekel, Gert Goch and Heinz G. Walther

BIMAQ--Bremen Institute for Metrology, Automation and Quality Science, 28359 Bremen, Germany; kru@biba.uni-bremen.de