Developing 1-Dimensional Transient Heat Transfer Axi-Symmetric MM to Predict the Hardness, Determination LHP and to Study the Effect of Radius on E-LHP of Industrial Quenched Steel Bar

1. Introduction

A study of the constitution and structure of all irons and steel must first start with the iron-carbon equilibrium diagram and the steel part of the phase diagram as shown in Figure 1.

Iron: is a general word used to describe metals that have pure iron as their main constituent.

Most iron wares around us are not made of chemically pure iron but are alloys, the most important of which is Carbon. Carbon is a major factor in understanding the difference between Iron, Steel and Cast iron. Adding some carbon to chemically pure iron makes steel. If the quantity of carbon is increased, cast iron will be obtained.

Steel: is the most commonly and widely used metallic material in today’s society as shown in Figure 2. It is one of the most signification engineering material and can be classified plain carbon steel and alloy steel.

Steel is an alloy of iron and carbon. The amount of carbon dictates whether the steel is hard or tough. Adding Carbon makes the iron harder. Steels can be hardened by heat treatment.

Carbon in steel may be present up to 2.03 percent. Steels with carbon content from 0.025 percent to 0.8 percent are called hypo-eutectoid steel. Steel with a carbon content of 0.8 percent is known as eutectoid steel. Steels with carbon content greater than 0.8 percent are called hyper-eutectoid steel. There are three major categories of steel which are as follows:

Figure 1.

Fe-Fe₃C Phase Diagram, [23, 24].Many of the basic features of this system influence the behaviour of even the most complex alloy steels. The iron-carbon diagram provides a valuable foundation on which to build knowledge of both plain carbon and alloy steels.

2. Medium carbon steels (carbon from 0.3 to 0.7 percent)

3. High carbon steel (carbon more than 0.7 percent).

Alloy Steels:are basically carbon steels with certain chemical elements added to improve the properties of the metal for specific applications or end products. Alloying elements include carbon, copper, sulphur, manganese, phosphorus, nickel, molybdenum, boron and chromium, and the resulting material is called an alloy steel as in the sample AISI-SAE 8650H, which will be studied in detail here. The 1-D mathematical model was developed to:

1. Predict the hardness.

2. Determine E-LHP.

3. Study the effect of radius on E-LHP.

Figure 2.

Fe-Fe3C Phase Diagram [Steel Part], [23, 25].

TTT Diagram: T (Time) T (Temperature) T (Transformation) diagram is a plot of temperature versus the logarithm of time for a steel alloy of definite composition as shown in Figure 3. It is used to determine when transformations begin and end for an isothermal (constant temperature) heat treatment of a previously austenitized alloy. TTT diagram indicates when a specific transformation starts and ends and it also shows what percentage of transformation of austenite at a particular temperature is achieved. In Figure 3 the cooling rates A and B indicate two rapid cooling processes. In this case curve A (water cooled) will cause a higher hardness than the cooling rate B (sea water cooled). The end product of both cooling rates will be martensite. Cooling rate B is also known as the Critical Cooling Rate, which is represented by a cooling curve that is tangent to the nose of the TTT diagram. Critical Cooling Rate is defined as the lowest cooling rate which produces 100% Martensite while minimizing the internal stresses and distortions.

Figure 3.

Rapid Quench, [23].

Many of the basic features of this system influence the behaviour of even the most complex alloy steels. For example, the phases found in the simple binary Fe-C system persist in complex steels, but it is necessary to examine the effects of alloying elements have on the formation and properties of these phases. The iron-carbon diagram provides a valuable foundation on which to build knowledge of both plain carbon and alloy steels in their immense variety [1]. Quenching is physically one of the most complex processes in engineering, and very difficult to understand. Simulation of steel quenching is thus a complex problem [2, 9-23]. Quenching is a heat treatment usually employed in industrial processes in order to control mechanical properties of steels such as toughness and hardness [9-23]. 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 [9-23]. 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, 9-23]. Quenching of steels is a multi-physics process involving a complicated pattern of couplings among heat transfer, because of the complexity, coupled (thermal-mechanical-metallurgical) theory and non-linear 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) which will be used in this study [9-23]. A larger diameter rod quenched in a particular medium will obviously cool more slowly than a smaller diameter rod given with a similar heat treatment [4].

It will be more important to know the lowest hardness point E-LHP (node E) once the radius of the quenched steel bar increases, in other words the E-LHP will be lower than the hardness on the surface (node 4), thus the radius of the bar will be inversely proportional to the hardness at E-LHP. 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 [9-23]. The heat transfer analysis in this study will be carried out in 3-dimensions. The three dimensional analysis will be reduced into a 1-dimensional axisymmetric analysis to save cost and computing time [9-23]. This is achievable because in 1-dimensional axisymmetric conditions, there is no temperature variation in theta (Ɵ) and (z)-direction for 1-D as shown in Figure 4 (a)(b) and Figure 5, the temperature deviations is only in (r)-direction [9-23]. The Galerkin weighted residual technique will be used to derive the mathematical model to predict the hardness at any point (node) of the heat treated quenched steel bar. Therefore [E-LHP] can be calculated where it’s exactly at half the length at the centre of the bar as shown in Figure 1. Experimentally, measurement of E-LHP is an almost impossible task using manual calculation techniques and furthermore the earlier methods only used hardness calculated at the surface (node 5) as shown in Figure 4 (a)(b) and Figure 5. This surface hardness is higher than E-LHP and this has negative consequences which can lead to the deformation and failure. In this chapter, 1-D line (radius) element will be used.

