Abstract

The modelling of an axisymmetric industrial quenched chromium steel bar AISI-SAE 8650H based on finite element method has been produced to investigate the impact of process history on metallurgical and material properties. Mathematical modelling of 1-dimensional line (radius) element axisymmetric model has been adopted to predict temperature history and consequently the hardness of the quenched steel bar at any point (node). The lowest hardness point (LHP) is determined. In this paper hardness in specimen points was calculated by the conversion of calculated characteristic cooling time for phase transformation t8/5 to hardness. The model can be employed as a guideline to design cooling approach to achieve desired microstructure and mechanical properties such as hardness. The developed mathematical model is converted to a computer program. This program can be used independently or incorporated into a temperature history calculator to continuously calculate and display temperature history of the industrial quenched steel bar and thereby calculate LHP. The developed program from the mathematical model has been verified and validated by comparing its hardness results with commercial finite element software results.

1. Introduction

Quenching is a heat treatment usually employed in industrial processes in order to control mechanical properties of steels such as hardness [1]. The process consists of raising the steel temperature above a certain critical value, holding it at that temperature for a specified time and then rapidly cooling it in a suitable medium to room temperature [2]. The resulting microstructures formed from quenching (ferrite, cementite, pearlite, upper bainite, lower bainite, and martensite) depend on cooling rate and on chemical composition of the steel [3].

Quenching of steels is a multiphysics process involving a complicated pattern of couplings among heat transfer, because of the complexity, coupled (thermal-mechanical-metallurgical) theory, and nonlinear nature of the problem, no analytical solution exists; however, numerical solution is possible by finite difference method, finite volume method, and the most popular one-finite element method (FEM) [4]. During the quenching process of the steel bar, the heat transfer is in an unsteady state as there is a variation of temperature with time [5]. The heat transfer analysis in this paper will be carried out in 3-dimensions. The three dimensional analysis will be reduced into a 1-dimensional axisymmetric analysis to save cost and computer time [4, 6, 16]. This is achievable because in axisymmetric conditions, there is no temperature variation in the theta (𝜽) direction and in (z) direction, the temperature deviations are only in (r). The Galerkin weighted residual technique is used to derive the mathematical model. In this paper, 1-dimensional line (radius) element will be developed.

2. Mathematical Model

The temperature history of the quenched cylindrical steel bar at any point would be calculated; three-dimensional heat transfer can be analyzed using one dimensional axisymmetric elements as shown in Figures 1, 2, and 3.

The linear temperature distribution for an element (radius) line, 𝑇 is given by:𝑇(𝑅)=𝑎1+𝑎2𝑅,(1) where 𝑇(𝑅) = nodal temperature as the function of 𝑅, 𝑎1 and 𝑎2 are constants, 𝑅 is any point on the (radius) line element.

2.1. Shape Function of the Axisymmetric Triangular Element

The shape functions were to represent the variation of the field variable over the element. The shape function of axisymmetric 1-dimensional line (radius) element expressed in terms of the 𝑟 coordinate and its coordinate as shown in Figure 4; which are derived to obtain the following shape functions as shown: 𝑆𝑖=𝑅𝑗𝑅𝑅𝑗𝑅𝑖=𝑅𝑗𝑅𝐿,𝑆𝑗=𝑅𝑅𝑖𝑅𝑗𝑅𝑖=𝑅𝑅𝑖𝐿.(2) Thus the temperature distribution of 1D radius for an element in terms of the shape function can be written as: 𝑇(𝑅)=𝑆𝑖𝑇𝑖+𝑆𝑗𝑇𝑗=𝑆(𝑟){𝑇},(3) where [S(r)] = [𝑆𝑖𝑆𝑗] is a row vector matrix and{𝑇}=𝑇𝑖𝑇𝑗(4) is a column vector of nodal temperature of the element.

Equation (3) can also be expressed in matrix form as: 𝑇(𝑅)=𝑆𝑖𝑆𝑗𝑇𝑖𝑇𝑗.(5) Thus for 1-dimensional element we can write in general: Ψ(𝑒)=𝑆𝑖𝑆𝑗𝜓𝑖𝜓𝑗,(6) where Ψ𝑖 and Ψ𝑗 represent the nodal values of the unknown variable such as in our case temperature also the unknown can be deflection, velocity, and so forth.

