Development and validation of an inverse method for identification of thermal characteristics of a short laser pulse

The paper presents development and validation of an inverse method for the identification of thermal characteristics of a short single laser pulse which stroke in a metal sample. The inverse method was applied to find unknown power of the laser pulse, the dimensionless shape parameter of the super-Gaussian function which describes the beam spatial profile as well as beginning and end times of the exposition of the metal sample to the laser pulse. The proposed inverse algorithm was based on the Levenberg-Marquardt technique as well as on temporal and spatial distributions of temperature on the rear surface of the sample, i.e., the opposite to the irradiated one, measured using the experimental stand. The performed investigations showed that the problem was ill-posed but good accuracy was obtained. The low sensitivity of registered temperature to changes in both power and duration of the pulse affected the retrieving accuracy most significantly. Moreover, the dependence of solution of the inverse problem on the initial guess was observed. The accuracy was also affected by low temporal resolution (500 Hz, with the exposure time from 0.2 to 1 ms) of the IR camera. This resolution affected the temporal sampling of measured temperatures. Despite these problems, the method was able to retrieve unknown pulse parameters with 20-25% accuracy.


Introduction
The high-power lasers (1-100 kW) are applied in many engineering fields, e.g., in material processing for cutting, drilling, heat treating, etching, machining and welding [1], in electro-optical systems for optical pumps for solid state lasers, fibre amplifiers, LIDAR, target illuminators and designators [2] as well as in optical data storage and in the free-space optical communication [2]. They are also utilized in the medicine [3] (e.g., in dermatology, ophthalmology, stomatology and surgery) and cosmetics [3]. The defence industry is another area of wide applications of the high-power lasers [4]. They may also serve as heat sources in materials characterization [5]. The interactions of high-power beams with matter are very complex and depend on the laser power and duration of the exposition [6]. The process of absorption of the laser radiation by matter involves many phenomena which are difficult to describe and model [7]. Heating, melting, evaporation, formation of clouds of the evaporated material, shock waves generation, ionization and formation of plasma are examples of these phenomena [7,8].
In various high-power laser beam applications the beam intensity distribution and spot size are formed in an optic system. These two parameters as well as the laser power are very important because they determine how the laser beam will interact with the target. However, due to imperfections, manufacturing tolerances and other effects like input laser distortions, optical distortion, heating, overall instabilities and nonlinear effects the beam profile might deviate from the assumed one [9]. These processes affect the beam profile changing it from the desired Gaussian one to super-Gaussian or even more complex form. Therefore, monitoring, diagnostic and optimization of the laser beam power, temporal and spatial distribution of beam intensity as well as spot size are very important for repeatable operation of many devices. Among the facilities for beam profile analysis are the scanning aperture profilers which use slits, knife-edges, or pinholes with single large area detectors, or camerabased profilers using CCD's or photodiode arrays. Most of the available measuring devices cannot be applied for direct analysis of the high-power lasers in industrial environment and their response time and mounting is inadequate for fast or on-line measurements [9]. For example, the high sensitivity of the camera profilers requires the laser power to be reduced by many orders of magnitude using beam sampling or optical attenuation [10]. Recently a new profiling technique which based on the Rayleigh scatter from the beam have been proposed [10]. This method overcomes the power obstacle and allows dynamic measurement and diagnostic of the beam parameters as well as determination of M 2 parameters of laser beams with power from 1kW to 100 kW. However, most of instrumentation for the beam monitoring is very expensive, fragile and requires great care in its usage. Recently, Kujawińska et al. [11] proposed a new approach to measure laser beam characteristics which is based on the maps of temperature and displacements acquired at the rear surface of the heated metal sample with the aid of the high-speed infrared thermography and the 3D Digital Image Correlation (DIC) method [12], respectively. The thermographic method registers the temporal and spatial variations of temperature while the 3D DIC method allows to track sample deformation resulting from thermal stresses induced by significant sample heating by the laser pulse. Moreover, the proposed approach utilizes inverse thermal analysis to find unknown laser beam parameters [13] in which the surface heat flux resulting from the laser interaction with the sample is retrieved. The developed method is very simple and therefore may be easily applied in the industrial or field conditions with a relatively low cost compared to the other methods.  In this paper, development and validation of an inverse thermal method for the identification of thermal characteristics of a short single super-Gaussian laser pulse which incident on a metal sample is presented. The method was applied to find four unknown laser pulse parameters, i.e., its power, the dimensionless shape parameter as well as beginning and end times of the exposition. The inverse algorithm based on the Levenberg-Marquardt (LM) technique [14] which used the temporal and spatial distributions of temperature on the rear surface of the sample, i.e., the opposite to the irradiated one. The temperature profiles were recorded by a high-speed IR camera. The developed method is a first step in the implementation of the new and more general approach to measure laser beam characteristics proposed in [11]. The second step will be related to the incorporation into the recovering procedure the maps of displacements registered by the 3D DIC method.

