Numerical modelling of the Water-Quenching Process validated through experiments with IN718 Nickel-Based Superalloy

Quenching is often a necessary step during manufacturing to tailor microstructures of metals and alloys for desired applications. Although employed extensively throughout human history, it relies heavily on experience and often trial and error. Therefore, the availability of computational tools, capable of describing the physics of the process during quenching, would help effective design, lower manufacturing costs, and make the process more environmentally friendly. This work analyses and validates a numerical procedure for immersion water quenching using computational fluid dynamics. The methodology employs the partitioned approach. Hence, the heat transfer information is exchanged at the interface of the part and the surrounding environment (i.e., water). An energy equation and an Eulerian two-fluid model describe the solid and fluid regions, respectively. The validation experiment includes water quenching of an IN718 nickel-based superalloy solid cylinder, in axial direction, with embedded thermo-couples for the measurements of cooling curves during the process. As such, no phase transformation within the solid material is assumed. The computational procedure focuses on determining the heat transfer coefficient due to fluid dynamics and the quenchant (i.e., water) phase change during cooling. The cylinder experiences effects of various heat transfer regimes, including film, transition and nucleate boiling, and natural convection intensified by buoyancy forces due to the presence of two fluid phases. A good agreement with validation data is achieved presenting an innovative numerical tool capable of predicting with high accuracy the temperatures and cooling rates of metal alloys during water quenching. Some complications are observed and reported where vapour movement is obstructed. Therefore we highlight further steps to alleviate any issues.


Introduction
Quenching is a heat treatment procedure used to modify the microstructure of most engineering material alloys (i.e., both ferrous and non-ferrous) to improve their properties for desired applications [1,2].Depending on the material and requirements, a metallic component is typically mechanically processed (e.g., forging, hot-rolling) or heat treated at a specific temperature, followed by rapid cooling in a fluid, (i.e., the quenchant) which can be water, polymer solution or oil depending on the required cooling rate (i.e., water imposing the fastest).We consider immersion cooling in water in the current work.
Due to its complex multi-physics nature, this metallurgical process continues to pose significant challenges.The development of computational tools and analysis methodologies improving the predictability would (i) decrease product wastage rate [1]; (ii) shorten manufacturing have used the commercial software AVL FIRE.Wang et al. [11,12] were the first to implement immersion quenching computational methodology within that specific software.They solved a conjugate heat transfer problem between solid and fluid regions by assessing the interface temperature and heat transfer coefficient.An energy equation governed the solid region.Meanwhile, the Eulerian-Eulerian two-fluid model (  ) described the fluid behaviour using a separate set of continuity and momentum equations for each phase completed with a mixture energy equation.Film and nucleate boiling were the only heat transfer regimes (HTR) considered.Nevertheless, the nucleate boiling mass source was estimated only by multiplying the film boiling mass source.
The majority of subsequent studies extended the work within AVL FIRE.Srinivasan et al. [13,14] introduced a new approach for modelling the solid-fluid interface, implemented transition boiling HTR, neglected nucleate boiling HTR and modified the drag force.Greif et al. [15] continued the work and combined the simulation procedure with structural mechanics analysis in ABAQUS.Jan et al. [16] then used the software to predict immersion quenching behaviour of an engine cylinder head.Later work [17,18] considered variable Leidenfrost temperature, lift and wall lubrication interfacial momentum forces.Further, Zhang et al. [19] considered separate energy equations for gas and liquid using thermal homogeneous assumption.The last two research papers using AVL FIRE discussed the influence and calibration of various heat transfer framework coefficients, which usually need to be chosen and tuned to get valid results [20,21].
Excluding the work done on AVL FIRE, there have been four other largely relevant attempts to simulate the metallurgical process numerically, to the best of the authors' knowledge.The first used the mixture model and bubble crowding in Fluent [22].Other utilised the Finite Element -Galerkin method with monolithic coupling, tracking the liquid-vapour interface using the Level Set method and anisotropic mesh adaptation [23,24].Finally, [25] used drift-flux mixture-model in Comsol assuming constant vapour temperature.
Further available quenching and boiling research covers a large body of numerical modelling considering CFD [26][27][28][29], focusing on phase transformation [30][31][32][33] and experimentation [34][35][36].That research is very valuable and aids our understanding of quenching.This knowledge serves as inspiration to the authors and various aspects of the presented models might be further used to improve our modelling capabilities.Though, there might often be substantial differences comparing to our research due to the fact that the existing literature focused on a particular boiling regime, or did not include phase change, or used different coolants or added momentum to the coolant due to agitation.
Immersion quenching is a complex process, and even though the computational research in the area is extensive, the success varies.Furthermore, computational results are often difficult or impossible to reproduce due to a lack of experimental and computational data.This article shows, validates and discusses experimental and computational results produced with OpenFOAM using the methodology as described in [10,37] and demonstrates its usability across the whole boiling curve.Hence, convective heat transfer, nucleate boiling, transition boiling (partial film boiling), and film boiling stages have all been considered [38].
This article consists of six main sections.After the introduction and a brief description of the computational methodology, a section about experimental work used for the validation is presented.Only then we show the computational setup.Finally, the numerical results are introduced and discussed, followed by a summary with conclusions and future work.