2.2. Natural Area Coordinate

Using the natural length coordinates and their relationship with the shape function by simplification of the integral of Galerkin solution.

The two length natural coordinates 𝐿1 and 𝐿2 at any point P inside the element are shown in Figure 5 from which we can write:𝐿1=𝑅𝑗𝑅𝑅𝑗𝑅𝑖=𝑙1𝐿,𝐿2=𝑅𝑅𝑖𝑅𝑗𝑅𝑖=𝑙2𝐿.(7) Since it is a one-dimensional element, there should be only one independent coordinate to define any point P. This is true even with natural coordinates as the two natural coordinates 𝐿1 and 𝐿2 are not independent, but are related as:𝐿1+𝐿2=1or𝐿1+𝐿2=𝑙1𝐿+𝑙2𝐿=1.(8) The natural coordinates 𝐿1 and 𝐿2 are also the shape functions for the line element, thus:𝑆𝑖=𝑅𝑗𝑅𝑅𝑗𝑅𝑖=𝑅𝑗𝑅𝐿=𝐿1,𝑆𝑗=𝑅𝑅𝑖𝑅𝑗𝑅𝑖=𝑅𝑅𝑖𝐿=𝐿2,𝑆𝑖=𝐿1,𝑆𝑗=𝐿2,𝑅=𝑅𝑗𝐿2+𝑅𝑖𝐿1=𝑅𝑖𝑆𝑖+𝑅𝑗𝑆𝑗,𝜕[𝑆]𝑖𝜕𝑟=𝜕𝐿1𝜕𝑟=1𝑅𝑗𝑅𝑖=1𝐿,𝜕[𝑆]𝑗𝜕𝑟=𝜕𝐿2𝜕𝑟=1𝑅𝑗𝑅𝑖=1𝐿.(9)

2.3. Develop Equation for All Elements of the Domain

Derivation of equation of heat transfer in axisymmetric one-dimensional line (radius) elements by pplying the conservation of energy to a differential volume cylindrical segment as shown in Figure 6,𝐸in𝐸out+𝐸generated=𝐸stored.(10) The transient heat transfer within the component during quenching can mathematically be described by simplifying the differential volume term; the heat conduction equation is derived and given by1𝑟𝑑𝑑𝑟𝐾𝑟𝑟𝑑𝑇𝑑𝑟+1𝑟2𝑑𝑑𝜃𝐾𝜃𝑑𝑇𝑑𝜃+𝑑𝑑𝑧𝐾𝑧𝑑𝑇𝑑𝑧+𝑞=𝜌𝑐𝑑𝑇𝑑𝑡,(11)𝑘𝑟 = heat conductivity coefficient in r-direction, W/m°C, 𝑘𝜃 = heat conductivity coefficient in θ-direction, W/m°C, 𝑘𝑧 = heat conductivity coefficient in z-direction, W/m°C, 𝑇 = temperature, °C, 𝑞 = heat generation, W/m3, 𝜌 = mass density, kg/m3, c = specific heat of the medium, J/kgK, t = time, s.

The assumptions made in this problem were:

(i) for axisymmetric situations one dimensional line (radius) element, there is no variation of temperature in the Z-direction as shown in Figures 1, 2, and 3. Because we already assumed that in steel quenching and cooling process of the steel bar is insulated from convection at the cross section of the front and back.

It means that we have convection and radiation at one node only which is on the surface (node 5), in our research we focus to calculate LHP which is at (node 1), where it is the last point will be cooled, this give the maximum advantage to make our assumption more safe, because it is the last point which will effect by convection and radiation from the front and back cross-section of the steel bar therefore we can write, (𝜕𝑇/𝜕𝑧=0),

(ii) for axisymmetric situations, there is no variation of temperature in the 𝜽-direction, because it is clear from Figures 1, 2, and 3 that the temperature distribution along the radius will be the same if the radius move with angle θ, 360° therefore, (𝜕𝑇/𝜕𝜃=0),