Quenching of steels in general has been and continues to be an important commercial manufacturing process for steel components. It is a commonly used heat treatment process employed to control the mechanical properties of steels. In this chapter hardness in specimen points was calculated by the conversion of calculated characteristic cooling time for phase transformation t_8/5 to hardness [5-8, 9-23]. Temperature histories must be performed to obtain more accurate transformation kinetics; an adequate tool has been produced for investigating the impact of process history on metallurgy and material properties. Mathematical modelling of axisymmetric industrial heat treated quenched steel bar on the finite element method has been developed to predict temperature history then the hardness at any point (node) even inside the bar can be determined and also the effect of the radius on the temperature history subsequently on LHP can be studied. The temperature history needs to be properly understood in order to efficiently produce high quality components.

It is clear that the first point (node) will be completely cooled after quenching (surface node) because it is located on the surface in contact with the cooling medium, then the other points (nodes) on the radial axis to the centre accordingly will be cooled and the last point will be completely cooled after quenching (centre node)[15, 34].

Thus the maximum hardness will be on the surface node subjected to fast cooling, then the hardness will decrease from the surface node on the radial axis to the centre node of the quenched steel bar. Consequently the lowest hardness point of the quenched steel bar will be detected at the centre node [15, 34].

The lowest hardness point (E-LHP) should be expected inside the heat treated quenched steel bar at the half of the length of the bar centre (centre node). To prove this experimentally is an almost impossible task using manual calculation techniques.

It will be more important to know the E-LHP (centre node) when the radius of the quenched steel bar increases because the lowest hardness point will be lower than the hardness on the surface (surface node). This means that increasing the radius of the bar is inversely proportional to LHP (centre node), while the hardness at the surface (surface node) will be the same [15, 34].

No published information are available till date on this aspect. This work represents a contribution towards understanding of steel behaviour at elevated temperature during quenching at the LHP (centre node) of the steel bar. The results of this study may prove useful to obtain the hardness of the lowest point of the steel bar in order to reach the maximum benefit against the deformation and failure of the component [15, 34].

2. Mathematical model

Three dimensional heat transfers can be analyzed using one dimensional axisymmetric element as shown in Figure 4 (a)(b) and Figure 5 [9-22].

Figure 4.

a) Theaxisymmetric one dimensional line (radius) elementfromthe domain, on the cylindrical steel bar which had been heated and then submerged in water. (b) The1-dimensional elementfromthe domain on the axisymmetric rectangular cross section when the radii equal 12.5mm, theselected 4 elements with 5 nodes and the boundary at node j [5] for an element 4.

Figure 5.

Theaxisymmetric one dimensional line (radius) elementfromthe domain, the selected 4 elements with 5 nodes and the boundary at node j [5] for an element 4.

2.1. Methodology of building the F. E. Model in details

The temperature distribution inside the cylindrical steel bar at thermal equilibrium will be calculated. These are special classes of three-dimensional heat transfer problem:

1. Geometrically axisymmetric.

2. Each thermal load is symmetrical about an axis.

This three-dimensional heat transfer problem may be analyzed using one-dimensional axisymmetric elements as shown in Figure 4 (a)(b) and Figure 5 [9-22].

The finite element method is applied to the one-dimensional cylindrical coordinates heat transfer problem.

The finite element formulation is developed with the Garlekin Weighted-residual method. The appropriate working expressions of the conductance matrix, capacitance matrix and thermal load matrix are derived in details.

The time dependent solution is obtained by applying the Backward Difference Scheme.

2.1.1. Meshing the engineering problem of the domain

Since the modelling work is on one-dimensional axisymmetric elements then line element has been selected in this study. Let us consider a cylindrical chromium steel bar as shown in Figure 4 (a) which had been heated and then submerged in water.

The linear temperature distribution for an element (radius) line, T is given by:

$${\text{T}}^{(\text{R})}={\text{ a}}_{\text{1}}+{\text{ a}}_{\text{2}}\text{R}$$E1

Where,

T^(R) = nodal temperature as the function of R

a₁ and a₂ are constants.

R is any point on the (radius) line element.

2.1.2. Shape function of 1-D axisymmetric element

The shape functions were to represent the variation of the variable field over the element. The shape function of axisymmetric 1-Dimensional line (radius) is element expressed in terms of the r coordinate and its coordinate are shown in Figure 6;

Figure 6.

one-dimensional linear temperature distributions foran element (radius) line in global coordinate system.

