Simulation of Temperature Field in Steel Billets during Reheating in Pusher-Type Furnace by Meshless Method

: Using a meshless method, a simulation of steel billets in a pusher-type reheating furnace is carried out for the first time. The simulation represents an affordable way to replace the measurements. The heat transfer from the billets with convection and radiation is considered. Inside each of the billets, the heat diffusion equation is solved on a two-dimensional central slice of the billet. The diffusion equation is solved in a strong form by the Local Radial Basis Function Collocation Method (LRBFCM) with explicit time-stepping. The ray tracing procedure solves the radiation, where the view factors are computed with the Monte Carlo method. The changing number of billets in the furnace at the start and the end of the loading and unloading of the furnace is considered. A sensitivity study on billets’ temperature evolution is performed as a function of a different number of rays used in the Monte Carlo method, different stopping times of the billets in the furnace, and different spacing between the billets. The temperature field simulation is also essential for automatically optimizing the furnace’s productivity, energy consumption, and the billet’s quality. For the first time, the LRBFCM is successfully demonstrated for solving such a complex industrial problem.


Introduction
Reheating steel billets after continuous casting is necessary to temper the billets for hot rolling.Steel billets are usually heated up to 1200 • C. Five types of continuous industrial reheating furnaces are known for moving steel billets inside the furnace [1], of which the walking beam and pusher industrial reheating furnaces are the most widely used types in the steel industry.This work is focused on the simulation of temperature conditions in the gas-fired pusher-type furnace [2].These furnaces operate mostly in a continuous mode.The steel billets are charged at the furnace entry, heated while stepwise moving inside the furnace, and released at the furnace exit [3].The continuous mode of this process is interrupted at the start of charging until the furnace is filled up and at the end of the process, when the billets are only discharged without charging.
Heat is transferred to each steel billet by moving inside the furnace from the furnace entry to the furnace exit.The heat transfer mechanisms include radiation from the gases of the burners, radiation from the furnace walls, and convection [4].The heat transfer by thermal radiation accounts for more than 90% of the total heat flux impinging on the billet surface [5].The required temperature of the billets at the furnace outlet is determined by the hot rolling requirements, which consider the size of the billets, rolling speed, and steel composition [6,7].The process should ensure the required uniform temperature distribution of the billet at discharge from the furnace.
The control of the furnace, natural gas consumption, and productivity mainly depend on the knowledge of the temperature field of the billets.Measuring the temperature inside the billets with thermocouples is possible for experimental purposes only.However, performing such measurements in a standard way during production is impossible.Computer simulation thus represents an affordable way to replace the mentioned measurements.The temperature field simulation is also essential for automatically optimizing the furnace's productivity, energy consumption, and the billet's quality.Automatic optimization is performed by calculating the related object functions and using the evolutionary optimization algorithm to minimize them.Several simulation models have already been developed to predict the reheating process of steel billets in a reheating furnace.The first published model on the heating process in regenerative furnaces considering the moving billets was introduced in [8].A simplified furnace geometry was considered, and the billets were modeled as a single body located on the furnace floor.Later, Kim and Huh [9] modeled the thermal behavior of the furnace system as a steady-state heating process.The temperature field of the billets has been calculated using the finite difference technique.Hsieh et al. [10] applied the model introduced in [9] to model the turbulent flow of the combustion process using complex real-world furnace geometry.The details are provided in [1].
The solution of the temperature field in the middle of each of the billets during the reheating process is, in the present work, obtained based on the meshless Local Radial Basis Function Collocation Method (LRBFCM) with explicit time-stepping for the first time.This method has already been demonstrated for solving different steel processing steps such as continuous casting [11], hot rolling [12], cooling at the cooling bed after rolling [13], and microstructure evolution [14].

Heat Transfer Equation
The temperature field in the steel billets is governed by the heat diffusion equation where ρ, h, T, and k stand for density, mass-specific enthalpy, temperature, and thermal conductivity, respectively.It is assumed that the relation with c p (T) standing for the, in general, temperature-dependent specific heat at constant pressure.