Method for the identification of the short laser pulse parameters
In the developed method, temporal and spatial distributions of temperature on the rear side (opposite to irradiated one) of the metal plate heated by the laser beam are utilized [11] in the identification of the laser pulse parameters. The schematic of the method is presented in Figure 1. A short pulse is emitted by the laser diode and hits the front surface of a metal plate. The thermal response of the sample at the rear surface serves as the input for the inverse algorithm which allows to retrieve the beam parameters.
In order to increase the thermal response on the rear side of the plate the thermal resistivity of the sample along its thickness should be as low as possible [11]. The temperature in the region of the sample irradiated by the laser beam might increase substantially, therefore high thermal stresses might also appear. Considering both thermal and structural aspects of the problem, 1 mm thick aluminum sheet metal (AW2017A T4) was selected among others available on the market [11]. Moreover, the thickness of the sample in the irradiated region was reduced to 0.5 mm by circular milling as shown in Figure 2 to lower its thermal resistivity.

Model of the forward problem
The non-destructive character of the proposed identification method was assumed. The duration of the laser beam was very short and expected maximal temperatures in the sample were below the melting point of aluminium (Tm = 883.0 K)see [11]. Therefore, only heat conduction in the plate was considered, as described by the following equation: with boundary conditions given by following formulas (see Figure 2): -for the circular hollow: -for side walls: -for other walls: and with the following initial condition: where: c is the specific heat, hheat transfer coefficient, n normal vector, rsurface reflectivity, ttime, Ttemperature, wwall, thermal conductivity, density,  -Stefan-Boltzmann constant and surroundings. The profile of the heat flux, generated at the specimen surface, was assumed to be rectangular and super-Gaussian with respect to time and space, respectively, and was given by the following relationship: with the following normalization constant: where: p is the dimensionless shape coefficient of the super-Gauss function, Plaser power, Rdistance from the beam centre to a given location, R0distance over which the Gaussian profile (p = 2) decreases to e -2 of its axial value, ts and teinitial and end time of the laser pulse, respectively, and  -Euler gamma function.
The developed model was implemented applying the ANSYS Workbench environment. The geometry and mesh were generated in the ANSYS DesignModeler and ANSYS Meshing, respectively. The grid consisted of almost 198 thousand of elements. The mesh skewness and orthogonal quality were below 0.7 and above 0.3, respectively. These parameters were within acceptable range. The forward problem was solved applying the control volume method [15] based solver ANSYS Fluent and its customization interface, i.e., user-defined functions.

Experimental measurements
The experimental stand is ready for more complex measurements which include monitoring of both temperature and 3D displacements at the investigated samples (see Figure 3). In tests the sample was illuminated by single pulse or train of 5 pulses of the Nd:YAG laser of the following parameters: wavelength of 1.06 µm, frequency of 10 Hz and impulse duration in the range of 0.2-1 ms. The energy of a single pulse varied from 2 up to 5 J and was focused on the sample by a lens into a circular area of diameter equal to 12 mm. The effects of the interaction between the laser beam and the solid sample have been observed in the form of both the temperature and displacement maps provided by the IR camera and 3D DIC [12] systems, respectively, captured at the rear surface of the sample. The setup ( Figure 3A) consisted of one infrared camera (FLIR SC 7500, spectral range of 1.5-5.1 µm, temperature range of 5-300C, accuracy of 1C) used for monitoring of the temperature and two fast CMOS cameras (Optronis CR30002) equipped with 100 mm focal length lenses used for the 3D DIC measurements. The sampling rate of the IR camera was 500 Hz. To ensure sufficient amount of illumination for the 3D DIC measurements and to avoid additional heating of the sample, two 100 W LED reflectors had been used for illumination of the field of view (see Figure 3B). The IR camera was mounted at a heavy free-standing tripod and directed at a slight angle towards the illumination axis to avoid reflections from the sample ( Figure 3B). The samples have been mounted in a sample holder made of constructional aluminum profiles. The holder fixed shorter edges of the sample while the longer ones were left unconstrained. It has been bolted to the optical table to ensure its stability.