Theyare derived to obtain the following shape functions [9-22];

$$\begin{array}{l}{\text{S}}_{\text{i}}=\left(\frac{R_j-R}{R_j-R_i}\right)=\left(\frac{R_j-R}{L}\right)\text{(a)}\\ {\text{S}}_{\text{j}}=\left(\frac{R-R_i}{R_j-R_i}\right)=\left(\frac{R-R_i}{L}\right)\text{(b)}\end{array}$$E2

Thus the temperature distribution of 1-D radius for an element in terms of the shape function can be written as:

$${\text{T}}^{(\text{R})}={\text{ S}}_{\text{i}}{\text{T}}_{\text{i}}+{\text{ S}}_{\text{j}}{\text{T}}_{\text{j}}={\text{ S}}^{(\text{r})}\left\{\text{T}\right\}$$E3

Where [S^(r)] = [S_i S_j] is a row vector matrix and $\left\{T\right\}=\left\{\begin{array}{c}{T}_i\\ T_j\end{array}\right\}$is a column vector of nodal temperature of the element.

Eq. (3) can also be expressed in matrix form as:

$${\text{T}}^{(\text{R})}=\left[{\text{S}}_{\text{i}} {\text{S}}_{\text{j}}\right]$$E4

Thus for 1-dimensional element we can write in general:

$${\Psi }^{(\text{e})}=\left[{\text{S}}_{\text{i}} {\text{S}}_{\text{j}}\right]\left\{\begin{array}{c}{\psi }_i\\ {\psi }_j\end{array}\right\}$$E5

Where Ψ_i and Ψ_j represent the nodal values of the unknown variable which in our case is temperature. The unknown can also be deflection, or velocity etc.

2.1.3. 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 L₁ and L₂ at any point p inside the element are shown in Figure 7 [9, 10, 19]. From which we can write:

$$\begin{array}{l}{L}_1=\frac{R_j-R}{R_j-R_i}=\frac{l_1}{L}\text{(a)}\\ L_2=\frac{R-R_i}{R_j-R_i}=\frac{l_2}{L}\text{(b)}\end{array}$$E6

Figure 7.

Two-node line element showing interior point pand the two naturals coordinates L₁ and L₂.

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 L₁ and L₂ are not independent, but are related as:

$$L_1+L_2=1 \text{or} L_1+L_2=\frac{l_1}{L}+\frac{l_2}{L}=1$$E7

The natural coordinates L₁ and L₂ are also the shape functions for the line element, thus:

$$\begin{array}{l}{S}_i= \left(\frac{R_j-R}{R_j-R_i}\right) = \left(\frac{R_j-R}{L}\right) =L_1\\ S_j= \left(\frac{R-R_i}{R_j-R_i}\right)= \left(\frac{R-R_i}{L}\right) =L_2\\ S_i=L_1,S_j=L_2,\end{array}$$E8

$$R=R_j{L}_2+R_i{L}_1=R_i{S}_i+R_j{S}_j$$E9

$$\frac{\partial {\left[S\right]}_i}{\partial r}=\frac{\partial L_1}{\partial r}=\frac{-1}{R_j-R_i}=-\frac{1}{L}$$E10

$$\frac{\partial {\left[S\right]}_j}{\partial r}=\frac{\partial L_2}{\partial r}=\frac{1}{R_j-R_i}=\frac{1}{L}$$E11

2.2. Develop equation for all elements of the domain

Derivation of equation of heat transfer in axisymmetric 1-dimensionalline (radius) elements [9-22].By applying the conservation of energy to a differential volume cylindrical segment has been done. As shown in Figure 8;

Figure 8.

Axisymmetric element from an axisymmetric body.

$$E_{in}–E_{out}+E_{generated}=E_{stored}$$E12

(12)

The transient heat transfer within the component during quenching can mathematically be described by simplifying the differential volume term [4, 12]; the heat conduction equation is derived and given by;

$$\frac{1}{r}\frac{d}{dr}\left(K_{r}r\frac{dT}{dr}\right)+\frac{1}{r^2}\frac{d}{d\theta }\left(K_{\theta }\frac{dT}{d\theta }\right)+\frac{d}{dz}\left(K_z\frac{dT}{dz}\right)+q=\rho c\frac{dT}{dt}$$E13

Where;

k_r = heat conductivity coefficient in r-direction, W/m∙°C.

k_θ = heat conductivity coefficient in θ-direction, W/m∙°C.

kz = heat conductivity coefficient in z-direction, W/m∙°C.

T = temperature, °C.

q = heat generation, W/m³.

ρ = mass density, kg/m³.

c = specific heat of the medium, J/kg∙K.

t = time, s.

2.3. The assumption made in this problem was

For axisymmetric situations one dimensional line (radius) element, there is no variation of temperature in the Z-direction as shown in Figure 4 (a)(b)and Figure 5. This is because we have 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 E-LHP which is at [node 1], where it is the last point that will be cooled. This gives the maximum advantage to make our assumption safer, because it is the last point that will be affected by convection and radiation from the front and back cross section of the steel bar.

