2.2. Numerical model
The present study establishes co-simulation of two commercial software, along with secondary adaptation scripts to perform numerical simulations of mobile induction hardening. Electromagnetic simulation was performed using the eddy current solver in ANSYS Maxwell and the thermal-mechanical-metallurgical simulation in the induction hardening process is executed in QFORM UK 3D. The effect of scanning speed is taken account in the calculation by evaluating iteration in each position of the inductor coil. The computation data exchanges of each solver are illustrated in Fig. 6. There are still several physics phenomenon in the real induction hardening process which are not modeled in this co-solver as the limitations in this study, such as additional heat generation due to plastic deformation, stress induced phase transformation, and influence of microstructure on electromagnetic field.
2.2.1 Electromagnetic-temperature model
Electromagnetic solver in this model is based on the Faraday’s law of electromagnetic with a modification to consider the effect of moving inductor as proposed by P. Karban and M. Donátová [11] in Eq. (1).
$$\nabla \text{x}\left(\frac{1}{{\mu }}\nabla \text{x} \text{A}\right)+{\sigma }\frac{\partial \text{A}}{\partial \text{t}}-{\sigma }\text{v} \text{x} \text{A}={\text{J}}_{\text{i}\text{n}\text{d}}$$
1
The first term represents the behavior of magnetic vector potential A in the certain medium with permeability µ. The second term expresses the effect of changing magnetic field inducing eddy currents along material with conductivity σ due to the time derivative of A. The third term is the advection of magnetic field due to the motion of inductor coil with speed v. Meanwhile Jind is the current density of induction coil. Solving this equation with input data and boundary conditions presented in Table 3 for the computation of the magnetic vector potential, then Eddy currents in each node is simply calculated as Eq. (2). where j is imaginary unit and f can be considered as the alternating current frequency.
$${\text{J}}_{\text{e}\text{d}\text{d}\text{y}}=\text{j} \text{x} \left(2{\pi }\text{f}\right){\sigma }\text{A}$$
2
The important result in this study is the heat generation induced by electromagnetic solver to be exported to the temperature model solver to calculate the temperature history in each node. Total of heat generation due to electromagnetic excitation is based on the calculation of ohmic loss as shown in Eq. (3) on the volume of gear tooth model.
$${\text{Q}}_{\text{e}}=\frac{{\left|{\text{J}}_{\text{e}\text{d}\text{d}\text{y}}\right|}^{2}}{{\sigma }}$$
3
During heating process, this amount of heat generation is exported to the temperature model and solve the equation that follows Fourier’s heat transfer in Eq. (4) to obtain temperature history in each node.
$$\nabla . \left({\lambda } \nabla \text{T}\right)={\rho }\text{c}\left(\frac{\partial \text{T}}{\partial \text{t}}+\text{v} . \nabla \text{T}\right)-{\text{Q}}_{\text{e}}$$
4
where \({\lambda }\), ρ, and c are thermal conductivity, density, and specific heat capacity of material, respectively. The temperature history of each node is exported as a feedback to the electromagnetic solver to recalculate heat generation again after considering the effect of temperatures on the gear tooth model where the electrical properties of material are temperature-dependent as shown in Fig. 9c adopted from reference [24]. This looping iterations are terminated when \({\text{T}}_{\text{i}+1}<{\text{T}}_{\text{i}}\) in the heating simulation and in the cooling simulation is decided by the time of stop condition.
Table 3. Boundary conditions of electromagnetic simulation
During cooling or quenching simulation, the convectional heat transfer equation as shown in Eq. (5) is incorporated to the Eq. (4).
$$-{\lambda } . \frac{\text{d}{\text{T}}_{\text{s}}}{\text{d}\widehat{\text{n}}}={\alpha }({\text{T}}_{\text{s}}-{\text{T}}_{\text{m}})$$
5
The convectional heat transfer comes from two sources, the first is the air-medium room which is set at 200 C, 0.6, and 30 W/(m2K) for temperature, emissivity, and heat transfer coefficient, respectively. The second source is a set of water sprayer that moves following inductor coil which consists of eight single nozzles with 5 mm-vertical and horizontal pitch distance, 2 mm-nozzle diameter, 150-nozzle angle, 80 mm-height, 20-attenuation coefficient of the jet, and 400 C-temperature as shown in Fig. 7. The temperature-dependent heat transfer coefficient data of water medium and the implementation to the Eq. (4–5) are adopted from our previous study [23]. To create smooth and accurate numerical calculations in the both electromagnetic and temperature solvers, the tetrahedral finite elements fining in the surfaces are executed for each iteration. In the electromagnetic solver, the remeshing is executed on the certain depth of surface that follows skin effect equation as shown in Eq. (6) [8], then this mesh refinement broadens to areas with massive temperature changes as shown in Fig. 8.
$$\delta =\sqrt{\frac{1}{\pi f\mu \sigma }}$$
6
2.2.2 Temperature-phase transformation model
As shown in Fig. 6, there are two kind of computation data exchanges in this multi physics co-solver. The first is exporting the historical temperature (temperature vs time) data of each node in the temperature solver to the phase transformation solver to then calculate the percentage of each phase transformed and the second is the additional heat generation in Eq. (4) obtained from the phase transformations process. The temperature history of each element determines the percentage of transformed child phase from mother phase based on the transformation kinetic model. For austenitization (initial material-to-austenite phase transformation) during heating process, the numerical computation is modeled by an assumption where the real process was the isothermal precipitation kinetics from a homogeneous supersaturated solid solution, hence the Johnson–Mehl–Avrami–Kolmogorov (JMAK) equation [25] as expressed in Eq. (7) is implemented.
$${{\xi }}_{\text{A}\text{u}\text{s}}\left(\text{t},\text{T}\right)=1-\text{e}\text{x}\text{p}\left\{-{\left[{\text{D}}_{0}.\text{e}\text{x}\text{p}\left(\frac{-\text{Q}}{\text{R}{\text{T}}_{\text{a}\text{b}\text{s}}}\right)\right]}^{\text{n}}\right\}$$
7
where D0, Q, R, Tabs, and n are rate constant, (effective) activation energy, gas constant, absolute temperature, and Avrami exponent, respectively. The values of each parameter were adopted from references where D0 = 2.84 x 1020 s-1, Q = 453811 J/mol and n = 1.525 [26]. During quenching simulation, the percentage of transformed austenite stored in each node decomposites to the percentage of child phases, such as ferrite, pearlite, bainite, and martensite based on the cooling rate occurs in the node. The ferrite, pearlite, and bainite transformations are determined by temperature-time transformation (TTT) data which are collected from interpolating experimental data from Saunders et al. [27] based on the current chemical composition of material. However, the numerical computation still need the data of the start-transformation temperatures for each phases and they are adopted from interpolating numerous experimental testing in the alloy heat treatment by Trzaska [28]. Meanwhile the percentage of martensite transformation is numerically calculated follows a model from Koistinen and Marburger [29] as shown in Eq. (8).
$${{\xi }}_{\text{m}\text{a}\text{r}\text{t}}\left(\text{T}\right)=1-\text{e}\text{x}\text{p}\left(-{{\alpha }}_{\text{K}\text{M}}({\text{T}}_{\text{m}\text{a}\text{r}\text{t}}-\text{T}\right)$$
8
where \({{\alpha }}_{\text{K}\text{M}}\) and \({\text{T}}_{\text{m}\text{a}\text{r}\text{t}}\) are Koistinen-Marburger parameter and martensite start temperature, respectively. The values of those parameters are obtained from reference [30] which \({{\alpha }}_{\text{K}\text{M}}\) and \({\text{T}}_{\text{m}\text{a}\text{r}\text{t}}\) are 0.017 K-1 and 351.110C, respectively. Further illustration figures in kinetics transformation and properties of each 100%-phases are provided in our previous study [23]. However, there are some updates obtained from recent studies [30, 31] regarding the properties of SCM440 material as presented in Fig. 9a-b, while the hardness values of each phase as shown in Fig. 9d are obtained from another Trzaska’s study [32] with considering the effect of temperature.
During phase transformation simulation, mixed phases with different percentage and properties values are possibly stored in each node. The node’s properties are determined using a proportional mixing rule. For example, determining the hardness distribution of each node which is the main focus in this study are simply expressed as in Eq. (9).
$$\text{d}\text{H}{\text{V}}_{\sum }=\sum _{\text{j}=1}^{{\text{m}}_{{\Phi }}}\text{d}{{\xi }}_{\text{j}}\left(\text{t}\right).{\text{H}\text{V}}_{\text{j}}\left(\text{T}\right)$$
9
where \(H{V}_{\sum }\) and \({HV}_{j}\left(T\right)\) are hardness of a mixture phase and j-th temperature-dependent phase respectively.
For heat generation due to phase transformation process in this model follows an equation provided in the reference [30] for latent heat in the C-Mn steels. Therefore, by inputting the weight percentage of carbon value of 0.415, then Eq. (4) obtains an additional heat generation of + 168.07 kJ/kg and − 168.07 kJ/kg during heating and cooling simulation, respectively.
After quenching process, the procedural standard requires an additional heat treatment process, namely tempering process with tempering temperature and time of 1500C and 1 hour, respectively. This additional process aims to improve the ductility of material and reduce the distortion. However, this process generates a new phase presence called tempering martensite and surely changes the hardness of material. Therefore, this process should be simulated to compare the material’s hardness with the real final processed specimen. The diffusion transformation kinetics model as expressed in Eq. (7) can be implemented to predict the percentage of transformed tempering martensite. The constant parameters are obtained from interpolating the material’s chemical composition with numerous experimental results in reference [33]; D0 = 2.75 x 107 s-1, Q = 122212 J/mol, and n = 0.1. The hardness of 100%-tempering martensite also can be extracted in that study as 248.66 HV.
2.2.3. Temperature-stress/strain model
Massive temperature changes in the workpiece will cause volume expansion or shrinkage in the workpiece and cause changes in strain, this mechanism is called the strains due to temperature changes (\({\dot{{\epsilon }}}_{\text{t}\text{h}}^{\text{i}\text{j}}\)). This is one form of strain in the heat treatment process, besides that there are strains due to elastic behavior (\({\dot{{\epsilon }}}_{\text{e}}^{\text{i}\text{j}}\)), strains due to plastic behavior (\({\dot{{\epsilon }}}_{\text{p}}^{\text{i}\text{j}}\)), strains due to phase transformation (\({\dot{{\epsilon }}}_{\text{t}\text{r}}^{\text{i}\text{j}}\)), and strains due to phase transformation induced plasticity (\({\dot{{\epsilon }}}_{\text{t}\text{p}}^{\text{i}\text{j}}\)), thus the total increment of the strain tensor is expressed as Eq. (10).
$${\dot{{\epsilon }}}_{\sum }^{\text{i}\text{j}}={\dot{{\epsilon }}}_{\text{t}\text{h}}^{\text{i}\text{j}}+{\dot{{\epsilon }}}_{\text{e}}^{\text{i}\text{j}}+{\dot{{\epsilon }}}_{\text{p}}^{\text{i}\text{j}}+{\dot{{\epsilon }}}_{\text{t}\text{r}}^{\text{i}\text{j}}+{\dot{{\epsilon }}}_{\text{t}\text{p}}^{\text{i}\text{j}}$$
10
The strains tensor due to temperature changes \({\dot{{\epsilon }}}_{\text{t}\text{h}}^{\text{i}\text{j}}\) are the accumulation of volume expansion/shrinkage of each phase percentage as formulated in Eq. (11).
$${\dot{{\epsilon }}}_{\text{t}\text{h}}^{\text{i}\text{j}}=\sum _{\text{i}=1}^{\text{n}}{{\xi }}_{\text{i}}{{\alpha }}_{\text{i}}\text{d}\text{T}$$
11
where \({{\xi }}_{\text{i}}\) and \({{\alpha }}_{\text{i}}\) are volume fraction and thermal expansion coefficient of i-th phase respectively with n is the number of transformed phases in the certain volume. During heating simulation, there are two numbers of transformed phase; retained initial material and austenite phase, while in quenching simulation, retained austenite phase and the non-austenite phase are the transformed phases. The thermal expansion coefficients for these two types of phase are adopted from reference [30].
Qform software implements model from Prandtl-Reuss [34, 35] in Eq. (12) to solve the elastic and plastic stress-strains problem.
$${\dot{{\epsilon }}}_{\text{e}}^{\text{i}\text{j}}+{\dot{{\epsilon }}}_{\text{p}}^{\text{i}\text{j}}=\frac{\text{d}\left\{{\sigma }\right\}}{2\text{G}}+\frac{3{\stackrel{-}{\text{d}{\epsilon }}}^{\text{p}}}{2\stackrel{-}{{\sigma }}}.\left\{{\sigma }\right\}$$
12
where \(\left\{{\sigma }\right\}\), \(\stackrel{-}{{\sigma }}\), \({\stackrel{-}{\text{d}{\epsilon }}}^{\text{p}}\), and G are stress deviator increment, effective strain or strain intensity, plastic strain increment intensity, and modulus respectively.
2.2.4. Phase transformation-stress/strain model
There are two types of strain that are influenced by phase transformations; \({\dot{{\epsilon }}}_{\text{t}\text{r}}^{\text{i}\text{j}}\) and \({\dot{{\epsilon }}}_{\text{t}\text{p}}^{\text{i}\text{j}}\). The volume strains due to phase transformation are dependent on the temperature and the volume fraction of each initial and child phase as formulated in Eq. (13–14).
$${\delta }{\text{V}}^{\left(\text{i}\right)}=\frac{{\text{V}}_{\text{C}\text{P}}^{\left(\text{i}\right)}-{\text{V}}_{\text{I}\text{P}}^{\left(\text{i}\right)}}{{\text{V}}_{\text{I}\text{P}}^{\left(\text{i}\right)}}$$
13
$${\dot{{\epsilon }}}_{\text{t}\text{r}}^{\text{i}\text{j}}=\sum _{\text{i}=1}^{\text{n}{\Pi }}\left\{{\delta }{\text{V}}^{\left(\text{i}\right)}.{\int }_{0}^{\varDelta \text{t}}{{\xi }}_{\text{C}\text{P}}^{\left(\text{i}\right)}\text{d}\text{t}\right\}$$
14
where CP and IP refer to child phase and initial phase. Meanwhile, the strains due to phase transformation induced plasticity during austenitization and martensite transformation are calculated as Eq. (15), where the transformation induced plasticity (TRIP) constants K are 10− 4 and 6.5x10− 5 [31] during austenitization and martensite transformation, respectively.
$${\dot{{\epsilon }}}_{\text{t}\text{r}}^{\text{i}\text{j}}=\frac{3}{2}\text{K}.{{\sigma }}^{\text{i}\text{j}}$$
15