Boundary Conditions
Neumann boundary condition is employed as a magnitude of heat flux q in the normal direction n on the boundary where the heat flux inside the billet is defined by Fourier's law q = −k∇T.The total heat flux is composed of the convective heat flux q conv , radiative heat flux from the surface of the furnace and billets inside the furnace q radS , radiative heat flux from the burning gas q radBG , and radiative heat flux from the inner gas q radG .The constants c conv , c BG , c G depend on the shape, position of the billets, and the position of the furnace's inner walls.
The billets are positioned close to the bottom of the furnace.So, the bottom of the billets only receives the radiation from the bottom of the furnace.The burning gas influences from the top of the furnace.The top of the billets receives radiation from all the contributors.The radiation from the burning gas to the billet is geometrically shadowed by the neighboring billets situated at the left or right side of the considered billet.The constants c conv , c BG , c G for the surface of billets are shown in Table 1, where d is the distance between two neighboring billets; w b is the square dimension of the billets.Since there is no obstruction from other billets, the constants are c conv = 1, c BG = 1/2, c G = 1 on the left side of the first billet, the right side of the last billet, and the edge of all billets, respectively.

Top
Left-Right Bottom The furnace has dimensions L F × W F × H F (length, width, and height).A symmetry plane is considered at W F /2.The billets have dimensions L B × W B × H B , and there is a space between the billets, denoted by S B .The distance between the bottom of the billets and furnace bottom is denoted by D FB .The distance of the center of the first billet from the start surface of the furnace is D LBe , and the center position of the last billet from the end surface of the furnace is denoted by D LBo .The billets are moving from the start to the end surface of the furnace.
The radiative heat flux from the surface of the furnace and billets is calculated in three dimensions and used to calculate the temperature in the two-dimensional middle slices of each of the billets.The surface of the furnace and the billets is divided into finite surface planes (FSPs), depicted in Figure 1. received by the i-th FSP [15].A ray tracing algorithm is used in the simulations.The number of rays is defined per square meter of the surfaces involved.The rays are randomly emitted and tracked until they are absorbed.The net radiative heat flux of the i-th FSP q radS,i is given as a difference between the emitted heat flux from the i-th surface and received heat flux from all surrounding FSPs written as where T j is the average temperature of the j-th FSP, and the sum goes over all J FSPs, where A j represents the area of the j-th FSP.The visibility between FSPs is given by a view factor F j →i ∈ [0, 1].It defines the fraction of the radiation emitted by the j-th FSP and received by the i-th FSP [15].A ray tracing algorithm is used in the simulations.The number of rays is defined per square meter of the surfaces involved.The rays are randomly emitted and tracked until they are absorbed.The temperatures T Bs , s = 1, 2, 3, 4 of the four head and four tail planes of each of the billets are calculated by solving Equation (1) four times in one dimension by assuming the Dirichlet boundary condition at the center of each quarter of the middle slice from the 2D calculation and the Neumann boundary condition in Equation (2) of each plane at the head of the billet (Figure 2a).The one-dimensional equations are calculated by the finite difference method The one-dimensional equations are calculated by the finite difference method where subscripts 0 and l indicate the initial time-step and the position of the point.Half of the billet length is represented by 201 equidistant nodes.T Cs = T 1 is the average temperature value of a quarter of the middle slice (Figure 2d) computational domain of the billet.T Bs = T 201 (in Figure 2d) is used as the average temperature of the head plane of the billet in Equation ( 6).The temperature T m of plane m (in Figure 2c) is the average temperature value of the boundary nodes belonging both to the 2D slice and plane m.The temperature T n of plane n is calculated in the same way as the temperature of plane m.The temperature of the planes i and j (in Figure 2b) is This procedure is repeated four times in order to cover all four sides of the billet.The temperature of the billets is the same at the head as at the tail because of the symmetry.The radiative heat flux from the burning gas q radBG and the radiative heat flux from the inner gas q radG are described as the difference between the emitted heat flux from the gas and the heat flux emitted from the billets where σ = 5.6704 × 10 −8 Wm −2 K −4 is the Stefan-Boltzmann constant, α is the absorption of the gas, ε is the emissivity of the gas [16], and T BG and T G stand for the temperature of the burning gas and inner gas, respectively.T BG is obtained from the balance of the heat flux at the central point of the front furnace wall when the furnace is empty.The temperature of the inner gas is considered the same as the temperature of the nearest furnace wall (see [17]).