Therefore we can write,$$\left(\frac{\partial T}{\partial z}=0\right)$$

For axisymmetric situations, there is no variation of temperature in the θ-direction, because it is clear from Figure 4 (a)(b) and Figure 5 that the temperature distribution along the radius will be the same if the radius moves with an angle θ, 360^o.

Therefore,$$\left(\frac{\partial T}{\partial \theta }=0\right)$$.

The thermal energy generation rate$$\dot{q}$$represents the rate of the conversion of energy from electrical, chemical, nuclear, or electromagnetic forms to thermal energy within the volume of the system. The conversion of the electric field will be studied with details in the 2^nd part of our research, however in this manuscript no heat generation has been taken into accounted.

Therefore, $$\dot{q}$$= 0

After simplifying, the Eqn. (13) becomes;-

$$\frac{k}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)={\rho }_{}c\frac{\partial T}{\partial t}$$E14

And also known as residual or partial differential equation

$$\left\{\Re \right\}=\frac{k}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)-{\rho }_{}c\frac{\partial T}{\partial t}=0$$E15

2.4. Galerkinweighted residual method formulation

From the derived heat conduction equation, the Galerkin residual for 1-dimensional line (radius) element in an unsteady state heat transfer can be obtained by integration the transpose of the shape functions times the residual which minimize the residual to zero becomes;

$${\displaystyle {\int }_V{[S]}^T{\left\{\Re \right\}}^{(e)}dv=0}$$E16

Where, [S]^T = the transpose of the shape function matrix

$${\left\{\Re \right\}}^{\left(\text{e}\right)}$$= the residual contributed by element (e) to the final system of equations.

$$\frac{k}{r}{\displaystyle \int \stackrel{1}{\overbrace{{\left[S\right]}^T\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)dv}}-{\displaystyle \int \stackrel{2}{\overbrace{{\left[S\right]}^T\rho c\frac{\partial T}{\partial t}dv}}}}=0$$E17

2.5. Chain rule

The term 1 and 2 of Eqn. (17) can be re-arranged using the chain rule which states that;

$${\left(fg\right)}^{-}=f{g}^{-}+g{f}^{-}$$E18

Therefore,$$f{g}^{-}={\left(fg\right)}^{-}-f^{-}g$$ then$$\frac{\partial }{\partial r}\left({\left[S\right]}^{T}r\frac{\partial T}{\partial r}\right)$$=$${[S]}^T\left\{\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)\right\}$$+$$r\frac{\partial T}{\partial r}\frac{\partial {\left[S\right]}^T}{\partial r}$$

Term 1 of Eqn. 17 is rearranged thus;

$${[S]}^T\left\{\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)\right\}=\frac{\partial }{\partial r}\left({\left[S\right]}^{T}r\frac{\partial T}{\partial r}\right)-r\frac{\partial T}{\partial r}\frac{\partial {\left[S\right]}^T}{\partial r}$$E19

By substitute Eqn. (18) in to Eqn. (17), get

$$=\stackrel{A}{\overbrace{\frac{k}{r}{\displaystyle \int \left\{\frac{\partial }{\partial r}\left({\left[S\right]}^{T}r\frac{\partial T}{\partial r}\right)\right\}dv}}}-\stackrel{B}{\overbrace{\frac{k}{r}{\displaystyle \int \left\{\frac{\partial {\left[S\right]}^T}{\partial r}r\frac{\partial T}{\partial r}\right\}dv}}}-\underset{C}{\underbrace{{\displaystyle \int {\left[S\right]}^T\left\{{\rho }_{}{c}_{}\frac{\partial T}{\partial t}\right\}dv}}}$$E20

Term A is the heat convection terms and contributes to the conductance and thermal load matrix. Term B is the heat conduction terms and contributes to the conductance matrix. Term C is the transient equation and contributes to the capacitance matrix.

Where,

$$\begin{array}{l}\stackrel{A}{\overbrace{\frac{k}{r}{\displaystyle \int \frac{\partial }{\partial r}\left({\left[S\right]}^{T}r\frac{\partial T}{\partial r}\right)}dV}}=\\ \underset{(1)}{\underbrace{-2\pi hrz{\left[S\right]}^T\left[S\right]{\left\{T\right\}}^e}}+\underset{(2)}{\underbrace{2\pi hrz{\left[S\right]}^T{T}_f}}-\underset{(3)}{\underbrace{2\pi rz{\epsilon }_s\sigma {\left[S\right]}^T\left[S\right]{\left\{T\right\}}^e{\left(\left[S\right]{\left\{T\right\}}^e\right)}^3}}+\underset{(4)}{\underbrace{2\pi rz{\epsilon }_s\sigma {\left[S\right]}^T{T}_f^4}}\end{array}$$E21

