Advances in the numerical investigation of the immersion quenching process

A numerical investigation of the immersion quenching process is presented in this paper. Immersion quenching is recognized as one of the common ways to achieve the desirable microstructure, and to improve the mechanical properties after thermal treatment. Furthermore it is important to prevent distortion and cracking of the cast parts. Accurate prediction of all three boiling regimes and the heat transfer inside the structure during quenching are important to finally evaluate the residual stresses and deformations of thermally treated parts. Numerical details focus on the handling of the enthalpy with variable specific heat capacity in the solid. For two application cases, comparison between measured and simulated temperatures at different monitoring positions shows very good agreement. The study demonstrates the capability of the present model to overcome the numerical challenges occurring during immersion quenching and it is capable of predicting the complex physics with good accuracy.


Introduction
To satisfy ecological requirements, modern automotive industry is using various means to make vehicles become more efficient and more ecological. To further reduce fuel consumption and emissions, heavier metals have been increasingly replaced by those high strength-to-weight ratio and low-cost aluminium alloys [1]. Temperature sensitive microstructural defects can lead to early failure of working parts at the service temperature. Furthermore, Residual stress created during heat treatment can lead to stress corrosion cracking, fatigue cracking and distortion in components. Among many heat treatment techniques, optimized immersion quenching process has been extensively applied to improve the mechanical properties and to obtain the desirable microstructure, mechanical and metallurgical characteristics of final metallic alloy components during production. In addition, immersion quenching could produce a relatively even temperature distribution during the rapid cooling process, thus reduces the residual stress level and consequently lowers the risk of distortion and cracking of the cast parts [2].
In the immersion quenching process, rapid heat transfer happens immediately after the metal piece with a microstructure desirable temperature is dipped into the sub-cooled liquid pool. Water, oil or polymers are widely used as liquid quenchant [3]. In most cases, this microstructure dependent temperature is much higher than the normal Leidenfrost point, so film boiling, transition boiling and nucleate boiling are found during the immersion quenching cooling process [4]. Optimal  uniformity is essential for distortion control and prevention of cracking problems. This demands the heat transfer between the quenchant and hot workpiece as uniform as possible throughout all three boiling regimes. Although experimental researches [5,6] have showed invaluable insight into the fluid mechanics of the quenching process, they are generally incapable of continuously monitoring fluid flow in commercial quench tanks during heat treatment process. CFD modeling [7][8][9][10][11][12] is being extensively used to investigate fluid flow and heat transfer in quenching process and has clearly demonstrated the enormous prospect of promoting quenching performance with the objective of retrofitting for improved performance or designing new and innovative quench systems. Detailed description of the immersion quenching process and numerical aspects are described in the publications of Kopun et al. [9][10], Greif et al. [11] and Srinivasan et al. [12]. Previous works [9][10][11][12][13] focus mainly on constant thermal properties, while in many cases, thermal properties of the solid work pieces, i.e. specific heat capacity, C p and thermal conductivity  depend on temperature. The specific heat capacity, C p is particularly a very critical parameter during the numerical simulation of the heat transfer. When the specific heat capacity C p is constant over a temperature difference of T, the static enthalpy change can be easily taken as h = C p T; whereas when C p is temperature dependent, the static enthalpy change is not C p T. In this work, a correct treatment of the spatially and temporally changing C p is developed and a more general expression for temperature dependence of thermal enthalpy change is presented.
The paper is organized as follows. Section 2 presents the mathematical model. Simulation set-up and details of two different industrial quenching processes are briefly introduced in section 3. The comparisons of the simulated results from two industrial quenching processes with the experimental data are presented and discussed in section 4. Conclusions and recommendations are given in section 5.

Mathematical model
In the framework of the Eulerian multi-fluid model, each phase is taken as interpenetrating continua coexisting in the flow domain with inter-phase transfer terms accounting for phase interactions. Twosets of ensemble averaged balance equations for mass, momentum and energy [14] take the following form:

Continuity equation
The compatibility condition is exerted on the continuity equation: where the indices k refers to the phase (l for liquid and v for vapour). The volume fraction of each phase is denoted as stands for density and is the velocity vector. The inter-phase mass change rate (in this particular case, boiling) is  k .

Momentum equation
where p is pressure; is gravity vector,  is stress tensor and interfacial velocity is denoted as .
where the terms on the right hand side represents forces due to drag, F D , wall lubrication, F wl and lift F L .

