A V S Oliveira

J Teixeira

V Schick

D Maréchal

M Gradeck

S Denis

Using a linear inverse heat conduction model to estimate the boundary heat flux with a material undergoing phase transformation

Keywords: Inverse method, Cooling, Steel, Heat transfer coefficient, Quenching, Latent heat Nomenclature

Modeling heat conduction during steel-cooling processes is challenging because solid-solid phase transformations can take place, introducing highly non-linear phenomena to the heat equation. Using a linear model for the inverse problem can be helpful for estimating the dissipated heat flux, especially because complex phase transformations models that use empirical parameters are not always available. However, the performance and applicability of such a simplified inverse model for a strongly non-linear process is still unclear. This study presents several numerical simulations to evaluate the accuracy and limitations of solving a linear inverse heat conduction problem (i.e. constant thermophysical properties and no internal heat source) to estimate the boundary heat flux during cooling of a material undergoing phase transformations. Preliminary simulations with stable materials but with temperature-dependent thermophysical properties showed that the linear model performs well to estimate the boundary heat flux, except when the material has a highly temperature-dependent specific heat, like pure iron near its Curie temperature. Then, we performed several simulations for 42CrMo4 steel, which undergoes phase transformations and, hence, has phaseand temperature-dependent thermophysical properties and latent heat of phase transformations as an internal heat source. The latent heat has a very small effect on the heat flux estimation for fast cooling conditions; however, it can lead to an interpretation bias in medium and slow cooling conditions due to temporary underestimates of the heat flux, reaching errors up to 100%. Even when the estimated heat fluxes seem accurate (average errors smaller than 10%), further estimates of the temperature evolution or phase transformations kinetics are inaccurate (range of uncertainties as large as 200 • C and 20%, respectively) because of the phase-dependent thermophysical properties.

Hence, using a linear inverse heat conduction problem for a material undergoing phase transformations is acceptable only for fast cooling conditions (h > 1500 Wm -2 K -1 ) and exclusively to estimate the boundary heat flux without further simulations.

Abbreviations for phases

Introduction

In many steel heat-treating processes, the sample is submitted to an austenitization treatment at high temperature to obtain a material completely formed by austenite. Then, the controlled and continuous cooling of the steel leads to phase transformations in order to get the final microstructures to ensure the desired mechanical properties at the end of the process. In addition, thermal gradients and phase transformations generated during cooling lead to stresses and deformations of the piece that must be mastered in order to obtain adequate final product with controlled residual stresses and distortions. The cooling can be performed using different methods, from natural convection in air, promoting a slow temperature decrease in the material, to quenching with water sprays or jets, reaching very high heat transfer coefficients. Slow cooling rates tends to enhance formation of ferrite and pearlite for low-alloyed steels, while faster cooling enhances martensite formation. The products of austenite decomposition also depend on the steel composition and initial austenitic microstructure (grain size, presence of crystallographic defects, and others).

In any case, phase transformations make the heat conduction analysis and boundary heat flux estimation more complicated because of two main points. First, phase transformations change the material thermophysical properties in time and space as these properties are both temperature-and phase-dependent. Second, phase transformations are exothermic processes during cooling, so there is a latent heat of solid-solid phase transformations that acts as an internal heat source. This latent heat can be very significant as it may reach values as high as a few hundreds of MJ/m 3 , which has motivated, for example, its use for thermal energy storage in a solid-state phase-change material module [START_REF] Sharar | High-capacity high-power thermal energy storage using solidsolid martensitic transformations[END_REF].

This problem is already well-known by metallurgical researchers, who performed direct simulations of the heat conduction (i.e. by imposing the boundary condition at the surface in order to predict the temperature evolution inside the material). In 1985, Fernandes et al. [START_REF] Fernandes | Mathematical model coupling phase transformation and temperature evolution during quenching of steels[END_REF] performed numerical simulations of steel cooling processes, which were validated with experimental data, and showed that the internal heat source can be sufficiently high to promote a recalescence in the material, which is a temporary increase in the sample temperature. Furthermore, it was shown that internal stresses that develop during cooling in a metallic piece cannot be disregarded in the prediction of transformation kinetics as well as in the prediction of the cooling laws [START_REF] Denis | Influence of stresses on the kinetics of pearlitic transformation during continuous cooling[END_REF] 1 . Chen et al. [START_REF] Chen | A study on heat conduction with variable phase transformation composition during quench hardening[END_REF], with a numerical simulation study, observed that neglecting the latent heat of phase transformations and internal surface radiation affects the temperature evolution calculation during the cooling of a hollow cylinder. Although they claim that neglecting both parameters introduce errors, the latent heat had a much stronger effect than the radiative heat transfer on the internal surface. Edalatpour et al. [START_REF] Edalatpour | Effect of phase transformation latent heat on prediction accuracy of strip laminar cooling[END_REF] also observed this important effect of the latent heat with simulations of a metal strip cooled by water jets. Moreover, they analyzed the effect of thermophysical properties changes during the cooling, which has some influence but is still minor compared to the latent heat. Fu et al. [START_REF] Fu | Latent heat model of a medium plate in case of jet quenching[END_REF] also analyzed the effect of the latent heat of phase transformations after modeling an industrial-scale quenching experiment. In 2018, Huang et al. [START_REF] Huang | Effect of controlled cold air distribution on temperature profile and phase transformation of wire loops in the stelmor air-cooling process[END_REF] performed fluid flow and heat transfer simulations of a Stelmor air-cooling process of wire loops and they included a phase-transformation model to take into account the internal heat source due to the latent heat of phase transformations.

We find as well several studies using inverse methods to estimate the boundary heat flux considering the phase and temperature dependence of thermophysical properties and the presence of the latent heat of phase transformation. For example, Azim et al. [START_REF] Azim | Experimental validation of inverse heat conduction method: quenching of steels xc42 and xc80[END_REF] solved numerically the inverse problem coupled with phase transformation calculations using the finite difference method with Beck's function specification method, including the effect of internal stresses in the transformation kinetics. They showed that the heat fluxes can be predicted quite well for martensitic transformation during cooling. On the other hand, for pearlitic transformation that leads to recalescence, the prediction of heat fluxes needs to include as well the effect of internal stresses on the phase transformation kinetics. After integrating a modified pattern search method into a finite element model, Wang et al. [START_REF] Wang | Estimation of heat transfer coefficient and phase transformation latent heat by modified pattern search method[END_REF] estimated simultaneously the heat transfer coefficient and the latent heat of phase transformations. They showed that neglecting the latent heat of phase transformations could induce a bad estimate of the temperature evolution in the material. Telejko [START_REF] Telejko | Analysis of an inverse method of simultaneous determination of thermal conductivity and heat of phase transformation in steels[END_REF] proposed an inverse method based on the finite element method to estimate simultaneously the material latent heat of transformation and thermal conductivity. Although the total generated heat was well estimated by his model, estimating the rate of generation was more difficult because it depended, for example, on the starting point chosen in the inverse calculations.

The problem is that necessary data to model phase transformations are often not accessible for thermal engineers who only desire to estimate the dissipated heat of their cooling system. These consist of thermodynamic data, thermophysical properties for each phase and as a function of temperature, latent heat of phase transformations, and finally, the knowledge of the phase transformations kinetics, including the effects of stresses or initial plastic deformations. For this reason, we find many fast cooling studies where the researchers used samples made of materials that do not undergo phase transformations, like Ni201 [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF][START_REF] Labergue | Combined three-color lif-pda measurements and infrared thermography applied to the study of the spray impingement on a heated surface above the leidenfrost regime[END_REF], which is >99% nickel, aluminum alloys [START_REF] Sako | Quenching of a heated wall with spatial temperature gradient using a liquid film through oblique jet impingement[END_REF] and, more commonly, austenitic stainless steels such as AISI 304 [START_REF] Leocadio | Experimental investigation of heat transfer characteristics during water jet impingement cooling of a high-temperature steel surface[END_REF][START_REF] Gomez | Film boiling in quench cooling with hightemperature jets[END_REF][START_REF] Tsukamoto | On the onset of quench during spray cooling: The significance of oxide layers[END_REF][START_REF] Leocadio | Rewetting and boiling in jet impingement on high temperature steel surface[END_REF][START_REF] Karwa | Experimental investigation of free-surface jet impingement quenching process[END_REF][START_REF] Lee | Quasi-steady front in quench subcooled-jet impingement boiling: Experiment and analysis[END_REF] or AISI 309 [START_REF] Kotrbacek | Influence of water temperature on spray cooling at high surface temperatures[END_REF]. Most of the studies with stainless steel considered the material's temperature-dependent thermophysical properties in their inverse method to estimate the boundary heat flux. In turn, several researchers set constant values for the thermophysical properties in their inverse models, like Karwa and Stephan [START_REF] Karwa | Experimental investigation of free-surface jet impingement quenching process[END_REF] for AISI 304, Sako et al. [START_REF] Sako | Quenching of a heated wall with spatial temperature gradient using a liquid film through oblique jet impingement[END_REF] for aluminum alloy A5052, and the researchers who used Ni201 samples. Kita et al. [START_REF] Kita | Quenching mechanism of spray cooling and the effect of system pressure[END_REF] study is a recent example that used a test section made of pure iron in their spray cooling experiments, which is a material with a highly temperature-dependent specific heat. They used the lumped system analysis to estimate the dissipated heat by spray cooling. This is a strong simplification with several limitations that were recognized by the authors in their analysis.

