Transient thermal response with nonlocal radiation of a blast furnace main trough

Applied Mathematical Modelling 105 (2022) 197–225 Contents lists available at ScienceDirect Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm Transient thermal response with nonlocal radiation of a blast furnace main trough P. Barral a , b , c , L.J. Pérez-Pérez a , c , P. Quintela a , b , c , ∗ a Department of Applied Mathematics, Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain b Technological Institute for Industrial Mathematics (ITMATI), Santiago de Compostela 15782, Spain c Instituto de Matemáticas (IMAT), Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain a r t i c l e i n f o Article history: Received 19 July 2021 Revised 14 December 2021 Accepted 19 December 2021 Available online 29 December 2021 Keywords: Blast furnace trough Heat transfer Transient numerical simulation Nonlocal radiation Adaptive time-stepping Manufactured solution test a b s t r a c t A mathematical model for the transient thermal behaviour of the main trough of a blast furnace (BF) is presented. The proposed model consists of the transient heat equation with mixed radiation-convection boundary conditions to model the cooling process. The heat equation is coupled with an integral equation posed on the inner boundary, which models the radiative heat exchange on the internal cavity formed by the trough and the refractory cover placed over the trough. The main scope of this work is to address the evolution of the temperature field during a full BF tapping. A reliable algorithm, capable of simulating entire trough campaigns, is presented. The open-source computing platform FEniCS is used to numerically solve the model using a finite element method. A manufactured solution test for the heat diffusion coupled with 2D nonlocal radiation is defined with the purpose of verifying the implementation, comparing the performance of different time discretization schemes and the adaptive time stepping algorithm. Concerning the BF tapping problem, the results show that during the time interval corresponding to a single tapping, the temperature in the radiation enclosure swiftly reaches the steady state value. Nevertheless, to obtain a steady state in the bulk of the solids, much longer time scales are needed due to the large thermal inertia of the structure. ©2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction The blast furnace (BF) is a high temperature, counter current flow reactor wherein the iron ore is smelted and reduced to produce hot metal. The molten hot metal and the slag, formed as a by-product of the process, accumulate in the hearth of the BF and are regularly drained as their level becomes sufficiently high. Their discharge takes place on a BF main trough or runner, a long receptacle lined with castable refractories designed to transport and separate slag and hot metal. Increasing BF trough availability and avoiding breakouts is a major concern for the metallurgical industry, as the trough refractory linings suffer a high degree of wear [1] . Usually, trough campaign life is capped when a certain tonnage of hot metal is produced or critical damage is detected, point at which it is relined. To improve the ironmaking process, it is of ∗Corresponding author at: Department of Applied Mathematics, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain. E-mail addresses: [email protected] (P. Barral), luisjavier[email protected] (L.J. Pérez-Pérez), peregrina.quint[email protected] (P. Quintela). https://doi.org/10.1016/j.apm.2021.12.029 0307-904X/© 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Nomenclature Greek symbols αAbsorptivity tTime step [s] εEmissivity ηDiscretized radiosity unknown vector [W/m 2 ] λParameter on Dirichlet temperature profile [K/s] Blocking factor Boundary of computational domain R Radiation cavity  R Intersection between the radiation cavity and the computational domain D R Slag upper surface Computational domain ωKernel of nonlocal radiation integral equation ρDensity [kg/m 3 ] σStefan-Boltzmann constant [W/(m 2 ·K 4 )] ςReflectivity τh Mesh of the computational domain ξDiscretized temperature unknown vector [K] Latin symbols c p Specific heat at constant pressure [J/(kg ·K)] c n Step size change ratio at n th time step ˆ c n Limited step size change ratio at n th time step e n Vector of local error estimate at n th time step G Irradiation [W/m 2 ] h Heat transfer coefficient [W/(m 2 ·K)] IIdentity operator k Thermal conductivity [W/(m ·K)] k in Thermal conductivity at the insulation lining [W/(m ·K)] KNonlocal radiation integral operator n Outward-pointing unit normal vector N h Dimension of V h q rad Heat flux due to radiation [W/m 2 ] Q Integral operator for radiative heat flux r n Norm of the local error estimate at the n th time step R Radiosity [W/m 2 ] tTime [s] t end Final time of the problem [s] t p Internal time node corresponding to the pth RK stage T Time interval of the problem T Temperature [K] T ext,c,n c Convection external temperature on n c C [K] T ext,r,n c Radiation external temperature on n c C [K] T ext,c,R Convection external temperature on R [K] T 0 Initial temperature [K] V h Vector space for spatial discretization x Spatial coordinate [m] X p Runge-Kutta method pth stage Subscripts C External boundaries D Dirichlet n c Index on convection boundaries R Radiation Superscripts a Academic problem G Restriction to the mesh of R 198 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 n n th time step n c Index on convection boundaries interest to assess and reduce the wear of refractory materials, which could allow operators to increase operation safety and trough productivity by reducing the frequency of repairs. Nonetheless, wear and degradation of the refractory linings is a complex phenomenon, with several factors involved. Some researchers argue that the main driving factor of trough refractory wear is mechanical erosion due to the fluid flow [2] . This is supported by the observation that the areas with higher turbulence in the liquids are those where increased wear rates are measured, both inside BF hearth [3] and in the main trough itself [4] . However, other authors identify chemical corrosion as the main cause responsible for refractory wear [5] , since the computed shear stress is low. The turbulence in the hot metal can enhance corrosion in the trough, which explains the non-uniform measured degradation patterns. Furthermore, the wear profiles are most intense at the locations where the hot metal-slag interface is located, which often dictates the trough campaign life [6] . In addition, as refractories are porous materials, the fluids infiltrate the lining. This may lead to severe spalling, as cycles of heating and cooling occur with the thermal expansion in the infiltrated refractory being different from that within the non-penetrated refractory [7] . Most of the available research has been focused on the different phenomena that occur during the various stages of the smelting process within the BF, such as the behaviour of the flow inside the BF hearth during its drainage [8,9] or the hearth refractory lining erosion [10] . In [11] , the state of the art concerning numerical simulation of processes inside the BF is reviewed. On the other hand, the research concerning phenomena related to the main trough is not as numerous. From the experimental standpoint, several works, such as [12–14] , have resorted to scale models to study the complex flow structures that are generated in the hot metal and slag flow in the BF trough, especially in the jet impingement region, where heavier turbulence is observed. Several authors have used CFD simulations to investigate these flow patterns, such as [15] and [16] , where the aim was to assess slag-metal separation and to reduce hot metal losses with slag. Moreover, in [5] , the structures generated by the jet impact in the hot metal pool in the trough were investigated. In [2] , the effect of parameters such as the taphole inclination angle on the wall shear were analysed by solving a similar model. In [17] , CFD was also used to study the main characteristics of the flow, comparing the numerical findings with experiments performed on a scale model. Regarding thermal modelling in the trough, [18] investigated the steady state thermal response of two different 2D cross-sections, accounting for radiative exchange with the slag surface. In the previous work [19] , the authors studied the 3D steady state of the temperature in the final part of the trough, assuming simplified flow conditions. In [20] , the hot metal flow in the trough was studied, including the conjugate heat transfer with the working lining. A thermomechanical model was solved to predict the thermal stress and the fatigue life of the trough. In the present work, we aim to investigate the evolution of the temperature in a cross-section of a BF trough. Specifically, we focus on assessing the temperature build-up during a BF tapping of 1.8 h of duration. To this end, it is essential to account for the radiative heat exchange due to the high temperatures that are reached. To study its effect, a nonlocal radiation condition is used in the radiation enclosure defined by the trough and the cover, which involves the solution of an integral equation coupled to the heat equation to model heat conduction in the solids. To the best knowledge of the authors, there are no previous works specifically dealing with the temperature evolution during a BF tapping and the analysis of the numerical methodology. To assess the correctness of the numerical algorithm, a manufactured solution test is performed, also known as method of manufactured solutions (MMS). The method constructs an analytical solution for the set of governing equations, determining suitable initial and boundary conditions [21,22] . When only heat diffusion is present, the derivation of the analytical solution is straightforward. However, for more complicated physics, such as conjugate heat transfer, the use of the MMS becomes more difficult [23] . In the case of coupled heat diffusion and nonlocal radiation, the main difficulty is dealing with the integral operator. For 3D geometries, if the domain is chosen as a spherical cavity, the derivation of a test greatly simplifies, as shown in [24] . On the other hand, for the 2D case, such derivation is substantially more intricate. In [25] , a procedure to obtain an analytical solution for the radiosity integral equation on a circular cavity was proposed. This procedure was also applied in [26] to obtain a test including steady heat diffusion within an annulus surrounding the cavity. Nonetheless, the proposed analytical test has the shortcoming of not being smooth. In this work, we extend this approach to provide a strategy that allows to derive the analytical solution to the radiosity corresponding to a given temperature by solving a Cauchy problem defined on the circumference, which leads to a smooth temperature field. This paper is structured as follows. In Section 2 , we describe the main features of the physical problem. In Section 3 , we state the mathematical model that we propose to find the transient temperature field in the trough cross-section, including the nonlocal radiation model. In Section 4 , the numerical approach that we follow to solve the model is described, presenting both the spatial and time discretization, as well as the adaptive time stepping algorithm, which allows to significantly reduce the computational cost. In Section 5 , a manufactured solution test is presented and solved, assessing the performance of different time discretization schemes, validating the numerical methodology. In Section 6 , the computed numerical results for the BF trough case are presented and discussed. Lastly, in Section 7 , the conclusions derived from this work are summarized. 199 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 1. Schematic of a typical BF main trough. Fig. 2. Main trough cross-section. 2. Physical problem The smelted hot metal and slag are drained from the BF hearth by drilling open an orifice in the BF wall, the taphole, a process known as tapping. The resulting molten stream exiting the BF is collected in the main trough, where separation by density difference takes place. Slag, being substantially lighter, stays on top and, after reaching the skimmer, is directed towards the slag runner, which is used to transport the slag to the granulator, as depicted in the empty trough schematic shown in Fig. 1 . Hot metal flows through a passageway underneath the skimmer and, after overflowing the iron dam, is transported to the tilting runner through the smaller iron runner. The end of the tapping is indicated by gas bursting out of the taphole, at which point it is plugged. Modern tapping practices involve main troughs designed so that a pool of fluids is kept between casts, which not only helps mitigate the kinetic energy of the impact of the falling jet stream [6] , but also leads to better fluid separation and preserves heat in the trough [27] . Even though the design of the BF main trough widely varies among casthouses, their general features are similar. In Fig. 2 (a), the cross-section of the design employed in this work is depicted, which is composed of three different material layers. The outermost one, known as working lining or wear lining, is lined with a castable refractory, usually with an alumina-based formulation containing silicon carbide and carbon [28] . As it is directly in contact with hot metal and slag, it sustains harsh operation conditions and thus suffers most of the wear during the BF trough campaign life. The safety lining serves as a back-up layer that should withstand an eventual working lining perforation. Lastly, the insulation lining is made with a thermal insulating material that maintains the temperature within suitable values. The trough is built on a supporting structure, formed by thin steel shells placed on a set of beams. Several lateral fins are incorporated to the structure, which enhance cooling with the surrounding air. A system of removable refractory covers is placed on the trough. Its main function is to suction the pollutant fumes produced during the BF tapping process. To this end, the air suction is performed next to the taphole, which causes a forced air stream to develop under the cover, flowing in the opposite direction to the slag and the hot metal. The trough is completely covered by this system, although in Fig. 1 it has been cut for clarity purposes. The covers block thermal radiation from the slag surface, allowing maintenance of heat within the fluids in the trough. Since both slag and the 200 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Table 1 Properties of the materials. k [W/ ( m ·K ) ] c p [ J / ( kg ·K )] ρ[ kg / m 3 ] Cover refractory 2.5 1296 2750 Insulation lining k in (T ) 1050 800 Safety lining 1.8 1172 2900 Working lining 2.6 1212 2500 refractory materials can be assumed to be opaque to radiation, the internal part of the runner and the cover form a radiation enclosure. 3. Mathematical model 3.1. Computational domain A 2D section located approximately halfway within the main trough is considered. This cross-section, situated before the slag runner and depicted in Fig. 2 (a), corresponds to x 3 = 8 m, where the coordinate x 3 is the distance to the BF. The region that fills with a pool of hot metal and slag during the tapping is excluded from the computational domain. Instead, the thermal behaviour of the fluids computed in [19] is imported, as it can be assumed that they quickly reach a temperature steady state. Therefore, the computational domain corresponds to the solid materials in the cross-section, as depicted in Fig. 2 (b). The steel casing shown in Fig. 1 is not included in the computational domain, as it is not compatible with a 2D thermal model due to the presence of lateral fins and supporting beams underneath it. The boundary D R , depicted in Fig. 2 (b), corresponds to the air-slag free surface. Even though the temperature in the fluids is not computed by the model, radiation emission from the slag plays a key role, drastically increasing the temperature in the exposed zones of the working lining and the cover. Its position slightly varies during the tapping, but once it reaches the level of the slag runner, it remains approximately at a steady horizontal position as slag begins to drain from the main trough. Since the trough fills swiftly during the first tapping due to the high discharge rate from the BF, we assume that the position of D R corresponds to a fixed x 2 coordinate, x 2 = 1 . 5 m, which coincides with the height at which the slag runner is located. 3.2. Model equations In , we solve the transient heat equation to model the evolution of the temperature: ρc p ∂T ∂t −div (k (T ) ∇T ) = 0 , (1) for t ∈ T = (t 0 , t end ] and x = (x 1 , x 2 ) ∈ . The symbols k , c p and ρdenote the thermal conductivity, specific heat and density, respectively. The insulation lining is composed of a good thermal insulator with a low thermal conductivity. Therefore, large temperature gradients are expected within it, ranging from values reaching up to 10 0 0 K to room temperature in the external part. To model its thermal conductivity, the following temperature-dependent function is considered: k in (T ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 . 13 , if T ≤293 , 1 980 0 0 (9 T + 10103) , if 293 <T <1273 , 0 . 22 , if T ≥1273 , which was obtained fitting the available experimental measurements at different temperatures. To model the remaining material properties, as we lack experimental data, the constant values shown in Table 1 are used, according to the data supplied by the steelmaking company collaborating in the research project PID2019-105615RB-I00. The working lining is an Al 2 O 3 -SiC-C based refractory castable, with a thickness of around 0.5 m. Even though factors such as refractory composition or porosity may strongly influence refractory wear, these are neglected since the study is mainly concerned with the thermal modelling. 3.3. Boundary conditions The boundaries of the computational domain are decomposed as displayed in Fig. 2 (b). The outer boundaries 1 C , 2 C and 3 C correspond to the top, the lateral walls and the bottom of the main trough, respectively. At these boundaries, we consider mixed convection-radiation boundary conditions: −k (T ) ∂T ∂n = h n c (T −T ext,c,n c )+σε n c (T 4 −T 4 ext,r,n c ) , (2) 201 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Table 2 Heat transfer coefficients. h [W/(m 2 ·K)] h 1 7 h 2 8 h 3 4 h R 35 Fig. 3. Computed temperatures in the fluids (from Barral et al. [19] ) and proposed fitted temperature profile T F (solid lines) at different locations. on n c C , n c = 1 , 2 , 3 , and for all t ∈ T . The values for the heat transfer coefficients, h n c , are gathered in Table 2 . In particular, the values corresponding to the outer boundaries are computed using standard correlations for vertical and horizontal plates assuming natural convection with the surrounding air extracted from Bergman et al. [29] , Cengel [30] . Since the steel casing is not included in the computational domain, as the fins and ribs built-in on it have a substantial impact in the cooling of the main trough, the heat transfer coefficients are modified to take into account their presence. Specifically, we compute a modified coefficient h n c = ˆ h n c A real /A , where ˆ h n c is the coefficient obtained using the correlation for natural convection assuming that the boundary is a flat plate, A real is the surface area of the boundary considering the presence of fins and ribs and A is the corresponding surface area for the flat plate. The values for the ambient temperature and the emissivity are considered equal in all cases ( T ext,r,n c = T ext,c,n c = 293 K and ε n c = 0 . 7 ). In addition, σ= 5 . 67 e −8 W/(m 2 ·K 4 ) is the StefanBoltzmann constant. The boundary D corresponds to the interface between the solids and fluids. On it, a Dirichlet condition using a known temperature is set: T (x , t) = T D (x , t) , (3) for all x ∈ D and t ∈ T . In (3) , the temperature T D is calculated as indicated below. At the beginning of the tapping, the temperature in the fluids rises quickly due to the hot metal and slag coming from the BF hearth at a very high temperature. Using the average hot metal tonnage produced in each single tapping, the estimated mean velocity of hot metal in the trough is around 0.1 m/s, so its residence time is approximately 3 min. As each tapping lasts for around 60–90 min, the flow develops rapidly and the temperature reaches a quasi-steady state in the fluids. In the previous work of the authors [19] , the steady state temperature of the main trough was computed using a 3D thermo-hydrodynamic model, assuming that the drained liquids were at a constant temperature, T = 1773 K, when they were discharged from the BF. At the vertical centreline crossing the fluids in the 2D domain (i.e. for x 3 = 8 m and x 1 = 1 . 5 m), the temperature depicted in Fig. 3 (a) was obtained. With h B , h hm and h slag , we denote the positions at which the bottom of the trough and the hot metal and slag surfaces are located, respectively. Only the closest part to the slag upper surface ( x 2 ∈ [1 . 35 , 1 . 5] ) is shown, as for lower heights the temperature is almost homogeneous. To approximate the computed temperatures in [19] , we propose the following fitting function, also displayed in Fig. 3 (a): T F (x 2 ) = −1 . 52 e −9 x 63 . 46 2 + 1773 , being h B ≤x 2 ≤h slag , where h B = 1 . 04 m and h slag = 1 . 5 m. To model the temperature rise in the fluids during the tappings, we assume that the approximated steady profile T F is reached swiftly. Moreover, we neglect any variations along the x 1 coordinate, which were found to be much smaller in magnitude when compared to variations in the x 2 coordinate, as depicted in Fig. 3 (b), where the temperature along 202 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 4. Temperature profile T D (x , t) with homogeneous T 0 = 293 K. horizontal lines crossing the fluids at different height values are plotted. Specifically, we consider that the temperature of slag and hot metal is equal to T D (x , t) = ⎧ ⎨ ⎩ T 0 (x ) −3 λ2 (t −t 0 ) 2 T 0 (x ) −T F (x 2 ) −2 λ3 (t −t 0 ) 3 (T 0 (x ) −T F (x 2 )) 2 , if t 0 < t ≤t steady (x ) , T F (x 2 ) , otherwise , (4) where t 0 and T 0 denote the initial time and temperature values. A constant value T 0 = 293 K is selected. In addition, the time required to reach the steady temperature is assumed equal to t steady (x ) = t 0 + T F (x 2 ) −T 0 (x ) λ, (5) where the value λ= 3 / 2 [K/s] is used in (4) and (5) , which ensures that the steady state in the fluids is reached in approximately 15 min. The profile (4) is constructed such that it is continuous for all t ∈ [ t 0 , t end ] , strictly increasing if t ∈ (t 0 , t steady (x )) for x ∈ D , and with zero time derivative at the endpoints of the latter interval. In Fig. 4 , the resulting profile is depicted for all x 2 , assuming a homogeneous initial temperature. The dependence on the initial temperature allows the model to be readily used for subsequent tappings after the first one, since the temperature within the retained pool of fluids in the trough may deviate substantially from a homogeneous value. It is noteworthy that the temperature profile (4) implicitly assumes that slag is completely molten, as both slag and hot metal were considered in [19] fully liquid. Nevertheless, partial solidification of the slag surface can occur, which would result in a higher temperature gradient than that shown in Fig. 3 (a). Lastly, on the intersection of the radiation cavity with the computational domain,  R , the thermal problem is coupled with a nonlocal radiation term to account for the heat flux due to radiation, q rad . Since this radiation model involves an integral equation defined on the cavity boundary, it is described in detail in Section 3.4 . Moreover, we also consider a convection contribution to the heat flux on  R , which becomes: −k ∂T ∂n = h R (T −T ext,c,R ) + q rad , (6) on  R for all t ∈ T . The radiative heat flux, q rad , is obtained using the radiation model, whereas the temperature of the surroundings is T ext,c,R = 293 K. The heat transfer coefficient, h R , takes the value shown in Table 2 , and is estimated based on the air flow prediction, and the corresponding computed heat flux, using similar 3D numerical results that include the forced air circulation below the cover to those reported in [19] . 3.4. Nonlocal radiation model To model the radiation heat exchange in the cavity formed by the refractory cover, the slag and the part of the working lining not covered by fluids, a nonlocal term is applied on its boundary, analogous to that used in [31] for an axisymmetric cavity formed by black surfaces in a problem related to silicon purification. Also, a similar problem was studied in [24] , considering an axisymmetric grey enclosure to model radiation in the internal cavity of a workpiece within an induction furnace. In our setting, the same model is employed to obtain the radiation heat flux on the nonlocal radiation cavity, with the difference that our cavity is 2D. 3.4.1. Radiosity integral equation The complete radiation cavity is denoted as R , as shown in Fig. 2 (b) (see the dashed line). We assume that the air inside this cavity is transparent to radiation, while the solid refractories, slag and hot metal are assumed to be opaque to radiation. 203 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 It should be noted that the fumes produced during the BF tapping process may not be transparent to radiation. However, since the air suction is performed next to the taphole, the forced air circulation below the cover flows countercurrent to slag and hot metal, with air being suctioned mainly from the opposite side of the refractory cover. Thus, the air composition at the computational domain is expected to be similar to that of the casthouse environment air. Moreover, we assume that R behaves as a grey and diffuse surface. Hence, emitted or reflected radiation does not depend on wavelength or direction, and the radiation heat flux on R satisfies q rad (x ) = R (x ) −G (x ) , (7) for all x ∈ R and where R and G are functions defined on R , known as radiosity and irradiation, respectively. The first one corresponds to the outgoing radiative energy whereas the latter is the incoming radiative energy at each point. Note that in (7) and in the remainder of this section, the dependence on time of all fields is omitted to simplify the writing. Radiosity is the sum of emitted radiation and reflected irradiation. The emitted part is related to the temperature at the boundary by the Stefan-Boltzmann law. Thereby, on R : R (x ) = ε R (x ) σT 4 (x ) + ς R (x ) G (x ) , (8) where ς R and ε R denote the reflectivity and emissivity on R , respectively. We assume that there exists a constant ε 0 such that 0 < ε 0 ≤ε R (x ) ≤1 . (9) Provided that the reflectivity represents the fraction of incoming radiative energy that is reflected, we can instead rewrite (8) in terms of the absorptivity αR = 1 −ς R . At thermal equilibrium, Kirchhoff’s law for thermal radiation implies that the absorptivity satisfies ε R (x ) = αR (x ) . Then, (8) becomes R (x ) = ε R (x ) σT 4 (x ) + (1 −ε R (x )) G (x ) , on R . (10) For any x ∈ R , the irradiation is related to the radiosity on the remainder of the radiation boundary R \ { x } as follows: G (x ) =  R ω(x , y ) R (y ) ds y , (11) where ωis the integral kernel, sometimes known as differential view factor or simply view factor [32] . When R is the boundary of a two-dimensional domain, ωis given by the following expression: ω(x , y ) = n (x ) ·(y −x ) n (y ) ·(x −y ) 2 | x −y | 3 ( x ,y ) , with x  = y , (12) with n denoting the outward-pointing unit normal vector from . The factor accounts for parts of the cavity that may block the view between points and is defined as (x , y ) = 0 , if xy ∩  = ∅ , 1 , if xy ∩ = ∅ , (13) where xy is the segment connecting the points x and y . It can be shown that the view factor satisfies [33]  R ω(x , y ) ds y ≤1 , (14) for any x ∈ R \ S, Sbeing the set of non-smooth points of R . If the equality is attained in (14) , then is an enclosure , which intuitively means that the radiation from the cavity formed by R cannot escape to the surroundings. Let Kdenote the following linear integral operator: K(R )(x ) =  R ω(x , y ) R (y ) ds y , (15) for x ∈ R . Substituting (11) in (10) and using (15) , the following integral equation for the radiosity R is obtained: (I −(1 −ε R (x )) K)(R )(x ) = ε R (x ) σT 4 (x ) , on R . (16) Under sufficient regularity of R and if the emissivity satisfies (9) , it can be shown that the previous equation has a unique solution [34] . To evaluate the radiative heat flux on the radiation cavity, the operator Q is introduced: Q (R )(x ) = (I −K)(R )(x ) , (17) which by virtue of expression (7) and the definition of the integral operator K, introduced in (15) , directly relates the radiosity to the radiative heat flux on R : q rad (x ) = Q (R )(x ) , on R , where R is computed by solving (16) . The coupling with the thermal problem is done through the radiative heat flux, q rad , on the radiation cavity as a boundary condition. With this purpose, the convection contribution is added to q rad to obtain the boundary condition (6) . Note that the 204 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 boundary of the domain only includes the subset  R of the radiation boundary R (see Fig. 2 (b)). Despite the temperature being known at the slag surface, the corresponding radiosity on its surface is unknown, as it depends also on the radiosity on the rest of the cavity. Thus, (16) must be solved on R assuming that T | D R = T D , i.e., that the slag upper surface temperature is equal to the value discussed in Section 3.3 . In more general cases, when the medium is emitting, absorbing or scattering radiation, the full integro-differential radiative heat transfer equation must be solved instead of the integral Eq. (16) [32] , which is much more computationally costly and requires specific discretization methods. 3.5. Strong form Considering (1) subjected to the boundary conditions (2), (3) and (6) , as well as a homogeneous initial temperature T 0 = 293 K, we obtain the following strong formulation of the problem. Problem P. Find T (x , t) and R (x , t) such that: ρc p ∂T ∂t −div (k (T ) ∇T ) = 0 , in ×T , T (t 0 ) = T 0 in , T = T D , on D ×T , (18) −k (T ) ∂T ∂n = h n c (T −T ext,c, n c ) + σε n c (T 4 −T 4 ext,r,n c ) , on n c C ×T , (19) −k ∂T ∂n = h R (T −T ext,c,R )+q rad , on  R ×T , q rad = Q (R ) , on R ×T , (I −(1 −ε R ) K)(R ) = ε R σT 4 , on R ×T , where n c = 1 , 2 , 3 . 3.6. Variational formulation Let L 2 () denote the Lebesgue space of square-integrable functions in . Moreover, let H 1 () be the Sobolev space of functions in L 2 () such that their distributional derivative is also in L 2 () and let H 1 0 , D () be the subspace of H 1 () formed by those elements with zero trace on D . Hereafter, we assume that the data are smooth enough for the following operations to make sense. For each t ∈ T , we multiply (18) by a test function v = v (x ) . Then, we integrate in the corresponding domain, apply a Green’s formula, and use the boundary conditions, leading to the following variational formulation. Problem VP . For each t ∈ T , find T (t) ∈ H 1 () and R (t) ∈ L 2 (R ) such that ∂T (t) ∂t ∈ L 2 () , T (t) = T D (t) on D , and satisfying:   ρc p ∂T (t) ∂t v dA +   k (T (t)) ∇T (t) ·∇v dA +   R h R T (t) v ds + 3  n c =1  n c C h n c T (t) + σε n c T 4 (t) v ds (20) = 3  n c =1  n c C h n c T ext,c,n c + σε n c T 4 ext,r,n c v ds +   R (h R T ext,c,R −Q (R (t))) v ds, for all v ∈ H 1 0 , D () , and verifying: (I −(1 −ε R ) K)(R (t)) = ε R σT 4 (t) , on R , (21) for t ∈ T with T (0) = T 0 . 4. Numerical solution The model is solved numerically using the open-source computing platform FEniCS [35] . Its components automatically generate efficient low level code to assemble the corresponding linear system from a symbolic representation of the weak formulation in the UFL language [36] . A triangulation τh of the domain is constructed such that = ∪ K∈ τh Kand conforms to the different material layers in the BF trough. A method of lines is used to discretize Problem VP : first, the spatial 205 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 7. Numerical approximation of the temperature for Problem P a at t = 5 s, implicit Euler scheme, t/ 4 and Mesh 5. Fig. 8. Approximated radiosity for Problem P a , with implicit Euler and different step size values at t = 5 s with Mesh 5. 5.2. Fixed step size Instead of the time discretization proposed in Section 4.2 using an adaptive time stepping technique, firstly we consider a fixed step size to solve Problem P a . To this end, we compare the performance of an implicit Euler scheme and the ESDIRK 3/2a RK. In comparison, the former is much easier to implement and is fairly inexpensive numerically, which motivates its use in various industrial problems on different fields [44,45] and also on similar applications [31] due to its unconditional stability. Furthermore, the formal order in time of the Euler scheme is one, which matches the expected order of convergence of the radiosity approximation in space. In both cases, the reference value t = 0 . 5 s is set, which is successively divided by a factor of two to analyse the convergence in the considered error norms as the mesh and the step size are refined. Concerning the use of the implicit Euler scheme, in Fig. 7 , the contours of temperature are shown at time t = 5 s. Moreover, in Fig. 8 , the numerical approximations for the radiosity are compared with the exact solution for the different step sizes, also at time t = 5 s. Except for the smallest ones, the errors are very large. In Table 4 , the computed e t h,L 2 (T h ) relative error is gathered. Moreover, in Table 5 , the e t h,L 2 (R h ) relative error approximating the radiosity is displayed. We observe the expected formal asymptotic behaviour as both the time step and mesh size are refined, i.e., O (h 2 + t) convergence in the L 2 (0 , t a end ;L 2 (a )) norm for the temperature and O (h + t) convergence for the radiosity in the L 2 (0 , t a end ;L 2 (a R )) norm. It is also worth noticing that the relative error due to the time discretization is too large to see the optimal convergence rates as the mesh size is refined, given a constant t. For instance, in Table 4 , considering t/ 32 , we observe that the error is reducing with rate O (h 2 ) , which is the expected rate for P 1 finite elements in the L 2 (a ) norm. The behaviour continues up until reaching Mesh 3, where the convergence in space stagnates. This is also the case in Table 5 . When the mesh is sufficiently fine, such as Mesh 2 in Table 4 , successively considering refined step sizes yields an order one convergence rate of the error. 212 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Table 4 e t h,L 2 (T h ) relative error, implicit Euler scheme. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 t5.97e—02 6 . 44 e −02 6 . 53 e −02 6 . 56 e −02 6 . 56 e −02 6 . 57 e −02 t/ 2 2 . 96 e −02 3 . 34 e −02 3 . 41 e −02 3 . 43 e −02 3 . 44 e −02 3 . 44 e −02 t/ 4 1 . 39 e −02 1.67e—02 1 . 74 e −02 1 . 76 e −02 1 . 76 e −02 1 . 76 e −02 t/ 8 6 . 91 e −03 8 . 15 e −03 8 . 70 e −03 8 . 86 e −03 8 . 90 e −03 8 . 91 e −03 t/ 16 5 . 27 e −03 3 . 86 e −03 4.29e—03 4 . 43 e −03 4 . 47 e −03 4 . 48 e −03 t/ 32 5 . 63 e −03 1 . 88 e −03 2 . 07 e −03 2 . 20 e −03 2 . 24 e −03 2 . 24 e −03 Table 5 e t h,L 2 (R h ) relative error, implicit Euler scheme. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 t5.32e—01 5 . 51 e −01 5 . 55 e −01 5 . 56 e −01 5 . 57 e −01 5 . 57 e −01 t/ 2 3 . 24 e −01 3.45e—01 3 . 50 e −01 3 . 51 e −01 3 . 51 e −01 3 . 51 e −01 t/ 4 1 . 73 e −01 1 . 93 e −01 1.97e—01 1 . 98 e −01 1 . 99 e −01 1 . 99 e −01 t/ 8 8 . 63 e −02 1 . 01 e −01 1 . 05 e −01 1.06e—01 1 . 06 e −01 1 . 06 e −01 t/ 16 4 . 97 e −02 5 . 26 e −02 5 . 40 e −02 5 . 45 e −02 5.47e—02 5 . 47 e −02 t/ 32 4 . 42 e −02 3 . 01 e −02 2 . 81 e −02 2 . 79 e −02 2 . 78 e −02 2.78e—02 Table 6 e t h,L 2 (T h ) relative error, ESDIRK 3/2a. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 t6 . 54 e −03 2 . 07 e −03 1 . 77 e −03 1 . 79 e −03 1 . 80 e −03 1 . 80 e −03 t/ 2 6 . 66 e −03 1 . 43 e −03 4 . 57 e −04 3 . 30 e −04 3 . 28 e −04 3 . 30 e −04 t/ 4 6 . 70 e −03 1 . 44 e −03 3 . 65 e −04 1 . 01 e −04 5 . 74 e −05 5 . 51 e −05 t/ 8 6 . 71 e −03 1 . 44 e −03 3 . 66 e −04 9 . 13 e −05 2 . 33 e −05 9 . 78 e −06 t/ 16 6 . 71 e −03 1 . 44 e −03 3 . 67 e −04 9 . 17 e −05 2 . 25 e −05 5 . 66 e −06 t/ 32 6 . 71 e −03 1 . 44 e −03 3 . 67 e −04 9 . 18 e −05 2 . 26 e −05 5 . 63 e −06 Fig. 9. Approximated radiosity for Problem P a , with ESDIRK 3/2a and different meshes at t = 5 s, t/ 32 . Instead, if the ESDIRK 3/2a RK method is used to solve the manufactured solution test P a with the same reference step size, the accuracy of the computed solution improves significantly. In Fig. 9 , the numerical approximation for the radiosity with this method is displayed, for the different mesh sizes and t/ 32 . From Mesh 2 to Mesh 6, the approximation essentially overlaps the exact solution, in contrast with the results shown in Fig. 8 for the implicit Euler scheme. In Table 6 , the e t h,L 2 (T h ) relative error is shown whereas the same error for the radiosity is displayed in Table 7 . For the coarsest meshes, reducing the time step size yields no error reduction, as the space discretization error dominates. On the other hand, for the finest meshes, a fast decay of the error is observed. Comparing the results in Table 4 with those in Table 6 , we observe that the RK scheme yields a significant error reduction on the temperature approximation, which reaches up to three orders of magnitude for the t/ 32 step size and Mesh 6. 213 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Table 7 e t h,L 2 (R h ) relative error, ESDIRK 3/2a. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 t4 . 82 e −02 2 . 39 e −02 1 . 79 e −02 1 . 63 e −02 1 . 59 e −02 1 . 58 e −02 t/ 2 5 . 29 e −02 2 . 18 e −02 1 . 08 e −02 5 . 79 e −03 3 . 67 e −03 2 . 91 e −03 t/ 4 5 . 38 e −02 2 . 22 e −02 1 . 08 e −02 5 . 33 e −03 2 . 68 e −03 1 . 38 e −03 t/ 8 5 . 40 e −02 2 . 23 e −02 1 . 08 e −02 5 . 34 e −03 2 . 66 e −03 1 . 33 e −03 t/ 16 5 . 40 e −02 2 . 23 e −02 1 . 08 e −02 5 . 34 e −03 2 . 67 e −03 1 . 33 e −03 t/ 32 5 . 40 e −02 2 . 23 e −02 1 . 08 e −02 5 . 34 e −03 2 . 67 e −03 1 . 33 e −03 Fig. 10. Relative errors in log-log scale. Table 8 Maximum step size, number of accepted steps and rejected steps and acceptance rate using the H211b controller. Max. step size [s] Accepted steps Rejected steps Acceptance rate tol 3 . 32 e −01 44 14 7 . 59 e −01 tol/ 2 3 . 17 e −01 57 18 7 . 60 e −01 tol/ 4 2 . 38 e −01 72 16 8 . 18 e −01 tol/ 8 1 . 98 e −01 91 13 8 . 75 e −01 tol/ 16 1 . 62 e −01 117 13 9 . 00 e −01 tol/ 32 1 . 30 e −01 150 13 9 . 20 e −01 In Fig. 10 (a), a log-log scale plot of the e t h,L 2 (T h ) relative error considering Mesh 6 is shown. In the case of using the implicit Euler scheme, the expected theoretical rate O (t) is clearly observed whereas the ESDIRK 3/2a results in a convergence rate between O (t 3 ) and O (t 2 ) . Specifically, the estimated convergence rate using the RK scheme, computed by linear regression of the first 4 points, is 2.57. For the smallest step sizes, convergence rapidly deteriorates and stagnates around 5 e −6 . This behaviour owes to the space discretization error becoming dominant, as the O (h 2 ) rate in space does not allow it to fall below such value for the considered mesh sizes, as shown by the column corresponding to Mesh 6 in Table 6 , or in Fig. 10 (b), where the expected convergence rates in space are observed. There, the e t h,H 1 (T h ) is computed similarly as (40) . Hence, the error can only be further reduced by considering a finer mesh. Since the objective is to simulate long time intervals, adjusting the size of the time step as much as possible but keeping the error on the radiosity controlled, it can be concluded that the implicit Euler method is not the best option in this type of application. 5.3. Adaptive step size As described in Section 4.3 , when using the ESDIRK 3/2a scheme instead of considering a constant step size, the local error estimator provided by the method can be used to adapt the time step size. Hence, we choose a fixed tolerance value tol = 1 e −3 and consider atol = rtol = tol as the parameters in the error norm (32) . This tolerance value is successively reduced to solve Problem P a . In Fig. 11 , the computed step size values are displayed with respect to time for Mesh 6. The rejected steps are represented as squares. For all tolerance values, it is observed that the first step, supplied to the H211b controller as the reference value t = 1 s, is always rejected and adapted to a smaller step size. In Table 8 , the maximum step size, the number of accepted and rejected steps as well as the acceptance rate computed by the controller are displayed. The number of steps increases with decreasing tolerance. Conversely, the number of rejected 214 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 11. Computed step size values for the different tolerance values in Mesh 6. The coloured square symbols ( ) represent the rejected steps for each tolerance. Table 9 e t h,L 2 (T h ) relative error, ESDIRK 3/2a with H211b controller for adaptive step size. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 tol 6 . 71 e −03 1 . 46 e −03 4 . 30 e −04 2 . 50 e −04 2 . 36 e −04 2 . 36 e −04 tol/ 2 6 . 70 e −03 1 . 44 e −03 3 . 80 e −04 1 . 45 e −04 1 . 18 e −04 1 . 18 e −04 tol/ 4 6 . 70 e −03 1 . 44 e −03 3 . 68 e −04 1 . 07 e −04 6 . 33 e −05 6 . 04 e −05 tol/ 8 6 . 69 e −03 1 . 44 e −03 3 . 66 e −04 9 . 58 e −05 3 . 93 e −05 3 . 37 e −05 tol/ 16 6 . 69 e −03 1 . 44 e −03 3 . 65 e −04 9 . 18 e −05 2 . 71 e −05 1 . 77 e −05 tol/ 32 6 . 69 e −03 1 . 44 e −03 3 . 65 e −04 9 . 11 e −05 2 . 33 e −05 9 . 84 e −06 Table 10 e t h,L 2 (R h ) relative error, ESDIRK 3/2a with H211b controller for adaptive step size. Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 Mesh 6 tol 5 . 33 e −02 2 . 19 e −02 1 . 07 e −02 5 . 47 e −03 3 . 07 e −03 2 . 08 e −03 tol/ 2 5 . 36 e −02 2 . 21 e −02 1 . 07 e −02 5 . 35 e −03 2 . 77 e −03 1 . 56 e −03 tol/ 4 5 . 38 e −02 2 . 22 e −02 1 . 07 e −02 5 . 33 e −03 2 . 69 e −03 1 . 40 e −03 tol/ 8 5 . 39 e −02 2 . 22 e −02 1 . 08 e −02 5 . 33 e −03 2 . 67 e −03 1 . 35 e −03 tol/ 16 5 . 40 e −02 2 . 22 e −02 1 . 08 e −02 5 . 33 e −03 2 . 66 e −03 1 . 34 e −03 tol/ 32 5 . 39 e −02 2 . 23 e −02 1 . 08 e −02 5 . 33 e −03 2 . 66 e −03 1 . 33 e −03 steps does not show a clear increasing trend. Hence, the acceptance rate is approximately constant for the larger tolerance values and starts to increase from tol/ 4 onwards. In Table 9 , the e t h,L 2 (T h ) relative errors computed for the different tolerance values and meshes are shown. For the coarser meshes, the committed error in the space discretization dominates, and O (h 2 ) order is observed as the mesh is refined for a constant time step, up until Mesh 5. The convergence rate is slightly below order 1 with respect to the tolerance, O (tol) , as evidenced in Fig. 12 , where the e t h,L 2 (T h ) relative error corresponding to Mesh 6 is plotted for all the tolerance values. The slope of the error is approximately equal to 0.94. As in the previous discussion when using the ESDIRK 3/2a scheme without adapting the step size, the space discretization error is dominant, and the theoretical convergence rate is observed as the mesh is refined. Similar remarks apply to e t h,L 2 (R h ) , displayed in Table 10 for the different tolerance values. The computed results for the proposed manufactured solution test show that the numerical algorithm performs within expectations. When the mesh size and the step size are both suitable, the theoretical convergence rates are reached. Nonetheless, it is observed that for a smooth problem such as the presented one, with small radiative heat fluxes, the computed errors for the implicit Euler scheme are very large, even for the smallest step sizes. The ESDIRK 3/2a scheme enables to greatly reduce the time discretization error. However, this comes at the cost of increasing the computational time by a factor of three approximately, as three stage solutions have to be computed at each time step. Since the method has an embedded computationally inexpensive error estimator, the step size may be readily adapted using an error controller. In such case, for the H211b controller, we obtain comparable errors if the same number of steps are used, meaning that no significant gain is achieved. Nevertheless, for the problem in the main trough, the advantage of using an adaptive step size is more evident, since for a fixed tolerance value, the difference between the minimum and maximum step sizes computed by the controller becomes much larger. In the proposed manufactured solution test, as the temperature variations are relatively smooth and periodic in time, the maximum step sizes computed by the H211b are at most twice the smallest step size for each tolerance value. 215 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 12. e t h,L 2 (T h ) with Mesh 6 for various tolerance values, using ESDIRK 3/2a with H211b controller. Table 11 Meshes used for the mesh sensitivity test, maximum element size, minimum element size, number of triangle elements in the grid of the domain and number of edge elements on the nonlocal radiation boundary. h max h min # τh # τG h Mesh 1 1 . 14 e −1 5 . 37 e −2 3802 100 Mesh 2 6 . 34 e −2 2 . 69 e −2 15208 200 Mesh 3 3 . 17 e −2 1 . 34 e −2 60832 400 Mesh 4 1 . 58 e −2 6 . 73 e −3 243328 800 Fig. 13. BF trough cross-section Mesh 2. 6. Numerical results In this section, the numerical results obtained solving Problem Pare presented. The discretization provided in Section 4 and Algorithm 1 are used. The H211b controller is employed to adapt the time step size according to the error estimator embedded in the ESDIRK 3/2a RK scheme. To ensure that the considered mesh of the computational domain is suitable, in Section 6.1 , we compare the obtained solution for different mesh sizes. In addition, the sensitivity of the solution with respect to the tolerance parameter selected in the the H211b controller is analysed. Subsequently, in Section 6.2 , the obtained results for the selected tolerance and mesh size are discussed. 6.1. Sensitivity test To obtain the solution of Problem P, we consider the meshes gathered in Table 11 , where the total number of triangles in the mesh and the number of edge elements in the nonlocal radiation boundary R mesh are presented, as well as the minimum element size h min and the maximum element size h max . Mesh 2 is depicted in Fig. 13 . 216 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 14. Computed numerical solution for the considered mesh sizes. Fig. 15. Step size values. Table 12 Maximum step size, number of accepted steps, number of rejected steps and acceptance rate using the H211b controller with tol = 1 e −4 for the different mesh sizes used in the grid sensitivity test. Max. step size [s] Accepted steps Rejected steps Acceptance rate Mesh 1 5 . 39 e + 02 39 1 9 . 75 e −01 Mesh 2 7 . 04 e + 02 38 0 1 . 00 e + 00 Mesh 3 7 . 11 e + 02 37 1 9 . 74 e −01 Mesh 4 6 . 75 e + 02 36 1 9 . 73 e −01 A fixed tolerance value for the step size controller t ol = at ol = rt ol = 1 e −4 is considered, with a starting step size t = 1 s. Moreover, in the remainder of this section, we set μ=0 . 8 and tol fp =1 e −4 in Algorithm 1 . In Fig. 14 (a), the radiosity on R is displayed as afunctionof arc length on the boundary at t = t end = 1 . 8 h, for all mesh sizes. The starting point of the curve corresponds to the left corner of the slag surface ( x = (1 . 041 , 1 . 5) m) with R parametrized anticlockwise. Meshes 2, 3 and 4 produce almost identical results, whereas the computed radiosity with Mesh 1 shows small differences in comparison with the remaining meshes. In Fig. 14 (b), the computed temperature at point x = (1 . 5 , 2 . 5) m, located in the inner part of the middle of the cover, is displayed over time. For Mesh 1, the computed temperature is above the remaining mesh sizes, whereas Mesh 2 is slightly below the results for Meshes 3 and 4. In Fig. 15 (a), the accepted step size values computed by the controller are plotted for the considered grids. The rejected steps are also shown, depicted as squares. For all mesh sizes, the same trend is observed, as the step size rapidly grows upon reaching around 100 s. At this point, the temperature on the slag surface, D R , has increased enough for the radiative heat fluxes to become relevant. Thus, the step size is adjusted accordingly by decreasing its size. Then, once the radiosity on the nonlocal radiation boundary stabilizes, the step sizes swiftly rise back, increasing until reaching the end time t end . Provided that the same tolerance value is selected for all meshes, the differences among the computed step sizes are small. In Table 12 , the maximum step size, the number of accepted and rejected steps as well as the acceptance rate of steps are gathered, showing similar results for the considered meshes. 217 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Table 13 Maximum step size, number of accepted steps, number of rejected steps and acceptance rate using the H211b controller for various tolerance values. tol Max. step size [s] Accepted steps Rejected steps Acceptance rate 1 e −3 9 . 80 e + 02 22 1 9 . 56 e −01 1 e −4 7 . 12 e + 02 37 1 9 . 74 e −01 1 e −5 3 . 82 e + 02 69 0 1 . 00 e + 00 1 e −6 1 . 97 e + 02 140 0 1 . 00 e + 00 Fig. 16. Computed fixed point algorithm iterations for the second RK stage at each time step for different tolerance values. In view of the presented results, Mesh 3 is considered fine enough to obtain an accurate approximation of the solution for the proposed problem. Hereinafter, in all the displayed results, Mesh 3 is used to solve the problem. For this mesh size, if different tolerance values for the controller are selected, a similar trend for the time step size values is observed, as displayed in Fig. 15 (b). For tolerances tol = 1 e −5 and tol = 1 e −6 , the parameter tol fp in Algorithm 1 is adjusted so that the tolerance of the fixed point algorithm is as restrictive as that of the step size controller, i.e., t ol fp = t ol. In Table 13 , the maximum step size, the number of accepted and rejected steps as well as the acceptance rate computed by the controller are gathered for the selected tolerance values. The maximum step sizes are two orders of magnitude above the initial step size. For the selected tolerance values, both the computed temperature and radiosity are very similar. However, the choice of the tolerance significantly affects the performance of the fixed point algorithm. To illustrate this dependence, in Fig. 16 , the number of iterations computed to solve the second RK stage at each step is represented. The trends are similar, showing that the lower the tolerance, the lower the number of iterations required to converge, as the computed time step sizes are smaller. Furthermore, as the time step size increases, convergence of the algorithm can be compromised, as evidenced by the reported iteration count for tol = 1 e −3 towards the end time of the problem, which rises quickly and eventually shows an unconverged stage solution near t = 1 . 4 h, corresponding to the rejected step depicted in Fig. 15 (b) for that tolerance. This behaviour can also appear in the rest of the considered tolerances if larger t end are set, as the computed step sizes may also become too large for the fixed point algorithm to converge. In view of the presented results, tol = 1 e −4 is considered enough to obtain a sufficiently accurate approximation of the solution for the proposed problem. 6.2. Simulation of the BF tapping Using Mesh 3 and a tolerance parameter tol = 1 e −4 such that at ol = rt ol = t ol in the norm of the local error estimator (32) , Problem Pis solved, considering a time interval T = (t 0 , t end ] = (0 , 1 . 8] h. For these values and a starting step size t 0 = 1 s, as it was shown in Table 12 , only 37 time steps are required to solve the problem. However, as displayed in Fig. 15 (a), at the start of the tapping, small steps are needed to maintain the error estimator below the specified bounds, evidencing the convenience of using an adaptive step size controller for this problem. In Fig. 17 (a) and (b) the temperature contours at t = 1562 . 1 s and t = t end are depicted, respectively. High temperature values are observed only in the closest parts to the fluids, which heat rapidly following the profile determined by (4) . Most of the trough is unaffected due to its large thermal inertia, still remaining at the initial temperature T 0 = 293 K. However, those parts of the working lining and the refractory cover that are exposed to radiation emission from the slag surface heat quickly. In particular, temperature reaches up to 110 0 K in the middle of the cover. To better showcase the evolution of the temperature on the radiation enclosure, in Fig. 18 (a), the temperature on R is shown for different time values, as a function of arc length. The hottest part of the cavity corresponds to the slag surface, where the temperature rises swiftly. Therefore, radiation emission and radiosity increase, which leads to a quick temperature build-up in the remainder of R . In Fig. 18 (b), 218 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 17. Temperature contours at various time values. Fig. 18. Temperature and radiosity on R at various time values. the radiosity on R is displayed for the same time values. As in Fig. 14 (a), the highest radiosity value is observed at the slag surface, which drives heat transfer towards the rest of the radiation cavity. In contrast, if the end time value was selected as t end = 15 days, a steady state of the temperature in the cross section of the trough is almost reached. In Fig. 19 (a), the temperature contours in such situation are depicted. Fig. 19 (b) shows the critical isotherms (1423 K, 1473 K and 1523 K). These isotherms indicate an increase in the degree of wear due to the onset of chemical composition changes or phase transitions in the fluids or the solid refractories. In particular, the 1423 K isotherm corresponds approximately to the solidification of the hot metal. Keeping it inside the working lining ensures that any hot metal leaked in the safety lining solidifies, forming a protective layer and reducing the possibility of a major trough breakout. The 1423 K isotherm lies close to the safety lining in the bottom of the trough after 15 days of continuous cast, indicating risk of hot metal infiltration within the safety lining. Also, as already was suggested by Fig. 18 (a) and (b), the temperature and radiosity on R stabilize much quicker than the temperature in the remaining of the trough. This is clearly shown in Fig. 20 , where the evolution in time of the temperature computed at several points in the domain is depicted. Specifically, the temperature at points p 1 = (1 . 5 , 2 . 55) , p 2 = (1 . 5 , 0 . 54) and p 3 = (1 . 5 , 0 . 18) m is shown. These points are located along the vertical symmetry line of the domain at the middle of the cover, the boundary among the working and safety linings, and the boundary between the safety and insulation linings, respectively (see Fig. 2 (a)). The temperature at point p 4 = (0 . 65 , 2) m is also displayed, which has the particularity of not being directly exposed to incoming radiation from the slag surface. This results in a slower temperature increase when compared to the parts that are exposed, such as p 1 , where the temperature computed at t = 1 . 8 h was already very close to the steady state temperature. However, at t end = 15 days, the temperature on the insulation and safety linings is still not completely stable. In particular, a delay of almost a day for the temperature build-up to start is observed in the insulation lining, as evidenced by the depicted evolution of the temperature at p 3 . 219 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 Fig. 19. Continuous cast with t end = 15 days. Fig. 20. Computed temperature over time choosing t end = 15 days at selected points in . We remark that the proposed setting with t end = 15 days is far from realistic, since the BF draining process is not continuous, as described in Section 2 . The effect of the process stops to control the liquid level in the BF hearth may have a significant effect in the BF trough temperatures. Also, even though the trough that we used for the geometrical model has installed thermocouples to obtain temperature measurements, the devices are embedded in the insulating layer to protect them from any hot metal infiltration, which means that the large thermal inertia of the structure prevents the devices from recording any temperature change during the time period corresponding to a single BF tapping. Further work should be aimed to obtain an experimental validation of the study and to incorporate the effect that the BF process stops have in the temperature field. 7. Conclusion A coupled mathematical model of the transient thermal heating of a 2D cross-section of a BF main trough during the BF tapping was proposed, accounting for the different layers of refractory materials and the thermal radiation using a nonlocal model. The ESDIRK 3/2a RK method was used altogether with the H211b step size controller to modify the step size, allowing to keep the time discretization error below specified bounds. To couple the nonlocal radiation and the heat equation in the solids, a fixed point algorithm was proposed. The numerical algorithm was validated with a manufactured solution test, specifically designed to address the coupled heat conduction and radiation. To this end, the methodology devised in [25] was extended to derive a methodology to obtain the corresponding radiosity to a given temperature by the solution of a Cauchy problem. A benchmark case was proposed and its analytical solution was given. For this case, the performance of the ESDIRK 3/2a was compared to that of an implicit Euler method, considering both an adaptive step size and a constant step size. The proposed algorithm was found to perform well in adequately approximating the radiosity, which is a critical phenomenon in the temperature build-up, and also retains a good approximation of the temperature over long time scales. The algorithm was successfully applied to the BF cast problem setting, which features large thermal gradients and long time scales. The obtained results show that during a standard tapping of 1.8 h, both the radiosity and temperature swiftly reach their steady state values on the radiation enclosure, as the slag surface is the main heat driving factor towards the 220 P. Barral, L.J. Pérez-Pérez and P. Quintela Applied Mathematical Modelling 105 (2022) 197–225 exposed areas in the cavity. In those parts of the enclosure that are not directly exposed to thermal radiation from the slag, the temperature build-up is substantially slower. On the other hand, if a continuous cast of 15 days is assumed, the bulk of the solids has not fully reached a steady state due to the large thermal inertia of the structure. In particular, the insulation lining remains at the initial temperature for almost one day. After 15 days, the critical isotherms are found to be close to the safety lining, which indicate high risk of a hot metal leak into this back-up lining. The numerical findings suggest that the BF tapping process stops may have a significant effect in the temperature field in the refractories, preventing a steady state from being reached. The authors are also working in incorporating these stops to the proposed mathematical model, which would allow validation of the model with experimental measurements. Acknowledgements This work was partially supported by ERDF and Xunta de Galicia funds under the ED431C 2017/60 and the ED431C 2021/15 grants, by the Ministerio de Ciencia, Innovación y Universidades through the Plan Nacional de I+D+i (MTM201568275-R) and the grant BES-2016-077228, and by the Agencia Estatal de Investigación through project [PID2019-105615RBI00 / AEI / 10.13039/501100011033]. The authors also wish to thank the reviewers and the Subject Editor for their careful reading of the manuscript and their helpful comments, which have greatly improved its quality. Appendix A. Manufactured solution test derivation The aim of this appendix is, for a given temperature field, to find the corresponding radiosity such that the pair, for suitable data, is the solution of the coupled transient heat conduction in a 2D domain with a nonlocal radiation boundary. In 3D cavities, the derivation of this type of manufactured solution test can be done considering a spherical surface as the radiation cavity, which has the advantage of simplifying the integral kernel, becoming constant [24] . In the 2D case, it is possible to consider a circumference as cavity, which also leads to an integral kernel that is much easier to manipulate analytically than the general expression [25,26] . Nevertheless, the kernel is still not constant, significantly increasing the difficulty of the test derivation. In [25] , a manufactured solution test was proposed for this setting, without including the heat equation in the domain surrounding the radiation cavity. In [26] , the same procedure was applied to obtain a test including the coupled steady heat conduction in the surroundings. However, the insertion of a cut line in the domain was required to mitigate the effects of a singularity in the temperature gradient. Here, we define a new manufactured solution test which overcomes this limitation and incorporates dependence on time. To this end, a simpler approach to obtain the corresponding radiosity in the nonlocal radiation boundary to a given temperature is proposed, which relies on the periodicity of the data on the circumference. Specifically, we look for a temperature that solves the Problem P a , defined in (36) –(38) . The computational domain a , sketched in Fig. 6 (a), is selected as the region bounded by two circumferences of radius r 1 = 1 m and r 2 = 2 m. A nonlocal radiation boundary condition is applied on the inner circumference, denoted as a R . On the exterior circumference, labelled as a D , a Dirichlet boundary condition is set. In Section A.1 , a derivation of the simplified expression of the nonlocal radiation integral kernel on the circumference is presented. This expression is used in Section A.2 to propose a Cauchy problem whose solution for a given temperature yields a radiosity that satisfies the original integral equation on the cavity. Lastly, in Section A.3 , from a set temperature field, the corresponding radiosity is found by solving the defined Cauchy problem, resulting in a new manufactured solution test. A1. Simplified integral kernel We recall the expression (12) , which defines the integral kernel of the operator Kin the 2D setting. In the case of the cavity bounded by a R , we find that the blocking factor is (x , y ) = 1 . To derive the manufactured solution test, we follow a similar approach to the one detailed in [25] . First, we note that the proposed domain enables to derive a simplified expression of the integral kernel ω, denoted as ω a . Bearing in mind that any diameter perpendicular to a chord of the circumference divides it into two equal parts, the relation 2 r 1 cos β2 = 2 r 1 cos β1 = | y −x | holds for any x , y ∈ a R , being β1 and β2 the angles depicted in Fig. A.21 . Using this property, the kernel (12) can be expressed as ω a (x , y ) = cos 2 β1 2 | x −y | = | y −x | 8 , (A.1) since the inner radius is r 1 = 1 . Using a polar coordinate system and that r = r 1 = 1 , we establish that x 1 = cos θ, y 1 = cos ψ, x 2 = sin θ, y 2 = sin ψ, for some θ, ψ ∈ [0 , 2 π) . Furthermore, it is satisfied that | y −x | =  2 −2( cos ψ cos θ+ sin ψ sin θ) . 221