Drag force.
The drag force of the liquid phase is given as where C D stands for the drag coefficient, and d b is the vapour bubble diameter.

Lift force.
Lift force FL exerted on the liquid phase is described as: where C L is the lift coefficient.

Wall lubrication force.
Wall lubrication force assures the experimentally found zero void condition near vertical walls and it is described by where w n  is the unit wall normal vector pointing from the wall into the fluid, and C wl stands for the wall lubrication coefficient. For the above mentioned momentum interfacial forces, detailed models are referred to Kopun [9][10], Tomiyama [15] and Frank et al. [16].

Energy equation
In our previous works, vapour bubble condensation in fluid domain is assumed to be negligible, the phasic static enthalpy equation is then described by Thermal conductivity is presented as  T , the inter-phase heat transfer coefficient between liquid phase and vapour phase is denoted as C H . For simplicity, other terms are lumped into the source term

Boiling model
Based on the previous works of Kopun [9,10], the inter-phase mass change rate  k due to boiling is evaluated as: where C m is the closure coefficient due to uncertainty, C boil stands for boiling correction coefficient and and H fg denotes the vaporization latent heat. h boil is the boiling heat transfer coefficient. Clearly, where vapor bubble diameter d b stands for the length scale, the vapor thermal conductivity is denoted as  v and C p,v is the vapour specific heat.

Transition boiling mode.
In transition boiling regime, heat transfer coefficient, h boil is calculated as: where the transition boiling heat flux q TB is given as follows In Eq. (13), q CHF denotes the Critical Heat Flux, and q MHF stands for the Minimum Heat Flux. The critical heat flux is modeled as where  is the liquid phase surface tension. The Minimum Heat Flux correlation for the quenching system in Eq. (13) is K MHF serves as a model parameter and its typical value varies from 1 to 15. Further details about AVL quenching model can be obtained from AVL FIRE ® Multi-fluid model solver theory guide [18] and/or previous works [9][10][11][12][13].

Variable specific heat capacity
Specific heat capacity, C p is used to define enthalpy and to re-calculate temperature from enthalpy. When C p varies with temperature, C p becomes a critical parameter for the numerical simulation of the quenching process, in which temperature changes hundreds of Kelvin. In this case, the original definition of static enthalpy has to be used, During numerical simulation, the temperature dependent specific heat capacity is usually provided within a specified application range, [T 0 , T u ]. T 0 is the lower bound where Cp(T) begins to be valid and T u is the upper bound. Represented as polynomial function the specific heat capacity can be modeled by the following equation, For one aluminum work piece the polynomial coefficients are a 0 = -100.3, a 1 = 0.07751, a 2 = 0.0065, a 3 = -4×10 -6 , T 0 = 283.15 and T u = 1100 K. In practice, T 0 is lower than the initial quenchant pool temperature and T u is higher than the initial hot solid piece temperature. As from Eq. (16), enthalpy is defined from a reference temperature of 0 K, so specific heat capacity C p is needed in the range of 0 K till highest temperature of the hot piece. Figure 1a shows how the specific heat capacity varies with temperature if Eq. (17) is applied in the range of 0<T<1200. Clearly it is found in Figure  1a that negative C p appears, and this will introduce errors in the simulation and spoils the simulation results. Hence, a special treatment for variable C p is needed. Therefore Eq. (16) is re-written as where S(T) is an unknown function. Since the solid temperature change during quenching takes place within T 0 and T u , the numerical simulation focuses only on this interested temperature range of T 0 < T<T u . Furthermore, only the difference pf the enthalpy is needed during the simulation, so, an exact expression of S(T) in Eq. (18) is not strictly demanded. To keep the simulation procedure simple, here we set it as constant With the help of Eq. (18), temperature dependent C p,new in the hot solid part is defined as

Thermal homogeneous assumption
In previous works, the dispersed bubbly phase and the continuous phase are assumed to be in thermal equilibrium and a homogeneous enthalpy equation was solved. Rather to solve this homogeneous enthalpy equation, in this work, both enthalpy equations of bubble phase and continuous phase are solved. In order to achieve a thermal homogeneous field, a very high value is set as inter-phase heat transfer coefficient, C H in Eq. (8).