Inverse algorithm
In the presented version of the inverse algorithm single laser pulses were investigated. Their parameters were retrieved based on the temporal and spatial variations of temperature on the rear surface of the irradiated sample (thermal problem only, see Figure 4). These distributions were acquired applying the high-speed IR camera. The unknown pulse parameters were: the power P, dimensionless shape coefficient p as well as initial and end times of the pulse ts and te, respectively. These unknown parameters were represented by the following vector: q = [q1, q2, q3, q4], where q1 refers to the laser power P, q2 to the dimensionless shape coefficient p, q3 to the initial time of the pulse ts and q4 to the end time te. The inverse problem referred to minimization of the following last squares norm [14]: T  S r r (8) where the residual vector r was defined as a difference between temperatures obtained from the forward problem model (T m ) and from the experimental measurements (T e ), i.e.,   me  r T q T (9) The norm S given by eq. (8) depended on the unknown vector q. Therefore, the task of the inverse algorithm was to find vector q which minimized this norm. The LM technique [14] was applied to solve the defined inverse problem. This method uses the adaptive damping factor for regularization. Preliminary, the calibration of the method was necessary to find the initial value of this factor and additionally to select the stopping criteria. These parameters are case-dependent and were found by extensive testing of the developed inversed algorithm. The assumed values of the damping factor and stopping criteria are discussed later in this paper. Moreover, before the inverse method was applied to find unknown parameters of the laser pulse a thorough sensitivity analysis of the proposed method was performed and described by Pietrak et al. [13]. The investigations carried out in [13] answered the question about both the spatial locations on the sample surface and time moments in which the measured quantities (temperatures) were most sensitive to the changes in the unknown laser beam parameters.
The conclusion drawn from the investigation carried out in [13] was that the temperature acquisition should last long enough to include the cool down period at the rear surface of the sample. This was related to the shifting of the maximum temperature at the rear surface of the sample with respect to the maximum temperature at the front surface. Therefore, in the developed inverse algorithm the temperature history in the point C at the middle of the rear surface of the sample (Figure 4) was applied. These data were denoted by vector Tt. Moreover, the analysis carried out in [13] showed that the inverse procedure should also use the spatial temperature distribution data along the lines Lx and Ly shown in Figure 4 (denoted by vectors Tx and Ty, respectively) from at least one time instant (frame). These data are required to allow proper identification of the p-parameter, which describes the flatness of the super-Gaussian profile. The utilized frame should be collected just before the laser pulse ended when the temperature at the rear surface of the sample approaches maximum value. If necessary, additional frames might be also collected with a delay of a few milliseconds from the first frame.
After thorough testing of the developed inverse algorithm applying experimental data, a single frame collected close to the moment of maximum temperature at the rear surface of the sample was used in the procedure as suggested by Pietrak et al. [13]. Moreover, in order to improve convergence of the power and temporal parameters (i.e., P, ts and te) the information about time variation of temperature at the point C

T T T T T T T T T
The temperature history data from point C (denoted as Tt) were repeated several times in the residual vector in order to remove disparity between temporal and spatial resolution of the experimental data applied in the inverse algorithm. The former resolution was significantly lower than the latter one, i.e., 1 ms vis 2020 pixel, where 1 pixel corresponded to 0.41 mm, while the laser pulse duration varied from 0.2 to 1 ms. As mentioned earlier, the spatial temperature distribution data along lines Lx and Ly at the rear surface of the sample (vectors Tx and Ty) were taken from the thermogram map for a single time frame. The instant of the collection was close to the time in which maximal temperature at the rear surfaces in point C was attained, i.e., approximately 5 ms after the beginning of the laser impulse (the exact delays are not known because the sampling in the camera was not synchronized with the laser triggering).
In the adopted LM technique [14], the unknown vector q has been found iteratively according to the following formula: where superscript j denotes j-th iteration and the symbol J is the sensitivity (Jacobian) matrix, adaptive regularization coefficient and diagonal amplification matrix. The iterative loop started from the supplied initial guess q 0 . The elements of the Jacobian matrix were approximated by the following finite differences formula: (13) where the subscript i is the index of the element of temperature vectors T e and T m , k14index of the unknown parameter, qkperturbation of k-th unknown parameter. In practice, single computation of the Jacobian matrix required five runs of the numerical 4 MATEC Web of Conferences 240, 01021 (2018) https://doi.org/10.1051/matecconf/201824001021 ICCHMT 2018 model of direct problem, i.e., one for the current values of unknown parameters and four times with perturbed value of each unknown parameter separately. The perturbation vector q was taken as 1/10 of the initial guess vector q 0 . The matrix  was assumed after Marquardt [16] in the following form: The starting value of  was estimated equal to  0 = 50. Then in each iteration of the LM algorithm this coefficient was modified in the following way. If the norm S j+1 for the current iteration was smaller or equal to the norm S j for the previous iteration the regularization coefficient was decreased by the factor of 0.3. Otherwise  was increased by 0.3 -1 .
The following three stopping criteria were applied in the developed inverse algorithm: where the symbol . denotes Euclidean norm. The values of constants 1, 2 and 3 were set equal to 1 = 5, 2 = 5 and 3 = 0.002, determined for power in kW, temporal parameters in ms and least squares norm S defined by eq. (8). After the one of these stopping criteria was attained the unknown parameters were assumed found. 5. Implementation of the inverse procedure (q f denotes vector of estimated parameters, wvector of known system parameters).