(iii) the thermal energy generation rate ̇𝑞 represent the rate of the conversion of energy from electrical, chemical, nuclear, or electromagnetic forms to thermal energy within the volume of the system. Such as conversion is the electric field which will be studied with details in the 2nd part of our research, however in this manuscript no heat generation therefore,̇𝑞=0.(12) After simplifying, (11) becomes𝑘𝑟𝜕𝜕𝑟𝑟𝜕𝑇𝜕𝑟=𝜌𝑐𝜕𝑇𝜕𝑡,(13) and also known as residual or partial differential equation{}=𝑘𝑟𝜕𝜕𝑟𝑟𝜕𝑇𝜕𝑟𝜌𝑐𝜕𝑇𝜕𝑡=0.(14)

2.4. Galerkin Weighted Residual Method Formulation

From the derived heat conduction equation, the Galerkin residual for 1-dimensional line (radius) element in an unsteady state heat transfer by integration the shape functions times the residual which minimize the residual to zero becomes𝑉[𝑆]𝑇{}(𝑒)𝑑𝑣=0,(15) where [𝑆]𝑇 = the transpose of the shape function matrix, {}(𝑒) = the residual contributed by element (𝑒) to the final system of equations𝑘𝑟[𝑆]𝑇𝜕𝜕𝑟𝑟𝜕𝑇𝜕𝑟𝑑𝑣[𝑆]𝑇𝜌𝑐𝜕𝑇𝜕𝑡𝑑𝑣=0.(16)

2.5. Chain Rule

The terms 1 and 2 of (16) can be rearranged using the chain rule which states that;(𝑓𝑔)=𝑓𝑔+𝑔𝑓.(17) Therefore, 𝑓𝑔=(𝑓𝑔)𝑓𝑔,𝜕𝜕𝑟[𝑆]𝑇𝑟𝜕𝑇𝜕𝑟=[𝑆]𝑇𝜕𝜕𝑟𝑟𝜕𝑇𝜕𝑟+𝑟𝜕𝑇𝜕𝑟𝜕[𝑆]𝑇𝜕𝑟.(18) Term 1 of (16) is rearranged, thus[𝑆]𝑇𝜕𝜕𝑟𝑟𝜕𝑇𝜕𝑟=𝜕𝜕𝑟[𝑆]𝑇𝑟𝜕𝑇𝜕𝑟𝑟𝜕𝑇𝜕𝑟𝜕[𝑆]𝑇𝜕𝑟.(19) By substituting (19) in to (16), get=𝐴𝑘𝑟𝜕𝜕𝑟[𝑆]𝑇𝑟𝜕𝑇𝜕𝑟𝑑𝑣𝐵𝑘𝑟𝜕[𝑆]𝑇𝜕𝑟𝑟𝜕𝑇𝜕𝑟𝑑𝑣[𝑆]𝑇𝜌𝑐𝜕𝑇𝜕𝑡𝑑𝑣𝐶.(20) Term 𝐴 is the heat convection terms and contributes to the conductance and thermal load matrix. Term 𝐵 is the heat conduction terms and contributes to the conductance matrix. Term 𝐶 is the transient equation and contributes to the capacitance matrix, where 𝐴𝑘𝑟𝜕𝜕𝑟[𝑆]𝑇𝑟𝜕𝑇𝜕𝑟𝑑𝑉=2𝜋𝑟𝑧[𝑆]𝑇[𝑆]{𝑇}𝑒1+2𝜋𝑟𝑧[𝑆]𝑇𝑇𝑓22𝜋𝑟𝑧𝜀𝑠𝜎[𝑆]𝑇[𝑆]{𝑇}𝑒([𝑆]{𝑇}𝑒)33+2𝜋𝑟𝑧𝜀𝑠𝜎[𝑆]𝑇𝑇4𝑓4.(21) Note that term (1) and term (3) contributed to the conductance matrix since they contains the unknown temperature {𝑇}. Terms (2) and (4) contributed to the thermal load matrix as 𝑇𝑓 is the known fluid temperature. Term (3) and term (4) heat radiation are very important if our heat treatment is annealing (cooling in the furnace) or normalizing (cooling in air or jet air), but can be ignore (neglect) if the cooling is quenching in water as in our paper.