Note: term (1) and term (3) contributed to the conductance matrix since they contain the unknown temperature {T}. Term (2) and (4) contributed to thermal load matrix as T_f is the known fluid temperature. Term (3) and term (4) heat radiation very important if our heat treatment is Annealling [cooling in the furnace] or Normalizing [cooling in air or jet air], but it can be ignored if the cooling is quenching in water, sea water or oil as in this work.

From earlier explanations, derivation and after simplification we can formulate the conductance matrix in the r-direction for term B finally we get:

Term B (the conduction term) contributes to the conductance matrix

$$\underset{{\text{K}}_{\text{c}}}{\underbrace{\stackrel{}{\overbrace{\frac{k}{L}\left(R_j+R_i\right)\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left\{\begin{array}{c}Ti\\ \\ Tj\end{array}\right\}}}}}$$E22

Similarly, term C, the unsteady state (transient) which contributes to the Capacitance Matrix becomes:

Term C (heat stored) contributes to the Capacitance Matrix

$$\underset{\text{C}}{\underbrace{\stackrel{}{\overbrace{\frac{L\rho c}{6}\left[\begin{array}{cc}\left(3R_i+R_j\right)& \left(R_i+R_j\right)\\ \left(R_i+R_j\right)& \left(R_i+3R_j\right)\end{array}\right]\left\{\begin{array}{c}{\dot{T}}_i\\ \\ {\dot{T}}_j\end{array}\right\}}}}}$$E23

Term A (Heat Convection)

Term A₁ -Contributes to Conductance Matrix

Term A₁ (the convection term) contributes totheconductance matrix

$$\stackrel{}{\overbrace{\underset{{\text{K}}_{\text{h}}}{\underbrace{-2h{R}_j\left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right]\left\{\begin{array}{c}Ti\\ \\ Tj\end{array}\right\}}}}}$$E24

Term A₂ -Contributes to thermal load Matrix

Term A₂(the convection term) contributes to thermal loadmatrix

$$\stackrel{}{\overbrace{\underset{{\text{f}}_{\text{h}}}{\underbrace{2h{R}_j{T}_w\left[\begin{array}{c}0\\ 1\end{array}\right]}}}}$$E25

2.6. Construct the element Matrices to the Global Matrix

The global, conductance, capacitance and thermal load matrices and the global of the unknown temperature matrix for all the elements in the domain are assembled i.e. the element's conductance; capacitance and thermal load matrices have been derived. Assembling these elements are necessary in all finite element analysis [9-22].

Constructing these elements will result into the following finite element equation:

$${\left[K\right]}^{\left(G\right)}{\left\{T\right\}}^{\left(G\right)}+{\left[C\right]}^{\left(G\right)}{\left\{\dot{T}\right\}}^{\left(G\right)}={\left\{F\right\}}^{\left(G\right)}$$E26

Where:

$$\left[K\right]=\left[K_c\right]+\left[K_h\right]$$: is conductance matrix due to Conduction (Elements 1 to 4) and heat lossthrough convection at the element’s boundary (element 4 node 5) as shown in Figure 1, Figure 2 and Figure 3.$$\left\{T\right\}$$: is temperature value at each node, °C.$$\left[C\right]$$: is capacitance matrix, due to transient equation (heat stored)$$\left\{\dot{T}\right\}$$: is temperature rate for each node, °C/s.$$\left\{F\right\}=\left\{F_h\right\}+\left\{F_{\dot{q}}\right\}$$: is 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 points recurrence formulas will allow us to compute the nodal temperatures as a function of time. In this paper, Euler’s method which is known as the backward difference scheme (FDS) will be used to determine the rate of change in temperature, the temperature history at any point (node) of the steel bar [27, 28, 32].

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 (∆t ≠ 0), we have;

$$\left\{{\left[K\right]}^{\left(G\right)}\right\}{\left\{T\left(t\right)\right\}}^{\left(G\right)}+{\left[C\right]}^{\left(G\right)}{\left\{\frac{T\left(t\right)-T\left(t-\Delta t\right)}{\Delta t}\right\}}^{\left(G\right)}={\left\{F\left(t\right)\right\}}^{\left(G\right)}$$E27

With;

Ṫ = temperature rate (°C/s); T (t) = temperature at t s (°C); T (t-∆t) = temperature at (t-∆t) s, (°C)

∆t = selected time step (s) and t = time (s) (at starting time,t = 0)

By substituting the value of {Ṫ} into the finite element global equation, we have that;

$${\left[K\right]}^{\left(G\right)}{\left\{T\left(t\right)\right\}}^{\left(G\right)}+{\left[C\right]}^{\left(G\right)}{\left\{\frac{T\left(t\right)-T\left(t-\Delta t\right)}{\Delta t}\right\}}^{\left(G\right)}={\left\{F\left(t\right)\right\}}^{\left(G\right)}$$E28

Finally, the matrices become;

$$\left[{\left[K\right]}^{\left(G\right)}\Delta t+{\left[C\right]}^{\left(G\right)}\right]{\left\{T\right\}}_{i+1}^{\left(G\right)}={\left[C\right]}^{\left(G\right)}{\left\{T\right\}}_i^{\left(G\right)}+{\left\{F\right\}}_{i+1}^{\left(G\right)}\Delta t$$E29

