Multi-field coupling simulation of vacuum carburizing and high-pressure gas quenching process

ABSTRACT Vacuum carburizing and high-pressure gas quenching process is being widely used in the heat treatment industry. There is an urgent need for a mature technology to predict the influence of carbon content and quenching gas flow on the microstructure and mechanical performance of workpieces. In this paper, a coupling simulation system of the flow field, temperature field, phase field, and stress–strain field was established. And simulations and experiments were carried out for an 18CrNiMo cylindrical sample and a 20CrMnTi gear hub. The results prove that this system can realize the synchronous coupling simulation of the above physical fields, and the simulation results are accurate and applicable to complex shape parts, therefore can guide the design of the process control scheme of the vacuum carburizing and high-pressure gas quenching process.


Introduction
Vacuum carburizing and high-pressure gas quenching process is a new heat treatment technology which has several advantages including good product surface quality, small shape distortion, low pollution, and it potentially replaces the traditional chemical carburizing and oil quenching process [1,2]. However, the low heat exchange efficiency of gas limits the new technology's applications. Therefore, it is necessary to use numerical simulation approach to study the heat treatment process.
In vacuum carburizing and high-pressure gas quenching process, the vacuum carburizing step mainly involves mass transfer, while the high-pressure gas quenching step involves multi-physical behaviors, including flow of gas, heat transfer between gas and workpiece surface, heat conduction, phase transformation and deformation in workpieces. In most of the previous studies, a single physical phenomenon was simulated. For example, Elkatatny et al. [3] just considered the flow and heat transfer phenomena, and simulated the high-pressure gas quenching process with the isothermal stage using a computational fluid dynamics software. Schmidt et al. [4] and Wang et al. [5,6] carried out similar simulation research. In other studies, the influence of gas flow is simplified as the thermal boundary condition on the workpiece's surface, and only the temperature, structure, and stress-strain field in the workpiece are simulated. For example, Jung et al. [7] and Miao et al. [8] respectively used this method to simulate the high-pressure gas quenching process of H13 tool steel block and 20MnCrS5 helical gear, and analysed the influence of process conditions on the mechanical properties of workpieces.
In recent years, with the improvement of computer performance and the development of simulation technology, researchers tried to simulate the high-pressure gas quenching process with the fluid-solid coupling method. Zhang et al. [9] combined the CFD software Fluent and FEM software Marc, with the help of coupling simulation platform MpCCI, and developed a packet named Thermal Prophet for heat treatment, then accomplished the calculation of phase transformation and strain to simulate the high-pressure gas quenching process of 20CrMo cylinder. The cooling curve and phase structure obtained are roughly consistent with the experimental results. However, the configuration of this coupling simulation system is cumbersome and inefficient, so it is not suitable for simulating the actual workpieces with complex structures and large scales.
At present, most of the simulation studies of vacuum carburizing and high-pressure gas quenching process still use simple models either considering fluid domain or solid domain and low accuracy material parameters, and not closely combined with industrial production. Based on the coupling simulation platform CoSim developed by MSC company and Thermal Prophet, a complete and efficient coupling system was established, including flow field, temperature field, phase field, and stress-strain field. The effectiveness of the system has been verified by the 18CrNiMo cylindrical sample. On this basis, the microstructure and mechanical performance of a 20CrMnTi gear hub processed by the vacuum carburization and highpressure gas quenching were predicted, further proving that the proposed approach is suitable for a complex shape workpiece.