From earlier explanations, derivation and after simplifying we can formulate the conductance matrix in the 𝑟-direction for Term 𝐵 finally we get:

Term B (the conduction term) contributes to the conductance matrix:𝑘𝐿(𝑅𝑗+𝑅𝑖)1111𝑇𝑖𝑇𝑗𝑘𝑐.(22) Similarly, Term 𝐶 is the unsteady state (transient) which contributes to the Capacitance Matrix,

Term C (heat stored) contributes to the Capacitance Matrix: 𝐿𝜌𝑐63𝑅𝑗+𝑅𝑖𝑅𝑖+𝑅𝑗𝑅𝑖+𝑅𝑗𝑅𝑖+3𝑅𝑗𝑇𝑖𝑇𝑗𝐶.(23)

Term  A (heat convection): (i)Term 𝐴1: contributes to conductance matrix,

Term 𝐴1 (the convection term) contributes to the conductance matrix:2𝑅𝑗0001𝑇𝑖𝑇𝑗𝑘,(24)(ii)Term 𝐴2: contributes to thermal load Matrix:

Term 𝐴2 (the convection term) contributes to thermal load matrix 2𝑅𝑗𝑇𝑤01𝑓.(25)

2.6. Construct the Element Matrices to the Global Matrix

Assemble the global, conductance, capacitance, and thermal load matrices and the global of the unknown temperature matrix for all the elements in the domain, that is, the element’s conductance, capacitance and thermal load matrices have been derived. Assembling these elements is necessary in all finite element analysis.

Constructing these elements will result into the following finite element equation:[𝐊](𝐆){𝐓}(𝐆)+[𝐂](𝐆)̇𝐓(𝐆)={𝐅}(𝐆),(26) with [𝐊]=[𝐊𝐜]+[𝐊𝐡]. Conductance matrix due to conduction (elements 1 to 4) and heat loss through convection at the element’s boundary (element 4 node 5) as shown in Figures 1 and 2:{𝑇}: Temperature value at each node, °C,[𝐶]: Capacitance matrix, due to transient equation (heat stored),{̇𝑇}: Temperature rate for each node, °C/s,{𝐹}={𝐹}+{𝐹̇𝑞}: heat load due to heat loss through convection at the element’s boundary (element 4 node 5) and internal heat generation (element 4 node 5).

2.7. Euler’s Method

Two point recurrence formulas will allow us to compute the nodal temperatures as a function of time. In this paper, Euler’s method which known as the backward difference scheme (BDS) will be used to determine the rate of change in temperature, the temperature history at any point (node) of the steel bar [3].

If the derivative of T with respect to time t is written in the backward direction and if the time step is not equal to zero (Δ𝑡0), we have that,̇𝑇𝑇(𝑡)𝑇(𝑡Δ𝑡)Δ𝑡,(27) with ̇𝑇 = temperature rate (°C/s); 𝑇(𝑡) = temperature at 𝑡 s(°C); 𝑇(𝑡Δ𝑡) = temperature at (𝑡Δ𝑡) s, (°C) Δ𝑡 = selected time step (s) and t = time (s) (at starting time, 𝑡 = 0).

By substituting the value of {̇𝑇} into the finite element global equation, we have that[𝐾](𝐺){𝑇(𝑡)}(𝐺)+[𝐶](𝐺)𝑇(𝑡)𝑇(𝑡Δ𝑡)Δ𝑡(𝐺)={𝐹(𝑡)}(𝐺).(28) Finally, the matrices become[𝐾](𝐺)Δ𝑡+[𝐶](𝐺){𝑇}(𝐺)𝑖+1=[𝐶](𝐺){𝑇}(𝐺)𝑖+{𝐹}(𝐺)𝑖+1Δ𝑡.(29) From (29) all the right hand side is completely known at time 𝑡, including 𝑡=0 for which the initial condition apply.