Boundary condition and simulation set up
To evaluate the performance of the above mentioned simulation methodology, two immersion quenching processes, a simple step plate geometry and a real case aluminum cylinder head, are nummerically simulated here.
As its names says, the thickness of the step plate work piece varies along its length, and it is displayed in Figure 2a. Temperature history of three different positions, named as T 1 , T 2 and T 3 in the middle of the test piece were measured in the experiments.
The aluminium cylinder head of the second test case and its 9 measuring points are illustrated in Figure 2b. Figure 3a displays the polyhedral computational mesh and the boundary conditions of the step plate test case. The total computational domain is subdivided into 245830 grid cell, in which the solid step plate takes approximately 61000 cells. Fixed atmospheric pressure (1 atm) is applied at the outlet boundary on the top (light red plane), and a velocity boundary condition is exerted on the inlet boundary (at the bottom, light pink plane) to model the dipping process. The dipping velocity is used to determine the inlet velocity rising up from the initial water level (purple part) toward top. In this simulations, inlet velocity at the first 0.5 seconds takes a value of 0, which implies that the step plate is air cooled for the first 0.5 seconds, and then a constant velocity of 0.14 m/s lasts for 2 seconds, so that the final submerging depth has been reached; then inlet velocity is 0 again. The domain walls at the sides were treated adiabatic. Initial water pool temperature is 353.15K and step plate has an initial temperature of 774.15K.

Step plate
The step plate test piece was immersed into a water pool with bulk temperatures of 353 K. The alloy's specific heat capacity, C p varies from 884.8 J/(Kg.K) at 300 K to 1427.3 J/(Kg.K) at 1100 K. Leidenfrost temperature ranges from 560 K to 650 K, detailed variable Leidenfrost temperature model is referred to Kopun [10]. Different to the constant C p used in the previous works [10,13], the variable specific heat capacity as described Eq. (20) is used in this work. Figure 4 illustrates continuous phase volume fraction distribution at different times. At the initial stage, i.e., from 1 to 3s, film boiling occurs as it is seen that water volume fraction surrounding the solid is still higher than 0.85. With time proceeding, more vapour bubbles are generated and travelling to the top due to buoyancy, which implies rapid transit boiling or nucleate boiling appears between 3 and 15 s. At 20 s, boiling happens only at the upper part of the work piece; whereas no bubbles are observed in the lower part of the plate, indicating that in this region, the solid surface temperature is already below the saturation temperature. It is noted in Figure 4 that, during quenching, in the lower part of the work piece, locally generated vapour bubbles escape the surface and rise up; in the upper part of the solid workpiece (anti-gravity direction), locally generated vapour bubbles and the upward travelling bubbles serve like a vapour film, and this blanket effect leads to the longer staying of film boiling regime in the upper part of the work piece. Therefore, during immersion quenching, Leidenfrost temperature increases in gravity direction.
Temperature snapshots in thickness and width directions are presented in Figure 5. It can be seen that intensive cooling is achieved, and across the entire plane, a highly non-uniform temperature distribution is established after only 3 seconds. The temperature drops substantially in the thinner part of the solid. After 20 seconds, the upper part, as it is thicker and heat transfer is still taking place. After 25 seconds, the metal is cooled down to 372K and natural convection begins, so, till 30 seconds, the temperature further drops to 367 K. Figure 6 provides a quantitative comparison of the simulated and measured cooling curves at three monitoring points. In general, excellent agreement is achieved. It is obvious from the measured cooling curves in Figure 6 that, the Leidenfrost temperature at T3 is almost 740K and it is about 610 K at T1. As shown in Figure 2, monitoring point T3 locates in the lower tip of the plate and T1 is on the upper part of the plate, this again confirms that Leidenfrost temperature increases in the gravity direction. At monitoring point T3, the nucleate boiling regime is reached very early, almost immediately after immersion, and therefore it loses heat very fast; together with the fact that its 7 5th Global Conference on Materials Science and Engineering IOP Publishing IOP Conf. Series: Materials Science and Engineering 164 (2017) 012004 doi:10.1088/1757-899X/164/1/012004 thickness is thinner, it reaches the surrounding bulk water temperature very soon. The discrepancy between the measurements and simulated cooling curve at monitoring point T3 may come from the different Leidenfrost temperature. The simulated cooling curve first experiences a film boiling and then enters nucleate boiling so there is a lag between the measurements and simulated curve. The difference between the experiment and simulation at T2 could be partly contribute to measurement errors or the Leidenfrost temperature is not accurately predicted in the simulation. At monitoring point T1, simulated cooling curve perfectly fits with the measured ones and because it experiences longer film boiling and the workpiece thickness at this point is 16 mm, consequently, after about 25 seconds it reaches the surrounding bulk liquid temperature. t= 1.0 s t=3.0 s t=5.0 s t=10 s t=15 s t=20 s