Calculation of carburizing process in solid domain
Isothermal vacuum carburizing mainly involves mass transfer in the solid domain, which obeys the solid mass transfer differential equation. For the carburizing process, the equation can be simplified as: where t is the time, D C is the diffusion coefficient of carbon, and C is the mass fraction of carbon, also known as the carbon content. The boundary condition considered in vacuum carburizing is the carbon transfer coefficient between gas and solid surface.
Coupling calculation of temperature field, phase field, and stress-strain field in solid domain As shown in Figure 1, the temperature field, phase field, and stress-strain field in the solid domain obey certain physical laws respectively and affect each other, resulting in coupling phenomena. All of these physical quantities are affected by the element composition, which is mainly represented in the simulation as the influence of carbon content on material parameters. The solid heat conduction obeys solid heat conduction differential equation: where r is the density, T is the temperature, l is the thermal conductivity, c p is the specific heat capacity, andḞ q is the heat source power per unit volume, which usually comes from latent heat of phase transformation during heat treatment. Pearlite formation and bainite transformation belong to diffusional transformation, whose isothermal transformation kinetics can be described by Avrami equation [10]: where f is the volume fraction of new phase, t is the isothermal time, t 0 is the incubation period, and b, n are kinetic parameters. By using the Scheil superposition principle [11], this equation can be applied to continuous cooling transformation. Martensitic transformation belongs to displacive transformation, whose athermal transformation kinetics is described by Koistinen-Marburger equation [12] commonly: where M s is the starting temperature of martensitic transformation, and a A M is the kinetic parameter.
In the quenching process, the strain of material includes elastoplastic strain 1 ep , thermal strain 1 th , transformation strain 1 tr , and transformation induced plastic strain 1 tp . The increment of total strain is: Elastoplastic strain increment d1 ep and stress increment ds obey the elastoplastic mechanics constitutive equation: where elastoplastic stiffness tensor C ep , additional Figure 1. Coupling phenomena between temperature field, phase field, and stress-strain field.
strain term d1 0 , and additional stress term ds 0 are related to elastic stiffness tensor, yield strength, deviator stress, and plastic modulus. Thermal strain reflects the effect of temperature on atomic equilibrium spacing, and transformation strain reflects the effect of microstructure on atomic arrangement. These two phenomena can be considered together: where f k is the volume fraction of some phase, and b k is the linear expansion coefficient of this phase, which is a function of temperature. When phase transformation occurs in undercooled austenite, if a small tensile load is applied to the material, it will have additional plastic deformation with the phase transformation progress, which is called transformation induced plasticity. Transformation induced plastic strain is described by Desalos model [13] commonly: where K is the coefficient of transformation induced plastic strain, and s ′ is the deviator stress.
Hardness is often used to characterize the mechanical performance of workpieces. In the simulation, the influence of microstructure is generally regarded as a simple superposition of each phase, and the hardness of each phase at room temperature is determined by transformation temperature and transformation cooling rate, so the incremental form for calculating Vickers hardness is: Calculation of flow field in fluid domain The flow and heat transfer of fluid obey the conservation law of mass, momentum, and energy. The general control equation of fluid is: where u is the flow velocity, f is a general variable, G is the generalized diffusion coefficient corresponding to f, S is the generalized source term corresponding to f. Turbulence flow occurs in the highpressure gas quenching, which is difficult to calculate directly due to its complex structure and drastic changes. Therefore, the concept of timeaverage is introduced, that is, the instantaneous value of physical quantity f can be regarded as the sum of time-averaged value f and fluctuating value f ′ . Combining this concept with the momentum conservation equation, Reynolds averaged equations were raised: where p is the pressure, m is the diffusion coefficient. Based on this equation, several turbulence models were established [14,15].
Coupling simulation of fluid-solid heat transfer The boundary coupling method was used to deal with the fluid-solid coupling heat transfer, whose basic idea is establishing equations for the solid domain and the fluid domain respectively, and calculating the two domains alternately in each time step. When calculating one domain, the influence of the other domain is regarded as boundary conditions. In each time step, firstly, the interface temperature of the solid in the previous time step is obtained and transferred to the fluid domain through interpolation, used as boundary condition to calculate the flow field and temperature; then the interface temperature and equivalent heat transfer coefficient of fluid are obtained and transferred to the solid domain, used as boundary condition to calculate the temperature field, phase field, and stress-strain field.
With this method, the synchronous coupling calculation of fluid and solid was realized, and the influence of the flow field on the cooling effect of workpieces was revealed.
Multi-field coupling simulation system Figure 2 shows the structure of the multi-field coupling simulation system. In the fluid domain, the CFD software scFLOW was used to calculate the flow field and temperature; in the solid domain, the nonlinear FEM software Marc and packet Thermal Prophet were used to calculate the temperature field, phase field and stress-strain field; between the two domains, the coupling simulation platform CoSim was used to realize coupling heat transfer. In addition, user subroutines were programed to secondary develop scFLOW and Marc to meet specific needs.

Gas quenching simulation and verification of a sample
Vacuum carburizing and high-pressure gas quenching experiment These experiments and simulations aim to prove the multi-field coupling simulation system's function, and verify the accuracy of simulation results. The carburizing and gas quenching furnace used in this experiment is provided by Beijing Research Institute of Mechanical and Electrical Technology, and its geometric structure model is shown in Figure 3(a). The size of the furnace is 1150×710×770 mm, and its wall is made of graphite felt insulation material. There are two cavities at the upper and lower parts, which contain 20 circular gas inlets at the bottom respectively, and gas outlets at the bottom of the side wall. There are heating tubes near the furnace wall and a pair of graphite supports at the bottom of the furnace. The sample is a cylinder of 18CrNiMo, whose length is 100 mm and diameter is 15 mm, installed on the frame with deflectors around it.
The sample is vacuum carburized in the furnace, including carburization stage at 920°C for 42 min and diffusion stage for 145 min. Then, nitrogen is introduced from the upper gas inlets for high-pressure gas quenching. The flow velocity is 33 m/s and the gas pressure is 0.6 MPa.