From Eqn. (27) all the right hand side is completely known at time t, including t = 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 low carbon steel bar can be obtained such as hardness and strength.

3. Application

3.1. Calculation of 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 chromium steel AISI 8650H. The cylindrical chromium steel bar has been heated to 850°C. Then being quenched in water with T_{sea water} = 32°C and the convection heat transfer coefficient, h_{sea water} = 1250 W/m2∙°C.

The temperature history for the selected nodes of the cylindrical chromium steel AISI 8650H after quenching is being identified in Figure 9 and Figure 10.

The cylindrical bar was made from chromium steel AISI 8650H, with properties as seen below [26].

Thermal capacity, ρc (J/m³∙°C)

$$0\le T\le 650°\text{C}$$, $$\rho c=\left(0.004T+3.3\right)\times {10}^6$$, $$650<T\le 725°\text{C}$$,

$$\rho c=\left(0.068T-38.3\right)\times {10}^6$$E30

$$725<T\le 800°\text{C}$$,$$\rho c=\left(-0.086T+73.55\right)\times {10}^6$$ , $$T>800°\text{C}$$,

$$\rho c=7.55\times {10}^6$$E31

Thermal conductivity, k (W/m∙°C)

$$0\le T\le 900°\text{C}$$, $$k=-0.022T+48$$, $$T>900°\text{C}$$,

$$k=28.2$$E32

Where in our case the global conductance matrix [K]^(G), the global capacitance matrix [C]^(G) and the global thermal load matrix {F}^(G) can be computed easily as follow:

$${\left[K\right]}^{\left(G\right)}={\left[K_c\right]}^{\left(1\right)}+{\left[K_c\right]}^{\left(2\right)}+{\left[K_c\right]}^{\left(3\right)}+{\left[K_c\right]}^{\left(4\right)}+{\left[K_h\right]}^{\left(4\right)}$$E33

Where [K_c]⁽¹⁾, [K_c]⁽²⁾, [K_c]⁽³⁾, [K_c]⁽⁴⁾ are the conductance matrices due to conduction in 1-D element for the 1st element, the 2nd, the 3rd, and the 4th element respectively, while [K_h]⁽⁴⁾ because we note that there is convection in element 4 at node j(5) only as shown clearly in Figure 5 and Figure 9.

$${\left[C\right]}^{\left(G\right)}={\left[C\right]}^{\left(1\right)}+{\left[C\right]}^{\left(2\right)}+{\left[C\right]}^{\left(3\right)}+{\left[C\right]}^{\left(4\right)}$$E34

Where [C]⁽¹⁾, [C]⁽²⁾, [C]⁽³⁾, [C]⁽⁴⁾ are the capacitance matrices due to transient [unsteady state] in 1-dimensional line (radius) element.

$${\left\{F\right\}}^{\left(G\right)}={\left\{F_h\right\}}^{\left(4\right)}$$E35

We have convection in element 4 at node j⁽⁵⁾ only as shown clearly in Figure 5 and Figure 9.

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 an example the transient state temperature distribution results of the selected five nodes from the center [W₁] to the surface [W₅] of the quenched chromium steel AISI 8650H which were computed as shown in Figure 10 with dimension as in Figure 9. [26].

Figure 9.

Theaxisymmetric one dimensional line (radius) elementfromthe domain when the radius equal 12.5 mm, theselected 4 elements with 5 nodes and the boundary at node j [5] for an element 4 of chromium steel AISI 8650H.

Figure 10.

Graph of temperature history along WW cross-section when the radius = 12.5 mm fromMATLAB program.

3.2. LHP Calculation

3.2.1. Calculating the cooling time required

In this study, we choose to calculate the cooling time between 800^oC and 500^oC [5-8, 27-31]. Where, the characteristic cooling time, relevant for phase transformation in most structural steels is the time of cooling from 800 to 500°C (time t_8/5) [5-8, 27-31].

$$t_c=t_{800}-t_{500}$$E36

Interpolation Method:

From Figure 10 we can determine the time taken for node W₅ to reach 800^oC,

By interpolation method as the following:

Node W₅:t = 0.6 s when T = 800.151⁰C

Nolde W₅:t = t₈₀₀⁰_Cwhen T = 800⁰C

Node W₅:t = 1 s when T = 781.6454⁰C

t1	T1
t800	800
t3	T3

Foloving for t₈₀₀

$$\frac{t_{800}-t_1}{t_3-t_1}=\frac{{{}^{‵}T}_{800}-T_1}{T_3-T_1}$$E37

Thus,

$$t_{800}=\frac{(T_{800}-T_1)}{(T_3-T_1)}(t_3-t_1)+t_1=\frac{(800-800.151)}{(781.6454-800.151)}(1-0.6)+0.6=0.603\sec$$E38

$${\text{t}}_{\text{800}°\text{C}}=0.603\sec$$E39