Convective Heat Flux
The convective heat flux q conv is specified as where h and T Gas stand for the heat transfer coefficient and inner gas temperature, respectively.

Solution Procedure
The temperature field is solved by the LRBFCM [18] on a 2D slice, representing the billet cross-section, described in the Cartesian coordinate system.Instead of a grid generation over the computational domain, such as in FEM, collocation nodes without any special geometrical connectivity are uniformly distributed over the domain and at the boundaries.An unknown scalar or vector component field θ(p) at position p = p x i x + p y i y (p x , p y are Cartesian components and i x , i y are orthonormal Cartesian base vectors) is interpolated with multi-quadric RBF ψ(p) together with the augmented first-order N p = 3 polynomi- als ϑ(p) over the neighboring (N = 7) nodes (see Figure 3) Math.Comput.Appl.2024, 29, 30 6 of 16 The first and second partial derivatives are obtained from the interpolation as follows The explicit solution of the governing equation of the thermal model is considered The temperature field ( ) The details of the solution procedure are provided in [19].

Numerical Examples
This section shows the verification of the calculation of view factors (Example 1) and the thermal model (Example 2).Afterwards, the simulation results of an actual reheating x max and y max are the maximum horizontal and vertical distances between the considered nodes.p x0 , p y0 are the coordinates of the central node inside the influence domain, and c is the multi-quadric shape parameter.The collocation thus reads The first and second partial derivatives are obtained from the interpolation as follows The explicit solution of the governing equation of the thermal model is considered The temperature field T(p) in the thermal model is interpolated in the sense of Equation ( 9) The details of the solution procedure are provided in [19].

Numerical Examples
This section shows the verification of the calculation of view factors (Example 1) and the thermal model (Example 2).Afterwards, the simulation results of an actual reheating furnace are demonstrated with different process parameter configurations.Case 1 is calculated with three different ray densities.Case 2 is calculated with two different stopping times.Case 3 is calculated with different spacing between billets.In all the simulations, emissivity ε and absorption coefficients α are 0.8, and the time-step ∆t for the explicit solution is 0.1 s.The shape parameter c = 10 −3 is used in the meshless method (MM).Since the interpolation function is based on the scaled RBF, the results are not sensitive to c.A suitable shape parameter is obtained during the comparison with the reference solution and used in our related similar simulations.The details of the choice of c are shown in [19].Results in Sections 4.1 and 4.2 are compared with Finite Volume Method (FVM).

