{'fraction_non_alphanumeric': 0.025333600882913615, 'fraction_numerical': 0.006772348750877897, 'mean_word_length': 4.689212328767123, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heat transfer model for quenching by submerging

In quenching by submerging the workpiece is cooled due to vaporization, convective flow and interaction of both mechanisms. The dynamics of these phenomena is very complex and the corresponding heat fluxes are strongly dependent on local flow variables such as velocity of fluid and vapor fraction. This local dependence may produce very different cooling rates along the piece, responsible for inappropriate metallurgical transformations, variability of material properties and residual stresses. In order to obtain an accurate description of cooling during quenching, a mathematical model of heat transfer is presented here. The model is based on the drift-flux mixture-model for multiphase flows, including an equation of conservation of energy for the liquid phase and specific boundary conditions that account for evaporation and presence of vapor phase on the surface of the piece. The model was implemented on Comsol Multiphysics software. Generation of appropriate initial and boundary conditions, as well as numerical resolution details, is briefly discussed. To test the model, a simple flow condition was analyzed. The effect of vapor fraction on heat transfer is assessed. The presence of the typical vapor blanket and its collapse can be recovered by the model, and its effect on the cooling rates on different parts of the piece is analyzed. Comparisons between numerical results and data from literature are made.


Introduction
Quenching is a technological process that involves cooling of a heated workpiece in order to modify its mechanical properties. Depending on the material of the workpiece, the cooling rate needs to be very fast or can be mild. In particular, for low to medium alloy steels the cooling rates should be as fast as possible.
Focusing on steels, quenching produces metallurgical transformations that are responsible for modifications in mechanical properties (hardness, yield and ultimate stresses and toughness). This feature allow us to tailor the steel properties in a broad range that is not easy to replicate in common engineering materials.
In addition to metallurgical transformations that lead to mechanical modifications, geometrical distortions and generation of residual stresses can be developed during the process. In order to get a successful quenching, the final metallurgical phases and mechanical properties should be the desired ones, and distortions and residual stresses should be kept at a minimum.
As any other process of technological interest, several efforts were made in order to describe it through mathematical models. In this area it is important to remark that quenching of steels can be divided in at least three sub-problems. One major problem is the modeling of metallurgical transformations. This part is well known and the available models are robust and reliable [1].
On another side we have the generation of residual stresses and geometrical distortions. This complex phenomena are also well covered, even though new refined models can be found in specialized literature. Finally, these two sub-problems strongly rely on an accurate description of the temperature evolution inside the piece. Several attempts to put together a thorough description of the quenching process can be found in literature, see refs. [2]- [6]. These and other works take this approach, where each sub-problem is solved using the most adequate model. But, the thermal problem is usually solved taken a very rough description of heat exchange with the fluid through the boundaries conditions. It is important to remark that all the complexity and accurateness of each sub-model is lost if the thermal problem is over-simplified. In this work the analysis of heat transfer in the workpiece-fluid interface is tackled. The microstructure and residual stresses problems will be coupled in further works.
In quenching by immersion the piece transfers heat to the quenching media through vaporization of the fluid, convective flow and interaction of both mechanisms. Generalities about subcooled boiling stages can be found in [7]- [9]. The description of this phenomenon has to take into account the dynamic of a multiphase vapor-liquid flow and specific heat transfer contributions due to vaporization. The problem of multiphase flows with heat exchange is described in general terms in works such as [10]- [12]. Based on these models and including the heat partition model for boiling, a heat transfer model for quenching by submerging is developed here. The model includes the presence of vapor phase and how it affects the heat transfer from the piece to the quenching media.

Model
The analysis of multiphase flows presented here was based on the general methodology presented in [10]- [12]. Before the definition of the particular multiphase model it was necessary to analyze coupling parameters between phases, characteristic times and finally characteristic nondimensional numbers of the specific problem.
Different characteristic times were considered, such as those related to: • Mass transfer, τ M , the ratio between the latent heat of condensation of a bubble and the rate of heat transfered from it to the fluid. • Momentum transfer, τ V , which is the characteristic response time of a bubble.
• Energy transfer, τ T , the ratio between the sensible heat of a bubble and the rate of heat transfered from it to the fluid.
Those characteristic exchange times were compared to the characteristic residence time of the flow, τ r . These comparisons allow us to assess if each quantity is rapidly transfered from one phase to another, therefore reaching equilibrium, or if the dynamic of the transference needs to be described. Under our assumptions of dispersed vapor (say a fraction of liquid, α l , higher than 0.55) and with a bubble diameter (D b ) ranging from 10 −4 to 10 −3 m, all the characteristic times are remarkably smaller than the mean flow residence time. Based on the characteristic transfer times of each quantity, the corresponding coupling parameters were analyzed according to [10]. The obtained parameters indicated that mass and energy transfered due to phase change is highly coupled between phases, while momentum and energy transfered due to temperature differences have a negligible coupling.
Finally, the non-dimensional analysis showed that the drift number (N D ) and the phase change number (N ph ) are in the same order of magnitude, and the Eckert (N Ec ) number very small (see appendix for a detailed explanation of the meaning of these characteristic numbers). These results imply that the description of velocities for both phases is needed, and that the interfacial terms in the equation of energy are negligible.
Based on these results, for our particular problem, the main conclusions that were drawn are: • Liquid and vapor phases are in mechanical equilibrium, therefore the velocities of both phases are strongly correlated. • The heat delivered from vapor to liquid corresponds mainly to latent heat of condensation.
Heat corresponding to cooling of vapor from temperatures higher than saturation temperature are negligible. • The characteristic time of cooling of vapor phase is very small. Based on this conclusion and the previous one, the vapor temperature can be assumed constant at saturation value. • Both phases are strongly coupled through exchanges of mass and energy due to phase change.
Based on these conclusions, the most convenient multiphase model selected was the drift-flux mixture-model. This model is easier to solve than the fully detailed Euler-Euler model, and captures all the characteristics of the two-phase flow. Due to the constant vapor temperature assumption, only the conservation of energy of the liquid phase has to be described. The system of partial differential equations (PDE's) that describe this problem is presented in section 2.1. Subindexes m, v and l correspond to mixture, vapor and liquid respectively. The fraction of each phase is represented by α and the and symbols correspond to phase average and density-weighted average respectively.

Conservation equations
The multi-phase problem derived in previous section can be modeled by means of four main conservation equations, which are [12]- [13]: This model needs the description of the drift velocity (V vj ) and the mass transfer (Γ v ) between the phases. The selected expression for V vj obtained from [12] is presented in Eq. (5). Is important to note that the drift velocity only acts in the same direction of the gravity.
The rest of the closure relationships in order to get a well-posed problem were taken from [12]- [22]. In addition, the k − turbulence model for the liquid phase (taking into account the presence of disperse vapor) was used to describe the Reynolds stress tensor T T . The evolution of the temperature in the solid is solved by a separate heat conduction equation.

Boundary conditions
The most common ways to estimate the heat flux at the boundaries of the workpiece are: • by means of a heat transfer coefficient (h) depending on the global temperatures of the piece and cooling fluid, as it is presented in [2] - [5] or, • calculating the heat flux emitted by the piece (q w ) using the heat flux partition model, which is a function of global values of piece and fluid temperatures and fluid velocity, in a similar way as [17] - [22].
Those approaches are completely insensitive to local variations of fluid temperature and velocity or the presence of vapor in the fluid. In this work we propose to include local dependences in the formulation of the boundary conditions. We take the heat partition model as a base and expand it in the way explained in following paragraphs. It is important to remark that, up to this point, is not our aim to describe early stages such as vapor blanket formation. That specific phenomenon will be included in further studies. The heat transfered during boiling was described in a similar way than it is presented in [15]- [17], [19]- [22]. It is important to remark that due to the fact that our cooling fluid is oil, the functions presented for water were replaced for expressions based on physical properties of the fluid, rather than fitting constants [21]. The total heat flux extracted from the solid (q w ) result as the addition of three different mechanisms, as it is presented in its original form in Eq. (6), q 1f corresponds to the heat flux devoted to heat the fluid, while the others represent boiling.
The coupling between the temperature of the solid, the temperature of the liquid and the vapor fraction generated due to boiling is made in the interface solid-mixture. The heat fluxes related to evaporation (q e , q q ) are responsible for the flux of vapor generated, which in this case is considered injected through the boundary layer on the surface of the solid to the quenching flow. This balance of energy is represented in Eq. (7).
Due to the local characteristic of the formulations used in this work, the heat flux associated to the heating of the liquid (q 1f ) is obtained by means of a wall-law formulation instead of the global correlations typically used. The wall temperature was modified in order to capture the presence of vapor. The wall-law and the modified temperature are presented in Eqs. (8) and (9).
The relationship presented in Eq. (7) is exact but do not capture the effect of thermal insulation that produces a vapor blanket surrounding the piece. In order to damp the vapor generation in the presence of a considerable vapor fraction on the piece, the following modification on the injected vapor flux is proposed in Eq. (10).
According to this, the vapor injection is stopped if there is no liquid available. Therefore, the heat extraction due to q q and q e is also stopped. This physical aspect of the boiling phenomena usually is not accounted for in the current models. The final form of the modified heat partition model is presented in Eq. (11), taking into account that q 1f is obtained through wall-law formulations. The damping function presented here is a very basic first proposal and further analysis on it will be performed on future work.

Case studied
The case of a sphere submerged in an upward fluid flow was used to test the model developed here. The upward flow intends to represent the forced-flow condition of an industrial bath. The flow velocity was set at 1 m/s, based on characteristic flow velocities obtained during the specific analysis of the model.

Materials and geometry
The fluid was quenching oil, whose properties were obtained from [23], and the material of the sphere was Inconel 600. This alloy was selected in order to make experimental results form [23] and numerical simulations comparable. Summary of oil, its vapor, and Inconel properties is presented in Table 1. The initial temperatures of fluid and piece were 313 and 1150 K respectively. The domain was described on the r − z plane taking axial symmetry along the z-axis. A sketch of the domain, showing the spherical piece and the surrounding fluid is presented in Fig. 1. The sphere was 30 mm of diameter and the flow domain was a 200 mm long and 200 mm diameter cylinder. The sphere was placed in the centerline, at 60 mm from the bottom. The fluid domain extends in the back of the piece in order to allow an appropriate development of the wake.

Numerical solution
The PDE system presented in Eqs. (1 -4) can be assimilated to a regular fluid mechanics problem, with some extra sources, plus two additional conservation equations, one for the vapor phase and another for the liquid phase energy. Comsol Multiphysics allows the resolution of  [12]. For the conservation of energy on liquid phase, the extra terms due to vapor phase were included as user-defined sources. The resolution of the whole problem was made using the direct solver UMFPACK [24] and the time integration was made along 10 seconds using an implicit scheme of second order. The boundaries of the domain marked in Fig. 1 correspond to: (1) inlet of oil liquid at constant temperature, (2) sliding wall, (3) outlet, (4) symmetry axis, (5) dynamic and thermal wall-laws, and vapor injection due to evaporation. The mesh was generated in Comsol using triangular elements. The maximum length of the elements at each boundary was: 1 × 10 −3 for (1-4) and 3 × 10 −4 for (5). The maximum length in the domain was 3 × 10 −3 , with a growth ratio of 1. Numerical diffusion was needed to stabilize the finite element scheme due to the highly convective behavior of the conservation equations of vapor phase and energy of liquid. Isotropic diffusion was used in turbulence, dispersed phase and energy of the liquid; while streamline diffusion (anisotropic in the streamline direction for the dynamics of the flow and SUPG for thermal problems) was used in the equations. Diffusion coefficients were kept at a minimum (δ < 0.5, where δ is a tuneable parameter that controls the amount of artificial diffusion).
The dynamics of the initial boiling events are represented by characteristic times different than the ones studied here. Because of this reason, the model cannot cope with an initial condition of zero vapor and develop the whole boiling process. As first approximation to solve the problem, an initial distribution of vapor phase corresponding to the initial temperature jump was considered. The description of the initial stages of boiling and the development of the vapor  blanket is a much more complicated process that will be studied in further works.

Vapor blanket and temperature fields
The evolution of vapor around the sphere and its temperature field is presented through snapshots at three different times. In Fig. 2 the early stages of the cooling are presented.
Most of the surface of the piece is densely covered by vapor. In Fig. 3 the partial collapse of the vapor blanket and its thermal insulation effect is observed. In Fig. 4 the last stages of presence of vapor is depicted. Except for a small region in the wake that is still boiling, the temperature looks radially distributed.

Cooling curves, heat transfer coefficient and heat fluxes
In order to assess the local effect of the heat transfer model four different points of the sphere were analyzed. They are located: one at the center of the sphere (center ); and three at 1 mm below the surface, two along the symmetry axis (upper and lower ), and one in the horizontal plane (lateral ). The evolution of the temperature and the corresponding cooling rate at each point is presented in Fig. 5. There it is observed that lower and lateral zones had similar evolutions, while the upper point is clearly delayed. This difference is even more evident when the cooling velocities curves (T vs. dT /dt) are considered. These results highlight the effect of vapor as thermal insulator and how the heat transfer is locally affected along the surface of the piece. Time and spatial variation of the heat transfer coefficient (h) is depicted in Fig. 6, where the lower part of the circumference is the left side of the axis. This behavior is unlikely to be captured by the typical correlations of h depending on mean parameters of the flow. The evolution of heat fluxes on the lower, lateral and upper points as a function of the overheating temperature of the wall was also made (see Fig. 7). It is observed that for lower and lateral points, the heat flux recovers the trend of boiling heat transfer curve. Differences in the shape of the curves may be due to the geometry of the piece (sphere vs. plane surfaces/wires) and also the need of fine tunning of the model (damping of heat transfer due to the presence of vapor in the boundary layer of the piece).

Comparison to experimental data
Results for the center point are compared to data extracted from experimental data presented in [23]. Both geometries are different in shape but with similar volumes, and oil and alloy properties are fairly the same. This allowed us just to compare the order of magnitude of the numerical results. In Fig. 8 the evolution of temperatures and cooling velocities curves are presented. The shapes of the curves are not exactly the same, but the matching of the order of magnitude obtained for both curves is highly satisfactory. Differences at high temperature are attributed to the heat extracted due to latent heat of evaporation during early stages of boiling, that part of the process could still not be captured, as it was stated in the initial conditions discussion.

Summary and conclusions
Heat transfer during quenching should be described including boiling phenomena; this should be done through an appropriate multiphase flow model and adequate local description of heat fluxes exchanged with the treated workpiece. The reduction of the specific model from a general frame depends on characteristic values of the problem. A model that couples vapor dynamics, energy conservation of the liquid phase and heat transferred from the quenched piece was presented. Our model is based on the mixture model for multiphase flows and a heat partition model for heat transfer. Vapor fraction evolution and heat fluxes are obtained as main results.
The effect of thermal insulation of the vapor phase is clearly pointed out. Local analysis shows that cooling velocities are strongly correlated to the presence of vapor, and the resulting heat transfer coefficients are quite complex to be described using regular heat transfer correlations.
The developed model adequately captures the order of magnitude of the physical problem. • Eckert Number (N Ec ): ratio between the characteristic kinetic energy of the mixture vs. the characteristic enthalpy jump between the phases.
∆h lv0 • Phase change Number (N ph ): ratio between the characteristic momentum due to change of phase vs. the characteristic momentun of the vapor phase in the mixture.