Thermo-mechanical and metallurgical model of leaf spring industrial quenching process

. This study is a numerical analysis of the industrial quenching process for leaf springs developed by the CAVEO company. The leaf chosen for this study is of a parabolic profile made of EN-51CrV4 steel (AISI 6150). The aim of this study is to set up a numerical model to predict thermal, metallurgical, and mechanical behavior of a leaf spring from exit of the heating furnace to exit of the quenching bath going through a cambering operation. This study would therefore allow the company to switch from a development scheme based on experiments using physical prototypes tested on the production line to a new scheme based on virtual prototypes using numerical simulation. The development of the numerical model using the finite element method is carried out using the ABAQUS/Implicit solver coupled with two user subroutines Phase and UMAT. The first one have been developed to compute microstructure evolution and the second one to define the constitutive law taking into account phase transformations. This model helps us to follow the spatio-temporal evolutions of temperature and microstructure in the leaf, as well as the variation of the leaf deflection during the process. The proposed numerical model is supported by an experimental protocol based on infrared thermographic images, Rockwell-C hardness measurements, metallographic observations, and deflection measurements. Indeed, the results of the proposed thermo-mechanical and metallurgical model are closed to experimental results.


Introduction
During the manufacturing process, the leaf is heated until its austenitization temperature. After austenitization, it passes through a cambering operation to confer it the required deflection. Then a quenching operation is carried out in an oil bath to obtain a martensitic microstructure which offer better mechanical properties for the leaf.
The main objective of this study is to develop a reliable numerical model to predict the spatio-temporal temperature distribution, the structural-phase state and the deflection of the leaf during the quenching process. This numerical model will therefore allow the company to initially monitor and thus optimize the different parameters of this process, namely leaf austenitization temperature, oil temperature, oil bath holding time, etc. It will then facilitate the development of new leaf springs ranges and thus, resolve prototyping costs as well as customer response time.
This multi-physical problem is solved by considering the interaction between thermal, metallurgical, and mechanical phenomena. Numerous experiments have been carried out, considering structural transformations, to model the quenching (Belhadj et al. 2008), (Moumni et al. 2011) and (Esfahani et al. 2021). Two methods are used to model structural transformations, the first relates to transformations with diffusion depending on both temperature and time and the second relates to transformations without diffusion depending only on temperature. The works of (Avrami 1939), (Avrami 1940) and (Koisstinen and Marburger 1959) are the most used.
Several studies on the numerical simulation of various industrial hot forming processes followed by quenching of mechanical components (Karbasian and Tekkaya, 2010), (Pipard et al. 2013), and automotive components (Eşiyok et al. 2016), (Xuejun, et al., 2020) have been conducted. Most of these works are carried out on high elastic limit steels by integrating hot stamping and quenching operations in a single step (Naderi et al. 2008), (Xing et al. 2009), (Taylor et al. 2014) and (Merklein et al. 2016). The proposed model relates to a particular manufacturing process whose parameters differ from the published studies in terms of geometry, material, forming and quenching conditions. This article presents a numerical model of an industrial leaf quenching process. This multi-physical model takes into consideration the effect of quenching on the structural phase state of the leaf as well as their influences on the elasto-viscoplastic behavior.
This paper is structured as follows: the first part is dedicated to an experimental investigation of cambering and quenching operations to define its parameters. The leaf is subjected to temperature measurements, Rockwell-C hardness measurements, micrographic observations and deflection measurements. The results obtained from this experimental study will be used as input variables to the numerical model and to validate its results. The second part concerns the development of a numerical model using the finite element method on ABAQUS/Implicit software by means of user subroutines UMAT and Phase.

Experimental study
The leaf is heated until the austenization temperature, followed by a fast oil cooling. A hot cambering operation on a hydraulic press between heating and cooling is performed on the leaf, to give it the expected deflection. Two robots perform transfer operations of the leaf between furnace, press and oil bath.
The leaf chosen for this study has a parabolic profile ( Figure .1), measuring 1300 mm in length and 70 mm in width. Its thickness in zone 1 is 16 mm and ranges from 16 mm to 9 mm in zone 2. The leaf presents two symmetry planes (P1 (O, y, z) and P2 (O, x, y)).