Therefore, the nodal temperature can be obtained for a subsequent time given the temperature for the preceding time.

Once the temperature history is known the important mechanical properties of the steel bar can be obtained such as hardness and strength.

3. Application

3.1. Calculation the Temperature History

The present developed mathematical model is programmed using MATLAB to simulate the results of the temperature distribution with respect to time in transient state heat transfer of the industrial quenched steel bar. The cylindrical chromium steel bar has been heated to 850°C. Then being quenched in water with 𝑇water = 32°C, the water convection heat transfer coefficient, water = 5000 W/m2°C.

The temperature history at any point (node) of the cylindrical steel bar after quenched is being identified in Figures 6 and 7. The cylindrical bar was made from chromium steel 8650H, with properties as mentioned below.

Thermal capacity, 𝜌c (J/m3·°C):0𝑇650°C, 𝜌𝑐=(0.004𝑇+3.3)×106,650<𝑇725°C, 𝜌𝑐=(0.068𝑇38.3)×106,725<𝑇800°C, 𝜌𝑐=(0.086𝑇+73.55)×106,𝑇>800°C, 𝜌𝑐=7.55×106.

Thermal conductivity, k (W/m°C)0𝑇900°C, 𝑘=0.022𝑇+48,𝑇>900°C, 𝑘=28.2,

where in our case the global conductance matrix [𝐾](𝐺), the global capacitance matrix [𝐶](𝐺) and the global thermal load matrix {𝐹}(𝐺) can be computed easily as the following:[𝐾](𝐺)=𝐾𝑐(1)+𝐾𝑐(2)+𝐾𝑐(3)+𝐾𝑐(4)+𝐾(4),(30a) where [𝐾c](1), [𝐾c](2), [𝐾c](3), [𝐾c](4) are the conductance matrices due to conduction in 1-dimensional element for the 1st element, the 2nd, the 3rd, and the 4th element respectively while [𝐾](4) because we note that there is convection in element 4 at node 𝑗(5) only as shown clearly in Figures 3 and 7:[𝐶](𝐺)=[𝐶](1)+[𝐶](2)+[𝐶](3)+[𝐶](4),(30b) where [𝐶](1), [𝐶](2), [𝐶](3), [𝐶](4) are the capacitance matrices due to transient [unsteady state] in 1-dimensional line (radius) element:{𝐹}(𝐺)=𝐹(4).(30c)

We have convection in element 4 at node 𝑗(5) only as shown clearly in Figures 3 and 7.

With the input data and boundary conditions provided, a sensitivity analysis is carried out with the developed program to obtain the temperature distribution at any point (node) of the quenched steel bar, as example, is the transient state temperature distribution results of the selected five nodes from the center [𝑊1] to the surface [𝑊5] of the quenched steed bar which were computed as shown in Figures 7 and 8.

3.2. LHP Calculation
3.2.1. Calculating the Cooling Time Required

In this study, we choose to calculate the cooling time between 800°C and 500°C [3, 79], where the characteristic cooling time, relevant for phase transformation in most structural steels is the time of cooling from 800 to 500°C (time 𝑡8/5) [1017]:𝑡𝑐=𝑡800𝑡500.(31) From Figure 7 we can determine the time taken for node 𝑊1 to reach 800°C,𝑡800C=4.830sec.(32) By the same way the time taken for node 𝑊1 to reach 500°C is 𝑡500C=13.075sec.

So the Cooling time 𝑡𝑐 for node 𝑊1; 𝑡𝑐=𝑡500C𝑡800C=13.0754.830=8.245sec.(33)

For nodes 𝑊2 to 𝑊5, the cooling time 𝑡𝑐 calculated by the same way, the final results shown in Table 1.

Table 1 shows the cooling time 𝑡𝑐 and the rate of cooling ROC.

3.2.2. Calculating the Jominy Distance from Standard Jominy Distance versus Cooling Time Curve

Cooling time, 𝑡𝑐, obtained will now be substituted into the Jominy distance versus cooling time curve to get the correspondent Jominy distance. Jominy distance can also be calculated by using polynomial expressions via polynomial regression via Microsoft Excel.