Aluminium cylinder head
Immersion quenching of the current cylinder head was first numerically investigated by Kopun et al. [2]. The aluminium cylinder head is immersed into a water pool with a bulk temperatures of 353 K. The alloy's specific heat capacity, C p varies from 880 J/(Kg.K) at 293.15 K to 1166 J/(Kg.K) at 775.15 K. Different to the work of Kopun et al. [2], where also temperature dependent specific heat capacity and thermal conductivity were used, in the present work the enhanced approach with the original definition of the static enthalpy, i.e. Eq. (16) is applied. Consequently, the specific heat capacity C p is set to a constant value of 880 J/(Kg.K) when temperature is below 293 K and then varies from 880 J/(Kg.K) to about 1166 J/(Kg.K)). Figure 7 shows the water (continuous phase) volume fraction in the mid-plane at different times. Immediately after the hot cylinder heat touches the sub-cooled water, boiling begins and a vapour film is formed surrounding the cylinder head. With time preceding, boiling becomes stronger and stronger. t = 1s t = 2s t = 5s t = 10s t = 20s t = 45s and therefore the heat transfer is weak during this period, and this behaviour should be avoided. Furthermore, it can be seen in Figure 7 that gas is captured in the hollow rooms. Gas occupies the upper part of the hollow room and cannot escape from these hollow rooms. After rapid boiling comes (at time = 5s), water is consumed and more vapour bubbles are generated and these vapour bubbles and air almost fill these hollow rooms, this significantly reduces the heat transfer. After 45 seconds, hardly no bubbles are generated. Figure 8 shows the cylinder head temperature distribution in the mid-plane at different times. It is obvious that the temperature at the open corners drops faster since boiling happens there at all sides. As the vapour bubbles accumulate along the bottom of the cylinder head, and this vapour film prevents rapid boiling, the energy removal in the bottom region is weaker and so temperature drops slower. Again, due to vapour occupation in the hollow rooms, boiling only happens in the lower part of the hollow rooms and it removes heat there, this leads to much higher temperature in the upper part of the hollow room and greater temperature gradient. This is not desired during quenching as it produces higher thermal stress.
The quantitative comparison of the simulated and measured cooling curves at the twelve monitoring points is presented in Figure 9. Except at monitoring points 1, 6 and 9, the simulated cooling curves agree very well with the measured ones. At MP1 and MP6, the difference between simulation and experiments is related to the changing Leidenfrost temperature. The predicted Leidenfrost points are lower than the measured ones, and therefore longer film boiling is experienced at these two points in the simulation which leads to lagging to the measurements about 4 seconds. Furthermore, we had assumed a uniform initial temperature of 495 ˚C for the whole domain, but the measured one at MP9 is 510 ˚C. Consequently difference appears for this measurement point. It is further found in Figure 9 that, rapid boiling happens between 10 and 25 s and after 40 seconds, the cylinder head temperature is cooled down to the water saturation temperature, which is also observed in Figure 8.

Conclusions
A numerical investigation of the heat and mass transfer during the immersion quenching process of two geometries is performed in the frame of Multi-fluid model using the CFD code AVL FIRE ® . The treatment of the temperature dependent specific heat capacity is described in detail. Important physical phenomenon, e.g., different boiling regimes, vapour blankets and vapour accumulation in relatively closed region have been correctly predicted. For the immersion quenching of step plate, accumulation of vapour film in the upper region (anti-gravity direction) not only hinders heat transfer but also decreases the Leidenfrost point, thicker part downwards orientation quenching would be recommended to achieve a faster cooling. For the quenching of cylinder head, with current orientation, trapping of the vapour phase in the hollow room aggravates boiling and leads to highly non-uniformity of temperature gradient in hollow room's walls, immersion quenching should be performed in a way that vapour phase could easily escape from the hollow rooms. For almost all measurement points a very good agreement between numerical and experimental data have been achieved for both cases. It could be concluded that the robust improved CFD methodology provides a solid basis for detailed understanding and exploring into the essence of the quenching process in the workshop by utilizing CFD analysis. This methodology can be used to analyze fluid flow and heat transfer in immersion quenching processes with the objective of optimizing quenching performance for distortion control and reduction of cracking problems.