A Novel (1-D) Mathematical Modeling to Determine (E-LHP) of Industrial Transient Heat Transfer Quenched Chromium Steel 5147H, Sea Water Cooled

The temperature history of the quenched cylindrical steel bar at any point would like to be calculated; 3-dimensional heat transfer can be analyzed using 1-dimensional axisymmetric elements as shown in Figure 1.

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

T_(R)=a₁+a₂R      (1)

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 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 r coordinate and its coordinate as shown in Figure 2;

Which are derived to obtain the following shape functions as shown;

${{S}_{i}}=\left( \frac{{{R}_{j}}-R}{{{R}_{j}}-{{R}_{i}}} \right)=\left( \frac{{{R}_{j}}-R}{L} \right)$     (2a)

1.png

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

T_(R) = S_iT_i + S_jT_j = S_(r)] {T}     (3)

where, [S_(r)] = [S_i      S_j] is a row vector matrix and $\{T\}=\left\{\begin{array}{l}T_{i} \\ T_{j}\end{array}\right\}$ is a column vector of nodal temperature of the element.

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

${{T}_{(R)}}=\left[ {{S}_{i}}\quad {{S}_{i}} \right]\left\{ \begin{array}{*{35}{l}}   {{T}_{i}}  \\   {{T}_{j}}  \\\end{array} \right\}$       (4)

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

${{\Psi }_{(e)}}=\left[ \begin{array}{*{35}{l}}   {{S}_{i}} & {{S}_{j}}  \\\end{array} \right]\left\{ \begin{array}{*{35}{l}}   {{\psi }_{i}}  \\   {{\psi }_{j}}  \\\end{array} \right\}$        (5)

where, Ψi and Ψj represent the nodal values of the unknown variable such as in our case temperature also the unknown can be deflection, or velocity etc. 

2.png

Figure 2. One-dimensional linear temperature distribution for an element (radius) line in Global Coordinate system

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 L₁ and L₂ at any point p inside the element are shown in Figure 3. from which we can write:

$L_{1}=\frac{R_{j}-R}{R_{j}-R_{i}}=\frac{l_{1}}{L}$      (6a)

$L_{2}=\frac{R-R_{i}}{R_{j}-R_{i}}=\frac{l_{2}}{L}$       (6b)

3.png

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

Since it is a 1-D 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$ or $L_{1}+L_{2}=\frac{l_{1}}{L}+\frac{l_{2}}{L}=1$      (7)

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

$S_{i}=\left(\frac{R_{j}-R}{R_{j}-R_{i}}\right)=\left(\frac{R_{j}-R}{L}\right)=L_{1}$ and

$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}$         (8)

$R=R_{j} L_{2}+R_{i} L_{1}=R_{i} S_{i}+R_{j} S_{j}$        (9)

$\frac{\partial[S]_{i}}{\partial r}=\frac{\partial L_{1}}{\partial r}=\frac{-1}{R_{j}-R_{i}}=-\frac{1}{L}$       (10)

$\frac{\partial[S]_{j}}{\partial r}=\frac{\partial L_{2}}{\partial r}=\frac{1}{R_{j}-R_{i}}=\frac{1}{L}$       (11)

2.3 Develop equation for all elements of the domain

Derivation of equation of heat transfer in axisymmetric 1-D line (radius) elements by applying the conservation of energy to a differential volume cylindrical segment as shown in Figure 4;

4.png

Figure 4. Axisymmetric element from an axisymmetric body

E_in – E_out + E_generated = E_stored             (12)

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 by;

$\frac{1}{r} \frac{d}{d r}\left(k_{r} r \frac{d T}{d r}\right)+\frac{1}{r^{2}} \frac{d}{d \theta}\left(k_{\theta} \frac{d T}{d \theta}\right)+\frac{d}{d z}\left(k_{z} \frac{d T}{d z}\right)+q=\rho c \frac{d T}{d t}$       (13)

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

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

k_z = 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.4 The assumption made in this problem was

(1) For axisymmetric situations 1-D line (radius) element, there is no variation of temperature in the Z-direction as shown in Figure 1, 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, $\left(\frac{\partial T}{\partial z}=0\right)$.

(2) For axisymmetric situations, there is no variation of temperature in the θ-direction, because it is clear from Figure 1 that the temperature distribution along the radius will be the same if the radius move with angle θ, 360^o therefore, $\left(\frac{\partial T}{\partial \theta}=0\right)$.

(3) 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, $\mid \dot{q}=0$

After simplifying, the Eq. (13) become;

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

And also known as residual or partial differential equation

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

2.5 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;

$\int_{V}[S]^{T}\{\Re\}^{(e)} d v=0$      (16)

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

$\{\mathfrak{R}\}_{(e)}$= the residual contributed by element (e) to the final system of equations.

$\frac{k}{r} \int[S]^{T} \frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r}\right) d v-\int[S]^{T} \rho c \frac{\partial T}{\partial t} d v=0$    (17)

2.6 Chain rule

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

$(f g)^{-}=f g^{-}+g f^{-}$