Figure1: Geometry of a parabolic leaf before and after cambering
The leaf is made of EN-51CrV4 steel (AISI-6150). Its chemical composition follows standard EN 10089 (Table  1). The experimental investigation concerns measurements and observations, which are carried out, either on the production line or in the laboratory on samples taken from the leaf.

Infrared temperature measurements.
A calibrated infrared thermography camera takes the temperature measurements at three times of the process o At time t0: at the exit of the leaf from the heating furnace. o At time t1: after cambering. o At time t2: at the exit of the leaf from the oil bath.
Temperature evolutions obtained from infrared images (Figure 2a, 3a and 4a) by drawing lines (Li1 and Li2) on the surface of leafs. Leafs are transferred in pairs from heating furnace to quenching bath.
• At time t0: The camera displays an almost uniform temperature of 900°C at leafs surfaces (Fig 2b). This measured temperature shows that at t0 leafs are in the austenitic range according to the CCT diagram of EN 51CrV4 steel ( Figure 5).
• At time t1: After cambering the temperature displayed by the camera (figure 3) ranges between 780°C and 855°C in zone 2. This temperature drop is mainly due to the contact between the press and the leaf. The thermal losses undergone by the leaf until t1 does not cause a diffusive transformation of austenite according to the CCT diagram of EN 51CrV4 steel ( Figure 5). The leaf is therefore ready to the quenching operation.
• At time t2: After quenching in the oil bath at 53°C, figure 4 shows that the temperature in zone 1 is about 70°C and varies between 55°C and 70 °C in zone 2.

Rockwell-C Hardness measurements
At the end of the process when leaf reaches room temperature, Rockwell-C Hardness (HRC) measurements are performed on the surface along the leaf. The profile of Rockwell-C Hardness is shown in Figure 6 and it varies appreciably around 60 HRC.
Referring to the CCT diagram of the steel ( Figure 5), this Rockwell-C hardness corresponds to a cooling rate higher than the critical quenching rate. It suggests that austenite has been totally transformed into martensite.

Optical micrograph
Microstructural observations were performed on samples taken from three zones of the leaf: the central axis, the extremity, and between them. These samples reveal the same lath martensite structure in the optical microscope ( Figure 7). This morphology is in accordance with the Rockwell-C hardness measurements.

Deflection of the leaf
The deflection conferred by the cambering operation (figure 8) is modified during this industrial quenching process. Indeed the deflection measurements show a variation from 154 ± 3 mm during cambering to 130 ± after quenching which corresponds to a spring back of 24 mm. All experimental results will serve as input variables for the numerical simulation in first time and to validate it in second time.
Title of your paper 7

Thermo-mechanical and metallurgical model
Several physical phenomena, thermal, metallurgical, and mechanical are involved in the leaf quenching process. It is a multi-physical problem generating couplings between these different phenomena ( Figure 9). Thermo-metallurgical coupling reflects the effect of thermal on metallurgical phenomenon and vice versa. Depending on the cooling rate, the EN-51CrV4 steel initially in the austenitic state can give rise to various structures during the cooling of the leaf. These phase transformations generate latent heat that in turn vary the thermal history of the leaf. For this study, we only consider the transformation latent heat of austenite to martensite through the heat transfer coefficient during quenching that was determined in previous work (Slama et al. 2018). Indeed, this coefficient is obtained by quenching, under the same conditions as the leaf, a standard probe made of the same steel (EN-51CrV4).
Thermo-mechanical coupling defines the effects of mechanical on thermal phenomenon and vice versa. The quenching process generates thermal gradients within the leaf that cause strain and stress fields through thermal expansion as well as the variation of the mechanical properties of the steel. Moreover plastic deformation generates a heat source. Compared to convection and conduction heat flows, this heat source was neglected in this study.
Metallurgical-mechanical coupling presents the effects of mechanical on metallurgical phenomenon and vice versa. Transformation plasticity shows that stress state affects the phase transformation kinetics. Conversely microstructure transformations result in volume change that cause deformations. In addition the mechanical properties depend on the microstructure phases of the leaf. The transformation plasticity is not considered in this work.

Thermal model
During this quenching process, the thermal problem is governed by the following nonlinear transient heat equation:

