Validation of the Cooling Model for TMCP Processing of Steel Sheets with Oxide Scale Using Industrial Experiment Data

: To verify the mathematical model of the water-jet cooling of steel plates developed by the authors, previously performed experimental studies of the temperature of the test plates in a roller-quenching machine (RQM) were used. The calculated temperature change in the metal as it moved in the RQM was compared with the readings of thermocouples installed at the center of the test plate and near its surface. The basis of the model is the dependence of the temperatures of the ﬁlm, transition and nucleate boiling regimes on the thickness of the oxide scale layer on the cooled surface. It was found that the model correctly accounts for the oxide scale on the sheet surface, the ﬂow rates and combinations of the RQM banks used, the water temperature, and other factors. For all tests, the calculated metal temperature corresponded well with the measured one. In the experiments with interrupted cooling, the calculated temperature plots repeated the characteristic changes in the experimental curves. The main uncertainty in the modeling of cooling over a wide temperature range can be attributed to the random nature of changes in the oxide scale thickness during water cooling. In this regard, the estimated thickness of the oxide scale layer should be considered the main parameter for adapting the sheet temperature-control process. The data obtained conﬁrm the possibility of effective application of the model in the ACS of industrial TMCP (Thermo-Mechanical Controlled Process) systems.


Introduction
The current development of the global rolled steel market is characterized by a constant increase in the share of products produced by the Thermo-Mechanical Controlled Process (TMCP) [1][2][3]. The essence of this process lies in the combination of controlled hot rolling with controlled cooling in the temperature range of microstructural transformations in steel [4]. TMCP has the potential to impart such a combination of properties to rolled products (e.g., strength, ductility, toughness, cold resistance, weldability, etc.), which is not the case for other methods [5][6][7][8][9]. This feature is due to a special mechanism for the formation of fine-grained steel microstructures during phase transformations under rapid cooling conditions in combination with the deformed structure of the initial phase [4,10]. It is extremely difficult to achieve such a microstructure using other technological solutions (for instance, by microalloying) without controlled cooling. Moreover, in-line equipment for controlled cooling makes it possible to increase the efficiency of other technical and technological tools for imparting the required set of properties to the finished rolled products, including microalloying [11,12].
In practice, however, to benefit from TMCP, one must face the problems of ensuring precise control and uniformity of temperature and cooling rate, as well as flatness of the TMCP products [13]. When steel products are cooled with water, the biggest problems occur