Computational methodology
In this section, we present a summary of the computational methodology.Water immersion quenching is inherently a multiphase problem whereby solid(i.e. the sample), liquid and vapour occur concurrently and at different ratios depending on the local temperature.At the same time, air is treated as vapour, as it was done in previous studies without a lack of generality or accuracy [11][12][13][14][15][16][17][18][19][20][21].Similarly to AVL Fire, we have used the partitioned approach.Thus, solid and fluids are represented by discrete regions interacting via temperature boundary condition affected by thermal turbulent diffusivity, allowing for four heat transfer regimes: (i) convective heat transfer, (ii) nucleate boiling, (iii) transition boiling, and (iv) film boiling.The solid sample is described by the energy equation assuming no phase/microstructure transformation [32,33], while fluids are governed using Eulerian-Eulerian two-phase methodology with boiling and condensation capabilities.Undoubtedly, the prediction of the boiling curve is a remarkably challenging problem even more though because quenching is a transient process.It does not imitate exact results from quasi-steady-state experimental work except for film boiling and natural convection [3,10,37].A more detailed description is provided in [10,37].

Fluid region
The fluids are described with six volume averaged transport equations: phase-intensive conservation of mass, momentum and energy.These are supplemented with turbulence and phase continuity equations to close the system.It should be underlined here that each fluid phase is represented by its dedicated governing equations and variables (phase-intensive) and not by a single unified equation, increasing the accuracy, flexibility and capabilities of the solver.
The mass conservation is given by where , ,  and  are volume fraction, density, time, and velocity, respectively, and the subscripts  and  stand either for  −  or  − , but  ≠ .Moreover, ṁ is the mass source accounting for the phase change at the regions interface [10,37].
Similarly, the conservation of momentum is described by where  denotes the pressure field, common for both fluid phases. stands for the acceleration due to gravity, and    is the effective viscosity    , =  , +  , , where   is the kinematic viscosity, and   is the turbulent viscosity.  represents interfacial forces dictated by the fluid phases momentum interaction.We account for drag, lift, wall lubrication, virtual mass and turbulent dispersion [10,37].Following on, the conservation of energy in the form of specific enthalpy can be written as where ℎ represents specific enthalpy;   = 1 2 |  | 2 , and    , =  , +  , is the effective thermal diffusivity with   ,   standing for molecular and turbulent thermal diffusivity, respectively. is the sensible interfacial heat transfer coefficient evaluated as a product of the interfacial area concentration and a heat transfer coefficient calculated at the fluid phase interfaces.Finally,   is the liquid-vapour interface temperature.
Turbulence is considered through the use of the fluid phases mixture − model by Behzadi et al. [39].The only differences are the considerations of turbulence generation due to the presence of bubbles [40] and the usage of an effective density for averaging.The phase continuity equation obeying   +   = 1 is based on the work by Weller [41] and utilised in the form introduced by Rusche [42] to account for turbulent dispersion.

Solid region
The sole equation to solve in the solid region is energy conservation.No metallurgical or mechanical fields are taken into account.In reality, the system may sometimes be more complex, and phase transformation and plastic deformation generate heat [1].Nonetheless, this work is concentrated mainly on conjugate heat transfer and fluid dynamics as the primary research areas.
The solid energy equation is presented in the form of specific enthalpy: where subscript  stands for the solid region. is thermal diffusivity defined as   =    , , while  represents thermal conductivity and   specific heat capacity.