Fig.
The inverse algorithm was developed applying the GNU Octave and ANSYS Fluent software. The former was the master program and the latter was the slave, invoked with a journal file modified beforehand with the use of an in-house GNU Octave routine. The GNU Octave allowed for external solution of the forward problem (in the ANSYS Fluent), modification of input parameters to the thermal model (modification of input files which were then read by UDF macros in the ANSYS Fluent), utilization of the output parameters from the thermal model (reading of output files from ANSYS Fluent which were written by UDF macros), development of the inverse problem solution method by applying build-in libraries, iterative execution of the overall inverse algorithm and analysis of the obtained results. The general schematic of the solution methodology is presented in Figure 5.

Retrieving of the laser beam parameters
The aluminium sample used in the experiments had the following geometrical dimensions: 80.050.01.0 mm (widthlengththickness of the plate, respectively) and 2.00.5 mm (diameterthickness of the milling, respectively). Its thermophysical properties (aluminium alloy AW2017A T4) are given in Table 1 [11]. Four cases, e.g., two single and two multiple laser pulses were used in the investigations presented in this paper. During the experiments some of the laser beam parameters, i.e., the laser power, the beam spatial profile as well as beginning and end time of the exposition were known. Therefore, the measured data were used for the validation of both the forward problem and proposed inverse algorithm. The initial and boundary conditions as well as laser pulse parameters for considered cases are given in Table 2. The value of heat transfer coefficient at sample surfaces was assumed for natural convection in the case of horizontally oriented plates [17], while the value of surface reflectivity was taken for polished aluminium surface [17]. The initial and surroundings temperatures were measured during the experiments. Before application of the forward problem model in the inverse procedure its correctness and accuracy was proven. In the first step in the verification procedure convergence studies for time step and grid sizes were carried out. The forward problem was solved applying three different time step sizes, i.e., the base one (t = 10 -5 s), two times smaller and two times larger one as well as two additional meshes, i.e., the coarse and finer one with 10% less and 10% more elements, respectively. The spatial and temporal temperature distributions obtained for three time steps and three meshes were almost identical with the relative differences below 1%. Then in the second step, the validation using the experimental measurements conducted on the stand described above was carried out. The numerical and experimental results were found to agree very well. Detailed analysis showed that absolute differences between predicted and measured data were below 0.4 K and 0.1 K for the peak and cooling period, respectively. The comparison proved that the forward model gave physically-accurate results which matched the measurements for single and multiple pulses.
The results of the retrieving of four unknown laser beam parameters using experimental data are presented in Table 3 for cases defined in Table 2. For cases which involved multiple pulses (case 3 and 4), the data for the first pulse only were examined with the presented inverse algorithm. In Table 3 the obtained results are compared with the real laser pulse parameters which were known during experiments. The temporal variation of temperature in point C and spatial distribution of temperature along lines Lx and Ly for initial guess and final estimation as well as measured in the experiments for cases 2 and 4 are shown in Figure 6 and 7.
The tests carried out showed that the considered problem is ill-posed. However, relatively good retrieving accuracy of the unknown parameters of the laser pulse have been obtainedsee Table 3. The main reason of illconditioning is very low sensitivity of the measured temperature to the variation of the end time of the pulse. Inaccurate retrieval of the end time of the exposition to the laser pulse results in decrease in the accuracy of the laser power fitting. Decrease or increase of the laser pulse duration was related to decrease or increase of the amount of energy absorbed by the plate. The same effect was obtained by decreasing or increasing of the laser power. Therefore, the end time of the exposition and the laser power were mutually dependent parameters and there was no unique solution of the problem.
To correctly retrieve the end time of the pulse (pulse duration was on the order of 10 1 ms) as high as possibly temporal resolution of temperature registration was required (1 ms resolution was available in the measurements). The available temporal resolution the IR camera (much higher than duration of pulse) limited the accuracy of the exposition end time fitting. The accuracy of the retrieved parameters was found dependent on the initial guess. Good matches were obtained when the inverse algorithm started from the low laser power and short duration of the laser pulse. Moreover, the initial starting time of the exposition to the laser pulse should be closed to its real value. It should be noted that in practical applications the pulse duration is typically known and the main objective is to verify the spatial distribution of pulse energy. Moreover, for the specimen with a higher value of the thermal resistivity along its thickness (e.g., for thicker plate or sample made of material with the low value of thermal diffusivity) the conditioning of the inverse problem under consideration will be decreased. Such sample will have higher ability to damp thermal response on the rear surface of the sample. Fig. 6. A) Temporal temperature variation in point C as well as B) and C) distribution of temperature along lines Lx and Ly, respectively for approximately 5 ms after heating started for case 2 and for initial guess, final estimation and experimental measured. Despite the aforementioned problems and limitations, the accuracy of the proposed retrieving method was on the level of 20-25%. The limitations of the quality of the experimental data were related to damping of the temperature response on the rear surface of the sample (effect of material type and thickness of the sample) and relatively low resolution of the IR camera (the frame rate of 500 Hz at 2020 pixels, while the exposure times were in the range of 0.2-1 ms). The method recognized shorter (case 1, 3 and 4) and longer impulse (case 2). Moreover, the spatial distribution of the laser pulse matched well to the real datasee Figure 6 and 7. The under-or over-prediction of the pulse duration resulted in the over-or under-prediction of the laser power, respectively. However, the amount of absorbed energy was predicted very well.
The proposed inverse procedure has also disadvantages related to its iterative character and required numerical computation of the sensitivity matrix at each iteration of the solution procedure. The Jacobian is non-linear and depends on the unknown parameters and should be updated in each iteration. Moreover, the determination of the Jacobian matrix applying finite difference method required N+1 model evaluations, where N was the number of unknowns. Therefore, computations are quite expensive especially in the case of large grids. Fortunately, this problem may be solved by applying parallel computation of the forward problem.