Verification of Calculation of View Factors
The computational model's accuracy strongly depends on the accuracy of the calculation of view factors.First, it has to be proven that the view factors are correctly calculated, and second, as the number of rays increases, the results should converge to the exact solution.For this purpose, we have created a sample test environment, as shown in Figure 4.A steel billet is put inside a very long furnace, and only a small section, 400 mm long, is considered.Both ends of this section are considered to have perfect reflectivity.Just like in the reheating furnace simulations, we have calculated the view factors in 3D.The heat flux results are calculated at the middle thickness.
To verify the view factors, we focused on calculating the radiative heat flux on the boundaries of the billet.We compared our heat flux results with the results obtained by a reference solution computed with [20].In this case, the only source of heat flux at the boundaries is due to radiation from the other surfaces.Each furnace wall is divided into 160 planes, and each of the billet's sides are divided into 18 planes.Each of these planes has a dimension of 10 mm × 400 mm (width × length).There is no need for division along the 400 mm length of the furnace since the results do not change on that axis.Since the wall temperatures of the furnace and billet are fixed, the only error would come from the calculation process of the view factors.The results are shown in Figure 5.The computational domain is considered to be in the middle cross-section of the billet with a size of 180 mm × 180 mm.Uniform node distribution with 25 nodes at each side of the billet is used, and three different ray densities (10 6 /m 2 , 10 7 /m 2 , 10 8 /m 2 ) are randomly sent from each square meter.It is evident from Figure 5 that the view factors are correctly calculated, and with the increase in the number of the rays, the results converge to the reference value.The reference solution is also obtained by the ray tracing algorithm, and both simulations use the same number of planes.lation of view factors.First, it has to be proven that the view factors are correctly calculated, and second, as the number of rays increases, the results should converge to the exact solution.For this purpose, we have created a sample test environment, as shown in Figure 4.A steel billet is put inside a very long furnace, and only a small section, 400 mm long, is considered.Both ends of this section are considered to have perfect reflectivity.Just like in the reheating furnace simulations, we have calculated the view factors in 3D.The heat flux results are calculated at the middle thickness.To verify the view factors, we focused on calculating the radiative heat flux on the boundaries of the billet.We compared our heat flux results with the results obtained by a reference solution computed with [20].In this case, the only source of heat flux at the boundaries is due to radiation from the other surfaces.Each furnace wall is divided into 160 planes, and each of the billet's sides are divided into 18 planes.Each of these planes has a dimension of 10 mm × 400 mm (width × length).There is no need for division along the 400 mm length of the furnace since the results do not change on that axis.Since the wall temperatures of the furnace and billet are fixed, the only error would come from the calculation process of the view factors.The results are shown in Figure 5.The computational domain is considered to be in the middle cross-section of the billet with a size of 180 × 180 mm.Uniform node distribution with 25 nodes at each side of the billet is used, and three different ray densities ( 6 10 / m , 10 / m , 10 / m ) are randomly sent from each square meter.It is evident from Figure 5 that the view factors are correctly calculated, and with the increase in the number of the rays, the results converge to the reference value.The reference solution is also obtained by the ray tracing algorithm, and both simulations use the same number of planes.

The Verification of the Simulation Model of the Reheating Furnace
Once the calculation of the view factors is verified, the results in terms of temperature should also be verified.This time, the steel billet is positioned just above the bottom, like in a real furnace.This is a time-dependent process, and simulations are run up to 8000 s. in both simulations based on five reference points that are chosen in the billet, as shown in Figure 6.Both the MM and FVM temperature results for those five points are shown in Figure 7 as a function of time.In Figure 8, the temperature fields of the steel billet at 8000 s are shown.As expected, good agreements have been once again obtained.All the parameters used in the verification processes are listed in Table 2.

The Verification of the Simulation Model of the Reheating Furnace
Once the calculation of the view factors is verified, the results in terms of temperature should also be verified.This time, the steel billet is positioned just above the bottom, like in a real furnace.This is a time-dependent process, and simulations are run up to 8000 s. in both simulations based on five reference points that are chosen in the billet, as shown in Figure 6.Both the MM and FVM temperature results for those five points are shown in Figure 7 as a function of time.In Figure 8, the temperature fields of the steel billet at 8000 s are shown.As expected, good agreements have been once again obtained.All the parameters used in the verification processes are listed in Table 2.

Sensitivity Studies of an Industrial Reheating Furnace Simulation
We focus on the central cross-section of the billets by assuming no heat transfer along the length direction, since the longitudinal dimensions of the billets are large compared to the transversal dimensions.The billets with a square cross-section of four planes), as in Figure 1b.The temperature of the surface of the furnace j T is the same as that used in [17].
All the billets move from the previous position to the next in a stopping time interval Δ t .One new billet is transported in the furnace from the start side, and one is transported out of it from the end side when the furnace is full.Natural gas burners are positioned at the top of the furnace.The billets are made of steel type 70MnVS4.The temperature-dependent material properties are considered as obtained from JmatPro [21].A total of 725 uniformly distributed points on the middle slice (see Figure 10) of each billet and subdomains with seven nodes are used in the LRBFCM.The room temperature of 25 °C is used as the initial temperature of the billets.