Determination of simulation parameters
Due to the small size of the sample, the influence of other parts on flow velocity and temperature cannot be ignored. Therefore, based on the research of Pan Jiansheng et al. [16], the expanded solution domain method was used to model the carburizing and quenching furnace. Considering the symmetry of the system, half of the quenching furnace was modeled. Figure 3(b,c) shows the mesh model of the sample and the furnace. To accurately simulate the carburizing process, the mesh near the sample surface was refined.
In terms of the fluid model, the Realizable k-1 model proposed by Shih et al. [15] was applied as the turbulence model, and data from scFLOW material library was used as the material properties of nitrogen. In terms of the solid model, to reflect the effect of the carburization on material properties, material properties of 18CrNiMo and 65CrNiMo was measured, including transformation kinetic parameters, specific heat capacity, thermal conductivity, tensile curve, Vickers hardness, and transformation strain. The transformation induced plastic strain coefficient came from the test results of Qin Shengwei et al. [17], and the other material properties came from the calculation results of the software. As for other parts in the furnace, only the thermal  performance was considered, whose data were from scFLOW material library.

Simulation results and analysis
The sampling path of carbon content, microstructure, and hardness described below is located on the section 10 mm away from an end face. Figure 4(a) shows the simulated and measured carbon content of the sample surface after vacuum carburizing, which are relatively consistent. The measured carburized layer depth is 0.86 mm and the simulated carburized layer depth estimates 0.83 mm. Figure 4(b) shows the simulated carbon content distribution. The carburization at other positions is uniform except for the edge of end faces where carbon content is high.
In high-pressure gas quenching process, the flow field changes little with time as the flow boundary conditions are constant. Figure 5 shows the simulated flow velocity and gas temperature distribution near the sample after quenching for 58s. The gas flow from the upper gas inlets bypasses the sample, resulting in boundary layer separation, so the flow velocity on both sides is the highest, and the flow velocity on the top surface is relatively low, while vortex is formed near the bottom surface, making the flow velocity the lowest. At the same time, the temperature of gas near the bottom surface is significantly higher than that in other regions due to the vortex. Figure 6(a) shows the simulated and measured temperature history of the test point in the sample. The simulation results are close to the measured results except for the error between 30-90 s. There are two reasons for this: the simulated gas temperature at inlets is different from the actual value; the temperature sensor in the experiment is in contact with the gas, making the measured temperature lower than real. Figure 6(b) shows the distribution of the simulated sample temperature after quenching for 58s. The upper edge of the sample cools the fastest and the lower part cools slowly, which is the result of surface flow velocity, surface gas temperature and heat conduction of the solid. Figure 7(a,b) shows the simulated microstructure. Due to the rapid cooling rate of the whole sample, the microstructure of each position has no significant relationship with the cooling rate, but has an obvious correlation with carbon content. The high carbon content, on the one hand, prolongs the incubation period of bainite and prevents the occurrence of bainite transformation; on the other hand, decreases the M s and affects the phase transformation sequence in different positions. Therefore, the microstructure of the internal non-carburized region is almost the same, composed of 93% martensite and 7% bainite, while the surface carburized layer is almost completely transformed into martensite. The measured surface metallographic photo (Figure 7(c)) proves this result. Figure 8 shows the components of strain caused by various factors. The distribution of mechanical plastic strain, thermal strain, transformation strain, and transformation induced plastic strain are related to carbon content and microstructure. These four strain components all decreasing from the surface to the inside. Thermal strain and transformation strain differs a little from the surface to the inside, which is mainly affected by the difference of microstructure. Mechanical plastic strain and phase transformation plastic strain are very different at the inside and the surface, which is the result that the phase transformation sequence is controlled by carbon content. In the quenching process, the non-carburized region first completes the martensite transformation, and the martensite transformation occurs on the surface layer after the stress in the internal region has been relieved. The expansion trend caused by the transformation strain is restricted by the subsurface material, making the surface layer pressurized and the sub-surface layer tensioned, resulting in large plastic deformation and transformation induced plasticity.
This effect can also explain the simulated residual stress distribution. As shown in Figure 9 (the deformation in (b) is enlarged by 10 times for the convenience of observation), there are high residual compressive stress in the carburized layer, and residual     tensile stress in the subsurface. Due to different cooling rate and the influence of stress coordination, surface residual stress fluctuates along the circumference of the cylinder. In terms of deformation, in addition to uniform cooling shrinkage, there is no other obvious distortion, expect both ends slight upward warpage. Figure 10 shows the simulated and measured Vickers hardness. The hardness is high in carburized layer due to the fast cooling rate and high martensite volume fraction, and decreases with the increase of depth. In the non-carburized region, the simulation results are close to the measured results, but in the carburized layer the simulation results are relatively low. This error is mainly caused by the simulation method of hardness, which does not count the influence of residual stress, while residual compressive stress can significantly improve the hardness. At present, there is no general theory about the influence of residual stress on hardness, and the hardness simulated here can be used as the lower limit of expected hardness in actual production.
According to the simulation and experiment results, the multi-field coupling simulation system can synchronously simulate various coupling phenomena in vacuum carburizing and high-pressure gas quenching process, and intuitively reflect the influence of carburizing and cooling conditions on the mechanical performance of the sample. Simulation results accord with the reality and can be applied to the prediction of workpiece performance.