Conclusions
The paper presents a method of identification of thermal characteristics of the short single super-Gaussian laser pulse. The proposed method is based on the thermal imaging technique and utilizes the LM algorithm for retrieval of four unknown parameters of the laser pulse, i.e., the laser power, dimensionless shape parameter of the super-Gaussian function which describes the beam spatial profile as well as beginning and end times of the exposition of the sample to the laser pulse. The method involves thin metal plate illuminated by the investigated laser source as well as temporal and spatial distributions of temperature on the rear surface of the sample, i.e., the opposite to the irradiated one. The solution of the inverse problem was carried out applying the commercial software GNU Octave and ANSYS Fluent. Even though the data from the 3D DIC method were not yet considered in the retrieval procedure, it was shown that accurate determination of four considered parameters is possible, with 20-25% of accuracy.
The tests carried out showed that the considered inverse problem is ill-posed. The main reason for this behavior was a very low sensitivity of the measured temperature to the variation of the end time of the exposition to the laser pulse. Moreover, this parameter was found to influence the measurement data in the same way as the laser power. Therefore, these two parameters were mutually dependent and contributed to the solution non-uniqueness. The conditions under which the influence of ill-conditioning may be minimized were identified and described in this paper. The first one was related to required high temporal resolution of temperature registration. The second one referred to the proper form of the input data (vectors T e and T m ). The third one concerned the proper selection of the initial guess. The last one was related to the thermal properties and thermal resistivity of the aluminum sample. As shown in this paper, consideration of these conditions allows to obtain satisfying results. The method correctly recognized differences between the shorter and longer impulse duration. Moreover, the predicted spatial distribution of the laser pulse was found to correspond well to the experimental data. The under-or overprediction of the pulse duration resulted in the over-or under-prediction of the laser power, respectively. However, the amount of absorbed energy was predicted correctly.