Regions' interface
A critical component of our computational methodology is the precise definition of boundary conditions at the solid-fluid interface.This aspect is paramount for ensuring energy conservation across the interface and maintaining numerical stability throughout the simulation process.To our knowledge, the integration of approaches from [14], Peltola et al. [43], Kurul & Podowski [44], and Jayatilleke [45] represents a novel contribution to the literature.This work adopts and synthesises these methodologies to advance the understanding of fluidsolid heat transfer exchange.
The top-level idea is to satisfy temperature and its gradient continuity across the region's interface accounting for all three existing phases using the temperature boundary condition.That has access to the turbulent thermal diffusivity boundary condition which considers the fluid boiling effects.The interface temperature then follows: where the phase-intensive Biot number is and the effective thermal conductivity is formulated as    , = ( , +  , ) , .Due to space discretisation,   represents the distance between the boundary cell centre (superscript C) and the boundary face centre at the solid-fluid interface (superscript I).Alternatively, it is practical to express the phase-intensive Biot numbers as the fluid Biot number, taking into account all solid, vapour and liquid phases, including the boiling effects: where  (7) .
The thermal turbulence diffusivity boundary condition then assumes that heat conducted from the solid region to the interface must be further transferred to the fluid phases according to their volumetric fraction in the boundary cell.
where TB [14] FB [47] The HTRs are determined according to the interface temperature in comparison to   ,   and   , as characteristic temperatures [38] (Table 1).The temperature of the departure from nucleate boiling   is calculated following Schroeder-Richter et al. [46], and the Leidenfrost temperature   is estimated using wetting temperature    equal to 987 • C, local characteristics and vapour fraction.   is our only fitting parameter to replicate the experiment.In reality, the   is a function of local parameters such as wall temperature, orientation, geometry, fluid properties, etc. [1].Its models and correlations exhibit substantial uncertainty and can often be inaccurate [3].Therefore, it should be further investigated through targeted experimental and simulation studies in the future.

Experiment description and analysis
Prior to the numerical simulations, quenching experiments were conducted and the cooling curves were measured at specified locations for further analyses.The material used for this study was IN718 nickelbased superalloy, supplied in a form of a billet with 300 mm diameter and 2500 mm length.Further information about the alloy's microstructure and the need for fast cooling after heat treatment is given in the following metallurgical books [48,49] and the publication by Rahimi et al. [50].A solid cylinder with 200 mm diameter and 200 mm length was machined from the as-received billet to reasonably fine (i.e., RA ∼ 3 μm) surface finish condition.For the acquisition of cooling curves during quenching, holes with 1.1 mm diameter were drilled to the specified locations using an electrical discharge machine (EDM) for drilling.The measurement locations are visualised in Fig. 1b (orange dots) and the coordinates are tabulated in Table 2. N-type thermocouples, with 1 mm diameter, were positioned at each location and tightened using a high temperature cement to prevent any potential movement during quenching.The sample was placed in a furnace and heated to 980 • C and kept for 1 hour.The heated cylinder was then taken using a semiautomatic robotic arm and quenched in water vertically.Two type 316 stainless steel holders used for holding and transferring the sample from the furnace to the quenching tank, which were mounted on the top face of the cylinder using M16 screws (i.e., stainless steel) to minimise the interaction with the fluid dynamics during quenching.(Fig. 1a).The temperature of the quenchant (i.e., water) was at 20 • C, while the average initial specimen temperature measured in the heating furnace using thermocouples placed in the cylinder was about 982 • C.
For reproducibility, the experiment for each cooling condition (e.g., sample orientation) of the solid cylinder was repeated at least three times.After bringing the cooling curves into the same time coordinate (i.e., to correct for slightly different transfer times which was less than 3 s difference for different repeats of the same test), the data were repeatable with less than 5% variation in the results recorded by the same thermocouples at the same locations.
Analysing the experimental results is crucial in providing a comprehensive understanding of the process to set up the computational framework.Videos and photos were taken during the quenching process (Fig. 2) for visual observation and to reveal a correlation with the cooling curve data.In Fig. 2a, the cylindrical specimen was located in the heating furnace, when the door has just opened, for a forklift to transport it into the quenching vessel.Once the hot specimen is