Boundaries of Water-Cooling Zones and Parameters of Water Flow
Boundaries of water-cooling zones and parameters of water flow are used in blocks 2 and 3 of the diagram in Figure 2.
To determine the boundaries of an impingement zone produced by a spray jet (Figure 3), the authors, using analytical geometry methods, obtained the following formula for the arbitrary radius of the impact spot [35]: where θ is a polar angle of the given radius such as the angle between this radius and the polar axis of the nozzle (the polar axis is parallel to the plane of the sheet and, in most cases, is aligned with the longitudinal axis of the collector with nozzles); r θ is the nozzle outlet radius  Based on Formula (1), the area of the impact spot A [m 2 ] is calculated as follows: = 4 ( maj + 2 cos tg maj ) ( min + 2 cos tg min ) maj min (2) where maj = √1 + (sin tg ) 2 1 − (sin tg tg maj ) 2 (3) min = √1 + (cos tg ) 2 1 − (cos tg tg min ) 2 (4) dmaj and dmin are the major and the minor diameters of the nozzle outlet, respectively; φmaj and φmin are the open angles at the major and minor diameters of the nozzle outlet, respectively; χ is the nozzle rotation angle such as the angle between the polar axis and the major diameter of the nozzle outlet. Based on Formula (1), the area of the impact spot A [m 2 ] is calculated as follows: cos γ tgϕ maj d min + 2H cos γ tgϕ min k maj k min (2) where k maj = 1 + (sin χtgγ) 2 1 − sin χtgγ tgϕ maj 2 (3) k min = 1 + (cos χ tgγ) 2 1 − (cos χ tgγ tgϕ min ) 2 (4) d maj and d min are the major and the minor diameters of the nozzle outlet, respectively; ϕ maj and ϕ min are the open angles at the major and minor diameters of the nozzle outlet, respectively; χ is the nozzle rotation angle such as the angle between the polar axis and the major diameter of the nozzle outlet.
Compared to the known solution [36], Formula (2) takes into account the rotation of the nozzle around its longitudinal axis and the dimensions of the nozzle outlet orifice. Under the actual conditions of accelerated cooling systems, these factors can increase the impact area of the spray jet by up to 30%.
To calculate the speed of the spray drops in the impingement zone, the authors derived the approximate analytical solution of the differential equation of drop motion in the field of gravity, taking into account the resistance of the vapor-gas medium [37]. The key assumptions are the constancy of the inclination angle, the drag coefficient and the radius of the drop along its trajectory. These assumptions are valid for the operating conditions of spray nozzles for the accelerated cooling and quenching of metal sheets, namely: drop diameter d = 0.5 × 10 −3 . .
where the "+" sign under the square root refers to the drop going downward (i.e., for the top jet), and the "−" sign refers to the drop going upward (i.e., for the bottom jet); u is the speed of the water drop when impinging on the sheet surface [m·s −1 ]; u 0 is the speed of the water drop as it leaves the nozzle [m·s −1 ]; k u is the dimensionless speed parameter of the model (see below); H is the distance from the nozzle outlet to the sheet plane [m]; g is gravitational acceleration [m·s −2 ]; S is the complex parameter with the unit of [m −1 ]: c is the drag coefficient of the drop moving in a vapor-air medium (according to the data from [38] (p. 694), for a sphere, when Re = 10 3 . . . 10 4 , c = 0.45); d is the diameter of the drop [m]; ρ and ρ a are the density of the drop and the medium, respectively [kg·m −3 ].
The analytical expressions for the speed parameter of the model k u are obtained via adaptation to the results of numerical integration of the differential equation of the motion for the top and the bottom drops separately: − for the drops going downward: − for the drops going upward: where β 0 is the angle of the drop departure from the nozzle relative to the horizontal plane; for the axis of the spray jet, this angle is related to the jet inclination angle (see Figure 3) as The model allows us to calculate the dimensions of the supercritical water flow zone on the top surface (i.e., before the start of the hydraulic jump) for both circular and flat jets.
The supercritical flow zone from the circular jet is simplistically divided into two regions-inviscid and viscous flow ( Figure 4)-with the assumption of constant water height in the viscous flow region. Based on these simplifications, solving the differential equation of gradually varied flow in an open segmental channel with no slope [39] (p. 119), together with the ordinary hydraulic jump equation [40] (p. 653), leads to the following results [41].
− The boundary post-jump height, above which the jump becomes submerged [m]: − The boundary post-jump height, above which the jump occurs inside the inviscid flow region [m]: n is the roughness coefficient between water flow and sheet surface (for the considered conditions, n values are between about 0.007 and 0.020, depending on the laminar or turbulent flow mode, the roughness of the sheet surface and the phase state of the water. Constants in Formula (11) are obtained via approximation of the numerical solution of the differential equation for the flow height).  (11) n is the roughness coefficient between water flow and sheet surface (for the considered conditions, n values are between about 0.007 and 0.020, depending on the laminar or turbulent flow mode, the roughness of the sheet surface and the phase state of the water. Constants in Formula (11) are obtained via approximation of the numerical solution of the differential equation for the flow height).
-Radius Rs of the circular hydraulic jump [m] (subscript "s" denotes "spot"): (a) if ℎ 2 ≤ ℎ 2vis (jump occurs inside the viscous region), then (b) if ℎ 2vis < ℎ 2 ≤ ℎ 2sub (jump occurs inside the inviscid region), then (expression (13) is known as the Rayleigh formula for an ideal fluid [42]).  In contrast to the known formulas [43], the above model allows us to calculate the circular jump radius without iterative procedures.
In order to obtain a simple engineering estimate of the location of a hydraulic jump in a plane-parallel flow ( Figure 5), the authors of this article proposed to approximate the exponent y in the following well-known Pavlovsky formula for the Chezy coefficient [44] (p. 91): where C is the Chezy coefficient [m 1/2 s −1 ]; n is the roughness coefficient (the physical meaning is the same as in Formula (11)  − Radius R s of the circular hydraulic jump [m] (subscript "s" denotes "spot"): (a) if h 2 ≤ h 2vis (jump occurs inside the viscous region), then (b) if h 2vis < h 2 ≤ h 2sub (jump occurs inside the inviscid region), then (expression (13) is known as the Rayleigh formula for an ideal fluid [42]).
In contrast to the known formulas [43], the above model allows us to calculate the circular jump radius without iterative procedures.
In order to obtain a simple engineering estimate of the location of a hydraulic jump in a plane-parallel flow ( Figure 5), the authors of this article proposed to approximate the exponent y in the following well-known Pavlovsky formula for the Chezy coefficient [44] (p. 91): where C is the Chezy coefficient [m 1/2 s −1 ]; n is the roughness coefficient (the physical meaning is the same as in Formula (11)); h is the liquid-flow height [m]; y is the exponent that, in its original form, depends on roughness coefficient n and flow height h. Under conditions typical for industrial accelerated cooling systems (h = 10 −3 . . . 10 −1 m; n = 0.007 . . . 0.02) the function y = y(n, h) with an error of no more than 6% can be approximated by the function of only one variable-the roughness coefficient [45]: According to (15), the exponent y, can be considered independent of the flow height h. With this assumption, the differential equation of gradually varied flow in an open wide rectangular channel with no slope can be solved in its explicit form [45]: where L s is the length of the supercritical zone such as the distance from the border of the jet impact zone to the hydraulic jump toe [m]; α ≈ 1,05 is the kinetic energy correction factor; Q f is the specific flow rate per unit width [m 2 ·s −1 ]; h 0 is the entry flow height such as the height at the border of the jet impingement zone [m]; and h 1 is the pre-jump height [m], which is related to the post-jump height by the known ratio [44] (p. 251): where Ls is the length of the supercritical zone such as the distance from the border of the jet impact zone to the hydraulic jump toe [m]; α ≈ 1,05 is the kinetic energy correction factor; Qf is the specific flow rate per unit width [m 2 ·s −1 ]; h0 is the entry flow height such as the height at the border of the jet impingement zone [m]; and h1 is the pre-jump height [m], which is related to the post-jump height by the known ratio [44] (p. 251): The model considers four types of subcritical flow zone, i.e., the zone that lies beyond the hydraulic jump ( Figure 6): (a) "clamped layer"-the water layer between closely spaced jets on a limited area of the surface. The height of such a layer may be so great that the jets are no longer able to overcome it (the hydraulic jump becomes submerged); (b) "bounded layer"-the water layer between adjacent rows of jets or between a row of jets and a pinch roll separated by a relatively large distance; (c) "open layer"-the water layer that spreads over the surface of the sheet without any obstruction; such a layer is formed, as a rule, before the first or after the last cooling bank in the absence of special devices for removing water from the surface; The model considers four types of subcritical flow zone, i.e., the zone that lies beyond the hydraulic jump ( Figure 6): (a) "clamped layer"-the water layer between closely spaced jets on a limited area of the surface. The height of such a layer may be so great that the jets are no longer able to overcome it (the hydraulic jump becomes submerged); (b) "bounded layer"-the water layer between adjacent rows of jets or between a row of jets and a pinch roll separated by a relatively large distance; (c) "open layer"-the water layer that spreads over the surface of the sheet without any obstruction; such a layer is formed, as a rule, before the first or after the last cooling bank in the absence of special devices for removing water from the surface; (d) "shifted layer"-the water layer which is removed from the sheet by special devices (hydro or pneumatic separators). For each of the above types of water layer, the model calculates the height and the velocity of the flow; in addition, for the open layer type, it calculates the length of the subcritical zone Lf [46,47]. The method developed is based on dividing the water layer into three regions (Figure 7): direct flow, oblique flow and stagnation. Assuming that the water flow varies gradually, the equations of motion are solved for the median streamline in the direct and oblique flow regions. As a boundary condition in the oblique flow region, the equality-to-unity of the Froude number on the lateral edge of the sheet is assumed (as on the weir threshold [48] (p. 391)). For each of the above types of water layer, the model calculates the height and the velocity of the flow; in addition, for the open layer type, it calculates the length of the subcritical zone L f [46,47]. The method developed is based on dividing the water layer into three regions (Figure 7): direct flow, oblique flow and stagnation. Assuming that the water flow varies gradually, the equations of motion are solved for the median streamline in the direct and oblique flow regions. As a boundary condition in the oblique flow region, the equality-to-unity of the Froude number on the lateral edge of the sheet is assumed (as on the weir threshold [48] (p. 391)).
In the direct-flow region, continuity of the layer height at the border with the oblique flow region is assumed as a boundary condition. At the same time, the discontinuity of velocity at this border is allowed.
As a result, the authors have proposed a procedure that allows for the calculation of the height and speed of the water layer at all reference points of the subcritical flow zone.
The key parameter of this procedure is the "angular runoff coefficient" k ϕ , which means the sinus of the runoff angle, i.e., k ϕ ≡ sin ϕ. In the direct-flow region, continuity of the layer height at the border with the oblique flow region is assumed as a boundary condition. At the same time, the discontinuity of velocity at this border is allowed.
As a result, the authors have proposed a procedure that allows for the calculation of the height and speed of the water layer at all reference points of the subcritical flow zone.
The key parameter of this procedure is the "angular runoff coefficient" kφ, which means the sinus of the runoff angle, i.e., ≡ sin . For a bounded layer, the calculation of kφ is performed using a chain of formulas: − auxiliary parameter Φ: where n is the roughness coefficient (the physical meaning is the same as in Formula (11)); g is gravitational acceleration [m·s −2 ]; Qf is the specific flow rate per unit width [m 2 ·s −1 ]; Lf is the length of the subcritical zone [m]; and y is the exponent (15); − Friction parameter Ω Φ : − Shape parameter δ: For a bounded layer, the calculation of k ϕ is performed using a chain of formulas: − auxiliary parameter Φ: where n is the roughness coefficient (the physical meaning is the same as in Formula (11) − Friction parameter Ω Φ : − Shape parameter δ: − Angular runoff coefficient where B p is the width of the sheet [m]; x = 0.5 when Ω Φ /δ < 1; and x = 0 when Ω Φ /δ ≥ 1.
With a known value of k ϕ , the height and speed of the flow at all key points of the water layer can be calculated-for example, where u 22 and h 22 are the speed [m·s −1 ] and the height [m] of the water at the lateral edge, respectively (expression (23) follows from the equality of the Froude number-to-unity); To calculate the length of the subcritical flow zone L f in the case of the open water layer, the iteration procedure is realized. The essence of this procedure consists of solving, with respect to the desired length L f , the system of algebraic equations obtained for the bounded layer (with the angular flow coefficient according to (21)), together with the following equation obtained for the open layer: where V fs is the volume of water that runs off from one lateral

Heat Flux in the Water-Cooling Zones
Heat flux in the water-cooling zones is used in blocks 7-9 of the diagram in Figure 2.
The calculation of heat flux in the water-cooling zones at any surface temperature is performed using the concept of reference points of the boiling curve [50], which allows us to avoid the breaking of solutions at critical points in the boiling regime change. According to this concept, on a boiling curve [51,52] The heat flux at any surface temperature is calculated via linear interpolation between the reference points of the boiling curve. For example, if the current surface temperature refers to the transient boiling section, the heat flux is calculated as follows (the interpolation straight line for this section is shown in Figure 8): where q(t w ) is the heat flux as a function of the surface temperature [W·m −2 ]; t w is the sheet surface temperature (subscript w from "wall") [  To implement the described approach, the coordinates (i.e., temperature and heat flux density) of all the reference points of the boiling curves for each water-cooling zone are calculated in the model.
The End of Film Boiling temperature tEFB [°C] is calculated by the formula proposed by the authors: where tint is the temperature at the liquid-solid interface at the EFB point [°C]: ts is the liquid saturation temperature [°C]; tf is the current liquid temperature away from the cooled surface [°C]; tf,lim is the practical limit of the maximum achievable liquidphase temperature when impulse heating [°C] (for water at atmospheric pressure, we accepted tf,lim = 300 °C according to [53]; and ξ is the subcooling effect coefficient. For water at atmospheric pressure: Δtsub is the value of liquid subcooling [°C]: ρf and cf are the density [kg·m −3 ] and true isobaric specific-mass-heat capacity [J·kg −1 K −1 ] of liquid at the temperature of tf, respectively; λw, ρw and cw are the thermal conductivity [W·m −1 ·K −1 ], density [kg·m −3 ] and true isobaric specific-mass-heat capacity [J·kg −1 K −1 ] of the subsurface layer of the cooled sheet at the temperature of tw, respectively (see below); Re is the Reynolds number in the form: To implement the described approach, the coordinates (i.e., temperature and heat flux density) of all the reference points of the boiling curves for each water-cooling zone are calculated in the model.
The End of Film Boiling temperature t EFB [ • C] is calculated by the formula proposed by the authors: where t int is the temperature at the liquid-solid interface at the EFB point [ • C]: t s is the liquid saturation temperature [ • C]; t f is the current liquid temperature away from the cooled surface [ • C]; t f,lim is the practical limit of the maximum achievable liquidphase temperature when impulse heating [ • C] (for water at atmospheric pressure, we accepted t f,lim = 300 • C according to [53]; and ξ is the subcooling effect coefficient. For water at atmospheric pressure: ∆t sub is the value of liquid subcooling [ • C]: ρ f and c f are the density [kg·m −3 ] and true isobaric specific-mass-heat capacity [J·kg −1 K −1 ] of liquid at the temperature of t f, respectively; λ w , ρ w and c w are the thermal conductivity [W·m −1 ·K −1 ], density [kg·m −3 ] and true isobaric specific-mass-heat capacity [J·kg −1 K −1 ] of the subsurface layer of the cooled sheet at the temperature of t w , respectively (see below); Re is the Reynolds number in the form: u f is the liquid speed [m·s −1 ]; ν f is the kinematic viscosity of the liquid [m 2 ·s −1 ]; σ fs is the surface-tension coefficient of liquid at the saturation temperature [N·m −1 ]; g is the gravitational acceleration [m·s −2 ]; ρ fs and ρ vs are the density [kg·m −3 ] of liquid and vapor at the saturation temperature, respectively; and k τ is the model parameter [W·m −1 ·K −1 ] that specifies the value of turbulent thermal conduction of the liquid depending on the Reynolds number; for the conditions in question, k τ = 0.016 W·m −1 ·K −1 .
Formula (27) follows the well-known solution for the temperature, which is set at the boundary of two semi-infinite rods at the moment of their contact [54] (p. 401), with the following assumptions: (1) the molecular thermal conductivity of a liquid is negligibly small compared to its turbulent thermal conductivity, and (2) the turbulent thermal conductivity is a power function of the Reynolds number: where λ τ is the turbulent thermal conductivity [W·m −1 ·K −1 ]. Expression (29) is a linear interpolation of the dependence of the subcooling effect coefficient ξ on water subcooling in a possible variation range from 0.4 to 0.8 (this range for ξ is justified by the authors by analyzing the frequency of the potential contacts of liquid with the solid surface during oscillations of the liquid-gas interface phases near the EFB temperature).
In the presence of oxide scale on the steel surface, the thermophysical properties of the cooled subsurface layer used in the Formula (27) where ρ met is the density of steel; ρ sc is the apparent (i.e., including pores) density of oxide scale; and ψ sc is the volume fraction of oxide scale (including pores) in the subsurface layer; − thermal conductivity [W·m −1 ·K −1 ]: where λ met is the thermal conductivity of steel, and λ sc is the thermal conductivity of oxide scale; − true isobaric specific-mass-heat capacity [J·kg −1 K −1 ]: where c met is the true isobaric specific-mass-heat capacity of steel, and ε sc is the mass fraction of oxide scale in the subsurface layer. For steel, well-known formulas approximating the dependence of thermophysical properties on temperature are used (for example, [55], see Appendix B for the properties of the test plates). For scale, the formulas described in by the authors in [56][57][58][59] are taken.
The End of Transient Boiling temperature t ETB [ • C] is calculated by the formula: where t s is the water saturation temperature [ • C], and ∆t 0 ETB is the overheating of the surface at the ETB point for a subcooled water at free convection, i.e., without taking into account the water speed [ • C]: s sub is the coefficient of the influence of water subcooling on the first critical heat flux at free convection [60] (p. 205): (38) c f is the average isobaric specific-mass-heat capacity of water in the range from t f to t s [J·kg −1 K −1 ]; r v is the latent heat of vaporization of water [J·kg −1 ]; ρ fs and ρ vs are the density of water and vapor at the saturation temperature, respectively [kg·m −3 ]; and k w is the coefficient of the influence of water flow velocity on ETB temperature: u f is the water speed relative to the sheet surface [m·s −1 ]. Formula (37) is based on the assumption that the heat-transfer coefficient at the first critical point (i.e., at the ETB point) when the liquid is supercooled remains the same as for a saturated liquid. The validity of this assumption is confirmed, for example, by the experimental data cited in [61,62]. Formula (39) is based on the assumption that the volume vapor content of the near-wall layer in the first critical state (at the ETB point) does not depend on the fluid flow speed.
The heat flux at the reference points is calculated on the basis of known methods: − At the End of Film Boiling (i.e., at the EFB-point), according to Wang-Shi [63]: where Nu EFB = q EFB x/(λ f ∆t sub ) and Re where c f0 is the coefficient of friction between the liquid and the surface (evaluated by [65] (p. 289)); ϕ * is the vapor content in the near-wall two-phase layer in the first critical state (we evaluated ϕ * (1 − ϕ * ) ≈ 0.17 based on [64] (p. 312)); r v is the latent heat of the vaporization of water [J·kg −1 ]; s sub is the water subcooling coefficient (see (38)); and ρ fs and ρ vs are the density of water and vapor at the saturation temperature, respectively [kg·m −3 ].

Microstructure
Microstructure is used in block 10 of the diagram in Figure 2. The evolution of the microstructure in the sheet during cooling is predicted from isothermal Time-Temperature-Transformation (TTT) diagrams [67][68][69] (p. 356), using a step-by-step calculation scheme based on the additivity rule [70][71][72].
In order to use the TTT diagrams originally presented in graphical form, the authors modified the method of approximation by reference points proposed in [73,74] (p. 43). The essence of such modification is that in addition to the two reference points-at the beginning and at the "nose" of the C-shaped curve-another reference point below the "nose" is introduced (Figure 9). This makes it possible to significantly improve the approximation accuracy, especially in the lower part of the C-curve, i.e., at temperatures predominantly related to the bainite transformation. The method of 3-point approximation proposed by the authors allows each of the five basic C-shaped curves of the austenite isothermal transformation (start of ferrite formation, start and finish of pearlite formation, start and finish of bainite transformation) to be described by a single-parameter function of the following form: where y is the function of the logarithm of isothermal holding time τ (as an independent coordinate of the TTT diagram), w is the function of temperature t (as a dependent coordinate) and c is the parameter to be calculated by the following formula: The functions of the coordinates in Formulas (42) and (43) are as follows: where S means the logarithm of time: and U means the inverse temperature: Subscript "0" refers to the upper reference point (at the beginning of the C-curve), "N" refers to the reference point at the C-curve "nose" and 1 refers to the lower reference point (see Figure 9).
Within the above scheme of using TTT diagrams to predict the steel microstructure during cooling, the kinetics of isothermal austenite decomposition into ferrite and pearlite where y is the function of the logarithm of isothermal holding time τ (as an independent coordinate of the TTT diagram), w is the function of temperature t (as a dependent coordinate) and c is the parameter to be calculated by the following formula: The functions of the coordinates in Formulas (42) and (43) are as follows: where S means the logarithm of time: and U means the inverse temperature: Subscript "0" refers to the upper reference point (at the beginning of the C-curve), "N" refers to the reference point at the C-curve "nose" and 1 refers to the lower reference point (see Figure 9).
Within the above scheme of using TTT diagrams to predict the steel microstructure during cooling, the kinetics of isothermal austenite decomposition into ferrite and pearlite is calculated using the Kolmogorov-Johnson-Mehl-Avrami (KJMA) equation [75] (p. 128, 496), the bainite transformation via the Austin-Rickett equation [76,77] and the martensite transformation via the Koistinen-Marburger equation [78].

Temperature Distribution across the Thickness of the Sheet
Temperature distribution across the thickness of the sheet is used in block 11 of the Diagram in Figure 2.
A procedure for numerically solving the one-dimensional unsteady thermal conductivity equation for a flat, metal plate with oxide scale on both wide surfaces, with boundary conditions of the third kind, has been implemented.
The design scheme is shown in Figure 10.  Figure 2.
A procedure for numerically solving the one-dimensional unsteady thermal conductivity equation for a flat, metal plate with oxide scale on both wide surfaces, with boundary conditions of the third kind, has been implemented.
The design scheme is shown in Figure 10. The designations in this figure are as follows: H is the sheet thickness, including oxide scale; hsb and hst are the thickness of the oxide scale layer on the bottom and top surface, respectively; hm is the thickness of the metal body without oxide scale; nsb ≥ 3 is the quantity of nodes of the computation grid inside the bottom oxide scale layer (the quantity of nodes inside the top oxide scale layer should also be not less than 3); nm is the quantity of nodes of the computation grid inside the metal body; n is the total quantity of nodes across the thickness of the sheet with oxide scale (the value of n is determined automatically, based on the thickness of the elementary layer Δx specified in the initial data); i is the number of the current node; δb and δt are the thickness of internal elementary layer of the bottom and top oxide scale, respectively; and Δx is the thickness of the internal elementary layer of the metal body. The finite-difference equations to be solved at any time step, except the initial one (i.e., at time τ > 0), are as follows: − Inside the metal body ( n sb + 1 ≤ i ≤ n sb + n m − 2): where t k i is the temperature in the i-th node at the k-th point in time [ • C]; ∆τ is the time step between the k-th and k + − Inside the bottom scale layer ( 2 ≤ i ≤ n sb − 1): − Inside the top scale layer ( n sb + n m ≤ i ≤ n − 1): − Conjugation conditions between the bottom scale layer and the metal (i = n sb ): − Conjugation conditions between the top scale layer and the metal (i = n − n m + 1): − Boundary conditions at the bottom scale surface (i = 1): where α k b is the heat-transfer coefficient at the bottom surface of the sheet at the k-th point in time [W·m −2 ·K −1 ], and t k+1 ab is the ambient temperature at the bottom surface at the k + 1-th point in time [ • C]; − Boundary conditions at the top scale surface (i = n): where α k t is the heat-transfer coefficient at the top surface at the k-th point in time [W·m −2 ·K −1 ], and t k+1 at is the ambient temperature at the top surface at the k+1-th point in time [ • C].
The system of n Equations (50)-(56) is solved using an implicit finite-difference scheme via the Thomas (tridiagonal matrix) algorithm [80] (p. 83). The step of the calculation grid over the thickness of the sheet ∆x is about 0.1 mm, and the basic time step ∆τ is about 0.1 s (additionally divided by the borders between the designed cooling zones on the top and bottom surfaces of the sheet).

Experimental Procedure
In order to check the adequacy and to adapt the model, the authors used the data from previously conducted temperature measurements across the thickness of the steel plate during processing in a roller-quenching machine (RQM) of NKMZ design [81].
The technique of the experimental studies is described in detail in [82,83]. It consisted of the application of a measuring complex in the form of a test plate with embedded thermocouples, connected to the data-collection and recording system ( Figure 11).
The experiments involved two test plates of structural carbon steel 45  thermocouples were 3 and 2.5 mm from the top and bottom surfaces, respectively. A total of two series of five tests each were carried out. The first plate was used in test Nos. 1-5, and the second plate in test Nos. 6-10. Figure 11. Schematic of the measuring complex. Figure 11. Schematic of the measuring complex.
The scheme of each test involved charging a cold plate with thermocouples into a heating furnace from the RQM side (i.e., through the furnace outlet window), heating it to a preset temperature, discharging it from the furnace at a constant speed and cooling it in the roller-quenching machine. The heating temperature and RQM operating regime in each test are given in Table 1.
The signals from the thermocouples were recorded continuously during the whole process of heating the test plate in the furnace, its transportation and subsequent cooling in the RQM. As an example, Figure 13 shows the plots of temperature measured at different control points of the plate during a single test (No. 2). The temperature measurements in the furnace (i.e., for the time from t 0 to t 1 ) were used to calculate the thickness of the oxide scale layer on the plate surface. The measurement data after the plate was discharged from the furnace (after t 1 ) was compared with the results of the cooling model calculation.  The scheme of each test involved charging a cold plate with thermocouples into a heating furnace from the RQM side (i.e., through the furnace outlet window), heating it to a preset temperature, discharging it from the furnace at a constant speed and cooling it  Differential Equation (57) is non-linear, because the oxidation rate constant changes in time with temperature change [86] (p. 59): where T is an absolute temperature in degrees Kelvin, and A and B are material parameters. For steel 45, the values of these parameters were taken according to [87] (p. 146) as A = 11.41 kg·m −2 s −0.5 and B = 8274 K. To numerically solve differential Equation (57), we used Euler's method via an explicit finite-difference scheme: where Yi and Yi+1 are the specific mass of oxide scale at the beginning and at the end of the i-th time interval [kg·m −2 ]; Ti is the surface temperature measured at the beginning of the i-th time interval, K; and ∆τ is a measurement cycle, which is equal to 0.07 s. The transition from specific mass to scale thickness was carried out as follows: where hi is the oxide scale thickness corresponding to its specific mass Yi [m]; ηsc is the porosity of oxide scale in fractions of one; and ρsc is its true density [kg·m −3 ]. For furnace scale, formed when heated to 850-980 °С, we assumed ηsc = 0.15 [88] and ρsc = 5500 kg·m −3 [56].
Using the described procedure, we calculated the thickness of oxide scale at the exit of the furnace before each test. In this case, given that visually, no coarse scale was observed on the top surface of the plates after treatment in RQM, the initial thickness of the oxide layer upon loading into the furnace was taken to be zero (more precisely, 1 micron

Estimation of Oxide Scale Thickness on the Surface of the Test Plates
We accepted the diffusion mechanism of oxide scale growth, according to which the rate of its mass increase is inversely proportional to the mass of already-formed scale [84] (p. 49): where Y is the scale mass per unit surface [kg·m −2 ]; τ is time [s]; and K is the oxidation rate constant [kg·m −2 s −0.5 ].
The integration of (57) gives a parabolic equation of scale growth [85]: where Y 0 is the scale mass per unit surface [kg·m −2 ] at τ = 0 s. Differential Equation (57) is non-linear, because the oxidation rate constant changes in time with temperature change [86] (p. 59): where T is an absolute temperature in degrees Kelvin, and A and B are material parameters. For steel 45, the values of these parameters were taken according to [87] (p. 146) as A = 11.41 kg·m −2 s −0.5 and B = 8274 K. To numerically solve differential Equation (57), we used Euler's method via an explicit finite-difference scheme: where Y i and Y i+1 are the specific mass of oxide scale at the beginning and at the end of the i-th time interval [kg·m −2 ]; T i is the surface temperature measured at the beginning of the i-th time interval, K; and ∆τ is a measurement cycle, which is equal to 0.07 s. The transition from specific mass to scale thickness was carried out as follows: where h i is the oxide scale thickness corresponding to its specific mass Y i [m]; η sc is the porosity of oxide scale in fractions of one; and ρ sc is its true density [kg·m −3 ]. For furnace scale, formed when heated to 850-980 • C, we assumed η sc = 0.15 [88] and ρ sc = 5500 kg·m −3 [56].
Using the described procedure, we calculated the thickness of oxide scale at the exit of the furnace before each test. In this case, given that visually, no coarse scale was observed on the top surface of the plates after treatment in RQM, the initial thickness of the oxide layer upon loading into the furnace was taken to be zero (more precisely, 1 micron in order to avoid dividing by zero in the Formula (60) for the first cycle of measurements). The calculated values of the scale thickness at the exit of the furnace are shown in Table 2.

Results
According to the authors' cooling model, calculations were performed for the conditions of each of the ten tests. The following values of the basic input data (in addition to the data presented in Tables 1 and 2 and Figure 12) were taken: − The design and layout characteristics of RQM, according to the technical documentation of the equipment manufacturer (basic characteristics are given in [31,81,89]); − The thermophysical properties of the test plate as a function of temperature, according to the formulas [55] for medium-carbon steel (see Appendix B); − The thermophysical properties of oxide scale as a function of temperature, according to the authors' formulas [56][57][58][59]; − The oxide scale thickness was assumed to be constant throughout the cooling period of each test and equal to the values given in Table 2.
The results of the calculations in comparison with the corresponding experimental data are presented in the form of graphs of temperature changes over time. For example, Figures 14-16 Figure 14. Graphs of cooling in test No. 3. Numbers without an asterisk denote experimental cu according to the thermocouple data: 1-at the center of the thickness in the middle of the p width; 2a-at the top in the middle of the width; 2b-at the top at the edge; 3a-at the bottom in middle of the width. Numbers with one asterisk indicate calculated curves for the correspond control points across the thickness of the plate: 1* and 1**-at the center of the sheet thickness w regard to oxide scale and without regard to oxide scale, respectively; 2* and 2**-at a depth of 3 from the top surface with regard to oxide scale and without regard to oxide scale, respectively and 3**-at a depth of 2.5 mm from the bottom surface with regard to oxide scale and without reg to oxide scale, respectively. Point A is the first inflection of the calculated temperature graph, co sponding to the beginning of water cooling from below, and point B is the same in the test (its for several seconds is probably due to the thermal inertia of the thermocouple). Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 1* and 1**-at the center of the sheet thickness with regard to oxide scale and without regard to oxide scale, respectively; 2* and 2**-at a depth of 3 mm from the top surface with regard to oxide scale and without regard to oxide scale, respectively; 3* and 3**-at a depth of 2.5 mm from the bottom surface with regard to oxide scale and without regard to oxide scale, respectively. Point A is the first inflection of the calculated temperature graph, corresponding to the beginning of water cooling from below, and point B is the same in the test (its lag for several seconds is probably due to the thermal inertia of the thermocouple).

Figure 14.
Graphs of cooling in test No. 3. Numbers without an asterisk denote experimental curves according to the thermocouple data: 1-at the center of the thickness in the middle of the plate width; 2a-at the top in the middle of the width; 2b-at the top at the edge; 3a-at the bottom in the middle of the width. Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 1* and 1**-at the center of the sheet thickness with regard to oxide scale and without regard to oxide scale, respectively; 2* and 2**-at a depth of 3 mm from the top surface with regard to oxide scale and without regard to oxide scale, respectively; 3* and 3**-at a depth of 2.5 mm from the bottom surface with regard to oxide scale and without regard to oxide scale, respectively. Point A is the first inflection of the calculated temperature graph, corresponding to the beginning of water cooling from below, and point B is the same in the test (its lag for several seconds is probably due to the thermal inertia of the thermocouple). Figure 15. Graphs of cooling in test No. 5. The notations are the same as in Figure 14. The analysis of the data obtained shows the following.
(1) For all tests, the calculated temperature at the RQM outlet corresponds to the measured temperature, with deviations not exceeding 10 °C. Therefore, the total heat loss of the sheet is taken into account correctly, which indirectly confirms the adequacy of the calculation of the first and second critical surface temperatures corresponding to the changes in the film, transition and nucleate boiling regimes (i.e., the EFB and ETB temperature in Figure 8).  The analysis of the data obtained shows the following.
(1) For all tests, the calculated temperature at the RQM outlet corresponds to the measured temperature, with deviations not exceeding 10 • C. Therefore, the total heat loss of the sheet is taken into account correctly, which indirectly confirms the adequacy of the calculation of the first and second critical surface temperatures corresponding to the changes in the film, transition and nucleate boiling regimes (i.e., the EFB and ETB temperature in Figure 8). mocouple readings below 400-450 • C, the calculated temperature is much higher than the experimental one. In our opinion, this is due to the assumption made in the simulation of constant thickness of the scale during the entire cooling period. This assumption is generally not true, because during accelerated cooling, oxide scale can crack and be removed (partially or completely) from the sheet surface by water jets. At the same time, the influence of scale depends on the water boiling regime: at high temperatures corresponding to film and transient boiling, oxide scale, as a rule, increases the intensity of heat transfer to the surface; at stable nucleate boiling, on the contrary, reduces it [90]. Therefore, if oxide scale is removed from the sheet surface after stable nucleate boiling is achieved, it is accompanied by an increase in the intensity of cooling. In spray cooling, stable nucleate boiling of water usually begins at 200-250 • C at the surface, which, at the corresponding heat flux values, approximately corresponds to a temperature of 400-450 • C at a depth of 2.5-3 mm (for a plate 30 mm thick). This explains the above-mentioned overestimation of the calculated temperature in the noted experiments with surface thermocouple readings below 400-450 • C. This effect is also confirmed by model calculations. For example, in Figure 14, it is seen that at indications of surface thermocouples above 400-450 • C, the slope angle of the experimental plots is close to the slope angle of the design graphs with scale, and below this boundary, to the slope angle of design graphs without scale. (6) In practically all cases, there is a "lag" for 1-3 s of the experimental curves from the calculated ones at the very beginning of intensive cooling (see, for example, points A and B in Figure 14). This, most likely, can be explained by the thermal inertia of thermocouples [91], which manifests, to the greatest extent, as a sharp change in metal temperature. (7) To quantify the "proximity" of the experimental and calculated graphs, the value of the average cooling rate at a certain temperature interval was used. Table 3 summarizes the average cooling rate in the three typical temperature ranges: 800-400 • C, 400-200 • C and 200-100 • C. The value in each cell of this table is obtained by averaging over all experiments. It can be seen that calculated cooling rates are, in general, somewhat lower than experimentally determined (on average, by 12-20% at different temperature intervals). Higher cooling rates in the experiments (than in the simulation) can be explained by the factors mentioned above: the thermal inertia of thermocouples (see item 6) and scale removal from the plate surface during cooling in RQM (see item 5). In this case, if the thermal inertia of thermocouples affects only the "apparent" cooling rate, the reduction in the scale layer affects the actual intensity of heat transfer.

Discussion
The obtained results confirm that the main factor that introduces uncertainty into the process of the accelerated cooling of metal in production conditions is oxide scale on its surface. And the point is not that scale, especially peeled scale, distorts the readings of the workshop pyrometers. Much more important is the fact that oxide scale changes the real intensity of cooling because it shifts the boiling curve. When cooling a surface with oxide scale, the boiling curve, in general, is a composition of two curves ( Figure 17) [90]: at the initial stage, it follows the boiling curve for the surface covered with a continuous layer of scale, and then, due to scale cracking, it shifts in the direction of the boiling curve for a clean surface. Therefore, in the presence of scale, in general, two more reference points are added to the boiling curve: SHX-Start sHift due to oXide layer and EHX-End sHift due to oXide layer. surface. And the point is not that scale, especially peeled scale, distorts the readings of the workshop pyrometers. Much more important is the fact that oxide scale changes the real intensity of cooling because it shifts the boiling curve. When cooling a surface with oxide scale, the boiling curve, in general, is a composition of two curves ( Figure 17) [90]: at the initial stage, it follows the boiling curve for the surface covered with a continuous layer of scale, and then, due to scale cracking, it shifts in the direction of the boiling curve for a clean surface. Therefore, in the presence of scale, in general, two more reference points are added to the boiling curve: SHX-Start sHift due to oXide layer and EHX-End sHift due to oXide layer. The moment at which this "jump" from one boiling curve to the other occurs (see Figure 17) will determine both the course of the cooling process and the final temperature of the metal. Unfortunately, the cracking and removal of oxide scale from the sheet surface is very difficult to predict because it is a random event, which creates inevitable uncertainty when modeling cooling over a wide temperature range. Moreover, in practice, the The moment at which this "jump" from one boiling curve to the other occurs (see Figure 17) will determine both the course of the cooling process and the final temperature of the metal. Unfortunately, the cracking and removal of oxide scale from the sheet surface is very difficult to predict because it is a random event, which creates inevitable uncertainty when modeling cooling over a wide temperature range. Moreover, in practice, the initial thickness of the surface scale layer is also usually unknown.
This predetermines two main directions of further research: (1) predicting the thickness of scale during TMCP processing, taking into account its initial thickness, the strength of adhesion with the base metal and the dynamics of oxide thickness changes in the cooling unit; and (2) using the thickness of the oxide scale layer as the main parameter for the on-line adaptation of the temperature model in the ACS of accelerated cooling and quenching units.

Conclusions
Comparison with experimental data shows that the new version of the authors' mathematical model of jet cooling adequately takes into account the main design and technological factors determining the sheet temperature; this includes water consumption by the sections, the sequence of their inclusion, the temperature of cooling water and the thickness of scale on the surface. This makes it possible to recommend a new version of the model for use in control systems for accelerated sheet-metal cooling in a wide temperature range, up to full hardening below 100 • C. It is shown that the main disturbing factor that introduces uncertainty into the simulation results under real conditions is the thickness of the oxide scale layer on the sheet surface and the random nature of its removal in the process of jet cooling. Therefore, the rated thickness of the scale layer should be considered the main parameter for adapting the process of controlling the sheet temperature.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. service of the cooling model, and to A.L. Ostapenko for helpful advice when discussing the results of this work.

Conflicts of Interest:
The authors declare no conflict of interest Appendix A Figure A1. Graphs of cooling in test No. 1. Numbers without an asterisk denote experimental curves according to the thermocouple data: 1-at the center of the thickness in the middle of the plate width; 2a-at the top in the middle of the width; 2b-at the top at the edge; 3a-at the bottom in the middle of the width. Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 1*-at the center of the sheet thickness; 2*-at a depth of 3 mm from the top surface; 3*-at a depth of 2.5 mm from the bottom surface.  Figure A1. Figure A2. Graphs of cooling in test No. 2. The notations are the same as in Figure A1. Figure A2. Graphs of cooling in test No. 2. The notations are the same as in Figure A1. Figure A3. Graphs of cooling in test No. 4. The notations are the same as in Figure A1. Figure A3. Graphs of cooling in test No. 4. The notations are the same as in Figure A1. Figure A2. Graphs of cooling in test No. 2. The notations are the same as in Figure A1. Figure A3. Graphs of cooling in test No. 4. The notations are the same as in Figure A1. Figure A4. Graphs of cooling in test No. 6. Numbers without an asterisk denote experimental curves according to the thermocouple data: 2a-at the top in the middle of the width; 2b-at the top at the edge, 3a-at the bottom in the middle of the width; 3b-at the bottom at the edge. Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 2*-at a depth of 3 mm from the top surface; 3*-at a depth of 2.5 mm from the bottom surface. Figure A4. Graphs of cooling in test No. 6. Numbers without an asterisk denote experimental curves according to the thermocouple data: 2a-at the top in the middle of the width; 2b-at the top at the edge, 3a-at the bottom in the middle of the width; 3b-at the bottom at the edge. Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 2*-at a depth of 3 mm from the top surface; 3*-at a depth of 2.5 mm from the bottom surface. Figure A5. Graphs of cooling in test No. 7. The notations are the same as in Figure A4.   Figure A4. Figure A4. Graphs of cooling in test No. 6. Numbers without an asterisk denote experimental curves according to the thermocouple data: 2a-at the top in the middle of the width; 2b-at the top at the edge, 3a-at the bottom in the middle of the width; 3b-at the bottom at the edge. Numbers with one asterisk indicate calculated curves for the corresponding control points across the thickness of the plate: 2*-at a depth of 3 mm from the top surface; 3*-at a depth of 2.5 mm from the bottom surface. Figure A5. Graphs of cooling in test No. 7. The notations are the same as in Figure A4. Figure A6. Graphs of cooling in test No. 8. The notations are the same as in Figure A4. Figure A6. Graphs of cooling in test No. 8. The notations are the same as in Figure A4.

Appendix B
The temperature dependence of the thermal physical properties for the material of the test plates (steel 45) according to the approximation formulae from [55] for mediumcarbon steel: Figure A7. Graphs of cooling in test No. 9. The notations are the same as in Figure A4.