Names of the authors
Where : : Density in Kg/m3, ! : Specific heat in J/Kg.K, : Temperature in °C, : Thermal conductivity in W/m.K, ̇: The rate of heat generation term W/m 3 , Equation (1) is solved by taking into account the following boundary conditions: ü Convection and radiation between the leaf and the environment: Where • ℎ conv : Convective heat transfer coefficient between the leaf and the environment.
• ℎ ray : Radiative heat transfer coefficient between the leaf and the environment ü Convection between the leaf and the quenching oil : Where • h oil : Convective heat transfer coefficient between the leaf and the quenching oil.
• oil : Measured temperature of the oil bath.
ü Conduction between the leaf and the press : Where • h cond : Conductive heat transfer coefficient between the leaf and the press. • S1 et S2 : Temperature contact surfaces.

Metallurgical model
Title of your paper 9 During the quenching process the cooling of the leaf can lead to the transformation of the austenitic structure into ferrite-pearlite, bainite or martensite. The two models considered for the diffusive and non-diffusive transformations are detailed in the following paragraphs For the diffusive transformation, it consists of the transformation of austenite into ferrite, pearlite and bainite. Many models, describing phase transformations during continuous cooling, rely on isothermal kinetics using the additive rule initiated by Scheil (Scheil 1935). This rule represents the anisothermal transformation as a succession of isothermal transformations over very short time intervals. Pumpherey and Jones (Pumphrey and Jones 1948) proposed a model based on this additive rule, which assumes that the final microstructure is the sum of elementary isothermal transformations, which come from the discretization of the thermal cycle. This model ensures the transformation continuity between different levels by introducing a fictitious time. This time represents the time required to obtain the ;<= phase proportion if the transformation is isothermal at the Ti temperature. Thus, the diffusive transformation is modelled by Avrami (Avrami 1939) using the fictitious times as follows: Where ; * is the fictitious time associated with ;<= : Where and are the phase proportions at step i and i-1.
The parameters ; ( ) and ; ( ) are obtained by solving the following system of equation (Shah et al. 2012) : Where : • ts: Start time of transformation, estimated at 1%.
• te: End time of transformation, estimated at 99%.

Names of the authors
This work is performed within the framework of the hypo-elastic formulation (De Souza et al. 2011) of large deformation materials. A Von Mises plasticity criterion is assumed. The hypo-elastoplastic assumption consists of an additive decomposition of the total strain into an elastic ̇e l , a visco-plastic ̇, a thermal ̇t h and a microstructure transformation plasticity ̇p t parts. In this work, the last term is postponed ̇p t for a later study.
Thermal strain rate ̇t h The thermal deformations of the different phases are calculated according to a linear mixing law as follows: ̇= ∑ ; ;; With T is the temperature, I is the identity tensor and ; is the thermal expansion coefficient of the i th phase.

Viscoplastic strain rate ̇
The viscoplastic strain rate computation is based on the method of elastic prediction and plastic correction (De Souza et al. 2011). In this step, the plasticity criterion based on a variation of the stress tensor is determined by assuming that the strain increase is strictly thermo-elastic (no visco-plastic strain): That is, ̇= . A test stress tensor is calculated by the following equation: Where the exponents 0 and 1, respectively indicate the previous instant and the present instant. Elastic stress tensor used to verify the plasticity criterion is written as a function of the deviatoric stress tensor and the yield stress as follows: Where ; and T; are respectively the phase proportion and the yield stress of the i th phase with i in {austenite, ferrite, bainite, martensite}. The yield stress T; of the phase "i" is governed by a multiplicative viscous hardening law: T; = n T; *̂!1. ; *̂̇!1. ; ( ) Where n T; (̂!) is the quasi-static yield stress, ; *̂̇!1 is a viscous term introducing the effect of the plastic strain rate and ; ( ) is a term introducing the temperature effect. Johnson-Cook constitutive law is adopted in this study, which is widely used for hot forming processes simulation (Nie et al. 2016), (Bin et al. 2017). Then equation (16) terms is calculated as follows n T; = * ; + ; ̂! ? ! 1 ; *̂̇!1 = *1 + ; ln̂̇! ; * 1 (18) Where, ; , ; , ; , ; and ; are parameters of i th phase and : ̂! is the equivalent plastic deformation rate.
• ̂Ȯ ; and 285 ; are respectively, the plastic deformation rate and the reference temperature.
Von Mises Yield criterion value define the behavior of the material: • If ( , T ) ≤ 0 then ̇= and thus = • If * , T 1 > 0 then the visco-plastic deformation rate ̇ root of * , T 1 is obtained as follows: With is the normal to the plasticity surface which is given by: The value of ̂! is numerically calculated by the return mapping algorithm (Simo and Hughes, 1987). The calculation scheme is explained in section 4.5 dedicated to the numerical implementation of the model.

