Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 © 2014 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2014.1705 Original Scientific Paper Received for review: 2014-01-29 Received revised form: 2014-05-09 Accepted for publication: 2014-05-21 Numerical Investigations of Quenching Cooling Processes for Different Cast Aluminum Parts Rok Kopun1* - Leopold Škerget2 - Matjaž Hriberšek2 - Dongsheng Zhang3 - Wilfried Edelbauer3 i AVL - AST d.o.o., Slovenia 2 University of Maribor, Faculty of Mechanical Engineering, Slovenia 3 AVL List GmBH, Austria This paper discusses a recently improved computational fluid dynamics (CFD) methodology for virtual experimental investigation of the heat treatment for cast aluminium parts. The immersion quenching process of the heated work piece in a sub-cooled liquid pool is handled by employing the Eulerian multi-fluid modeling approach, which is implemented within the commercial CFD code AVL FIRE®. The applied heat and mass transfer rate is modeled based on a different boiling regime, which is controlled by the Leidenfrost temperature. The objective of the presented research is to present an updated quenching model by applying variable Leidenfrost temperatures . Furthermore, simulation results are compared with available measurements for a wide variety of quenching scenarios involving immersion cooling of the step plate and real cylinder head with different solid parts' orientations. The temperature histories predicted by the presented model fit very well with the available experimental data at different monitoring locations. Keywords: multiphase flow, immersion quenching, cast aluminium parts, computational fluid dynamics 0 INTRODUCTION Modern powertrain development of the internal combustion engine (ICE) is driven in the direction of reducing the weight ratio by replacing heavier metals with light low-cost alloys. An accurate prediction and optimization of the heat transfer characteristics within automotive, aerospace and processing industries is one of the more important influential factors for reducing fuel consumption and emission values [1]. For work pieces of complex geometrical nature, like the internal combustion engine's cylinder head and the deformation associated with thermal stress are a challenge to structural design optimization, where the difficulty arises from the ability to predict and analyze quench cooling heat transfer [2]. The immersion quenching cooling process from among all other heat treatment techniques has been long identified as one of the more important ways of achieving the desirable microstructure and mechanical properties of the metal piece [3]. This process provides even temperature distribution during the cooling process, leads to reduction of the residual stress levels, and consequently prevents distortion and cracking of the cast parts. During the immersion quenching process, a metal piece is heated up to a microstructure-dependent temperature, stays there for a while, and is then immediately submerged into a sub-cooled liquid like water, oil or polymers. Therefore, all three boiling regimes appear during the immersion quenching cooling process [4]. The film boiling regime occurs immediately after the heated piece is dipped into the sub-cooled liquid domain. The heat transfer of the cooling rate in this regime is relatively small, since the film is stable and the vapor blanket acts like an insulator. As soon as the surface temperature of the heated piece drops below the minimum film collapse temperature, also known as the Leidenfrost temperature, transitiive boiling starts, and the hot surface is partly wetted due to the collapse of the vapor film [5]. As the time proceeds, the hot surface temperature drops and the heat flux increases. When the wall heat flux reaches the maximum heat flux or the so-called critical heat flux, the nucleate boiling regime is entered. There the vapour film becomes unstable and starts to disappear. This results in a higher heat transfer rate between the metal and the fluid, and leads to faster cooling. As soon as the surface temperature of the heated part drops below the boiling temperature of the fluid, the convective heat transfer without phase change occurs, and the cooling flow rate becomes stable and low. Over recent years many analytical approaches for immersion quenching have been investigated and phenomena like vapor pocket generation, stagnation region, flow dynamics etc. have been successfully resolved by several researchers [2] to [5]. Based on the complexities and phenomena of multi-phase flows during immersion quenching applications, heat and mass transfer empirical models were proposed and developed. One of them is the model for mass transfer prediction, as proposed by Wang et al. [6], who assumed that the mass transfer rate due to boiling was proportional to the heat transfer rate in the fluid system. This model was later improved by Srinivasan *Corr. Author's Address: AVL-AST d.o.o., Ul. kneza Koclja 22, 2000 Maribor, Slovenia, rok.kopun@gmail.com 571 Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 et al. [7], and it has already been implemented within the commercial computational fluid dynamics (CFD) code AVL FIRE® where it is applied for the numerical simulations of quenching processes [8]. Published works like Wang et al[6], Srinivasan et al [7] and Greif et al [8], consider two separate domains, water and solid, which require separate simulations for each of them. The two simulations can be numerically coupled at the contact interface and communicate via the so-called ACCI (AVL code coupling interface) method after each time-step. The present paper describes a newly-developed quenching application method, where water and solid domains are treated as one. The methodology is known as the multi-material (MMAT) approach, where the surface temperature and local heat transfer coefficients are exchanged after each iteration and no longer after each time- step, as in the previously applied ACCI coupling method. The boiling phase change process between the heated part and the sub-cooled liquid side is handled by utilizing the Eulerian multi-fluid modelling approach. While for the fluid domain, general equations are solved within the framework of the multi-fluid modeling approach, only the energy equation is solved in order to predict the thermal field within the solid region. Additional to the MMAT approach, a new approach for modeling the variable Leidenfrost temperature threshold is also presented during this work. As long as the surface temperature is above the Leidenfrost threshold, film boiling appears, while transition boiling occurs when the surface temperature drops below the Leidenfrost temperature. Thus, the Leidenfrost threshold is a very important parameter for the numerical simulation of the quenching process, since it determines which regime is responsible for the calculation of the heat transfer coefficient. In previously published works [6] to [8] a constant Leidenfrost temperature threshold is assumed over the entire domain. This assumption has only limited validity, as discussed in the experimental work of Lubben et al. [9]. According to Lubben et al [9], it could be summarized that the Leidenfrost temperature for the quenching process depends on several factors, like the pool temperature (sub-cooling effect) [10], dipping and wetting velocity of the heated work piece [11], geometry and shape of the work piece, system pressure [12], surface roughness [5], etc. From previous work [13] it can be observed that, due to vapor bubbles rising from the bottom to the top of the quenched work piece, there is a higher vapor concentration at the top. This leads to the fact that transition boiling occurs at a lower temperature level in the upper part of the work piece, and consequently the Leidenfrost temperature is lower at the upper part of the quenched surface. Therefore, a model capable of resolving the spatially changing Leidenfrost temperature is of essential importance for achieving better agreement between the predicted results and the available measurement data. The temperature distribution within the solid part, obtained from the CFD simulation, can serve as a realistic input for subsequent finite element analysis (FEA) of thermal stresses within the quenched solid part [14]. This paper is organized as follows. In section 1 the mathematical model is described based on the previous work of Srinivasan et al. [7]. Brief introductions of the simulation set-up and experimental description are given in sections 2 and 3. Two different quenching orientations with various solid parts have been conducted, where the comparisons of the simulated results of a real-time quenching process with the available measurement data are presented and discussed in section 5. Conclusions and remarks are made in section 6. 1 MATHEMATICAL MODEL The Eulerian multi-fluid model considers each phase as interpenetrating continua coexisting within the flow domain, with inter-phase transfer terms accounting for phase interactions where conservation laws apply [15]. From the theoretical work of Drew and Passman [16] the ensemble averaged continuity and momentum equations are presented as: 1.1 Continuity Equation +v-(Wk)= Ž rkl, (1) 01 l=l,l*k subject to the compatibility condition: £ ak = 1, k = 1,..., N, (2) k=1 where a, p and stand for volume fraction, density and velocity vectors. The phase change rate (in this particular case, boiling) is rk and the subscript k is the phase indicator (k = l or k = v). Subscripts v and l denote the vapor and liquid phases in the current work. 572 Kopun, R. - Škerget, L. - Hriberšek, M. - Zhang, D. - Edelbauer, W. Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 1.2 Momentum Equation dak Pk vk dt + ^\ak P kv k vk ) = -ak^P + +^akTk +akPkg + Z Mu + v int Z_j kl'- (3) with p, g, t and vint are respectively the pressure, gravity vector, stress tensor and interfacial velocity. The interfacial momentum term is given by: Mm = Fd + FWL + Fl , (4) where FD indicates the drag force, FWL the wall lubrication force and FL the lift force. Detailed information concerning modeling the interfacial forces is described in the work of Kopun et al. [13] 1.3 Energy Equation With the assumption that the heat transfer rate between the vapor and liquid phases is rapid, then they are in thermal equilibrium [6] and [7]. In this case, the mixture enthalpy equation is solved, as described by: dpmh dt + V-(pm • v• h) = V-q +pm-d + dp + Pm - S' V + V- T-V + am^~ + dt N N + Z hv1+K, x r !=1J*k !=1J*k (5) with the mixture properties defined as: Vm =^at Vk > Pm =^ak Pk > Km = ^akKk > (6) k=1 k=1 k=1 where q is the heat flux, 0 is the specific enthalpy source, and the interfacial energy exchange between phases v and l is denoted as Hvl. Dynamic viscosity and thermal conductivity are presented as ^ and k, respectively. The heat flux q is given by: K _ q = Vh. (7) Cm ■ C ■ h ■ An, ■AT H, ' (8) where Cm, Cb, hb, Aint, AT' and Hfg are the closure coefficient, the boiling correction coefficient, the boiling heat transfer coefficient, interfacial area density, wall superheated temperature, and the latent heat of vaporization, respectively. The following subsections describe the applied models for heat transfer coefficient hFB of the different boiling regimes 1.4.1 Film Boiling Model The Bromely's model [17], originally applied for a horizontal tube, is employed to predict the film boiling heat transfer coefficient hFB is described by: hFB = 0.62 K pv Apg(Hfg + ÇAT') dbVv AT' 1/4 (9) where db is the length scale (vapor bubble diameter), is the vapor thermal conductivity, Ap = pl - pv is the density and C1 = 0.4-Cpv stands for the specific heat of the vapor. 1.4.2 Transition Boiling Model The heat transfer coefficient, for the transition boiling regime is given by: hr£ = QCHF+ QmHF -(1.0-*) , T^f _ T^f ■ v ' w sat where the corresponding heat flux is computed as Qchf, which denotes the Critical Heat Flux, and QMHF which stands for the Minimum Heat Flux. Detailed information about the applied quenching model in more detail can be obtained from AVL FIRE® Multi-fluid model solver theory guide [15] and the references [6] to [8], [13],[18] and [19]. 1.5 Variable Leidenfrost Temperature Model r 1.4 Boiling Model Based on the assumption that the boiling heat transfer rate is proportional to the phase change rate, the mass transfer predominantly controls the heat transfer. Thus, the phase change rate due to boiling can be written as: The Leidenfrost temperature is used to distinguish between film boiling and transition boiling regimes. As described above, it can be observed that the Leidenfrost temperature varies along the gravitational direction. Due to the higher vapor concentration at the top of the quenched work piece, transition boiling occurs there at a lower temperature level, and 573 Numerical Investigations of Quenching Cooling Processes for Different Cast Aluminum Parts Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 consequently, the Leidenfrost temperature is lower within the upper region and higher in the lower region. Based on extended measurements, a physically meaningful variable Leidenfrost temperature model has been developed. The proposed Leidenfrost temperature as a function of the vertical position depends on the vapor distribution is written as: TLeid (Z) = T AT Leid, range Leid ,ref (Z - Zmean X (11) where ATLeid,range is the range of the Leidenfrost temperature, T^ ref is the mean Leidenfrost temperature, z is the coordinate in the gravitational direction, and zmean the position where TLeid,ref is valid. Thus, the term Zmax - Zmin describes the maximum extension of the quenched work piece in the gravitational direction. The model implies that the Leidenfrost temperature decreases linearly from TLeid,max, at the bottom, to TLeid,min, at the top. Detailed information concerning the variable Leidenfrost temperature treatment can be obtained from the references [18]. 2 SIMULATIONS SET UP The presented MMAT approach triggers data exchange between the solid and liquid domains over the interface region, which is similar to the ACCI approach. The applied MMAT method is however different as the surface temperature and local heat transfer coefficients are exchanged after each iteration and no longer after each time-step, see Fig. 1. In the case of the ACCI methodology two separate simulations are numerically coupled at the interface and communicate via the so-called ACCI interface, Fig. 1a, whereas in the MMAT approach only one simulation with conformed computational mesh is applied describing water and solid domains Fig. 1b. Fig. 2. Computational domain set up with boundary conditions; a) step plate and b) cylinder head Fig. 2 displays the physical domains and the boundary conditions applied in the CFD simulations. In the current study, two different cases involving immersion quenching of the step-plate and real cylinder head configurations were investigated. The total number of cells for the step-plate test piece (Fig. 2a) was about 800,000 whereas the solid domain consisted of approximately 70,000 cells. The simulation set-up for the real cylinder head consisted aj b) MMAT -1 joined mesh INTERFACE Q solid V r.nr ✓ MAT 1 MAT 1 ■ ■ i :: 1 : Solve Mass Solve Momentum Solve Energy Eqn Solve Energy Eqn Fig. 1. Data exchange between liquid and solid domains; a) ACCI interface and b) MMAT interface 574 Kopun, R. - Škerget, L. - Hriberšek, M. - Zhang, D. - Edelbauer, W. Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 of approximately 2,400,000 cells for the liquid domain and about 2,100,000 cells within the solid domain Fig. 2b. With reference to the liquid domain, the prescribed atmospheric pressure is exerted at the outlet boundary on the top (green surface), and the inlet section is the given velocity boundary condition at the bottom (black surface). The dipping velocity is used to define the inlet velocity rising up from the initial water level (yellow part at the side) toward the outlet region. In current simulations, a constant velocity of 0.14 m/s is used until the final submerging depth has been reached. The domain walls at the sides are treated as adiabatic. Boundary conditions, dipping velocity, meshes and model parameters between the quenching scenarios (different orientations) during the entire simulation stayed the same Temperature measurements along the presented pieces were performed at different positions, and they are referred to as T1, T2, T3, T4, T5, T6, T7, T8 and T9, as shown in Fig. 3a for step plate and Fig. 3b for cylinder head. The numerical simulations were performed with commercial CFD code in which the Finite Volume approach with discretization of the governing equation and a SIMPLE algorithm for multiphase flows were used. The normalized residual limit for mass, momentum, and volume fraction were set at, the energy equation convergence limit was lowered to the value of. Only the energy equation was solved within the solid domain, whereas in the liquid part additional mass and momentum equations were solved. 3 EXPERIMENTAL WORK For the presented study, experimental investigations of the immersion quenching process were performed by the Nemak Company. Aluminium alloys made of AlSi7MgCu0 with material properties for thermal conductivity 183 [W/mK], specific heat 1166 [J/kg K] and density 2591 [kg/m3] were used for quenching simulation as displayed on Fig. 4. Inside the work piece three thermocouples for step plate and nine thermocouples for the cylinder head, of type "K" (Ni-CrNi) connected to the multi-meter (DEWE-TRON 2000) were used to record the temperature profiles. Fig. 3 shows the dimensions of the sample and the positions of the thermocouples, illustrated by dots, and referred to as Tb T2, T3, T4, T5, T6, T7, T8 and T9. The total height of the rectangular step plate solid was 201 mm, the width 149 mm and the variable thickness along the length, where the cylinder head shared the dimensions of 412*167x137 mm (length, width and height). Step plate Cylinderhead 1 J a) M&sm b> Fig. 4. Aluminium alloy; a) step plate, and b) cylinder head An air-forced oven was used to heat up the work piece to a uniform temperature of 780 K. Then the sample was carried from the oven to the pool filled with pre-heated water and submerged using the crane at a constant dipping velocity of 0.14 m/s. The environmental temperature was 294 K. The temperature profiles during quenching were Numerical Investigations of Quenching Cooling Processes for Different Cast Aluminum Parts 575 Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 recorded and several experiments were performed. Measurements at the test station were performed using different operating loads (different orientations and different aluminium alloys), with the cooling time set at 100 s. Detailed information concerning the measurement can be found in Kopun et al. [18]. 4 SIMULATION RESULTS AND DISCUSSION The comparison was performed next of the measured temperatures and numerically predicted results. This focused on the temperature distribution and on the phase distribution within the water domain represented by volume fraction. The presented study utilises four different orientations to cover the wide variety of quenching scenarios, as depicted in Fig. 5. The first orientation, shown in Fig. 5a stands for the step plate with the thick part up position, whereas by changing the solid orientation by 180°, the second orientation with thin part up appears; see Fig. 5b. The third orientation, shown in Fig. 5c features a so-called horizontal orientation which is similar to submerging a cup in an upside-down fashion. The fourth orientation, as demonstrated in Fig. 5d, features the submerging of a sideways cup, is so-called vertical orientation. Immersion quenching has been performed inside preheated water of 353 K. b) c) Fig. 5. Aluminium orientation of the quenched pieces; a) step plate thick part up, b) step plate thin part up, c) cylinder head horizontal orientation, and d) cylinder head vertical orientation 4.1 Case A: Step-Plate with Different Orientations The first case presented in this article demonstrates the numerical simulation of a so-called step plate immersed in a pool with a water temperature of 353 K. In Figs. 6 and 7 the measured and numerically predicted temperatures were compared at different monitoring points, Tb T2 and T3, along the height of the aluminium test piece. Fig. 6. Comparison between numerically predicted and measured solid part temperatures for step plate thick part up orientation Fig. 7. Comparison between numerically predicted and measured solid part temperatures for step plate thin part up orientation It can be seen from Fig. 6 that, for the thick part up case, the film boiling regime lasted up to 19s for areas around the monitoring point T1, where a fast transition regime was present later on The Leidenfrost temperature predictions varied between 660 K for point T3 to 550 K for the monitoring location T1, where excellent agreement between the numerical and measured data with the aforementioned new Leidenfrost temperature assumption were presented. 576 Kopun, R. - Škerget, L. - Hriberšek, M. - Zhang, D. - Edelbauer, W. Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 Similarly to the thick part up orientation, the results in the thin part up case were also well predicted and showed an excellent agreement between the numerical and available measured data, as seen in Fig. 7. Approximately 5 s shorter cooling time was achieved with the second orientation, where the Leidenfrost temperature prediction stayed the same as in the aforementioned case (550 K for monitoring point T3 and 660 K for the monitoring point Tj). The maximum deviation between the numerical and measured data was less than 2 s for all monitoring points in both orientations. The liquid volume fraction and temperature distribution within the structure at different time instants are displayed in Fig. 8a) for the thick part up case and Fig. 8b) for the thin part up case. It can be seen that the vapour film around the entire solid piece was present after 5 s, where a different trend in the temperature distribution was quite evident. The step plate in the thick part up case was still incompletely cooled down after 20 s (see Fig. 8f), where the uniform temperature had already been establish for the thin part up position. Detailed comparison between different pool temperatures can be obtained from Kopun et al. [18]. 4.2 Case B: Cylinder Head with Different Orientations The model's capability of predicting quench rates was tested on the simplified cylinder head, where the comparison between the measured and numerically predicted temperature history at nine different monitoring locations along the solid piece are presented, see Figs 9 and 10. A wide range of cooling regions were measured and carried out within the aluminium solid region, where different submerging orientations had been observed. Fig. 9 shows the numerically predicted results compared with the average value of the corresponding measured cooling curves of monitoring points Tj to T9 for horizontal orientation. In regard to all monitoring points it can be seen that the simulated film and transition boiling regime were in very good agreement with the available measured data. This implied that the proposed variable Leidenfrost temperature model could reasonably re-produce the cooling history. The predicted Leidenfrost temperatures varied from 563 K at point T3 to 713 K at monitoring point T5. The maximum deviation between the simulation and experiment was less than 4 s at all monitoring points. It was further observed that the entire aluminium cylinder head cooled down to the pool temperature after approximately 40 s. By changing the cylinder head position to vertical orientation, it was discovered in Fig. 10 that the predicted boiling regimes, film and nucleate were well described, and again they agreed well with the measurement values. All the model's parameters were exactly the same as in previous case for the horizontal orientation. Due to the larger vertical extension, only the Leidenfrost temperature range had increased. A.Trange in Eq. (11) was set to 200 K for the vertical orientation instead of 150 K for the horizontal. The predicted Leidenfrost temperature in the vertical case varied from 543 K at monitoring point T5 to 743 K at monitoring point T9. The maximum deviation between the measured data and the numerical results for all monitoring points was less than 4 s. The liquid volume fraction and the surface temperature distribution at different time instants during the entire quenching cooling process for Numerical Investigations of Quenching Cooling Processes for Different Cast Aluminum Parts 577 Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 Time (s) Fig. 9. Comparison between numerically predicted and measured solid part temperatures for horizontal quenching orientation with water temperatures of353 K at monitoring points T1 to T9 Fig. 10. Comparison between numerically predicted and measured solid part temperatures for vertical quenching orientation with water temperatures of 353 K at monitoring points Ti to Ta horizontal orientation are presented in Fig. 11. It was found in Fig. 11a that rapid boiling occurred at a time of close to 5 s, where transition boiling was the dominant regime and the surface temperature was somewhere around 753 K. As soon as the surface temperature dropped below the Leidenfrost 578 Kopun, R. - Škerget, L. - Hriberšek, M. - Zhang, D. - Edelbauer, W. Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 Time: 5 s Time: 15 s Time: 25 s Time: 50 s (e) (f) (g) (h) Fig. 11. Liquid volume fraction and the surface temperature distribution at different time steps: a) 5 s, b) 15 s, c) 25 s and d) 50 s for horizontal and vertical orientations temperature, transition boiling and later nucleate boiling took place. At a time of 25 s, it was observed in Fig. 11c, that boiling occurred only at the upper parts of the cylinder head, where there were still some hot temperature areas, especially in those regions of higher material thickness. The temperature of the bottom region was already too low for phase change. Due to the special geometrical configuration it can be seen in Fig. 11d that the gas phase was trapped in the middle part of the cylinder head, and couldn't escape. This was caused by the low positions of the inlet/ outlet channels. Thus, quenching of this cylinder head in the horizontal position should not be recommended. By comparing the vertical orientation with the horizontal case, it can be seen that the surface temperature distribution was significantly different between the two cases. After 5 s, the bottom area reached a temperature of about 570 K, (see Fig. 11e), whereas for the horizontal orientation case, the surface temperature of the same region was close to 750 K. It can be seen in Fig. 11g that almost half of the workpiece had cooled down to the liquid pool temperature after a time of 25 s, whereas the phase change distribution was present only at the upper parts of the cylinder head. As demonstrated in Fig. 11h, similar to the case with horizontal orientation, a homogeneous temperature field was established after 50 s and no more vapour was produced. 5 CONCLUSIONS The applied CFD code AVL FIRE® is capable of predicting real time quenching effects for simple step plate and cylinder head applications with different solid parts' orientations. It can be seen that the implemented model's extensions for variable Leidenfrost temperature approximation in combination with the Multi-Material model showed good agreement between numerical results and measurements. The film and transition boiling regimes were well described and numerically predicted. The variable Leidenfrost temperature had a significant effect on the simulation. It has been demonstrated that by changing the quenched pieces' orientation, the solid parts cool with different trends. It was found that, based on geometrical configuration, the gas phase was trapped inside the structure regarding the horizontal case. By changing the orientation to vertical, this trapping phenomenon did not occur anymore. The cooling curves were well predicted, Numerical Investigations of Quenching Cooling Processes for Different Cast Aluminum Parts 579 Strojniški vestnik - Journal of Mechanical Engineering 60(2014)9, 571-580 and there was a very good agreement between the numerical and measured data. It can be concluded that the introduction of a Multi-Material model, which would enable exchange of simulation data during the iterative solution process would lead to further improvements in computational results after each iteration. This, together with the improved CFD model provides a solid basis for future research within the field of real-time immersion quenching applications. 6 FUNDING This research work was partly funded by the European Union, European Social Fund and SPIRIT Slovenia, Slovenian Public Agency for Entrepreneurship, Innovation, Development, Investment and Tourism. Parts of the results presented in this paper were obtained within the FFG funded project QUENCH-IT, project number 828697. 7 REFERENCES [1] Abdulhay, B., Bourouga, B., Dessain, C. (2011). Experimental and theoretical study of thermal aspects of the hot stamping process. Applied Thermal Engineering, vol. 31, no. 5, p. 674-685, DOI:10.1016/j. applthermaleng.2010.11.010. [2] Srinivasan, V., Greif, D., Basara, B. (2012). On the heat and mass transfer modeling to simulate quenching heat treatment process. 6th International Quenching and Control of Distortion Conference, Chicago, [3] Moravčik, R., Stefanikova M., Čička R., Čaplovič L., Kocurova, K., Šturm, R. (2012). Phase transformation in high alloy cold work tool steel. Strojniški vestnik -Journal of Mechanical Engineering, vol. 58, no. 12, p. 709-715, D0I:10.5545/sv-jme.2012.531. [4] Babu, K., Prasanna Kumar, T.S. (2009). Mathematical modelling of surface heat flux during quenching. The Minerals, Metals & Material Society and ASM International, vol. 41B, p. 214-224, D0I:10.1007/ s11663-009-9319-y. [5] Meduri, P.K., Gopinath W.R., Vijay D.K. (2009). Wall heat flux partitioning during subcooled forced flow film boiling of water on a vertical surface. International Journal of Heat and Mass Transfer, vol. 52, no. 15-16, p. 3534-3546, D0I:10.1016/j. ijheatmasstransfer.2009.02.040. [6] Wang, D.M, Alajbegovič, A., Su, X.M., Jan, J. (2003). Numerical simulation of water quenching process of an engine cylinder head. Proceedings of ASME FEDSM, 4th ASME JSME Joint Fluids Engineering Conference, Hawaii, D0I:10.1115/FEDSM2003-45538. [7] Srinivasan, V., Moon, K., Greif, D., Wang, D.M, Kim, M. (2010). Numerical simulation of immersion quench cooling process using an Eulerian multi-fluid approach. Applied Thermal Engineering, vol. 30, no. 5, p. 499509, D01:10.1016/j.applthermaleng.2009.10.012. [8] Greif, D., Kovacic, Z., Srinivasan, V., Wang, D.M., Suffa, M. (2009). Coupled numerical analysis of quenching process of internal combustion engine cylinder head. BHM Journal, vol. 154, no. 11, p. 509517, D0I:10.1007/s00501-009-0514-6. [9] Lubben, T., Frerichs, F., Zoch, H.-W. (2011). Rewetting behaviour during immersion quenching. Strojarstvo, vol. 53, no. 1, p. 45-52. [10] Dhir, V.K., Purohit, G.P. (1978): Subcooled film boiling heat transfer from spheres. Nuclear Engineering and Design, vol. 47, no. 1, p. 49-66, D0I:10.1016/0029-5493(78)90004-3. [11] Drucker, M., Dhir, V.K. (1981). Effects of high temperature and flow blockage on the reflood behavior of a 4-rod bundle. EPRI Report INIS, OSTI ID: 5714419. [12] Cheng, S.C., Lau, P.W.K., Poon, K.T. (1985). Measurements of true quench temperature of subcooled water under forced convective conditions. International Journal of Heat and Mass Transfer, vol. 28, no. 1, p. 235-243, D0I:10.1016/0017-9310(85)90025-0. [13] Kopun, R., Greif, D., Edelbauer, W., Zhang, D., Tatschl, R., Stauder, B. (2013). Advances in numerical investigation of immersion quenching at different pool temperatures. 22nd International Conference SAE, Sao Paulo, p. 369, D0I:10.4271/2013-36-0369. [14] Trzepiecinski, T., Lemu, H.G. (2014), Frictional conditions of AA5251Aluminium alloy sheets using drawbead simulator tests and numerical methods. Journal of Mechanical Engineering, vol. 60, no. 1, p. 51-60, D0I:10.5545/sv-jme.2013.1310. [15] AVL LIST GmbH (2013). FIRE CFD solver, Eulerian multi-fluid model. Solver Theory Guide, Graz. [16] Drew, D.A., Passman, S.L. (1999). Theory of Multicomponent Fluids. Springer, New York, D0I:10.1007/b97678. [17] Bromely, L.A. (1950). Heat transfer in stable film boiling. Chemical Engineering Progress, vol. 58, p. 6772. [18] Kopun, R., Skerget, L., Hribersek, M., Zhang, D., Stauder, B., Greif, D. (2014). Numerical simulation of immersion quenching process for cast aluminium part at different pool temperatures. Applied Thermal Engineering, vol. 65, no 1-2, p. 74-84, D0I:10.1016/j. applthermaleng.2013.12.058. [19] Kopun, R., Greif, D., Zhang, D., Stauder, B., Skerget, L., Hribersek, M. (2013). Numerical investigation for immersion quenching of single and clustered test piece configuration. JSAE Annual Congress, Nagoya. 580 Kopun, R. - Škerget, L. - Hriberšek, M. - Zhang, D. - Edelbauer, W.