Table 2
Coordinates of the locations of the thermocouples used for the measurements of cooling curves.The origin and directions of the coordinate system are highlighted in Fig. 1b.submerged in water (Fig. 2b), wetting fronts appear at the cylinder's bottom edge, where the first contact with coolant took place.Additionally, increased cooling is seen above the thermocouple wires.They disturb the thermal boundary layer and cause the vapour film to detach.Video footage also reveals that the location repeatedly cools and heats.Thus, the vapour film is assumably reintroduced.Nonetheless, this phenomenon is not apparent from thermocouple measurements, most likely due to its short period (measurement data are acquired every 0.2 s) or thermocouple locations.In the same photo, the wetting front also emerges at the top of the cylinder surface.The wetting fronts eventually meet, and the cylinder does not emit thermal energy in the visible spectrum anymore (Fig. 2c).This cooling stage is characterised by large bubbles being chaotically released from the cylinder bottom.The wetting front is usually defined by a sudden liquid appearance at a specimen surface.Unfortunately, the current experimental setup cannot precisely detect this event.Therefore, in describing the figures above, the wetting front term was used more vaguely as a situation when the visible thermal radiation halts.
The whole specimen cools for more than 2500 s because it stores an extensive amount of energy relative to its surface area.Nonetheless, the research will concentrate on a limited time range only.The crucial and most challenging period to predict is the initial 25 s within the coolant when film, transition and nucleate boiling are dominant.Although the simulation of the rest of the experiment is a matter of additional computational time and should not pose any significant challenge, some results of longer simulations (200 s within the water) will be presented and analysed in Section 5.
In this manuscript we provide temperature measurements at various time points and temperature-time derivatives against temperature and time, instead of boiling curves in terms of heat flux and temperatures.Temperature-time derivatives are directly proportional to the divergence of heat flux; hence, they are directly related to boiling curves, as shown here It should be mentioned that, from a manufacturing point of view, the cooling curve is of primary interest while local heat fluxes are of minor importance.This is because the cooling curves are those defining the microstructure of the metal alloy, and as a consequence its properties.
In addition, it is impossible to measure local heat flux for large samples during quenching with today's technologies, hence the approach taken here (i.e., inverse heat transfer analysis), which is also a common practice by leading industries (e.g., Rolls Royce and Aubert & Duval who have sponsored these experiments).
The cooling curve results recorded during the experiment at two thermocouple locations namely,  1 (cylinder bottom centre) and  8 (cylinder side), are chosen for discussion here.The data are shown in Fig. 3, and the time corresponds to Fig. 2. The immersion process begins at 31 s.The time before captures specimen heating and transfer from the furnace.The cooling curves (Fig. 3a) start with a rapid temperature drop at the point of immersion.That is significantly stronger for  1 comparing to  8 as its submerged first.Furthermore, the cooling rate is smaller for  8 because its affected by generated vapour and released air.Further on, inflexion points are notable.These correspond to the peaks of cooling rates (Fig. 3b) and represent the Peak Heat Flux (  ), hence, the Departure from Nucleate Boiling.Again, this phenomenon happens sooner for  1.Soon after this moment, the cooling behaviour becomes very distinctive.On one side, the vapour cannot detach freely from the bottom surface, creating vapour film and increasing temperature ( 1 at about 37 s).The cooling rate curve exposes this behaviour by showing negative values with Minimum and Peak Heat Fluxes ( ,   ′ ).Conversely, the cooling trend at  8 is more straightforward.It does not exhibit clear MHF neither multiple PHFs.Nucleate boiling and convective heat transfer are encountered, and temperature declines towards the surrounding quenchant.The reader should be further aware that these curves are measured within the cylinder domain at around 3 mm from the surface (for exact locations see Table 2), hence, the actual characteristic temperatures are at lower temperatures affected by thermal conduction between the measurement locations and the cylinder surface.
This analysis suggests that the cooling behaviour at different locations might vary substantially, complicating the heat transfer prediction.Secondly, multiple heat flux peaks (  ) can occur during  one quenching procedure at an identical location.Moreover, they can take place at very different temperatures and heat fluxes.Thirdly, the cooling does not have to begin with film boiling.It can be omitted,  8, or appear later due to the vapour build-up,  1.Thus, nature does not necessarily obey the simplified pool boiling example from literature [38].These observations align with the previous discussion on variability of Leidenfrost temperature uncertainty and extends to the departure from nucleate boiling.

Computational setup
This section describes the computational domain, its space discretisation, and the boundary and initial conditions.The intricacy of the mesh sensitivity is also discussed.

Geometry, initial and boundary conditions
The computational domain is described using Fig. 4. It takes advantage of the cylinder and quenching tank symmetry, allowing it to simulate only one quarter of the actual domain.The symmetry assumption holds when the forklift, collars and thermocouple wiring are neglected.This decision was made for better mesh control, allowing us to focus on its impact.A  boundary condition (BC) bounds the domain, and at the top is an .Depending on the region of interest, the regions' interface BCs are labelled  __ or __ .Further, air (vapour) is above the water, which reaches 1100 mm.The cylinder submerging effect is replicated using the bottom  as an  where only liquid enters at constant velocity until it reaches the required height.Subsequently, the inlet BC is switched to the  definition.Table 3 provides an overview of the BCs definitions.  are specified using no-slip BC for both momentum equations and constant temperature at 20   4).The same equation is used for the interface temperature by the solid region.The heat transfer coefficient and, hence, the impact of boiling is estimated using the procedure from Section 2.3, thus via the computation of the thermal turbulent diffusivity.Finally, symmetry is used for both regions to avoid excessive computational power needs.
The physical properties of all three phases are shown in Table 4.All characteristics are assumed constant except the vapour density, which follows the ideal gas law.Liquid density can be defined as constant because the impact of fluid phases density difference is considered far more significant.The cylinder material is Inconel 718, as stated in the description of the experiment previously.