Therefore, $f g^{-}=(f g)^{-}-f^{-} g$ then $\frac{\partial}{\partial r}\left([S]^{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[S]^{T}}{\partial r}$

Term 1 of Eq. (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([S]^{T} r \frac{\partial T}{\partial r}\right)-r \frac{\partial T}{\partial r} \frac{\partial[S]^{T}}{\partial r}$       (18)

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

$\begin{align}  & =\overbrace{\frac{k}{r}\int{\left\{ \frac{\partial }{\partial r}\left( {{\left[ S \right]}^{T}}r\frac{\partial T}{\partial r} \right) \right\}dv}}^{A}- \\  & \overbrace{\frac{k}{r}\int{\left\{ \frac{\partial {{\left[ S \right]}^{T}}}{\partial r}r\frac{\partial T}{\partial r} \right\}dv}}^{B}-\underbrace{\int{{{\left[ S \right]}^{T}}\left\{ {{\rho }_{{}}}{{c}_{{}}}\frac{\partial T}{\partial t} \right\}dv}}_{C} \\ \end{align}$     (19)

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{align}  & \overbrace{\frac{k}{r}\int{\left\{ \frac{\partial }{\partial r}\left( {{\left[ S \right]}^{T}}r\frac{\partial T}{\partial r} \right) \right\}dv}}^{A}= \\  & \underbrace{-2\pi hrz{{\left[ S \right]}^{T}}\left[ S \right]{{\left\{ T \right\}}^{e}}}_{(1)}+\underbrace{2\pi hrz{{\left[ S \right]}^{T}}{{T}_{f}}}_{(2)}- \\  & \underbrace{2\pi rz{{\varepsilon }_{s}}\sigma {{\left[ S \right]}^{T}}\left[ S \right]{{\left\{ T \right\}}^{e}}{{\left( \left[ S \right]{{\left\{ T \right\}}^{e}} \right)}^{3}}}_{(3)}+\underbrace{2\pi rz{{\varepsilon }_{s}}\sigma {{\left[ S \right]}^{T}}T_{f}^{4}}_{(4)} \\ \end{align}$

Note: term (1) and term (3) contributed to the conductance matrix since they contains the unknown temperature {T}. Term (2) and (4) contributed to the thermal load matrix as Tf 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 can be ignore if the cooling is quenching in sea water as in our 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

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

Term C (heat stored) contributes to the Capacitance Matrix

$\underbrace{\overbrace{\frac{k}{L}\left( {{R}_{j}}+{{R}_{i}} \right)\left[ \begin{matrix}   1 & -1  \\   -1 & 1  \\\end{matrix} \right]\left\{ \begin{matrix}   Ti  \\   {}  \\   Tj  \\\end{matrix} \right\}}^{{}}}_{{}}$       (20)

$\underbrace{\overbrace{\frac{L\rho c}{6}\left[ \begin{matrix}   \left( 3{{R}_{i}}+{{R}_{j}} \right) & \left( {{R}_{i}}+{{R}_{j}} \right)  \\   \left( {{R}_{i}}+{{R}_{j}} \right) & \left( {{R}_{i}}+3{{R}_{j}} \right)  \\\end{matrix} \right]\left\{ \begin{matrix}   {{{\dot{T}}}_{i}}  \\   {}  \\   {{{\dot{T}}}_{j}}  \\\end{matrix} \right\}}^{{}}}_{{}}$      (21)

5.png

2.7 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 is necessary in all finite element analysis [11-13].

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

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

where:

$[K]=\left[K_{c}\right]+\left[K_{h}\right]$: is 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 Figure 1.

{T}: is temperature value at each node, °C.

[C]: is capacitance matrix, due to transient equation (heat stored)

$\{ \dot{T}\}$: is temperature rate for each node, °C/s.

$\{F\}=\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.8 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 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 [5, 14-17].

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 (Dt¹0), we have;

$\{ \dot{T}\} \approx\left\{\frac{T(t)-T(t-\Delta t)}{\Delta t}\right\}$      (25)

With;

$\dot{T}=$ temperature rate (°C/s); T(t)= temperature at t s (°C); T(t-Dt)=temperature at (t-Dt) s, (°C)

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

By substituting the value of $\{\dot{T}\}$ into the finite element global equation, we have that;

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

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

From Eq. (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 chromium steel bar can be obtained such as hardness and strength.

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 8a and Figure 8b;

8a.png

Figure 8a. The temperature distribution just before the steel bar becomes completely cooled

8b.png

Figure 8b. The temperature distribution at the moment that the entire steel bar becomes completely cooled after 1511s

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

From the graphs shown in Figure 6 by mathematical model and Figure 9 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 7 the cooling time, Jominy-distance and consequently the hardness of the quenched chromium steel 5147H at any point (node), even the lowest hardness point (LHP) is determined by ANSYS too, the final results shown in Table 2 and Figure 10.

9.png

Figure 9. Temperature-time graph from ANSYS

Table 2. Cooling time, Cooling rate, Jominy distance and HRC for the nodes A₁ to A₅, sea water cooled by ANSYS

10.png

Figure 10. Hardness distribution by ANSYS along AA cross section for the nodes A₁ to A₅ from the centre to the surface respectively at half the length at the centre of the quenched steel bar

From our results we found that in the mathematical model for the 1st node with W₁ in the center, we found that HRC=48.360. While in ANSYS for the same node A₁, we found that HRC=46.071. And for the nodes on the surfaces W₅ and A₅, it was found that HRC=51.486 and 48.357 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 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, this has negative consequence which can result to the deformation and failure of the component.