Some researchers performed their experiments using steels in which phase transformations take place. One example is the study by Nobari et al. [START_REF] Nobari | Heat transfer of a stationary steel plate during water jet impingement cooling[END_REF], who performed jet cooling experiments using a sample made of HSLA steel. In their inverse problem model, they considered the material temperature-dependent properties but not the internal heat source due to phase transformation. They used the inverse method developed by Zhang [START_REF] Zhang | Study of boiling heat transfer on a stationary downward facing hot steel plate cooled by a circular water jet[END_REF], who recognized that neglecting the latent heat of phase transformations could result in underestimated values for the estimated heat flux. Kashyap et al. [START_REF] Kashyap | Transient bottom jet impingement cooling of steel[END_REF] also neglected the latent heat in their inverse problem model. Lee et al. [START_REF] Lee | Role of quenching method on cooling rate and microstructure of steels: Variations in coolant and its flow arrangement[END_REF] is a good example of the difficulty thermal engineers have when testing steels because of phase transformations. They used a low-alloy steel sample in some of their cooling experiments to analyze post-quenching mechanical properties and microstructures; however, to estimate the boundary heat flux in their inverse problem, they performed tests with an AISI 310S stainless steel sample, avoiding the strong non-linearity in the heat conduction modeling due to phase transformations.

Despite this large literature on phase transformations and cooling processes in the last decades, thermal engineers still attempt to find solutions to overcome the phase transformations complexity in inverse heat conduction problems to estimate the boundary heat flux. Using a simplified heat conduction model for steels, i.e. assuming constant thermophysical properties and neglecting the internal heat source due to phase transformations, is certainly the easiest solution. However, how well does it perform to estimate the boundary heat flux with inverse methods for a material that undergoes phase transformations? The literature still lacks quantitative information about the accuracy of using a linear model for the inverse problem when phase transformations takes place. Furthermore, if such a simplified model can be used, it is also important to understand its limitations and in which conditions it provides reliable estimates of the boundary heat flux.

In this study, we performed several numerical simulations to evaluate the accuracy and limitations of using a linear heat conduction model to estimate the boundary heat flux in the presence of phase transformations in the material. Different cooling conditions were tested for different materials, three that do not have phase transformation but have temperature-dependent thermophysical properties (AISI 304, Ni201 and Fe) and one that undergoes phase transformations (42CrMo4 steel). Consequently, the latter has not only temperature-and phase-dependent thermophysical properties but also an internal heat source due to the latent heat of phase transformations. We expect that this paper will guide future thermal engineering investigations to establish in which conditions non-linear phenomena can be neglected without compromising the heat flux estimates, which can be useful to compare different cooling conditions or quenching systems. We also discuss how acceptable errors in the heat flux estimates (at least for thermal engineers) can result in large uncertainties in further simulations of the material phase transformations kinetics, which would be unacceptable, for example, in metallurgical applications.

Methods and materials' properties

In this study, we solve both the direct and inverse one-dimensional heat conduction problems but using different methods for each case. In the direct problem, we impose a heat transfer coefficient at the boundary x = 0 of a planar wall and simulate thermocouple measurements T T C at x = x T C . We solve this direct problem using the finite difference method considering temperature and phase-dependent thermophysical properties and an internal heat source due to phase transformations. We introduce the kinetics of phase transformations and material properties calculation further in this section. For the inverse problem, we use the simulated temperature measurements T T C to estimate the heat flux at the wall φ w (t) using an inverse method based on a pseudo-analytical solution of the heat equation considering constant properties and no internal heat source -the aforementioned linear heat conduction model. Beck's function specification method [START_REF] Beck | Inverse Heat Conduction: Ill-Posed Problems[END_REF] was used for regularization. After presenting in detail the direct problem calculation with phase transformations and the inverse method to estimate the boundary heat flux, we explain at the end of this section how we processed the simulation data to analyze the results.

We should highlight that, although this study is based only on simulation results, the methods and models we adopted to solve the direct problem of heat conduction with phase transformation have already been validated with experimental data for different steels and published in the last thirty-five years by IJL research group in University of Lorraine [START_REF] Fernandes | Mathematical model coupling phase transformation and temperature evolution during quenching of steels[END_REF][START_REF] Azim | Experimental validation of inverse heat conduction method: quenching of steels xc42 and xc80[END_REF][START_REF] Devynck | Influence of boiling heat transfer and phase transformations on the deformation of a steel tube during quenching by impinging water jets[END_REF][START_REF] Jeyabalan | Effect of carbon and nitrogen enrichment in the austenitic field on the formation of microstructures and residual stresses in carburized and carbonitrided low-alloyed steel parts: experimental study and simulation[END_REF]. Therefore, the simulated thermocouple signals and phase transformation kinetics represent well the results that would be obtained experimentally if the boundary condition was the one imposed in the simulations.

With that said, using simulated results is an interesting approach for this study because we can separate the effect of each parameter and evaluate the performance of the linear inverse method for each cooling condition.

Furthermore, we limited these simulations to a one-dimensional geometry because phase transformations introduce three significant sources of non-linearity to the heat equation, whose effects are better assessed in a simpler domain.

Although the one-dimensional analysis is realistic, for example, in fast jet-cooling processes at the impact location [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF], extending the present analysis to a two or three-dimensional geometry would certainly be more realistic of a general steel-cooling process. Nevertheless, it would add variables, like the boundary condition profile and the number and positioning of thermocouples in the sample, that create a large number of scenarios to evaluate and would be chosen arbitrarily. We evaluated the effect of these two parameters in a previous instrumentation study [START_REF] Oliveira | Instrumentation study for aspersion cooling experiments aiming optimal estimation of the heat flux dissipation by 2D inverse heat conduction problems[END_REF] for a two-dimensional domain with constant thermophysical properties, and we showed that the inverse method performance was directly and concomitantly affected by the number of thermocouples and the sharpness of the boundary condition profile. Because of phase transformations, the material thermophysical properties are dependent on both the temperature and the phase composition (volume fraction of austenite, ferrite, pearlite, bainite and martensite). Therefore, the heat equation that describes our problem is:

Direct problem: numerical solution

∂ ∂x λ(y, T ) ∂T ∂x + q(∂y/∂t) = ρ(y, T )c p (y, T ) ∂T ∂t ( 1
)
where T is the temperature, y is the phase composition (volume fraction of each product), λ, ρ and c p are, respectively, the material thermal conductivity, specific mass and specific heat, and q is the internal heat source because of the latent heat of phase transformation. We can estimate a given thermophysical property Π (specific mass, thermal conductivity or specific heat) in any point of the material based on each phase fraction y i (x, t) and local temperature T (x, t) by using the following linear rule of mixtures:

Π (x, t) = 5 i=1 y i (x, t)Π i [T (x, t)] (2)
where Π i is the thermophysical property of the phase i among the five phases considered. The internal heat source q is proportional to the rate of transformation (∂y i /∂t) and to the latent heat of transformation per unit of volume ∆H i of the product i being formed, thus:

q = ∆H i ∂y i ∂t (3)
As the thermophysical properties depend on y and T , q depends on ∂y/∂t, and y depends on T and t (as we will present in section 2.2 when introducing the phase transformations kinetics simulation), it is evident that Eq. 1 can be highly non-linear. Thus, we solve this equation using the implicit finite difference method [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations[END_REF] of the following equation that is obtained from Eq. 1 after using the product property on the first derivative, dividing the entire equation by λ and assigning j for the element index and k for the time step:

1 λ k+1 j ∂λ ∂x k+1 j ∂T ∂x k+1 j + ∂ 2 T ∂x 2 k+1 j + qk+1 j λ k+1 j = 1 a k+1 j ∂T ∂t k+1 j (4)
where a = λρ -1 c -1 p is the thermal diffusivity. Moreover, we use the following boundary conditions at x = 0 (1 st element, energy balance by finite volume) and x = L (N x -th element, insulated surface):

     - h(T k+1 1 -T amb ) λ k+1 1 -∂T ∂x k+1 1 + qk+1 1 λ k+1 1 ∆x 2 = 1 a k+1 1 ∆x 2 ∂T ∂t k+1 1 for x = 0 T k+1 Nx = T k+1 Nx-1 for x = L (5)
This system of equations can be structured in a tridiagonal matrix form and can be solved using, for instance, LU decomposition to find the values for {T } k+1 (the curly brackets {•} indicate a vector composed of values for each element). However, because of phase transformations, we had to adopt an iterative calculation to converge the calculated temperature by updating both the thermophysical properties and the internal heat source at the time k + 1 (Fig. 2). For each time step, {T } k+1 is estimated using previous profiles for the thermophysical properties {Π} k and internal heat source { q} k . Then, all the phase transformations equations (presented in the next subsection) are solved using this calculated {T } k+1 , which implicates re-estimating the phase composition profile {y} k+1 , as well as {Π} k+1 and { q} k+1 . After, we compare the calculated {T } k+1 with the precedent vector {T ref }. If the summation of their difference is larger than a maximum acceptable deviation ϵ 0 (set as 0.1 • C in this study), convergence is still not achieved, so the finite difference method equations (Eqs. 4 and 5) are rewritten with the new values for {T } k+1 , {Π} k+1 and { q} k+1 and the iteration restarts. Otherwise, the calculation converged, the temperature T T C of the simulated thermocouple at the time step t k+1 was computed by interpolation at x = x T C and the calculation proceeded to the next time step.

If the material does not have phase transformation and only has temperature-dependent properties, the calculation method is the same but without the phase transformation calculations, hence without considering phase-dependent properties (thus, Π(T ) instead of Π(y, T )) nor the internal heat source ( q = 0).

Phase transformations: kinetics and modeling

The solution to the direct problem is directly linked to the phase transformations kinetics, so we must introduce this phenomenon and its related equations to close the heat conduction model. If a steel sample is at the austenitization temperature for sufficient time, its phase composition is entirely austenite. As the material is cooled down and its temperature decreases, phase transformations take place and form either ferrite, pearlite, bainite or martensite. The transformation products and their rates of formation depend on the evolution in the material temperature, which are typically represented in a temperature-time-transformation diagram. There are two major types of phase transformation: diffusive transformation, where the kinetics depends on time and temperature, typically found for ferrite, pearlite and bainite transformations; and displacive transformation, typically martensitic, which starts when the temperature is below the martensite transformation temperature M s , and its kinetics depends only on the temperature.

Diffusive transformations of a phase i in isothermal conditions can be described by the Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation [START_REF] Johnson | Reaction kinetics in processes of nucleation and growth[END_REF][START_REF] Avrami | Kinetics of phase change. i general theory[END_REF][START_REF] Avrami | Kinetics of phase change. ii transformation-time relations for random distribution of nuclei[END_REF][START_REF] Avrami | Granulation, phase change, and microstructure kinetics of phase change. iii[END_REF][START_REF] Jung | Finite-element simulation of quenching incorporating improved transformation kinetics in a plain medium-carbon steel[END_REF]:

y i (t) = y d y max i (T ) 1 -exp -k i (T ) [t -t i (T )] ni(T ) (6)
where y i (t) is the volume fraction of the phase i formed at the instant t, y d is the fraction of austenite still available for transformation in the material, y max i (T ) is the maximum volume fraction of phase i that can be formed at the temperature T , which, in some cases, corresponds to a thermodynamic limit, k i (T ) and n i (T ) are temperaturedependent coefficients of the sigmoid function, which are usually obtained experimentally, and t i (T ) is the incubation time for phase i, i.e. the time to start phase-i transformation. Nevertheless, cooling processes are non-isothermal, so Eq. 6 is no longer valid as presented. Therefore, we use the additivity hypothesis, which consists of calculating progressive isothermal transformations for a differential time dt, increasing the phase i fraction by dy i . As described by Jeyabalan et al.

[36], we can rewrite JMAK equation in the following form:

dy i = y max i (1 -ξ i ) dξ e i (7)
with:

dξ e i = n i k 1 n i i ξ e i 1-1 n i dt (8)
where ξ e i = -ln (1 -ξ i ) is the extended volume for the phase i and

ξ i = y i / (y d y max i
) is the transformation progression of the same phase, which means phase-i transformation is still possible up to ξ i = 1. Therefore, Eqs. 7 and 8 allow applying JMAK's equation to a non-isothermal cooling process.

As t i is also temperature-dependent, the additivity hypothesis is also used to predict the incubation period by using the sum of Scheil S [37], which is an incubation progression in non-isothermal processes that is defined by:

S(t n ) = f Nt n k=1 ∆t k t i (T k ) (9)
where ∆t k is the time step, T k is the temperature at the same time step and N tn is the number of time steps until t n . While in an isothermal processes phase transformations starts at t i , in non-isothermal processes it only begins when S = 1. The coefficient f is a correction factor that is necessary when the additivity hypothesis is not valid, which can happen when the temperature is below the bainite start temperature B s . For example, if f = 0, which is the case of 42CrMo4 steel, the one we used in our simulations, bainite transformation does not depend on the incubation period of transformations that would take place at higher temperatures (ferrite and pearlite), hence the sum of Scheil is restarted for the incubation of the bainite transformation.

In turn, the martensitic transformation, which is displacive, occurs when T < M s and is described by the Koistinen-Marburger law [START_REF] Koistinen | A general equation prescribing the extent of the austenite-martensite transformation in pure iron-carbon alloys and plain carbon steels[END_REF]:

y M (T ) = y d (M s ) {1 -exp [-α (M s -T )]} ( 10
)
where y M is the martensite volume fraction in the material at T , α is an empirical coefficient, M s is the martensite start temperature, and y d (M s ) is the available austenite fraction when martensitic transformation started.

In the course of phase transformations, the product formed at each time step is kept until the end of the cooling, so the hardness of each new product being formed contributes to the final material hardness. Because austenite decomposition products have hardnesses that depend on the temperature of formation (except martensite), the material hardness is also estimated using a rule of mixtures but considering the phase hardness at the temperature at which it is formed. In other words, the Vickers hardness HV calculation is performed progressively by additivity departing from austenite hardness HV A (initial condition), as follows for a given time t n :

HV (t n ) = HV A + Nt n k=i {HV i [T (x, t k )] -HV A } dy i ( 11
)
where HV i is the hardness at room temperature of the phase being formed at T (t k ). Note the summation in Eq. 11

is the additive part that is calculated progressively with the time step, increasing or reducing the material hardness with the transformation product when compared to austenite.

Inverse problem: pseudo-analytical solution

With the knowledge of one temperature measurement (in our case, the simulated measurement T T C (t) using our direct problem solution), we can estimate the heat flux evolution φ w (t) at x = 0 of our 1D heat conduction problem using an inverse method. As we explained in the introduction, we analyze in this study how a simplified model for the inverse method can perform with a material undergoing phase transformations. We used the same method as in previous studies [START_REF] Oliveira | Thermocouple response time estimation and temperature signal correction for an accurate heat flux calculation in inverse heat conduction problems[END_REF][START_REF] Oliveira | Inverse ARX (IARX) method for boundary specification in heat conduction problems[END_REF] based on the pseudo-analytical solution of the following transient 1D heat equation considering constant thermophysical properties and no internal heat source:

∂ 2 T ∂x 2 = 1 a ∂T ∂t (12)
The solution can be found using the thermal quadrupoles method [START_REF] Maillet | Thermal Quadrupoles: Solving the Heat Equation through Integral Transforms[END_REF] with Duhamel's theorem [START_REF] Õzişik | Heat Conduction[END_REF], which assumes the heat flux φ w (t) is constant at each time step. Therefore, we can write the temperature response at x T C as:

θ T C (t k+1 ) = T T C (t k+1 ) -T 0 = k m=0 X k-m φ(t m ) (13)
X k = - 1 λ t k+1 t k Z(x T C , τ )dτ (14)
Z = L -1 1 p α cosh p α L sinh p α L cosh p α x T C -sinh p α x T C ( 15
)
where T 0 is the initial temperature, θ T C = T T C -T 0 is a temperature difference, and p is the Laplace variable. The inverse Laplace transform of Z is performed using the Stehfest algorithm [START_REF] Stehfest | Algorithm 368: Numerical inversion of laplace transforms [d5[END_REF].

The inverse method we used is the function specification method proposed by Beck [START_REF] Beck | Inverse Heat Conduction: Ill-Posed Problems[END_REF], which consists in filtering the temperature measurement noise using a functional for the heat flux at N f ts future time steps. We used the simplest functional of assuming constant heat flux for these future time steps: considering we are estimating the heat flux at a time step k, we have φ(t k ) = φ(t k+1 ) = ... = φ(t k+N f ts ). Therefore, with few mathematical manipulations and assuming past heat fluxes were already estimated, we can write Eq. 13 in the following form for the future N f ts time steps:

                   (X 0 ) φ(t k ) = θ T C (t k+1 ) - k-1 m=0 X k-m+1 φ(t m ) for k (X 0 + X 1 ) φ(t k ) = θ T C (t k+2 ) - k-1 m=0 X k-m+2 φ(t m ) for k + 1 . . . N f ts -1 m=0 X m φ(t k ) = θ T C (t k+N f ts ) - k-1 m=0 X k-m+N f ts φ(t m ) for k + N f ts -1 (16)
This system of N f ts equations can also be written in the matrix form using a sensitivity vector X f ts composed of the terms multiplying φ(t k ) on the left-hand side and the result vector d f ts composed of the terms on the right-hand side of Eq. 16:

X f ts φ(t k ) = d f ts (17)
Finally, we use the least squares method to estimate the heat flux at the time step k:

φ(t k ) = X T f ts X f ts -1 X T f ts d f ts (18)

Simulation steps and results comparison

Figure 3 presents a flowchart of the calculation steps. In summary, we solved the direct problem using the finite differences method (section 2.1) imposing a constant h, resulting in an imposed heat flux at the boundary, and considering both the temperature and phase-dependent thermophysical properties and the internal heat source due to phase transformations (section 2.2) -when applicable. Therefore, we are simulating a steel cooling experiment, so the outputs are a simulated thermocouple measurement T T C at x T C and profiles of the material phase fraction and hardness at the end of the cooling. Afterwards, using the simulated thermocouple measurement, we estimated the heat flux at the boundary using the inverse method presented in the previous section assuming constant thermophysical properties and no internal heat source (section 2.3). This estimated heat flux is compared to the imposed one (comparison 1), which is our target. It is important to observe that this comparison is only possible in a numerical study, as the imposed boundary condition is unknown in a real experiment. Hence, we solved again the direct problem using the finite differences method (sections 2.1 and 2.2 again) but now using the estimated heat flux as a boundary condition at x = 0, resulting in an estimated thermocouple response and estimated phase composition and hardness profiles. These results can be compared, respectively, to the simulated thermocouple response (comparison 2) and final phase and hardness profiles (comparison 3) of the baseline simulation. This last step serves as a validation procedure of the estimated heat flux that could be performed in a real experiment, as we compared measurable and estimated outputs. Hence, these three comparisons allowed us to evaluate the quality of the heat flux estimate but also the accuracy of further simulations using the estimated heat flux.

For the inverse method, we must assign constant values for the thermophysical properties in the model. For this choice, the material thermophysical properties were taken at the mean temperature T of the initial and final values of the cooling process (the maximum and the minimum, respectively). Figure 4 presents an example of this property choice for the thermal conductivity, which was done for all the thermophysical properties. Although a single mean temperature value can be assigned for AISI 304, Ni-201 and pure iron, whose properties are only temperaturedependent, another hypothesis is necessary for the simulations with 42CrMo4 steel since its thermophysical properties are also phase-dependent. In this case, we solve the inverse problem using the two extreme conditions: either the material has the properties of the austenitic phase ("Prop. A" scenario) or those of the transformation products (ferritic phases, thus "Prop. F" scenario). Therefore, we should expect to have a lower and upper boundary for the estimated results (heat flux, thermocouple response and phase and hardness profiles). This "Prop. A" and "Prop.

F" results, as well as their range, will be clearer during the analysis of the results in section 3.2.

Simulated materials and properties

As mentioned in the introduction, we first performed the simulations with materials without phase transformation to evaluate separately the effect of temperature-dependent thermophysical properties. The simulated materials were AISI 304 and Ni201, two materials that are commonly used in the literature to replace steel in cooling experiments (AISI 304 and Ni201), and pure iron (Fe), a material with highly temperature-dependent thermophysical properties.

More precisely, its specific heat presents a substantial peak near the Curie temperature at about 770 • C, resulting in a thermal diffusivity that varies strongly with temperature as well. Figure 5 presents the thermophysical properties for these three materials. For the simulations with a material undergoing phase transformations (42CrMo4 steel), as a matter of simplification, we used the same thermophysical properties for the "ferritic phases" (i.e., for ferrite, pearlite, bainite and martensite) and the data for austenite shown in Fig. 5. The input data for the phase transformation calculation (the temperature-dependent JMAK parameters k i (T ) and n i (T ) for each product used in Eq. 8) are obtained from the IT (isothermal transformation) diagram of 42CrMo4 steel, which is presented in Fig. 6a. They were taken from Fernandes et al. [START_REF] Fernandes | Mathematical model coupling phase transformation and temperature evolution during quenching of steels[END_REF], as well as the Koistinen-Marburger coefficient (α = 1.1 • 10 -2 K -1 ) and the values of ∆H i , the volumetric latent heat of phase transformation, for each product formation: 5.9 • 10 8 J/m 3 for ferrite and pearlite, 2.4 • 10 8 J/m 3 for bainite, and 4.4 • 10 8 J/m 3 for martensite. Finally, Fig. 6b presents the hardness of each transformation product as a function of its formation temperature, which is necessary to calculate the material hardness with Eq. 11.

As the material's initial and final temperatures of the simulations are already known (starting at T 0 = 850 • C and ending at the ambient temperature T amb = 20 • C), the values of the thermophysical properties at T = 435 • C used in the inverse method can already be defined for each material. They are presented in Table 1. Table 1: Thermophysical properties at the mean temperature (435 • C) that were used in the inverse method for each material.

Material

Results and discussion

First, we present the results of cooling simulations for AISI 304, Ni201 and Fe to evaluate the effect of temperaturedependent properties on the boundary heat flux estimation. Two cooling rate conditions were analyzed, one with 290 a high heat transfer coefficient (h = 10, 000 Wm -2 K -1 ), typically found in quenching processes using water jets or sprays, and one with a low heat transfer coefficient (h = 100 Wm -2 K -1 ), which is of the same order as air-cooling processes. The simulated data acquisition rate of the virtual thermocouple measurement T T C is 50 Hz and 0.5 Hz for the simulation with the high and the low heat transfer coefficients, respectively. These values were chosen to be representative of a real experiment where the investigator tries to acquire as much data as possible to observe the temperature transients but excessive to respect the limit of the data acquisition system.

Then, using the same heat transfer coefficients (10,000 and 100 Wm -2 K -1 ), we present results with 42CrMo4 steel to evaluate the effect of latent heat of phase transformations and phase-dependent properties on the inverse method estimations and discuss the physical processes during the cooling. Finally, we analyze the results for different h to understand when the use of a linear heat conduction model is appropriate for estimating the heat flux with inverse methods when the material undergoes phase transformations. For all the results with 42CrMo4 steel, we also discuss the effect of the estimated heat flux on further simulations to estimate the material final phase composition and hardness profiles, which are important parameters in metallurgy.

Although we did not apply any noise in the virtual thermocouple measurements, we fixed N f ts = 3 future time steps for all the inverse method calculations to account for the effect of the regularization method, which is always necessary to process real experimental data.

Results for materials without phase transformation

In this section, we discuss the simulation results based on comparisons 1 and 2 from Fig. 3, i.e. the imposed and estimated heat flux, and the target and estimated temperature response at the thermocouple position. For the condition with a high heat transfer coefficient, we also present the temperature evolution at x = L because there is a high temperature gradient in the x-direction. Target values are represented in the graphs as filled lines, while estimated values are shown as dashed lines.

Figure 7 presents simulation results for h = 10, 000 Wm -2 K -1 and for AISI 304, Ni-201 and pure iron. For AISI 304 and Ni201, whose thermophysical properties are less temperature-dependent than pure iron, the estimated heat fluxes match well the target values, as well as the temperature evolution both at the thermocouple location and at the rear face -although it is slightly underestimated for AISI 304 and slightly overestimated for Ni201. In fact, the deviation from the estimated and target heat fluxes is only significant at the very beginning of the cooling process, where the peak heat flux is estimated to be approximately 6 MWm -2 while the target value is about 8 MWm -2 . This is partially caused by the inverse method, which assumes stationary conditions for t < 0 (zero heat flux until t = 0), affecting the first heat flux estimates, but also by the fast transient observed in the first milliseconds of the cooling process. This local underestimation would be less noticeable if the data acquisition rate were higher, as the hypothesis of constant heat flux used in the functional of Beck's regularization method would be closer to reality.

Nevertheless, the estimated heat flux is very close to the imposed one throughout the cooling, presenting average absolute errors of 5.8% and 4.5% in the first 20 s of cooling for AISI 304 and Ni201, respectively. These errors are smaller than those that could be introduced by model errors, like the thermocouple position or the values of thermophysical properties even if they were constant. For example, in a previous study [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF] in which we performed jet cooling experiments with a 20mm-thick Ni201 plate, we estimated that model errors could result in an uncertainty of 5-10% in the heat flux estimate.

For the pure iron simulations, a material whose thermophysical properties are highly temperature-dependent (especially the specific heat -see Fig. 5), not only the heat flux peak is highly underestimated -less than 5 MW m -2 while it should be 8 MW m -2 -but also the underestimation is persistent until t ≈ 40 s, resulting in a highly overestimated temperature evolution both at the thermocouple position and at the rear face. The mean absolute error of the heat flux estimate in the first 20 s of cooling is 12.5%, more than twice the errors observed with AISI 304 and Ni201. This discrepancy is a result of using a much lower specific heat in the inverse problem model (631 J kg -1 K -1 ) than the true values found between 600 • C and 800 • C (between 800 and 1300 J kg -1 K -1 ), i.e.

near the Curie temperature. When the thermocouple measurement is within this range during the cooling, the lower temperature variation caused by the high specific heat is interpreted by the inverse method as a lower heat flux at the boundary. When this estimated heat flux is applied to the system in the direct problem calculation for validation, the estimated temperature variation becomes highly overestimated because of this underestimated heat flux but also because, once again, the high specific heat in this 600-800 • C temperature range hinders the material cooling. This result complies with the study by Blackwell and Beck [START_REF] Blackwell | A technique for uncertainty analysis for inverse heat conduction problems[END_REF], in which they observed the importance of the material volumetric heat capacity (ρc p ) by heat conduction simulation of a planar wall with constant thermophysical properties. They concluded that it was necessary to be as precise as possible on this property value to improve the estimated heat flux accuracy.

This effect caused by the high specific heat of pure iron at the Curie temperature is even more noticeable in slow cooling conditions (Fig. 8, in which only the temperature at the thermocouple is presented because the temperature gradient in the x-direction is very low). As before, the slower temperature decrease near the Curie temperature is interpreted by the inverse method as a lower heat flux; however, the decrease in the estimated heat flux is very steep, which could lead to an interpretation bias. For example, in a real experiment where the target heat flux is unknown, the investigator could conclude that there was a problem in the cooling system in the first minutes of the experiment or imagine a phenomenological reason for this heat flux behavior that does not exist. Although the heat flux was also underestimated in the previous conditions of high heat transfer coefficient, there was no risk of interpretation bias as we see in this case. As a result, the estimated thermocouple signal is largely overestimated, similarly to the high heat dissipation scenario. For pure iron, we found a 19.5% mean absolute error for the estimated heat flux in the first 1000 s. For the cases with AISI 304 and Ni201, both the heat flux and temperature evolution estimates are very satisfactory, presenting mean absolute errors for the heat flux estimation in the first 1000 s of only 3.1% and 2.3%, respectively. Therefore, these materials are very appropriate for experiments to evaluate cooling systems either for both high and low heat dissipation processes, which justifies their use in many past studies [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF][START_REF] Labergue | Combined three-color lif-pda measurements and infrared thermography applied to the study of the spray impingement on a heated surface above the leidenfrost regime[END_REF][START_REF] Leocadio | Experimental investigation of heat transfer characteristics during water jet impingement cooling of a high-temperature steel surface[END_REF][START_REF] Gomez | Film boiling in quench cooling with hightemperature jets[END_REF][START_REF] Lee | Quasi-steady front in quench subcooled-jet impingement boiling: Experiment and analysis[END_REF][START_REF] Gradeck | Solution of an inverse problem in the hankel space -infrared thermography applied to estimation of a transient cooling flux[END_REF].

Results for 42CrMo4 steel: effect of phase transformations

In this section, we present the results of several simulations with 42CrMo4 steel, a material that undergoes phase transformations. Similar to the previous section, we analyzed the results with this material for a high and a low heat transfer coefficient (10,000 and 100 Wm -2 K -1 , respectively). While the first case is a fast cooling condition that results in the formation of martensite and bainite, the second one is a slow cooling process resulting in ferrite and pearlite as main products (some bainite was also formed). The objective is to discuss in detail these cases that have very different thermal behaviors and, hence, phase transformations kinetics, in order to evaluate the effect of the phase-dependent properties and latent heat due to phase transformation. As explained in section 2.4 and Fig. 4, the inverse method was solved using the mean thermophysical properties of austenitic and ferritic phases ("Prop.

A" and "Prop. F", respectively), so both results were compared with the target values. In each figure, we added a shaded area between the two inverse solutions estimates to aid the visualization of the results' range between these two extremes. In the end, we present and discuss results for different heat transfer coefficients of cooling to establish in which conditions a linear heat conduction model can be used to estimate the boundary heat flux when phase transformations take place.

Fast cooling condition

Figure 9 presents results for h = 10, 000 Wm -2 K -1 at the x = 0 boundary, where mainly martensite and a small amount of bainite are formed at the end of the cooling process. The results shown on the left-hand side of the figure (index "I" in the subfigure identification) are for direct problem solutions including all the phase transformations phenomena, i.e. temperature and phase-dependent properties and internal heat source due to the latent heat. In turn, the results on the right-hand side (index "II") are for direct problem solutions neglecting the latent heat of phase transformation in the model, allowing us to evaluate its effect on the linear inverse model estimates.

We see in Fig. 9a.I that the estimated heat flux is very close to the target when using the austenite thermophysical properties in the inverse method -which was expected since there is only austenite at the beginning of the cooling process. The inset figure shows that even the peak heat flux is fairly well estimated, being as accurate as those obtained previously with AISI 304 and Ni201 for the estimated heat flux. On the other hand, using ferritic thermophysical properties in the inverse method results in overestimated heat fluxes, even though the heat flux estimate is still relatively close to the target.

One could state that the estimated heat flux range between "Prop. A" and "Prop. F" results is narrow or that the error for the estimated heat flux using austenitic properties is very small. However, the effect of the heat flux estimate errors is visible when estimating the temperature evolution at the thermocouple position or at the rear face (Fig. 9b.I). Starting with "Prop. A" results, the estimated steady-state temperature of the body at the end of the cooling is much higher than the ambient temperature -an effect of the latent heat of phase transformations, as shown later. This non-realistic result means that small errors in the heat flux estimate have a significant effect on further estimations of the temperature evolution. Because the ranges of temperature estimates are large, the range of results

for the phase transformations kinetics is also large (Fig. 9c.I and d.I), especially for the martensite formation. In fact, the overestimated steady-state temperature in the "Prop. A" scenario resulted in a residual volume fraction of austenite of almost 20%, while the target volume fraction was only 3% (Fig. 9e.I). Moreover, in the same scenario, the volume fraction of bainite at the end of the cooling was overestimated in half of the domain (x > 10 mm). As a consequence, the final profile of material hardness (Fig. 9f.I) has a relatively wide range between "Prop. A" and "Prop. F" results for this half of the domain. Even if in this tested case the target values (either for heat flux and temperature and phase transformations evolution) are always inside the range using "Prop. A" and "Prop. F" results in the inverse method, such large ranges are unacceptable from a metallurgical point for predicting microstructures and hardnesses. These large errors are a significant issue for further predictions of residual stresses and deformations in metallic pieces, which are extremely sensitive to the temperature gradients and the transformation kinetics [START_REF] Denis | Prediction of residual stress and distortion of ferrous and non-ferrous metals: Current status and future developments[END_REF]. The graphs of the right-hand side of Fig. 9 present the same type of results as those discussed above but, as mentioned before, the latent heat of phase transformations is neglected in the direct problem model ( q = 0) to evaluate its effect on the accuracy of the linear inverse model. The latent heat has a very small effect, virtually negligible, in the estimation of the heat flux at the boundary (Fig. 9a.II), especially when estimating the peak heat flux. Nevertheless, the estimated temperature ranges (Fig. 9b.II) are larger and even converge to unrealistic temperature values (below room temperature, which is a consequence of imposing a heat flux at the boundary instead of the heat transfer coefficient) as compared with the previous case. The ranges for the phase transformations kinetics (Fig. 9c.II and d.II) are still large and the estimation of the final profiles of volume phase fractions (Fig. 9e.II) and material hardness (Fig. 9f.II) show a slightly lower range than when latent heat was present. Thus, we can see that, for fast cooling conditions, the presence of the latent heat of phase transformations practically does not affect the heat flux estimation with the simplified inverse model and the error of the estimated heat flux is mainly related the choice of the thermophysical properties. Therefore, the heat fluxes that were estimated by Nobari et al. [START_REF] Nobari | Heat transfer of a stationary steel plate during water jet impingement cooling[END_REF] and Kashyap et al. [START_REF] Kashyap | Transient bottom jet impingement cooling of steel[END_REF] in their jet cooling experiments using steel are probably correct even though they neglected the latent heat of phase transformations in their inverse problem model. However, even small deviations in heat flux estimated have large effects on further simulations of temperature evolution and, consequently, phase transformation kinetics.

Slow cooling condition

Figure 10 presents the same type of results as Figure 9, but now for a slow cooling condition where h = 100 The effect of the latent heat on the inverse method results is very clear in Fig. 10a.I when the estimated heat flux drops steeply between 200 s < t < 500 s and deviates significantly from the target value. This happens because the aforementioned plateau of the thermocouple measurement is interpreted by the inverse method as a loss of cooling, a similar behavior we observed in the results with pure iron because of its specific heat peak (section 3.1). Moreover, the large drop in the estimated heat flux can lead the investigator to the same interpretation bias as discussed in section 3.1 for the specific heat effect. The temperature evolution (Fig. 10b.I) is largely overestimated once phase transformations start at approximately t = 200 s because not only the dissipated heat flux is largely underestimated by the inverse method but also because the internal heat source takes place in the validation calculation where the estimated heat flux is used as boundary condition of the direct problem. This behavior of the estimated heat flux is strictly related to the presence of the internal heat source, since the heat flux would be correctly estimated if the latent heat did not exist (Fig. 10a.II). Therefore, for slow cooling conditions, the latent heat of phase transformations must be considered in the inverse problem model to estimate the boundary heat flux.

Wm -2 K -1 ,
The results regarding phase transformation are in accordance with the temperature evolutions: as the estimated It is important to notice that using ferritic properties in the inverse method ("Prop.

F" results) provided quite accurate results for the heat flux and temperature evolution estimates. This occurs because phase transformations start and finish relatively early, so there is no austenite in the material in most of the cooling process (for t > 800 s). As a consequence, without latent heat, the phase transformations kinetics would be correctly estimated (for the tested case) with a narrow range between "Prop. A" and "Prop. F" results (Fig. 10c.II, d.II, e.II and f.II). Therefore, for slow cooling conditions, the choice of the thermophysical properties for the linear inverse heat conduction model plays a much smaller role than neglecting the latent of phase transformations.

Influence of the cooling rate

As we have discussed the estimated results using the simplified inverse heat conduction model for two different cooling conditions, we evaluate in this section the performance of this model for different cooling rates. We varied the heat transfer coefficient at the boundary from 50 to 10,000 Wm -2 K -1 , hence considering different phase transformation sequences, from full ferrite-pearlite to full martensite at the end of the cooling simulation. To simplify the results' analysis, we solved the inverse problem only once in these simulations by using the average thermophysical properties between "Prop. A" and "Prop. F" at T . Figure 11 presents the mean errors for the estimated heat flux and average material hardness, represented respectively as ϵ φw and ϵ HV f , as a function of the average hardness at the end of the cooling HV f . The advantage of using HV f is that we can correlate the estimation errors to the phase transformations processes, the final fraction of each phase (presented above the graph), and the degree of heat dissipation -for the 42CrMo4 steel, a higher heat transfer coefficient leads to a final product with a higher hardness.

These three parameters are respectively defined as follows:

ϵ φw = 1 φ w,targ 1 N t Nt k=1 [φ w,est (t k ) -φ w,targ (t k )] 2 (19) ϵ HV f = 1 N x Nx j=1 [HV f,est (x j ) -HV f,targ (x j )] 2 (20)
HV f = 1 N x Nx j=1 HV f,targ (x j ) (21)
where φ w,targ is the average heat flux for 0 < t < t Nt , where N t is the number of time steps used to calculate ϵ φw , so:

φ w,targ = 1 N t Nt k=1 φ w,targ (t k ) (22)
Before discussing the results in Fig. 11, we must highlight that N t used in Eqs. 19 and 22 is the number of points from the initial condition, where the entire body is at 850 • C, until the instant when the surface temperature reaches 100 • C in the simulated experiment. This was done to concentrate the error analysis at higher heat fluxes, the range that interests most thermal engineers, and to avoid an excessive amount of data at very low heat fluxes, where the cooling duration is much longer. Furthermore, we preferred not using relative errors for the heat flux error calculation because they are significantly higher for lower heat fluxes, which are the denominator in the relative error calculation. Thus it could bias the analysis for the highest heat fluxes where we are more interested, especially the peak heat flux. For this reason, the normalization of the heat flux error was done with an average heat flux φ w,targ for the first N t values.

By analyzing the results in Fig. 11, we observe a systematic decrease in the heat flux estimation error with the increase in the heat transfer coefficient, which is indicated on the curve. This complies with what we discussed in the previous section: the effect of the latent heat of phase transformations becomes less and less important as the heat dissipation is intensified at the boundary. Nevertheless, the error in the final hardness estimation does not present any correlation with the heat transfer coefficient. This is because the temperature evolution estimates are highly affected both for fast and slow cooling processes, respectively by the choice of the thermophysical properties and by neglecting the latent heat of phase transformations. As a consequence, incorrect estimates of the temperature evolution can result incorrect estimates of the material final hardness profile, since the hardnesses of ferrite, pearlite and bainite are highly temperature-dependent (see Fig. 6). For example, slow cooling processes resulting in ferrite, pearlite and bainite formation are sensitive to error in the cooling rate estimate, as discussed in section 3.2.2. That was also the scenario for h = 150 Wm -2 K -1 , where the error in the material hardness estimate has a significant peak (Fig. 11). In this case, the target volume fraction of bainite in the piece was about 70% but the estimated value was between 5% and 18% along the thickness, while pearlite volume fraction was estimated to be more than 50% but the target was practically zero. This poor estimate of the phase transformation kinetics resulted in poor estimates of the material phase fraction composition and, hence, of the material hardness.

When the cooling process is fast enough and results in a predominant fraction of martensite (over 80%), the estimation error of the heat flux is less than 15%, which can be acceptable in studies aiming to compare different cooling systems as it is similar to uncertainties in studies using materials without phase transformation [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF][START_REF] Lee | Role of quenching method on cooling rate and microstructure of steels: Variations in coolant and its flow arrangement[END_REF].

Furthermore, the more martensite is formed during the cooling process, the smaller the error in the final hardness estimation. Indeed, as martensite fraction depends only on temperature and not on time, an error in the cooling rate estimate would not affect the estimation of the final material phase fraction nor the hardness at room temperature. Even when the phase fractions of martensite and bainite are similar at the end of the cooling, as for h = 1, 500 Wm -2 K -1 in Fig. 11, the simplified inverse model performs well in the estimation of the heat flux with no interpretation bias. In fact, this value of h seems to be the threshold that divides reliable estimations from interpretation bias. Figure 12 presents comparisons of the target and estimated heat fluxes in time for different heat transfer coefficients (normalized values). These results were also obtained by solving the linear inverse problem using the average properties of austenitic and ferritic phases at T . Although the estimated heat flux has a larger error for h = 1, 500 Wm -2 K -1 than for a higher h, the estimated value decreases constantly, with no abnormal behavior during the cooling. On the other hand, for lower heat transfer coefficients, the effect of the latent heat is visible as the heat flux is underestimated for a given period, which becomes more and more noticeable as the heat transfer coefficient decreases. Therefore, for the simulated conditions and material, there is no interpretation bias for h > 1, 500 Wm -2 K -1 . To illustrate what this value of heat transfer coefficient means in a real application, the peak heat flux for this h is about 1 MW m -2 , which is relatively low for jet or spray cooling with water in industrial conditions. It was, in fact, the lowest peak heat flux observed in our previous study [START_REF] Oliveira | Experimental study of the heat transfer of single-jet impingement cooling onto a large heated plate near industrial conditions[END_REF] for the lowest water jet flow rates and far from the impact location.

Hence, we can conclude that the use of a simplified heat conduction model for the inverse problem, in which it is assumed constant thermophysical properties and no internal heat source due to the latent heat of phase transformations, has an acceptable performance in fast cooling processes that results mostly in martensite formation and only for the heat flux estimation. On the other hand, it has to be noticed that this conclusion is no longer valid when predicting temperature evolutions and gradients in a piece, which control the progress of martensitic transformation and impacts largely residual stresses and deformations.

Hence, in real fast cooling experiments of steels like spray aspersion, jet impingement or bath immersion, the investigator may use the linear inverse model to estimate the boundary condition even if the material undergoes phase transformations. If the estimated heat transfer coefficient is higher than 1,500 Wm -2 K -1 , the results probably have acceptable errors -the same order that would be found using a material without phase transformation in the experiment -and with no risk of interpretation bias. However, if lower values are obtained, the results are probably affected by the latent heat of phase transformations and the estimates are not reliable, so the investigator risks misinterpreting the results. In this case, numerical inverse methods are necessary to account for non-linear phenomena due to phase transformations in the heat equation, which has a much higher computational cost than using the linear inverse method as it requires solving the direct problem. For example, in this one-dimensional study, the direct problem for h = 10, 000 Wm -2 K -1 took about 600 s to be solved, while the linear inverse method took only 2 s (Intel Core i5 12400 4.40 GHz, 6 cores, 32 GB RAM). Hence, the linear inverse method is an important tool for the investigator who desires to obtain fast data processing and estimates of the boundary condition.

Conclusions

Heat conduction in materials undergoing phase transformations is difficult to model because it is highly nonlinear and the models and the necessary data to take into account phase transformations are not always accessible.

Hence, in this study, several one-dimensional heat conduction simulations were performed to evaluate the accuracy and limitations of using a linear, simplified model for the inverse heat conduction problem to estimate the boundary heat flux when phase transformations take place. In this linear model, the thermophysical properties are assumed to be constant and the latent heat of phase transformations is neglected. Simulations with materials without phase transformation were also performed to separate the effect of temperature-dependent properties from other non-linear phenomena intrinsic of phase transformations, like phase-dependent properties and the existence of the latent heat of phase transformation.

The main conclusions of this study are:

The linear inverse model performs well for materials without phase transformation that show small variations of thermophysical properties with temperature, like AISI 304 and Ni201, which are two very commonly used materials replacing steel in fast cooling studies for determining heat fluxes.

On the other hand, it performs poorly for materials without phase transformation but with highly temperaturedependent properties like pure iron, which has a high peak of specific heat near its Curie temperature. In this case, the estimated heat flux is strongly underestimated and can lead to an interpretation bias of the results for slow cooling conditions.

For a material undergoing phase transformations like 42CrMo4, the steel whose properties were used in this study, the linear inverse model performs well for fast cooling conditions, as the latent heat of phase transformations has a practically negligible effect on the heat flux estimation.

However, for slow cooling conditions, the internal heat source due to the latent heat affects the linear inverse model results, as it can largely underestimate the boundary heat flux while phase transformations take place, reaching errors up to 100%. This momentary underestimated heat flux can also bias the results interpretation.

In any cooling condition, further simulations of temperature evolution or phase transformation kinetics using the estimated heat flux as boundary condition provide uncertain results for fast cooling conditions and incorrect estimates for slow cooling conditions, the latter because of the latent heat of phase transformations. For example, small mean errors in the heat flux estimate, as low as 10%, resulted in uncertainties in the temperature evolution and phase fraction estimates by 200 • C and 20%, respectively. Therefore, the use of the linear inverse model must be restricted only to estimate the boundary heat flux.

Simulation results for several cooling rates, i.e. heat transfer coefficients at the boundary, showed that the linear inverse model provides accurate heat flux estimates with no risk of interpretation bias only for fast cooling conditions where h > 1, 500 Wm -2 K -1 . Otherwise, a numerical inverse method including phase transformation calculation must be used, which is significantly more time consuming than using the linear inverse method.

Funding

This

  (K or • C) ϵ 0 maximal acceptable deviation (K or • C) λ thermal conductivity (Wm -1 K -1 )Π arbitrary material property (a.u) ρ specific mass (kgm -3 ) θ temperature difference (K or • C) phase transformation (Jm -3 ) q internal heat source (Wm -3)T mean temperature (K or • C) a thermal diffusivity (m 2 s -1 )c p specific heat (Jkg -1 K -1

Figure 1

 1 Figure 1 presents the transient one-dimensional heat conduction problem studied in this article, consisting of a planar wall with the x = 0 boundary submitted to a transient heat flux φ w (t) = h [T w (t) -T amb ], where h is the heat transfer coefficient (assumed constant in our simulations), T w (t) = T (x = 0, t) is the surface temperature and T amb is ambient temperature. The heat flux is considered positive if heat is extracted from the body. The x = L boundary is considered thermally insulated. The calculated temperature evolution T T C (t) at x = x T C simulates the temperature measurements of a virtual thermocouple, which is, for instance, the experimental data we would have in a real cooling experiment. In this study, we fixed some values of the problem geometry and the initial condition: L = 20 mm, x T C = 1 mm, T (x, t = 0) = 850 • C. The ambient temperature was also set as T amb = 20 • C.

Figure 1 :

 1 Figure 1: Illustration of the 1D heat conduction problem.

Figure 2 :

 2 Figure 2: Flowchart of the iterative calculation to include phase transformations in the heat equation solution, and algorithm of the phase transformation model.

Figure 3 :

 3 Figure 3: Flowchart of the calculation steps highlighting the parameters compared in the results analysis. White boxes and filled arrows are related to the first simulation with the imposed h and heat flux, generating the target values, while gray boxes and dashed arrows are related to estimated values to compare with the target ones. Legend: FD = finite difference; TDP = temperature-dependent properties; PT = phase transformation; PAS = pseudo-analytical solution; CP = constant properties.

Figure 4 :

 4 Figure 4: Choice of the thermophysical property at the mean temperature T for the linear inverse heat conduction model (example with the thermal conductivity).

Figure 5 :

 5 Figure 5: Thermophysical properties of austenite, ferritic phases (i.e. pure iron and all the phase transformations products during the cooling), AISI 304 and Ni-201 as a function of the temperature: (a) thermal conductivity; (b) specific heat; (c) specific mass; (d) thermal diffusivity.

Figure 6 :

 6 Figure 6: Characteristics of 42CrMo4 steel: (a) IT (isothermal transformation) diagram, where filled (---) and dotted (• • • • • • ) lines represent respectively the start and the end of the phase transformations and dashed lines (---) represent key phase transformations temperatures (A 3 , A 1 , Bs, Ms); (b) Vickers hardness for each phase, measured at room temperature, as a function of the formation temperature.

Figure 7 :

 7 Figure 7: Cooling simulation results with h = 10, 000 Wm -2 K -1 for AISI 304, Ni201 and pure iron (Fe): comparison of the target and estimated heat flux (comparison 1 in Fig. 3) and target and estimated temperature evolutions (comparison 2) at the thermocouple position x = x T C and at rear surface x = L. Filled (--) and dashed lines (---) are, respectively, target and estimated values.

Figure 8 :

 8 Figure 8: Cooling simulation results with h = 100 Wm -2 K -1 for AISI 304, Ni201 and pure iron (Fe): comparison of the target and estimated heat flux (comparison 1 in Fig. 3) and target and estimated temperature evolutions (comparison 2) at the thermocouple position x = x T C and at rear surface x = L. Filled (--) and dashed lines (---) are, respectively, target and estimated values.

Figure 9 :

 9 Figure 9: Cooling simulation results with h = 10, 000 Wm -2 K -1 (fast cooling) for 42CrMo4 steel, considering or neglecting the internal heat source due to latent heat of phase transformation (graphs on the left-hand and right-hand sides, respectively): (a) heat flux evolution; (b) temperature evolution at the thermocouple position and rear face; (c) product formation evolution at the thermocouple position and (d) at the rear face; (e) phase fraction profiles at the end of the cooling; (e) material hardness profiles at the end of the cooling. Filled lines (--) are the target, dashed lines (---) are estimated values using properties of the austenitic phase (Prop. A), and dashed-dotted lines (-• -• -) are estimated values using properties of the ferritic phases (Prop. F).

  forming mostly ferrite and pearlite at the end of the cooling process. The results presentation in the figure follows the same organization as previously: results taking into account the latent heat as an internal heat source in the direct problem solution are presented on the left-hand side, while results with the direct model neglecting the latent heat are shown on the graphs on the right-hand side. The presence of the latent heat is visible in Fig. 10a.I and b.I for 200 s < t < 400 s, where the target heat flux and the target thermocouple measurement become almost constant because the latent heat due to the formation of ferrite and pearlite (Fig. 10c.I and d.I) compensates the heat dissipated at the boundary. The material temperature decreases again after 400 s, and some bainite is formed (especially near the heat exchanging surface) until phase transformations are finished at about t = 850 s.

Figure 10 :

 10 Figure 10: Cooling simulation results with h = 100 Wm -2 K -1 (fast cooling) for 42CrMo4 steel, considering or neglecting the internal heat source due to latent heat of phase transformation (graphs on the left-hand and right-hand sides, respectively): (a) heat flux evolution; (b) temperature evolution at the thermocouple position and rear face; (c) product formation evolution at the thermocouple position and (d) at the rear face; (e) phase fraction profiles at the end of the cooling; (e) material hardness profiles at the end of the cooling. Filled lines (--) are the target, dashed lines (---) are estimated values using properties of the austenitic phase (Prop. A), and dashed-dotted lines (-• -• -) are estimated values using properties of the ferritic phases (Prop. F).

Figure 11 :

 11 Figure 11: Normalized mean error for the estimated heat flux and mean error for the estimated average hardness as a function of the average hardness after the cooling process. The numbers over some data are the heat transfer coefficients used in the simulation in [Wm -2 K -1 ]. The shaded areas divide the plot in ranges of phase fractions that are shown above the graph.

Figure 12 :

 12 Figure 12: Comparison between target and estimated heat fluxes in time for different heat transfer coefficients. The heat fluxes were normalized by the initial target heat flux (φ w,0 = φw,targ(0)), while the time was normalized by the time t N t at which the heat-exchanging surface reaches 100 • C in the direct problem simulation.

  work was performed in the frame of the research project RESEM 2020 managed by the Institut de Recherche Technologique Matériaux Métallurgie Procédés (IRT M2P) and financially supported by the French program Plan d'Investissement d'Avenir (PIA).

In the particular case when austenite undergoes hydrostatic compression stresses, phase transformations are slowed down, because the austenite is stabilized.


Cédric Orelle
email: cedric.orelle@ibcp.fr

Lutz Schmitt
email: lutz.schmitt@hhu.de

Jean-Michel Jault

Keywords: ATP-Binding Cassette, ABC transporters, multidrug efflux pump, ATPase activity, antibiotics

The ATP-Binding Cassette superfamily ABC -transporters (see Glossary) are one of the largest family of nanomachines found in all living cells. They are especially abundant in microorganisms and allow the uptake of a wide diversity of nutrients (e.g. sugars, phosphate, vitamins), or the export of macromolecules, from harmful xenobiotics (e.g. antibiotics or antifungal compounds) to large toxins (e.g. hemolysin A ~ 110 kDa) [1]. Although many ABC transporters display a high specificity for a transported molecule, or a class of molecules (e.g.

the maltose transporter, Figure 1), those having a protective role in tissues or micro-organisms are capable to expel many different molecules with a wide diversity of structures [2,3]. Yet, these poly-specific multidrug ABC exporters can, in an opportunistic way, confer resistance to therapeutic treatments in pathogenic microorganisms or malignant cells; they are therefore a matter of deep concern in public health [2,4]. The ABC transporters share a common architecture with two transmembrane domains (TMD) and two nucleotide-binding domains (NBD), which adopt different topologies ranging from four isolated subunits to a fusion of all four domains (Table 1). The availability of many 3D structures confirms that most ABC transporters function according to an alternating access mechanism, with inward-and outward-facing conformations [5], controlled by ATP binding and hydrolysis [4,[6][7][8][9]. For multidrug ABC exporters, the inward-facing (IF) conformation is prone to bind drugs while the outward-facing (OF) conformation permits drug release outside the cell [10].

Many ABC transporters have a basal ATPase activity

As inferred from their name, ABC transporters are energized by ATP hydrolysis and conserved motifs involved in ATP binding and/or hydrolysis are found in the NBDs. In fact, in the OF conformation, the two NBDs are engaged in a tight interaction with two composite ATPbinding sites located at their interface [2,4] 1 and Figure 1). In such cases, the addition of transported drugs does not (or only barely) increase the ATPase activity of these transporters. Other multidrug ABC exporters from eukaryotic origin, e.g. the human P-glycoprotein (also known as MDR1 or ABCB1) or its homolog from the unicellular red alga Cyanidioschyzon merolae have a moderate basal ATPase activity that can be stimulated several fold by small concentrations of drugs [15][16][17][18] (Table 1 and Figure 1).

Why multidrug ABC exporters have a high basal ATPase activity ?

The high basal ATPase activity of multidrug ABC exporters has been very puzzling for decades and it has been proposed to be caused by the co-purification of an endogenous substrate or the transport of lipids in reconstituted assays [19][20][21]. Although this explanation might be appropriate for multidrug ABC exporters that often bind and/or translocate lipids or detergents (e.g. [20]), this cannot explain the basal ATPase activity found in very specific ABC transporters such as the vitamin B12 importer (Table 1). A thorough thermodynamic study of human P-glycoprotein led to the alternative conclusion that this basal ATPase activity is an intrinsic property of this transporter [22], as initially suggested based on theoretical models , this scavenging effect will compensate for the weak interactions between drugs and the transporter to make sure that most drug molecules are efficiently captured. However, given the low intrinsic affinity of multidrug transporters for drugs, these compounds might be released from the drug-binding site before being expelled by the transporter. Therefore, the corollary of a weak affinity is that multidrug exporters must react immediately upon drug binding. This will prevent the dissociation of the drug from the transporter as this could lead to passive diffusion inside the cell thereby causing irreversible damages (Figure 2, panel B). It is noteworthy that in this model, the ATPase activity turnover number calculated for the most active multidrug ABC exporters is around few ATP hydrolyzed/sec (i.e. few drugs transported/sec).

This activity is not tremendously elevated and one can wonder if this is enough to cope with, and prevent, the penetration of toxic molecules, especially if the number of transporters per cell is limited. It is important to keep in mind that despite the relatively low affinity of these transporters for most of the drugs, these compounds share a hydrophobic nature that should increase their time of residency inside the membrane (i.e. step 5 should be rather slow in Figure 2, as indicated by a dashed arrow). Therefore, this should give enough time to multidrug ABC exporters to catch most, if not all, harmful compounds.

Drugs can help to reach the OF state for some multidrug ABC exporters

Another point worth mentioning is that due to the poly-specificity of multidrug ABC transporters, it is unlikely that the binding of each drug would generate, by itself, a conformational change sufficient to overcome the energy barrier required to reach the OF conformation (Figure 3). Accordingly, 3D structures of multidrug ABC transporters, or the related lipid A exporter MsbA, show that limited conformational modifications occur upon drug (or substrate) binding [8,27,28], and similar conclusions were drawn for multidrug transporters from the Resistance-Nodulation-Division family [29,30] or for soluble multidrug-binding proteins that control the expression of multidrug efflux pumps [31,32]. ATP hydrolysis occurs in the OF conformation for multidrug ABC exporters and, at the most, drugs could shift the equilibrium towards this state [11,[33][34][35]. Indeed, it has been shown by single particle cryoelectron microscopy (cryo-EM) that in the presence of ATP/Mg 2+ , whether or not an inhibitor of the ATPase activity is also added (i.e. vanadate), the mammalian P-glycoprotein remains almost exclusively in the IF conformation, similar to the situation found in the resting state (i.e. in the absence of nucleotides). Drugs need to be supplemented, in addition to ATP/Mg 2+ , to be able to populate the OF conformation [36]. In line with this, the transition of Pglycoprotein from IF to OF conformations requires the presence of both, nucleotides and drugs, indicating that the binding of each effector alone is not sufficient to stabilize the OF state [37,38]. Altogether, these results show that P-glycoprotein prevails in IF conformation under ATP hydrolysis conditions and the presence of drugs will help to switch the equilibrium towards the OF conformation. This is congruent with a stimulatory effect brought by low concentrations of drugs on the basal ATPase activity of this transporter (Figure 3). This IF conformation is presumably stabilized by the presence of several hydrogen bonds near the top of the drugbinding site of eukaryotic P-glycoproteins [25]. On the other hand, in the cases of bacterial homologues such as MsbA, the two states appear equally populated under turnover conditions and the addition of vanadate or the use of an ATPase deficient mutant strongly stabilize the OF conformation [36]. The propensity to reach a stable OF conformation was also demonstrated for the bacterial multidrug exporter, BmrA, by different biophysical approaches including solidstate NMR [39,40], small angle neutron scattering [6], nano-differential scanning fluorimetry [6,41], hydrogen/deuterium exchange coupled to mass spectrometry [6,42], and crystallography and cryo-EM [7]. Therefore, bacterial multidrug exporters appear to easily fluctuate from IF to OF states under turnover conditions, independently of the presence of the drugs. This is presumably because the energy required to reach the OF state is intrinsically lower for microbial multidrug ABC exporters than for their eukaryotic counterparts (Figure 3).

Concluding remarks and future perspectives

Due to this easiness to reach the OF state, microbial multidrug ABC exporters have acquired the ability to hydrolyze ATP at a fast pace, even in the absence of drugs. In fact, in yeast cells expressing major multidrug ABC transporters, specific inhibition of Pdr5, for example, increased cellular ATP levels indicating that these transporters constantly hydrolyze ATP in vivo [43] (see Outstanding Questions). But this seemingly waste of chemical energy give them the ability to promptly translocate any drug captured. Therefore, the high basal ATPase activity of microbial multidrug ABC exporters is most likely the price to pay to make sure that they readily preserve microorganisms from drug penetration and cell death. A similar conclusion could be applied for the basal ATPase activity of other ABC transporters where this could help the transporter to quickly respond to the presence of scarce substrates [44]. An ensuing unsolved question is the coupling efficiency of multidrug transporters (i.e. how many ATP molecules are hydrolyzed per drug transported). It is likely that in the presence of drugs, some transporters hydrolyze ATP in a non-productive manner while others are better coupled to expel drugs.

Studies on that matter are scarce, but a recent work with the tripartite system MacAB-TolC demonstrated that it has a low ATPase activity but it transports substrate with very high coupling efficiency [45]. Another major issue regarding drug transport is how the drug is released from multidrug ABC exporters? A twist and squeeze mechanism has been proposed to allow the direct expulsion of the drug from the P-glycoprotein pump into the external milieu [46]. Alternatively, a floppase mechanism has been proposed for Pdr5 where the drug will be moved by the transporter from the inner-to the outer-leaflet of the membrane [3]. Such a mechanism was also put forward for multidrug-like bacterial transporters, e.g. Sav1866 or

MsbA, based on structural variations found in the OF conformation between these transporters and P-glycoproteins [15].

Future works will certainly enlighten our understanding on how these fascinating transporters work at the molecular level.

Outstanding Questions

What is the energetic burden associated with the basal ATPase activity of ABC transporters for cells? This question relates to both the basal level of ATPase activity of each ABC transporter and its expression level in the cell.

In vivo, can ABC transporters use other energetic sources than ATP (GTP, protons etc.)? It has been shown that PatA/PatB has a strong preference for GTP as the energy source but other pieces of evidences suggest that this might not be an oddity among ABC transporters (e.g. Pdr5 in some instances).

What is the number of ATP hydrolyzed per molecule transported? For ABC transporters with a high basal ATPase activity (e.g. microbial multidrug exporters), this is a very delicate, yet essential, question to address. It is likely that in the presence of drugs, some transporters might hydrolyze ATP in a non-productive manner (i.e. if their drug-binding site is empty) while others would expel drugs. Therefore, the ratio between ATP hydrolyzed and drug transported can be overestimated. To overcome this issue, one might need to use single molecule approaches where both ATPase and drug transport could be monitored at once. What are the physiological substrate(s) of multidrug transporters and to which extent some transporters export drugs in an opportunistic way? In microbes, several putative multidrug exporters often co-exist in a single organism. Therefore, it is likely that some of them exports natural compounds that are presumably chemically related to certain drugs. What is the mechanism that allows the release of the drug? A direct release from a drug exit site has been proposed but the drug can also be delivered in the outer leaflet of the membrane according to a floppase mechanism. Can this mechanism vary according to the origin of the transporter and/or the nature of the drug translocated ?

A B C D

[

  23]. This activity might then appear futile leading to a seemingly useless consumption of ATP in vivo. In contrast, we propose that the high basal ATPase activity confers upon microbial multidrug transporters the ability to readily expel noxious molecules before they penetrate inside the cell(Figure 2, Key Figure). Some observations support this postulate. First, multidrug transporters have been tailored to recognize many different drugs with unrelated structures. This is presumably due to the intrinsic flexibility of the transporter in the IF conformation allowing a plasticity and re-shaping of the drug-binding pocket[24], which most likely involves an aromatic hydrophobic network on top of this site in ABCB1 transporters[25]. This promiscuity of the drug-binding site was presumably acquired at the expense of their affinity for drugs which remains intrinsically low[17, 26]. But because drugs are likely recognized by multidrug transporters directly from the membrane[21], which can concentrate a drug up to 10 4 depending on its lipid-water partition properties[26]

  While some specialized ABC transporters, such as the maltose importer, have an ATPase activity strictly dependent on the presence of the substrate-binding protein and maltose[12], many ABC transporters show a basal ATPase activity in the absence of any added molecule. This basal

	ATPase activity ranges from few dozen of nmoles to few moles ATP hydrolyzed/min/mg
	proteins depending on the ABC transporter [2] (Table 1). A high level of basal ATPase activity
	is a typical feature found in microbial multidrug ABC exporters, such as BmrA from Bacillus
	subtilis [7], PatA/PatB from Streptococcus pneumoniae [13], or Pdr5 from Saccharomyces
	cerevisiae [14] (Table

. This OF state is competent for ATP hydrolysis that will reset the transporter in the IF conformation prone for a new catalytic cycle

[4, 11]

.