Mesh sensitivity analysis
This section describes the computational mesh and elaborates on its impact on the simulation results, by analysing data from three locations in the quenched cylinder.These are  1,  4 and  8, cylinder bottom, top and side, respectively.
Meshes for both regions are purely hexahedral, block-structured, and growth ratios are frequently exploited (Fig. 5).The solid region is often more refined than fluid, and the Arbitrary Mesh Interface interpolation procedure is used to cope with regions' mesh non-conformity.Several meshes are used to converge each region (Table 5).The focus is on the boundary cells' perpendicular dimensions, yet the mesh is refined in all directions.Fig. 5 shows an example of regions' meshes used in these analyses.The particular mesh combination has 31,140 fluid and 10,440 solid elements.The fluid elements are visualised in light blue, while the solid are grey.During the discretisation, the mesh smoothness, especially alongside the interface, was emphasised to cope with high gradients in multiple variables.
Table 5 shows number of elements and perpendicular dimensions of boundary cells for both regions.Throughout this section, various mesh combinations are labelled using the total number of elements in each region.
Fig. 6 displays cooling curves and vapour volumetric fraction using a solid mesh with 58,800 elements and various fluid meshes.It is evident that the sensitivity varies across the locations.Relatively small temperature differences are recognised at the cylinder's top (Fig. 6b) and side (Fig. 6c).However, very distinct behaviour is observed at the cylinder's bottom, where significant discrepancies are manifested (Fig. 6a).The fluid region refinement causes an apparent increase in the temperature together with a considerable difference in the cooling rates.The maximum cooling rate drops with refinement.In comparison, utilising the most refined mesh results in a less smooth curve.Additionally, using Based on the authors' previous experience with the quenching of a horizontal thin plate [10], a rise of vapour volume fraction is the expected cause of the increasing temperature at location  1. Visualisation (Fig. 6d) indeed shows the build-up of vapour hand in hand with the mesh refinement (2-7 s).The vapour film then acts as insulation compared to when liquid water is in contact with the solid interface and is allowed to boil.The thin plate in our previous publication was less affected by this phenomenon of vapour build-up.Its thickness was minimal, so the stored heat could be redirected to the top surface if insulated from the bottom.On the other hand, the cylinder's cooling process is very different.The heat is conducted towards the interface, smoothing the solid temperature gradient when vapour is generated.
Further analysing Fig. 6d, it is apparent that there is no noticeable difference in the vapour volume fraction between mesh configurations with 31,140 and 59,008 elements in the time range between 7 to 20 s.After that period, a sharp drop is noticeable.Nonetheless, the finest mesh still causes the vapour to reach its physical maximum almost throughout the whole visualised period.On the cylinder's side and top, significant fluctuation with a relatively small impact on the cooling data occurs.
An analogous analysis is performed for the solid mesh refinement using the fluid mesh with 59,008 elements (Fig. 7).The cooling curves (Fig. 7a, b, c) show good behaviour for all locations and are largely mesh independent.The main discrepancy seems to be at the bottom between 3 to 7 s, where the flow is very violent.
Following the mesh sensitivity analysis, it can be concluded that the behaviour depends very much on the surface orientation concerning gravitational acceleration and the extent of vapour movement restriction.The results are found mesh independent on various surfaces yet not at horizontal surfaces facing downwards.There, only the solid mesh converges.Increasing the refinement within fluid leads to vapour volume fraction growth, hence to surface insulation.It is further presumed that the behaviour might be improved by accounting for bubble dynamics (break-up, coalescence and agglomeration).It should be mentioned that the boiling models applied at the solid wall as well as their implementation are predominantly developed and validated for vertical and forced flows [51][52][53][54][55][56][57][58][59][60].The attempt to model the complete quenching process without agitation highlights the need for further development.For the data presented in Section 5, the meshes with 59,008 fluid and 451,500 solid elements are chosen.This mesh combination is independent of the solid mesh refinement and provides a realistic amount of vapour at the cylinder bottom.