The standard Table (cooling rate at each Jominy distance (Chandler, 1999)) can be used too [18].

Then Jominy distance of nodes 𝑊1 to 𝑊5 will be calculated, the final results shown in Table 2, where the rate of cooling, ROC, was defined by;ROC=800C500C𝑡𝑐=800C500C𝑡500C𝑡800CC/s.(34)

3.2.3. Predict the Hardness of the Quenched Steel Bar

The HRC of AISI-SAE 8650H can be calculated by using the relation between the J-Distance and the HRC from the Practical date Handbook, the Timken Company 1835 Duebex Avenue SW Canton, Ohio 44706-2798 1-800-223, the final results are shown in Figure 9 and Table 3.

4. Mathematical Model Verification

The same data input for the steel properties and boundary condition used in the mathematical model is applied to the ANSYS software to verify the temperature simulation results. The temperature distribution from the ANSYS analysis is depicted figuratively as shown in Figures 10(a) and 10(b).

The temperature time graph from the ANSYS analysis is depicted as shown in Figure 11.

From the graphs shown in Figure 8 by mathematical model and Figure 11 by ANSYS, it can be clearly seen that the temperature history of the quenched steel bar have the same pattern. The heat transfer across the steel bar is uniform. From Figure 11 the cooling time, Jominy-distance, and consequently the hardness of the quenched chromium steel bar at any point (node), even the lowest hardness point (LHP) is determined by ANSYS too, the final results shown in Table 4 and Figure 12.

From our results we found that in the mathematical model for the 1st node with 𝑊1 in the center, we found that HRC = 59.297. While in ANSYS for the same node 𝐴1, we found that HRC = 58.492.

And for the nodes on the surfaces 𝑊5 and 𝐴5, it was found that HRC = 61.421 and 61.295 for the mathematical model and ANSYS, respectively. From the above, it can be seen that there is a strong agreement between both results. The difference between all the results of the mathematical model and the ANSYS simulations can be accounted due to the fact that the ANSYS software is commercial purpose, and thereby has some automated input data. But the developed mathematical model is precisely for a circular steel bar axisymmetric cross section. However, there is strong agreement between both results and thereby the result is validated, where the comparison indicated reliability of the proposed model.

Also the results showed that the node on the surface will be the 1st which completely cooled after quenching because it is in the contact with the cooling medium then the other nodes on the radial axis to the centre, respectively, and the last point will be completely cooled after quenching will be at half the length at the centre. Hence LHP will be at half the length at the centre of the quenched industrial chromium steel bar. It will be more important to know LHP once the radius of the quenched steel bar is large because LHP will be low, in other words, it will be lower than the hardness on the surface, that means increasing the radius of the bar inversely proportional with LHP.

LHP calculation experimentally is an almost impossible task using manual calculation techniques also the earlier methods only used hardness calculated at the surface, which is higher than LHP, which has negative consequence that can result to the deformation and failure of the component.

5. Conclusion

A mathematical model of steel quenching has been developed to compute LHP of the quenched chromium steel bar at any point (node) in a specimen with cylindrical geometry. The model is based on the finite element Galerkin residual method. The numerical simulation of quenching consisted of numerical simulation of temperature transient field of cooling process. This mathematical model was verified and validated by comparing the hardness results with ANSYS software simulations. From the mathematical model and ANSYS results, it is clear that the nodes on the surface (𝑊5 and 𝐴5), respectively, cools faster than the nodes on the center (𝑊1 and 𝐴1) because 𝑡𝑐𝑊5<𝑡𝑐𝑊1 and 𝑡𝑐𝐴5<𝑡𝑐𝐴1, this means that the mechanical properties will be different such as hardness where the hardness on the surface nodes (𝑊5 and 𝐴5) will be higher than the hardness on the center nodes (𝑊1 and 𝐴1).

Acknowledgments

The authors would like to thank the University Tun Hussein Onn Malaysia for supporting this research under the Science Fund Grant. The corresponding author grateful to the Postgraduate Centre of UTHM for their support of this research, where they accepted him under the university Ph.D. scholarship.