Sensitivity Studies of an Industrial Reheating Furnace Simulation
We focus on the central cross-section of the billets by assuming no heat transfer along the length direction, since the longitudinal dimensions of the billets are large compared to the transversal dimensions.The billets with a square cross-section of 0.18 m × 0.18 m (W B × H B ) and length of 4 m (L B ) were studied.The internal dimensions of the furnace are 31.6m × 4.8 m × 1.6 m (L F × W F × H F ).The distance of the center of the first billet from the start surface of the furnace is D LBe = 0.72 m There are 126 positions (pos in Figure 9) with a width of 0.24 m in the furnace.The distance between the bottom of the billets and furnace bottom is D FB = 0.002 m (see Figure 9).The billets are positioned centrally along the width of the furnace, so only 0.4 m space is left at both ends.In this work, the surface of the furnace is divided into 128 planes (N L = 30, N H = 4).The surface of each billet is divided into 32 planes (top, bottom, left, and right: six planes.Head and tail: four planes), as in Figure 1b.The temperature of the surface of the furnace T j is the same as that used in [17].All the billets move from the previous position to the next in a stopping time interval ∆t.One new billet is transported in the furnace from the start side, and one is transported out of it from the end side when the furnace is full.Natural gas burners are positioned at the top of the furnace.The billets are made of steel type 70MnVS4.The temperature-dependent material properties are considered as obtained from JmatPro [21].A total of 725 uniformly distributed points on the middle slice (see Figure 10) of each billet and subdomains with seven nodes are used in the LRBFCM.The room temperature of 25 • C is used as the initial temperature of the billets.In Case 1, three different view factors' results are used for different positions and compared.The Monte Carlo ray tracing method is probabilistic.Therefore, each time the view factors are calculated, the results fluctuate.The magnitude of fluctuation may be reduced by increasing the number of ray samples.When the results are summed up, parts of the errors (over/under fluctuations) are cancelled out.That is why the difference in the results in this case is so small.Moreover, the result of the Monte Carlo method for calculating the view factors depends on the number of random trials.It takes approximately 2 days to calculate the view factors when 63 billets are in the furnace by using Intel Xeon 6146 CPU.

Case 2: Different Stopping Times at Each Position
In Case 2, the space between the billets is S B = 0.27 m between each neighboring two billets.The cross-section of the furnace with the billets is shown in Figure 9.The billets stop at each position for 60 s and 90 s, respectively.Additionally, 10 6 /m 2 ray densities are used in the Monte Carlo method to calculate the view factors.The results are shown in Figure 12.We see that the temperatures for the 90 s stopping time are higher than those with the 60 s stopping time.However, both results reach 1200 • C at the furnace end, and the billet has almost a uniform temperature.We observe a slightly different behavior of the billet center at temperatures between 700 • C and 800 • C because of the changes in the material properties due to phase change.

Case 3: Different Spacing between Two Billets
In Case 3, we consider S B = 0.27 m and S B = 0.51 m, respectively.The cross-section of the furnace with two different spacings between the billets is shown in Figures 9 and 13.The billets stop at each position for 60 s.Additionally, 10 6 /m 2 ray densities are used in the Monte Carlo method to calculate the view factor.The results are shown in Figure 14.In this case, the temperatures with larger spacing are slightly higher than those with shorter spacing.This is because the billets received more energy from the walls of the furnace.

Conclusions
The steel billets in a reheating furnace are simulated by using the LRBFCM.The temperature fields depend on the radiative and convective heat fluxes on the boundaries at different positions in the furnace.The Monte Carlo method is employed to determine the radiative heat flux.To verify the presented model, two comparisons with the reference In this case, the temperatures with larger spacing are slightly higher than those with shorter spacing.This is because the billets received more energy from the walls of the furnace.