Results and discussion
This section compares the numerical and experimental results followed by a discussion on the physics of the process.The presented information are devoted to the locations of the cylinder from which   cooling curves are recorded during the experiment.The focus is on temperature and cooling rates as a function of time and temperature, respectively.Also, each presented figure has a schematic sketch of a cylinder with highlighted location of the thermocouple for clarity.At the end of the section, vapour volume fraction, velocity and temperature spatial and time distributions are discussed.
Fig. 8 shows plots of experimentally measured and numerically predicted cooling curves at locations  1,  3,  4,  6,  8, and  9.The shown data are limited to the initial 200 s from the start of cylinder's submergence, where all the interesting and complex boiling and heat transfer phenomena occur.There is no modelling restriction in our computational approach if someone wants to further extend these simulations and with additional computational time thermal equilibrium between the solid part and the surrounding fluid will be eventually established.The agreement with experimental results at most locations is near a perfect match.The location at the cylinder's bottom surface,  1, exhibits a considerable cooling delay during the initial 30 s, which is also captured by the numerical prediction (Fig. 8a).The side locations,  6,  8, and  9 (Fig. 8d, e, f), show a temperature development when the slowest cooling occurs at the cylinder mid-height.The cause is vapour build-up and detachment of bubbles as they move upwards.Fig. 8b presents temperature data at the cylinder centre where the numerical results match with the experiment and temperature starts decreasing at about 150 s.The numerical cooling rate seems to be a bit slower.This might result from several factors but most likely from the accumulated inaccuracies at other locations and especially at the top surface  4.
Fig. 8c shows the results at the centre of the cylinder's top surface, where the computed temperature drops at a substantially faster rate.It is believed that the reason is a substantial amount of vapour acting as insulator, resulting from the specimen holder and forklift arms [37], which were not included in the computational model for simplicity.
The rest of this discussion concentrates on the initial 25 s unless cooling rate versus temperature data are discussed, otherwise.The initial 25 s period is dominated by transition and nucleate boiling stages, yet occasional film boiling and natural convection also occurs.For instance, the latter emerges at the cylinder bottom edge at about 22 s.The data presented in the figures also visualise the HTR changes at the location of interest nearest boundary faces.However, notice that these changes are predefined by the wall boiling BC (Section 2.3).The HTR changes due to vapour behaviour are not highlighted visually.
The cylinder bottom is one of the most complex locations for numerical predictions.The secrete is in the vapour distribution and dynamics.At specific locations, like the cylinder bottom centre (Fig. 9a), the vapour volume fraction is very high so that the cooling rate is significantly restricted.This behaviour is down to film boiling.The computational prediction does an excellent job demonstrating this response for the initial 15 s, after which the predicted temperature decrease lags behind the measured data.Nevertheless, quite distinct behaviour is observed at 3.5 mm above the bottom surface (Fig. 9c), where the experimental temperature decreases considerably faster despite being further from the regions' interface.This phenomenon is likely caused by high cooling rates at more distant locations where the heat is redirected.An example is  10 (Fig. 9b), where a similar cooling pattern is recognised.This location is about 47 mm away in radial distance from  1, closer to the cylinder bottom edge.The numerical results on both locations,  2 and  10, pose a significant temperature mismatch compared to the experiment.The authors see the vapour/bubble dynamics as the most reasonable ground to explain this.
The cooling curve data at the cylinder bottom are complemented with cooling rates, which are more sensitive to changes and can reveal previously unseen effects (Fig. 10).Neither experimental nor numerical results manifest initial film boiling at the centre location  1, while both exhibit two peaks with film boiling-like behaviour between them (Fig. 10a).The initial PHF is well predicted in terms of both magnitude and trend.Also, a loop in data is indicated nicely at about 720 • C, representing the heat-up when the cooling rate is negative.Directing the focus to the location deeper in the cylinder,  2 (Fig. 10b), only one prominent peak is visible during the experiment compared to the numerical solution, which predicts two.Thus, as mentioned before, this location is far more influenced by the location below,  1, than in the case of the experiment.Lastly, the numerical solution for  10 (Fig. 10c) indicates the maximum cooling rate at a significantly lower temperature since it follows an initial period of minimal cooling rate.It seems that it is affected by a vapour presence to a much more considerable extent than the other two locations right from the beginning of submergence.Now, attention is turned to the cylinder's vertical surface.Here, better agreement with the experiment is expected because it resembles the usual flow conditions applied in the problems solved in the literature to validate parts of the methodology used.The numerical results are compared to the experiment at four locations,  6,  7,  8 and  9, where one of them is deeper in the cylinder (Fig. 11).It is concluded that an excellent agreement is reached for the temperature trends, yet the numerical results generally underpredict the temperature.The most significant discrepancy is found mid-way up the cylinder height, which might be either a misprediction of the vapour detachment point and/or an impact of thermocouples wires (Fig. 11b).
The cooling rates plotted against temperature (Fig. 12) provide additional information.A decent agreement with the experiment is observed.The estimated curves do not seem to be shifted in any direction, and the values at the descent from the maximum cooling rates fit well to the experiment.Also, it is noticeable that better agreement is reached deeper in the cylinder (Fig. 12c) compared to the other locations.Perhaps the most erroneous prediction is given at the location J8, corresponding to the finding from the previous temperature-time figure .The final cooling data concerns the top of the cylinder, thus locations  4 and  5.The assessment is started again with cooling curves (Fig. 13).A substantial difference is detected at the central thermocouple,  4 (Fig. 13a).The experiment demonstrates film boiling-like behaviour due to vapour insulation.In contrast, the numerical solution indicates much faster cooling.The immediate reasoning points to the holder structure, present and mentioned in Section 3, but more is discussed together with results shown in Fig. 15.
The second thermocouple,  5, placed radially from the centre (Fig. 13b), provides substantially better numerical prediction.A reason for the existing discrepancy may have much in common with the neighbouring location.If that was insulated, heat would be conducted to  5, hence the higher temperature.Notice that the experimental results at  4,  5 manifest significantly slower cooling, which might result from lack of vapour as will be seen later (Fig. 15).
The cooling rates results provide the very same information (Fig. 14).They further confirm that the highest computed cooling rates from the whole cylinder are present at the top surface centre (Fig. 14a).The value of about 225 • C∕s is far above any other investigated location.For example, at location J5, the maximum predicted cooling rate is almost half (Fig. 14b).Besides that, the simulations captures a phenomenon of reduced cooling rate at about 900 • C followed by an immediate increase.The event is forecast at a slightly lower temperature but qualitatively well.The final part is devoted to a visual investigation of vapour, total velocity and solid temperature's spatial and time distributions.Fig. 15 depicts the vapour volume fraction field near the hot cylinder at various times along the cooling process.It begins with submerging when liquid advances up the cylinder.The fluid phases move violently, the air is partially trapped under the cylinder, and vigorous boiling appears.At about 2 s, the cylinder is already fully immersed, yet large bubble crowds and vapour columns are repeatedly being developed and released from its surroundings.After a few seconds, a different flow pattern emerges.Only smaller bubble crowds are released at about 8.6 s.Later, at 22 s, vapour keeps developing and flowing alongside the cylinder.However, no patches of high vapour volume fraction are visible.The only exception is the cylinder bottom, where the surface still seems to be covered by the vapour layer.Finally, the vapour generation is relatively tiny at about 55 s, but nucleate boiling is still the dominant HTR.Concerning the misprediction at the cylinder top surface, a liquid phase is always visible at the centre location.The vapour is not capable of disturbing the central liquid column nor trapping the liquid, so it evaporates entirely.It is hypothesised that the collars and forklift generate extra vapour and obstruct fluid movements.An upward-facing horizontal surface unimpeded by obstacles would allow free vapour detachment.Indication of this phenomenon was recorded during the experiment.Fig. 16 shows the same view at the computational domain as Fig. 15, yet the vapour volume fraction is swapped for streamlines of total velocity   =     +     .Before the complete cylinder immersion, at 1.2 s, the flow pattern is dominated by vapour, which rises alongside the cylinder and creates a vortex above it.The energy is advected resulting from buoyancy forces in vapour and its relative velocity to the specimen.Later, hand in hand with liquid surrounding the cylinder surface, the vapour is restricted in channels and columns with diverging and converging cross-sections resulting in velocity variations.Due to the development of the vapour columns, the flow patterns change frequently, and a large variety of vortices appear, but most often, a    vortex above the cylinder is notable.Finally, the flow calms down because less and less vapour is generated, leading to smaller velocity gradients and the disappearance of vortices.