The time taken for node W₅ to reach 500^oC

$$=\frac{(500-502.9792)}{(483.933-502.9792)}(22-20)+20=20.313\sec$$E40

thus

$${\text{t}}_{\text{5}00°\text{C}}=20.313\sec$$E41

So the Cooling time t_c for node W₅;

$$t_c=t_{5\text{00 }°\text{C }}-t_{\text{800 }°\text{C }}=\text{ 2}0.\text{313 }–0.\text{6}0\text{3 }=\text{ 19}.\text{71 sec}$$E42

For nodes W₁ to W₄, the cooling time t_c calculated by the same way, the final results shown in Table 1.

3.2.2. Calculating the Jominy distance from Standard Jominy distance versus cooling time

Cooling time, t_c 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.

In this chapter the standard Table [Cooling rate at each Jominy distance (Chandler, H., 1998)] will be used. Then Jominy distance of nodes W₁ to W₅ will be calculated by using the data from [Cooling rate at each Jominy distance (Chandler, H., 1998)] [33]. The final results shown in Table 1, where the Rate of Cooling, ROC, was defined as;

$$\text{ROC}=\frac{\text{800}°\text{C-500}°\text{C}}{t_c}=\frac{\text{800}°\text{C-500}°\text{C}}{t_{\text{500}°\text{C}}-t_{\text{800}°\text{C}}}\left(\raisebox{1ex}{${}^{o}C$}\!\left/ \!\raisebox{-1ex}{$Sec$}\right.\right)$$E43

3.2.3. Predict the hardness of the quenched steel bar

The HRC of chromium steel AISI 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 shown in Figure 11& Table 1:

Table 1.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes W₁ toW₅, sea water cooled by 1-D mathematical model when the radius = 12.5 mm.

Figure 11.

Hardness distribution along WW cross section for thenodes W₁ to W₅ from the centre to the surface respectively at half the length at the centre of the quenched steel bar sea water cooled by 1-D mathematical model, when the radius = 12.5 mm.

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 Figure 12(a) and Figure 12(b); Figure 12(a) shows the temperature distribution just before the steel bar becomes completely cooled and Figure 12(b) shows the temperature distribution at the moment that the entire steel bar becomes completely cooled after 1175s.

The temperature time graph from the ANSYS analysis is depicted as shown in Figure 13;

Figure 12.

a)(b)

From the graphs shown in Figure 10 by mathematical model and Figure 13 by ANSYS, it can be clearly seen that the temperature history of the quenched steel bar has the same pattern. The heat transfer across the steel bar is uniform. From Figure 13 the cooling time, Jominy-distance and consequently the hardness of the quenched Chromium steel 8650H at any point (node), even the lowest hardness point (E-LHP) is determined by ANSYS too, the final results shown in Table 2 and Figure 14.

Figure 13.

Temperature-time graph from ANSYS.

Table 2.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes D₁₁ toD₅₅, sea water cooled by ANSYS.

Figure 14.

Hardness distribution by ANSYS along DD cross section for thenodes D₁₁ to D₅₅ from the centre to the surface respectively at the half length at the centre of the quenched chromium steel bar where sea water is cooled, when the radius = 12.5 mm.

From the above results it was found that in the mathematical model for the 1st node with W₁ in the center, we found that HRC = 51.688. While by ANSYS for the same node A₁, it was found that HRC = 50.359. And for the nodes on the surfaces W₅ and A₅, it was found that HRC = 54.506 and 54 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 E-LHP will be at half the length at the centre of the quenched industrial Chromium steel bar.

5. Effect of the radius on the temperature history

Where there is strong agreement between the mathematical model and the Ansys simulation results and thereby the result is validated where, the comparison indicated reliability of the proposed model thereby; we will apply the mathematical model to study the effect of radius on the temperature history and on E-LHP.By the same way the temperature history of the 5 selected nodes has been obtained when the radii 18.75, 25, 37.5 and 50 mm is determined. The final results are shown in Figure 15, Figure 16, Figure 17, and Figure 18 respectively.

Figure 15.

Graph of temperature history along WW cross-section when the radius = 25 mm fromMATLAB program.

Figure 16.

Graph of temperature history along WW cross-section when the radius = 25 mm fromMATLAB program.

Figure 17.

Graph of temperature history along WW cross-section when the radius = 37.5 mm fromMATLAB program.

It is clear from the results that the nodes on the surfaces cooled faster than the nodes at half the lengths at the centres in other words the nodes on the surfaces will be the 1st to be completely cooled after quenching because it is in the contact with the cooling medium then the other points (nodes) on the radial axis to the centre respectively while the last point that will be completely cooled after quenching will be at half the length at the centre because the cooling time t_C of nodes W₅, W_5-12, W₅₅, W_5-23 and W₅₅₅ less than t_C of nodes W₁, W_1-12, W₁₁, W_1-23 and W₁₁₁, respectively. A larger diameter rod quenched in a particular medium will obviously cool more slowly than a smaller diameter rod. where Figure 10, Figure 15, Figure 16, Figure 17, and Figure 18 showed that when the radii equals 12.5, 18.75, 25, 37.5 and 50 mm it required 154, 246, 354, 606 and 924 sec respectively to be cooled decreasingly from the austenitizing temperature [850°C] to the fluid temperature [32°C]. Based on the above results we expect that if the radius is 100 mm, 3000 sec is required to reach the fluid temperature.