Conclusions
The steel billets in a reheating furnace are simulated by using the LRBFCM.The temperature fields depend on the radiative and convective heat fluxes on the boundaries at different positions in the furnace.The Monte Carlo method is employed to determine the radiative heat flux.To verify the presented model, two comparisons with the reference solution have been made.In the first, a small part of the billet is placed in the middle of a furnace with uniform wall temperatures.The convergence of the results of the view factors with an increasing number of rays is observed.In this way, we were sure that the boundary conditions were correctly implemented.In the second case, the billet is placed close to the bottom, and realistic wall temperatures are applied.The temperature results of the billet in the simulation after 8000 s match perfectly with the reference result.
Three test cases were performed to show the capabilities of this simulation system while using a real reheating furnace environment.These test cases are as follows: a different number of ray densities emitted for calculating the view factors, different stopping time of the billets at each position, and different spacing between two billets, respectively.We obtained the temperature results of the billets for any position and time in the furnace.
Case 1 demonstrates that there is no obvious difference in temperature results when more than 10 5 /m 2 ray samples are used in the Monte Carlo method.In the simulations, view factors are calculated once for each position independently of time and we also know from example 5.1 that the magnitude of fluctuations decreases with increasing ray density.Hence, by using the same stopping time at each position, the effect of the ray density may not be so noticeable in the results throughout the simulation.Case 2 demonstrates that a suitable billet stopping time can be found using the described model.Case 3 demonstrates that when a longer distance between each two billets is used, the increase in the temperature is more rapid, especially at the beginning.This simulator can help the steel industry to save energy and optimize productivity.It is extremely complicated to also simulate the temperature of the furnace walls simultaneously.The temperatures of the furnace walls are kept constant in present simulations; however, they would change when a billet enters or leaves the furnace.The presented temperature model of the furnace represents a crucial first step towards optimizing the billets' quality, furnace productivity, energy, and consumption.

Figure 1 .
Figure 1.A scheme of the division of the involved surfaces into FFSs.(a) A scheme of the furnace planes, and (b) a scheme of the billet planes.The temperatures

Figure 1 .
Figure 1.A scheme of the division of the involved surfaces into FFSs.(a) A scheme of the furnace planes, and (b) a scheme of the billet planes.The furnace length L F is divided into N L equidistant segments with top, bottom, front, and back surfaces.The furnace's start and the end surfaces are divided into N H vertically positioned equidistant surfaces.The billets are in the longitudinal direction divided into three segments, of which only one and a half of the central segment is considered due to the symmetry.The top, bottom, left, and right surfaces of each of the three segments are divided into two equidistant surfaces (2 × 8 in total).The temperatures of the billets are calculated in the middle of the second longitudinal segment.Each billet's head and tail are divided into four equal surfaces.The net radiative heat flux of the i-th FSP q radS,i is given as a difference between the emitted heat flux from the i-th surface and received heat flux from all surrounding FSPs written as

Figure 1 .
Figure 1.A scheme of the division of the involved surfaces into FFSs.(a) A scheme of the furnace planes, and (b) a scheme of the billet planes.The temperatures

Figure 2 .
Figure 2. A scheme of the division of the billets to calculate the average temperature for each plane of the billets.(a) A whole billet with 32 planes and middle slice, (b) half of the billet showing two planes i and j from the first segment, (c) half of the billet showing two planes m and n from the second segment, and (d) half of the billet showing the connection between the middle and head planes.

Figure 2 .
Figure 2. A scheme of the division of the billets to calculate the average temperature for each plane of the billets.(a) A whole billet with 32 planes and middle slice, (b) half of the billet showing two planes i and j from the first segment, (c) half of the billet showing two planes m and n from the second segment, and (d) half of the billet showing the connection between the middle and head planes.
horizontal and vertical distances between the considare the coordinates of the central node inside the influence domain, and c is the multi-quadric shape parameter.The collocation thus reads

Figure 3 .
Figure 3. Meshless domain and boundary discretization scheme with local influence domains consisting of 7 nodes.Influence domains with boundary nodes (red doted lines) and with only internal nodes (black dotted lines) are shown.

Figure 3 .
Figure 3. Meshless domain and boundary discretization scheme with local influence domains consisting of 7 nodes.Influence domains with boundary nodes (red doted lines) and with only internal nodes (black dotted lines) are shown.

Figure 4 .
Figure 4.A section of the furnace with a steel billet in the middle, 400 mm long.

Figure 4 .
Figure 4.A section of the furnace with a steel billet in the middle, 400 mm long.