Finite element simulation
Numerical simulations are carried out using the finite element code Abaqus/Standard coupled with user subroutines Phase and UMAT.

Geometry and mesh
The geometrical model relies on three components ( Figure. 10): • The leaf: deformable 3D body.
• Upper and lower jaws of the press: isothermal rigid shell.
• Upper and lower fingers of the press: isothermal rigid shell.
This leaf presents two symmetry plans that allow us to reduce the geometrical model to the quarter in order to minimize computation time.

Figure 10: Geometrical model
The leaf is meshed using C3D8T Abaqus element. It is a fully integrated 8 node hexahedral linear element with a temperature degree of freedom. The leaf mesh is constituted of 3900 elements and 5252 nodes. The press components were meshed by rigid shell elements. The press mesh has been refined at contact areas.
• The parameters of Johnson Cook's law A, B, n, C, and m for the different phases are taken from (Biermann et al. 2010) and (Hor et al. 2013).

Initial and boundary conditions
Ø Initial conditions are as follows: ü Initial temperature of the leaf is homogeneous and fixed at 900 °C (Average value of experimental measurements). ü Initial microstructure of the leaf is considered homogeneous, ; 643$ = 1.
ü Temperature of the press is set at 25 ° C. ü Temperature of quenching oil is 53 ° C (Experimental measurements).

Ø Boundary and interaction conditions :
ü Radiation of the leaf towards the environment: ɛ (emissivity) = 0.85 (Bergman et al. 2017). ü Convection between the leaf and the environment: hconv = 7 W/m 2 .K (Naderi et al. 2008). ü Conduction between the press and the leaf: hcond = 1200W / m2.K (Merklein et al. 2009). Many works have been interested in the identification of this heat transfer coefficient (hcond) between parts and forming tools (Karbasian et al. 2008) and (Merklein et al. 2009). These studies show that hcond coefficient depends on temperature and pressure. Other works demonstrate that this coefficient is mainly dependent on pressure (Tekkaya et al. 2007) and (Chang et al. 2016).
14 Names of the authors The work of (Merklein et al. 2009) determine the evolution of this coefficient as a function of the contact pressure. In this work, we used hcond = 1200W / m2.K corresponding to a pressure of 40 bar.

Implementation of the metallurgical models
A phase subroutine recovers the leaf temperature evolution during its cooling then calculates phase proportions according to ( Figure 13) and communicates them to Abaqus through internal variables.

Implementation of the mechanical behaviour
The implementation of the constitutive law in the Abaqus code goes through the computation of the tangent stiffness fourth order tensor defined by : The second order stress tensor can be decomposed as follows: Where s is the deviatoric part and p is the hydrostatic pressure: Then equation (24) can be decomposed as follows The term # vol is as follows (De Souza et al. 2011) : Where = is the compressibility modulus (indices 0 and 1 respectively indicate the previous and the current instant respectively) which is written as a function of the Young modulus and Poisson ratio as follows: We derive the expression of the deviatoric stress tensor in relation to the tensor of the strains to obtain the deviatoric term. The calculation's result gives: 16 Names of the authors Where ijkl = = 9 * ik I jl + I il I jk 1 is the shear modulus of the material at the current time.
Where ̂+ is the increment of the equivalent plastic strain between the current time 1 and the previous time 0.
Where H is the strain-hardening rate calculated from Johnson Cook's law as follows: Where and Given the nonlinearity of the equations of this law, the method followed for its numerical implementation is based on the elastic prediction and plastic correction algorithm, which is often used to solve plasticity problems.
The first step, which is the elastic prediction, consists of determining whether there will be a plastic deformation. The plastic strain increment will be determined by the second step, which is the plastic correction, and thus the stress tensor value obtained in the elastic prediction stage will be corrected. The numerical algorithm (Table 2) describes the approach implemented in the UMAT subroutine called at each time increment and convergence iteration for each integration point. Variables value at the previous time and at the current time are indicated respectively by the indices n and n+1 (Table 2). Table 2. Elastic predictor plastic corrector algorithm
Compute strain tensor in the current referential of the element = with is the rotational tensor The return-mapping algorithm (Table 3) used to determine the new strain tensor and the equivalent plastic strain rate if plasticity occurs. Then, the plastic strain and the new stress tensors are deduced.
Compute plasticity criterion