Summary
This paper investigated water quenching of IN718 nickel-based superalloy, in the form of a solid cylinder, by both experiment and simulation.This is because quenching is a must have process to optimise the size and distribution of precipitates in this class of alloys to tailor the mechanical properties for demanding applications.It has served as a quenching methodology validation case following [10].Besides the validation, the work consisted of experimental analysis, mesh sensitivity study and discussions.The following conclusions can be drawn: The experimental analysis demonstrated that the thermodynamic complexities heavily depend on the hot surface orientation.The solid temperature near the interface may develop a distinct cooling pattern when vapour is present and not allowed to flow freely.Additionally, it has been shown that the testing facility design may greatly affect the cooling.The cylinder collars and thermocouple wiring act as examples.
The mesh sensitivity study confirmed mesh independence to solid refinement, yet also difficulty to converge the fluid mesh.The fluid mesh becomes independent of the mesh only at the surfaces where sufficient mixing appears, hence at the top and vertical cylinder surfaces.The difficulty arises at the bottom, where the vapour movement is restricted, and its volume fraction reaches considerable levels.Furthermore, the developed numerical method has been validated against the experiment.A good agreement was reached at the cylinder's vertical surface.At the bottom, the numerical results agreed well with the experiment, yet only at the cylinder axis close to the surface.Closer to the cylinder's bottom edge, rather large differences were observed.The last cylinder surface, the top, showed exactly the opposite behaviour, with relatively satisfying cooling results close to the cylinder edge but a considerable discrepancy at its axis due to reasons explained previously.
Finally, wall boiling mass source and bubble dynamics (breakup, crowding, coalescence) play important roles.The computational methodology utilised within the research has dominantly been used for forced and vertical flows.However, the circumstances are very different in water immersion quenching without agitation.In the authors' opinion, that leads to an even higher importance of bubble dynamics than at the previously mentioned conditions.Appropriate bubble dynamics implementation would result in more reasonable prediction of interfacial terms and, ideally, improvements at the cylinder bottom and top surfaces or better fluid mesh behaviour.Nevertheless, the wall boiling mass source term directly related to the boundary cell size always seems to be an issue and is one of the main points for future work.
In conclusion, our computational approach innovatively models the entire quenching process with high accuracy.This innovation stems from integrating models that were previously developed independently.This integration is crucial for ensuring energy conservation across solid-fluid interfaces and for maintaining numerical stability across the simulation process, thereby yielding results that are both accurate and realistic.Our future objective is to develop computational tools capable of accurately predicting the Leidenfrost temperature directly from first-principles simulations, eliminating the reliance on empirical data.Furthermore, we underscore the necessity for ongoing development and validation of boiling models, particularly under conditions of no agitation and minimal fluid velocity, to enhance the fidelity and applicability of quenching process simulations.