Figure 5 .
Figure 5. Heat flux on boundary with different rays emitted from each plane with Monte Carlo method at 0 t = s.

FVMFigure 5 .
Figure 5. Heat flux on boundary with different rays emitted from each plane with Monte Carlo method at t = 0 s.

16 Figure 6 .
Figure 6.A small cross-section of the furnace with a steel billet on the bottom, 400 mm long.

Figure 7 .
Figure 7.A comparison of the results from the simulation and reference solution.

Figure 6 . 16 Figure 6 .
Figure 6.A small cross-section of the furnace with a steel billet on the bottom, 400 mm long.

Figure 7 .
Figure 7.A comparison of the results from the simulation and reference solution.

Figure 7 .
Figure 7.A comparison of the results from the simulation and reference solution.

Figure 8 .
Figure 8.The temperature field of the steel billet at 8000 s, computed by (a) the LRBFCM and (b) FVM.
and length of 4 m ( B L ) were studied.The internal dimensions of the furnace are ).The distance of the center of the first billet from the start surface of the furnace is = 0.72 m LBe D There are 126 positions (pos in Figure 9) with a width of 0 .2 4 m in the furnace.The distance between the bottom of the billets and furnace bottom is = 0.002 m FB D(see Figure9).The billets are positioned centrally along the width of the furnace, so only 0.4 m space is left at both ends.In this work, the surface of the furnace is divided into 128 planes ( surface of each billet is divided into 32 planes (top, bottom, left, and right: six planes.Head and tail:

Figure 8 .
Figure 8.The temperature field of the steel billet at 8000 s, computed by (a) the LRBFCM and (b) FVM.

Figure 9 .
Figure 9.The cross-section of the furnace with the billets.One space between neighboring billets.There are 63 billets in the furnace.

Figure 9 .
Figure 9.The cross-section of the furnace with the billets.One space between neighboring billets.There are 63 billets in the furnace.

Figure 9 .Figure 10 .
Figure 9.The cross-section of the furnace with the billets.One space between neighboring billets.There are 63 billets in the furnace.
, 10 / m , 10 / m ray densities are used for the Monte Carlo method to calculate the view factor, respectively.The results are shown in Figure11.In Case 1, three different view factorsʹ results are used for different positions and compared.The Monte Carlo ray tracing method is probabilistic.Therefore, each time the view factors are calculated, the results fluctuate.The magnitude of fluctuation may be reduced by increasing the number of ray samples.When the results are summed up, parts of the errors (over/under fluctuations) are cancelled out.That is why the difference in the results in this case is so small.Moreover, the result of the Monte Carlo method for calculating the view factors depends on the number of random trials.It takes approximately 2 days to calculate the view factors when 63 billets are in the furnace by using Intel Xeon 6146 CPU.

Figure 10 .
Figure 10.Sketch of the middle slice with uniformly distributed points.Red points are on the boundaries and black points are in the domain.

4. 3 . 1 . 16 Figure 11 .Figure 11 .
Figure 11.The results of Case 1: the temperature field in 5 representative points of the billet as a function of different ray densities per m 2 .4.3.2.Case 2: Different Stopping Times at Each Position In Case 2, the space between the billets is

16 Figure 12 .
Figure 12.The results of Case 2: the temperature field in 5 representative points of the billets as a function of different stopping times at each position.4.3.3.Case 3: Different Spacing between Two Billets In Case 3, we consider

Figure 12 .
Figure 12.The results of Case 2: the temperature field in 5 representative points of the billets as a function of different stopping times at each position.

Figure 12 . 14 .Figure 13 .Figure 13 . 16 Figure 14 .
Figure 12.The results of Case 2: the temperature field in 5 representative points of the billets as a function of different stopping times at each position.4.3.3.Case 3: Different Spacing Two Billets In Case 3, we consider

Figure 14 .
Figure 14.The results of Case 3: the temperature field in 5 representative points of the billets as a function of one (1S) and two (2S) spacings between the billets.

Table 1 .
Coefficients of heat flux for different surfaces of billets.