Temperature distribution
Temperature distributions are presented below at times t1 and t2: • At time t1, after cambering when the upper fingers of the press are raised, figure 14 shows a temperature decrease mainly due to the contact between the leaf and the press. According to this figure the temperature varies between 800 and 855 °C in Zone2. These thermal losses are in accordance with those measured. Indeed the infrared camera displays at instant t1 a temperature that varies between 780 and 855 in Zone2. • At time t2, after quenching in the oil bath, the leaf temperature is almost homogenous showing 60°C in the extremities and 65°C in Zone 1 ( figure 15). These numerical results are consistent with experimental measurements were temperature ranges between 55 and 70°C.

Phase transformations
The metallurgical behavior simulation predicts the spatio-temporal evolution of phase proportions in the leaf during the quenching process. These phase proportions are given at three cooling times: After cambering at time t1, after quenching at time t2 and at time t12 between them (during the quenching).
ü At time t12 during the cooling of the leaf in the oil bath, austenite gradually turns into martensite (Figure. 17). Ferrite-pearlite (FV2) and bainite (FV3) do not appear in these calculations.
ü At time t2, when the leaf exit the quenching bath, numerical simulation shows in figure 18 almost 91% of martensite and 9% of untransformed austenite.  The hardness measurements which are 60 HRc all over the leaf, corresponds to a martensitic structure according to CCT diagram of the EN-51CrV4 steel. This hardness value can only be reached if the leaf structure at time t1 is totally austenitic (Figure 16). In addition martensite structure is revealed by micrographic observation (Figure 7).
At this level of comparison, numerical and experimental studies complement each other and assure a first degree of validity of the thermos-metallurgical model. Therefore, this thermo-metallurgical model is able to predict for each geometry the temperature history and the microstructure evolution during this industrial quenching process.

Mechanical results
This work focuses on the prediction of the leaf deflection at the end of the quenching process in order to estimate the spring back caused by the cambering and the quenching operations. This spring back is illustrated through curves of figure 19.
During the cambering operation, the leaf follows the shape of the fingers of the press, previously arranged, to give it a deflection of 154 mm (Y-cambering Figure 19). After quenching the deflection becomes 128mm (Y-quenching) so the spring back is rated at 26mm. This value is closed to the measured one which is 24mm.

Conclusion
In this study, we modeled and simulated the thermal, metallurgical and mechanical behavior of a parabolic leaf during its manufacturing throughout cambering and quenching operations. This model allowed us to predict the space-time distributions of the temperature, the phase proportions as well as the spring back of the leaf.
To validate the numerical model, temperature measurements as well as deflection measurements were carried out on the production line. Hardness measurements and metallographic observations were then performed on samples from the leaf.
After cambering and quenching, the proposed model shows that temperatures in Zone 1 and at the extremities are respectively of 65°C and 60 °C. These results are closed to those measured which are respectively 70°C and 55°C. In addition, Hardness and metallographic observations are in accordance with phase proportions delivered by the numerical model which is a mainly martensitic structure. Mechanical results in terms of deflection are also satisfactory since it gives a spring back of 28 mm against 26mm measured. Hence, the proposed numerical model give us results close to experimental once. Thus this model presents a reliable and fast tool to ensure optimal quenching process parameters and a continuous control of the process. This will allow the company to develop rapidly new spring ranges at a lower cost. Data availability All the data have been presented in the manuscript.

Declarations
Competing interests The authors declare that they have no conflict of interests.
Ethical approval Not applicable.

Consent to participate
The authors declare that they all consent to participate this research.

Consent to publish
The authors declare that they all consent to publishthe manuscript. Figure 1 Geometry of a parabolic leaf before and after cambering Figure 2 Leafs temperature at time t0.

Figure 3
Leafs temperature at time t1.

Figure 4
Leaf temperature at time t1. Rockwell-C hardness pro le on mid-length of the leaf Figure 7 Optical micrograph of the leaf after quenching Coupling of physical phenomena during leaf quenching process  Leaf phase proportions at time t1 (FV1=ξ_aust).

Figure 17
Leaf phase proportions at time t12.

Figure 18
Leaf phase proportions at time t2. Figure 19