Figure 18.

Graph of temperature history along WW cross-section when the radius = 50 mm fromMATLAB program.

6. Effect of the radius on E-LHP

As explained above and from Figure 15, Figure 16, Figure 17, and Figure 18 of the temperature history the hardness when the radii equals 18.75, 25, 37.5 and 50 mm is determined, the final results shown in Figures 19, Table 3, Figure 20, Table 4, Figures 21, Table 5, Figure 22 and Table 6.

Table 3.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes W_1-12 toW_5-12, sea water cooled by 1-D mathematical model, when the radius = 18.75 mm.

Table 4.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes W₁₁ toW₅₅, sea water cooled by 1-D mathematical model, when the radius = 25 mm.

Figure 19.

Hardness distribution along WW cross section for thenodes W_1-12 to W_5-12 from the centre to the surface respectively at half the length at the centre of the quenched steel bar sea water cooled by 1-D mathematical model, when the radius = 18.75 mm.

Figure 20.

Hardness distribution along WW cross section for thenodes W₁₁ to W₅₅ from the centre to the surface respectively at half the length at the centre of the quenched steel bar sea water cooled by 1-D mathematical model, when the radius = 25 mm.

Table 5.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes W_1-23 toW_5-23, sea water cooled by 1-D mathematical model, when the radius = 37.5 mm.

Figure 21.

Hardness distribution along WW cross section for thenodes W_1-23 to W_5-23 from the centre to the surface respectively at half the length at the centre of the quenched steel bar sea water cooled by 1-D mathematical model, when the radius = 37.5 mm.

Table 67.

Cooling time, Cooling rate, Jominy distance and HRC for the nodes W₁₁₁ toW₅₅₅, sea water cooled by 1-D mathematical model, when the radius = 50 mm.

Figure 22.

Hardness distribution along WW cross section for thenodes W₁₁₁ to W₅₅₅ from the centre to the surface respectively at half the length at the centre of the quenched steel bar sea water cooled by 1-D mathematical model, when the radius = 50 mm.

From the above results it was found that the hardness on the surfaces nodes will be higher than the hardness on the centres nodes as shown in Figures 11, 19, 20, 21 and 22 where the hardness on the surfaces at the nodes W₅, W_5-12,W₅₅, W_5-23 and W₅₅₅ equals 54.506, 51.579. 49.857, 48.700 and 48.236 respectively while the hardness at half the length at the centres W₁, W_1-12, W₁₁, W_1-23 and W₁₁₁ equals 51.688, 47.125, 43.614, 41.222 and 38.931 respectively.Hence E-LHP will be at half the length at the centre of the quenched industrial steel bar. Hence lowest hardness point [E-LHP] will be at centre of half length of the quenched industrial chromium steel bar.

It is clear from the results as shown in Figure 23 and Table 7 that increasing the radius of the bar inversely proportional with E-LHP, then it’s more important to know E-LHP once the radius of the quenched steel bar is large because this has negative consequence which can result of the deformation and failure of the component.

Table 8.

Effect of radius on E-LHP of transient heat transfer axi-symmetric industrial quenched chromium steel AISI 8650H, sea water cooled by 1-D mathematical model.

Figure 23.

Effect of the radius on E-LHP

From the above explanation we can say that the developed 1-Dimensional mathematical modelling in this chapter solved design problems in transient heat transfer axi-symmetric industrial quenched chromium steel 8650H bar, where experimental calculation of E-LHP is an almost impossible task using manual calculation techniques. Also the earlier methods only used hardness calculated at the surface, which is higher than E-LHP.

7. Conclusion

A mathematical model of steel quenching has been developed to compute E-LHP of the quenched Chromium steel 8650H at any point (node) in a specimen with cylindrical geometry and to study the effect of radius on E-LHP. 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 [W₅, W_5-12, W₅₅, W_5-23, W₅₅₅ and D₅₅] respectively cools faster than the nodes on the center [W₁, W_1-12, W₁₁, W_1-23, W₁₁₁ and D₁₁] because t_CW5, t_CW5-12, t_CW55, t_CW1-23, t_CW555 and t_CD55 less than t_CW1, t_CW1-12, t_CW11, t_CW1-23, t_CW111 and t_CD11, this means that the mechanical properties will be different such as hardness where the hardness on the surfaces nodes [W₅, W_5-12, W₅₅, W_5-23, W₅₅₅ and D₅₅] will be higher than the hardness on the center nodes [W₁, W_1-12, W₁₁, W_1-23, W₁₁₁ and D₁₁] respectively.