Process scheme and determination of simulation parameters
This simulation aims to further verify the applicability of the multi-field coupling simulation system to complex shape workpieces. The gear hub as the research object is made of 20CrMnTi, with a height of 24 mm. The distance from the axis to the internal tooth top is 25.4 mm, and the distance to the external tooth top is 73.15 mm. The carburizing and quenching furnace model in the previous case is still used, but the frame and deflector are removed, and the gear hub is placed directly on the center of the tray. This gear hub is first vacuum carburized in the furnace, including strongly carburized at 915°C for 42 min and diffused for 90 min. Then, nitrogen is introduced from the upper gas inlets for high-pressure gas quenching. The flow velocity is 33 m/s and the gas pressure is 0.6 MPa.
In the simulation, half of the quenching furnace and the gear hub were modeled, and the surface mesh of the gear hub was refined. Simulation parameters of the fluid model were the same as in the previous case. In terms of the solid model, the transformation kinetic parameters, specific heat capacity, thermal conductivity, tensile curve, Vickers hardness, linear thermal expansion coefficient, and transformation strain of 20CrMnTi were measured, and the corresponding material parameters of 80CrMnTi were calculated. Other material properties came from the calculation results of the software.

Simulation results and analysis
The sampling part of carbon content, phase constituents of the microstructure, and hardness described below is located at the central horizontal section of the gear hub, from the tooth surface of the internal and external teeth to the inside. Figure 11 shows the carbon content distribution after vacuum carburizing. The internal and external teeth is basically the same,  and the carburized layer depth is estimated to be 0.73 mm. The carburization at other positions is uniform except for the surface edge and tooth top where carbon content is high. Figure 12 shows the flow velocity on the surface and temperature distribution of the gear hub after quenching for 44s. As the gas inlets are located above the furnace, the flow velocity on the upper surface is higher than that on the bottom surface. Along flow direction, the flow velocity on the internal and external teeth is higher than that on the top and bottom surfaces. Higher flow speed brings markable cooling effect on the temperature field. Figure 13 shows the calculated fractions of the phase constituents of the microstructure. The microstructure is affected by cooling rate and composition. On one hand, the martensite volume fraction is high in the upper part and external teeth due to high cooling rate and carbon content, and low in the lower part and inner ring, while the bainite transformation is inhibited by high carbon content in the surface carburized layer. On the other hand, the non-carburized region in the workpiece is mainly composed of a mixture of martensite, pearlite, and bainite. Figure 14 shows the residual stress distribution in the gear hub and distortion. The residual stress is affected by cooling rate and microstructure. In general, the surface carburized layer has a high residual compressive stress, while the residual stress in the noncarburized region is low. In terms of the distortion, in addition to uniform cooling shrinkage, the external teeth of the gear hub warp upward. In actual production, it is suggested to use the up-down alternating blowing scheme to reduce the distortion.
The simulation results show that the multi-field coupling simulation solution can be effectively used   to study the vacuum carburizing and high-pressure gas quenching process of a complex shaped workpiece. It is capable of predicting the fractions of the phase constituents in the microstructure and mechanical performance of the workpiece under complicate cooling conditions, consequently, offering a useful guidance for the optimized scheme in the production process for distortion control.

Conclusions
(1) Based on the secondary development of existing software, a coupling simulation system of the flow field, temperature field, phase field, and stress-strain field was established, and the multifield coupling simulation of the vacuum carburizing and high-pressure gas quenching process was accomplished. (2) The simulation and experiment of an 18CrNiMo cylindrical sample shows that this system can synchronously simulate the gas flow and cooling, phase transformation and mechanical changes of the workpiece in the vacuum carburizing and high-pressure gas quenching process, reflecting the influence of carbon content and gas flow conditions. The simulated microstructure and surface hardness are close to the actual values, which have high reference values. (3) The simulation of a 20CrMnTi gear hub shows that this system is suitable for the simulation of complex workpieces, especially for predicting the deformation in the vacuum carburizing and highpressure gas quenching process. Therefore, this coupling simulation system has broad application prospects in industrial production.

Disclosure statement
No potential conflict of interest was reported by the author(s).