Fig. 1 .
Fig. 1.(a) A photograph of the quenched cylinder with two collars, and (b) the thermocouples' locations, highlighted with orange dots, used for measurements of cooling curves during quenching.

Fig. 2 .
Fig. 2. Photographs of the cylinder quenching experiment with approximate time related to the thermocouples measurements.

Fig. 3 .
Fig. 3. Analysis of the experimentally measured cooling curves acquired at locations  1 and  8. (a) Cooling curves at the period 31 to 100 s; (b) Cooling rates at the period 31 to 100 s;   -Peak Heat Flux;  -Minimum Heat Flux.

Fig. 4 .
Fig. 4. Computational domain with dimensions in mm and boundary conditions names.
• C. The  BC behaviour depends on the flow direction.When an outflow occurs, mass flux is positive, zero gradient BC is employed for all momentum and energy equations.While if the International Journal of Heat and Mass Transfer 235 (2024) 126158

Fig. 5 .
Fig. 5. Visualisation of the computational domain combining meshes with 31,140 fluid (coloured in cyan) and 10,440 solid elements (coloured in grey).Pictures are truncated and not in scale.(a) Side view (b) Top view.

Fig. 9 .
Fig. 9. Comparisons between experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling curves during the initial 25 s are visualised at locations (a)  1, (b)  2, (c)  10.

Fig. 10 .
Fig. 10.Comparisons between experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling rates as a function of temperature are visualised at locations (a)  1, (b)  2, (c)  10.

Fig. 11 .
Fig. 11.Comparison between experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling curves at the initial 25 s are visualised at locations (a)  6, (b)  8, (c)  7, (d)  9.

Fig. 12 .
Fig. 12.Comparison of experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling rates as a function of temperature are visualised at locations (a)  6, (b)  8, (c)  7, (d)  9.

Fig. 13 .
Fig. 13.Comparison of experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling curves at the initial 25 s are visualised at locations (a)  4, (b)  5.

Fig. 14 .
Fig. 14.Comparison of experimental and numerical results acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The cooling rates as a function of temperature are visualised at locations (a)  4, (b)  5.

Fig. 15 .
Fig. 15.Spatial and time distribution of vapour volume fraction and solid temperature acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The figures represent only a segment of the computational domain.

Fig. 16 .
Fig.16.Spatial and time distribution of total velocity   =     +     streamlines, coloured by its magnitude, and cylinder temperature acquired using meshes with 59,008 and 451,500 computational elements for fluid and solid, respectively.The figures represent only a segment of the computational domain.

Table 1
can be substituted with heat flux representing particular heat transfer regime   , ,    ,   or   , .Subscripts  ,  ,  and  label film boiling, transition boiling, nucleate boiling and convective HTR, respectively.  stands for   , or   , .Heat transfer regimes partitioning.Subscripts  ,  ,  and  label film boiling, transition boiling, nucleate boiling and convective HTR, respectively, while subscripts  , ,  represent saturation, departure from nucleate boiling and Leidenfrost temperatures, in the given order.Superscript  stands for the interface between the solid and fluids' regions.  <     ≤   <     ≤   <     ≤

Table 3
Cylinder boundary conditions.  ,, stands for mass flux through a boundary cell face.

Table 4
Physical properties of water, vapour and the cylindrical specimen.

Table 5
Fluid and solid regions' mesh characteristics, highlighting the total number of elements, and boundary cells interface perpendicular dimensions at the cylinder top, bottom, and side.