Numerical simulation of the transient heat transfer in a blast furnace main trough during its complete campaign cycle

International Journal of Thermal Sciences 173 (2022) 107349 Available online 18 November 2021 1290-0729/© 2021 The Author(s). Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts Numerical simulation of the transient heat transfer in a blast furnace main trough during its complete campaign cycle P. Barral a,b,c, L.J. Pérez-Péreza,c, P. Quintela a,b,c,∗ aDepartment of Applied Mathematics, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain bTechnological Institute for Industrial Mathematics (ITMATI), 15782 Santiago de Compostela, Spain cInstituto de Matemáticas (IMAT), Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain ARTICLE INFO Keywords: Steelmaking Blast furnace trough Large-scale transient simulation Radiative heat transfer Global optimization ABSTRACT To achieve higher blast furnace (BF) main trough availability and to minimize the frequency of reparations is a key concern in the steelmaking industry. For this purpose, strategies to assess refractory wear are required, which is heavily influenced by the temperature in the refractory linings. In this work, a mathematical model to assess the transient behaviour of the temperature in a cross-section of a BF main trough during a complete campaign is presented. The scope is to investigate the effect that the casting stops have on the temperature in the trough. A sequence of problems corresponding to each BF tapping and the subsequent stop is determined using process data of a BF. The open-source finite element computing platform FEniCS is employed to solve the model. The discretization and the numerical algorithm have been presented and validated with a manufactured solution test in a previous work. The numerical results show that the effect of the stops during these campaign cycles is nonnegligible, preventing the bulk of the solid layers from reaching a steady state. Qualitative agreement with temperature measurements obtained with thermocouples embedded in the trough is observed. Since there is a significant degree of uncertainty concerning the placement of the devices, a minimization problem to adjust their positions within the corresponding feasible regions is proposed. At the identified positions, good levels of fit between the measured and the computed temperatures are achieved. The agreement decreases towards the end of the campaign cycle, being suggestive of severe refractory wear, especially at the laterals of the trough. 1. Introduction In the context of the integrated iron and steelmaking process, it is a major concern for the industry to minimize the wear suffered by refractory linings. Frequently, in the experimental and numerical studies about such process, the focus is put on the blast furnace (BF), the metallurgical furnace where the iron ore is reduced and smelted to obtain hot metal. However, the refractory linings used at the main trough, whose purpose is to transport and separate hot metal from slag after its extraction from the BF through the taphole, are subject to frequent repairs and operation problems with the consequent costs for the companies. Degradation and wear of the refractories used in the steelmaking industry is a complex phenomenon that involves several factors. A brief summary of these was given in [1]. Due to the high-temperature environment posing a significant challenge to experimentation and its critical role in the steelmaking process, a significant body of research using numerical experiences to investigate the BF process is available. ∗Corresponding author at: Department of Applied Mathematics, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain. E-mail addresses: [email protected] (P. Barral), [email protected] (L.J. Pérez-Pérez), [email protected] (P. Quintela). For instance, in [2,3], the transient flow inside the BF hearth was investigated. In [4], an inverse heat transfer model was developed to estimate the wear profile in the hearth wall using the measured temperatures by several thermocouples. In [5], the flow through the BF taphole was modelled attempting to avoid splashing due to entrained air, which is related to higher degrees of wear. Moreover, in [6], several models designed to predict the flow in the BF hearth and taphole were presented, aiming to cover all the various taphole conditions. The interested reader is referred to [7] for a thorough review of the state of numerical simulation inside the BF. Similar multiphysics models, incorporating either fluid flow, thermal radiation and coupled heat transfer have been developed for aluminium furnaces (see e.g. [8,9]). Although the available literature related to the main trough is smaller than that related to the BF itself, some efforts have been made, especially in recent years. In [10], several strategies to increase the main trough productivity were detailed. The intricate flow patterns in the trough have also been a matter of research by several authors. Some https://doi.org/10.1016/j.ijthermalsci.2021.107349 Received 10 June 2021; Received in revised form 16 September 2021; Accepted 20 October 2021 International Journal of Thermal Sciences 173 (2022) 107349 2 P. Barral et al. Nomenclature Greek symbols 𝛼𝑏Backtracking line search parameter 𝛼𝑐Simulated annealing parameter 𝜀Emissivity 𝜆Parameter on transient temperature profile [K/s] 𝛩Blocking factor 𝛤Boundary of computational domain 𝛤𝑅Radiation cavity  𝛤𝑅Intersection between the solid subdomain and the radiation cavity 𝛤𝐷 𝑅Slag upper surface 𝛺(𝑡)Computational domain at time 𝑡 𝜔Kernel of nonlocal radiation integral equation 𝜈Step length in projected gradient descent method 𝜌Density [kg/m3] 𝜎Stefan–Boltzmann constant [W/(m2K4)] 𝜏ℎMesh of the computational domain Latin symbols 𝑐𝑝Specific heat at constant pressure [J/(kg K)] 𝑐𝑛Limited step size change ratio at 𝑛th time step Feasible region 𝐝Descent direction ℎHeat transfer coefficient [W/(m2K)] ℎℎ𝑚 Heigh of the hot metal and slag interface [m] ℎ𝑠𝑙𝑎𝑔 Height of the slag free surface [m] ℎ𝐵Height of the bottom of the trough channel bottom [m] Identity operator 𝐽Cost functional 𝑘Thermal conductivity [W/(m K)] 𝑘𝑖𝑛 Thermal conductivity at the insulation lining [W/(m K)] Nonlocal radiation integral operator 𝐿Thermocouple length [m] 𝑀Number of unit cycles 𝐧Outward-pointing unit normal vector 𝑁Geometric dimension 𝑁𝑐Simulated annealing maximum iterations 𝑞𝑟𝑎𝑑 Heat flux due to radiation [W/m2] Integral operator for radiative heat flux 𝑟𝑛Norm of the local error estimate at the 𝑛th time step 𝑅Radiosity [W/m2] publications have relied on physical scale models, while other authors have used CFD to solve the Navier–Stokes equations to characterize the multiphase flow in the trough. In [11], several experiments were carried out using a physical model to investigate metal separation efficiency. To emulate slag and hot metal, oil and water were used, respectively. The variation of several process variables, such as the BF 𝑡Time [s] 𝑡𝑠𝑡𝑒𝑎𝑑𝑦 Required time to reach the steady state [s]  𝑡𝑐Time at which 𝑐th cast or stop begins [s] 𝑡𝑒𝑛𝑑 Final time of problem [s] Time interval [s] 𝑇Temperature [K] 𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐Convection external temperature on 𝛤𝑛𝑐 𝐶[K] 𝑇𝑒𝑥𝑡,𝑟,𝑛𝑐Radiation external temperature on 𝛤𝑛𝑐 𝐶[K] 𝑇𝑒𝑥𝑡,𝑐,𝑅 Convection external temperature on 𝛤𝑅[K] 𝑇0Initial temperature [K] 𝑊Weight functional Subscripts 𝑐Index of times at which BF activity changes from cast to stop 𝐶External boundaries 𝑐𝑎𝑠𝑡 Cast problems 𝑑Index on thermocouples 𝐷Dirichlet 𝐹Fluid 𝑚Index on unit cycles 𝑛𝑐Index on convection boundaries 𝑅Nonlocal radiation boundary 𝑆Solid 𝑠𝑡𝑜𝑝 Stop problems Superscripts 𝑚Index on unit cycles 𝑛𝑛th time step 𝑛𝑐Index on convection boundaries discharge rate or trough geometry was investigated. Flow separation, as well as its main structures and characteristics, were also experimentally studied in [12] and [13], describing similar findings. Concerning the numerical simulation in the main trough, [14] solved the hot metal flow to investigate the patterns resulting from the jet stream impingement, obtaining qualitatively similar results to those reported in the physical models. In [15], the flow structures were studied by solving a similar mathematical model, also incorporating slag as a different flow phase. Hot metal losses with slag were assessed, finding that geometry changes, such as decreasing the trough bottom slope, have a substantial impact on the slag–metal separation. Moreover, in [16], the multiphase turbulent flow was solved with ANSYS CFX. To validate the CFD results, a 1 ∶ 10 scale model was used to undertake extensive experiments. Flow modifiers were introduced to contain the region where most of the mixing occurs, achieving plug-like flow in a greater region of the trough. In [17], the fluid flow in the trough was also investigated, but incorporating buoyancy due to thermal effects. It was concluded that controlling the height level between the slag runner and the iron dam is essential to improve separation. In addition, in [18], OpenFOAM was used to solve the transient Navier–Stokes equations, also including the conjugate heat transfer with the solid refractory. It was found that the shear stress in the trough wall increases with higher taphole inclination angles. Regarding the thermal response of the through, in [19], the heat transfer in a 2D cross-section was modelled, assuming steady state conditions and considering a simplified treatment for the radiative heat transfer. In a previous work [20], the temperature field in a 3D domain corresponding to the final part of the main trough was studied, including the flow of hot metal and slag, and the refractory cover, also in a steady state. Moreover, in [1], a transient thermal model for a International Journal of Thermal Sciences 173 (2022) 107349 3 P. Barral et al. Fig. 1. Schematic of a typical BF main trough. 2D cross-section of the trough was proposed, considering the time span of a single BF tapping. There, the effect of the thermal radiation was addressed by using a nonlocal radiation boundary condition. The purpose of this work is to study the influence of the casting process stops on the temperature in the main trough. Usually, these stops are kept short. However, various factors can lead to longer stops occurring. In such case, the bulk of the solids may cool down significantly, resulting in large temperature swings and enhanced risk of spalling of the infiltrated refractory layers. Special attention is devoted to the evolution of the position of the so-called critical isotherms in the refractories, which can be used to diagnose the wear profiles in the trough. Therefore, the model developed in [1] for the transient heating during a single BF tapping is extended for the complete campaign cycle. To this end, two different submodels are incorporated: one for those times when the BF is not being cast and a different one for times corresponding to BF tappings. Real BF operation data is used to determine the sequence of problems to be solved, comprising more than 300 BF tappings and the corresponding stops, which span for around two months of duration. The numerical results are compared with experimental measurements obtained with thermocouples embedded in the trough. Provided that there is a substantial degree of uncertainty regarding their placement, a minimization problem is posed and solved. As far as the authors are aware, there is no previous literature on the simulation of complete campaigns, which requires long timescales, nor verification and calibration of their results with experimental data. The outline of this work is as follows. In Section 2, we detail the fundamental characteristics of the physical process. In Section 3, we describe the mathematical model for the transient behaviour of the temperature in a 2D cross-section of the BF trough. In Section 4, we present the computed numerical results. Moreover, in Section 5, the methodology to adjust the position of the thermocouples within their feasible regions is proposed. Lastly, in Section 6, we highlight the main conclusions derived from this paper. 2. Physical problem Slag and hot metal are stratified inside the hearth of the BF. When the hot metal level is sufficiently high, its extraction is done through the taphole, an orifice in the BF wall. The BF draining process is known as tapping. While the process is active, the fluids fall into one of the BF main troughs, where they separate by density difference. In Fig. 1, a schematic of an empty typical main trough is shown before the start of the casting regime. During its operation, it fills with slag and hot metal. Above the trough, a system of removable refractory covers is installed, which spans for the complete length of the runner, even though in Fig. 1 they have been cut for clarity. The purpose of this system is two-fold: first, it allows a more effective aspiration of pollutant fumes emitted during the tapping process (see [21]); and second, it avoids Fig. 2. Temperature schematic at the contact area with fluids of the working lining during a campaign cycle. excessive cooling due to radiation emission from the slag upper surface, which is of key relevance given its high temperature. The fume aspiration, performed next to the taphole, generates a forced air circulation below the covers, which enhances the cooling of the refractories. The BF trough design greatly varies among different casthouses. In this work, we use a design that corresponds to a main trough operated by a steelmaking company, similar to that considered in [20]. As most trough designs, it has three distinct material layers, as shown in Fig. 5. These layers are constructed with different refractory materials that are prepared to perform specific functions. The BF tapping process is not fully continuous, as each cast lasts for about 60-90 minutes with its end being indicated by BF gas bursting out through the taphole, time at which it is plugged. Afterwards, there is a stop to allow the liquid level inside the BF hearth to rise back (see [22]). Depending on the size of the BF and its number of tapholes, a different taphole may be used for the next cast. Large modern BF designs, with up to four tapholes, can be operated with two tapholes being open simultaneously (see [21]). We define a unit cycle as a single cast and the subsequent stop. Due to the severe degradation that sustains, the working lining is replaced after several unit cycles are completed to prevent damage on the safety and insulation linings. We denote as campaign cycle the set of unit cycles that each working lining withstands before its replacement. Between campaign cycles, the fluids in the trough are drained and reparations are performed. Then, the main trough is heated, and a new campaign cycle begins. In Fig. 2, a schematic of the temperature values on the contact area with the fluids of the working lining is depicted during a campaign cycle. The time values  𝑡𝑐, with 𝑐= 0,…,2𝑀− 1, indicate the times at which tappings International Journal of Thermal Sciences 173 (2022) 107349 4 P. Barral et al. Fig. 3. Temperature measurements in the insulation lining during a runner lifetime cycle. are started and ended, 𝑀being the number of tappings or unit cycles. Standard campaign cycles last up to 2months. A runner lifetime cycle is defined as the set of campaign cycles and the drains and repairs between them. When it ends, the main trough is completely rebuilt. Typically, runner lifetime cycles are much longer than campaign cycles and can reach several years of duration. In Fig. 3, temperature measurements within the insulation lining are depicted during campaign cycles, comprising a complete runner lifetime cycle. If the unit cycle stops are kept short, experimental measurements show that the temperature in the bulk of the solids composing the runner reaches a quasi steady state. However, in real unit cycles the downtime of the trough may be far from such situation. This explains the large temperature swings that are observed in Fig. 3 for some campaigns. These larger oscillations correspond to longer stops that can last for up to several days. 3. Mathematical model To assess the cooling related to the longer stops that occurs in real operation conditions and how it may affect the refractory materials, we aim to obtain the temperature field corresponding to a full campaign cycle of a main trough. Therefore, a time interval of around two months of operation needs to be solved. Since it would be extremely difficult to solve a detailed 3D model including the fluid flow for such time scales, we follow a simplified approach. Hence, we consider only the 2D crosssection corresponding to the middle part of the main trough (see Fig. 5), situated before the slag runner at 𝑥3= 8 m, where the coordinate 𝑥3 denotes the distance to the BF. The proposed model is an extension of that described in [1], incorporating the model corresponding to stops. 3.1. Tapping cycles A time period = (𝑡0, 𝑡𝑒𝑛𝑑 ]is considered, which corresponds to a complete campaign cycle. This period can be split in two subsets, 𝑠𝑡𝑜𝑝 and 𝑐𝑎𝑠𝑡, which denote the time interval during which there is no fluid discharge from the BF hearth, and the total cast time interval, respectively. Consequently, we have that = (𝑡0, 𝑡𝑒𝑛𝑑 ] = 𝑠𝑡𝑜𝑝 ∪𝑐𝑎𝑠𝑡.(1) The sets 𝑐𝑎𝑠𝑡 and 𝑠𝑡𝑜𝑝 may be subdivided in many smaller time subintervals, which correspond to the BF casts and stops during each unit cycle. For instance, if we consider a campaign cycle composed of 𝑀 casts, we find: 𝑐𝑎𝑠𝑡 = 𝑀−1 ⋃ 𝑚=0 ( 𝑡2𝑚, 𝑡2𝑚+1] = 𝑀−1 ⋃ 𝑚=0 𝑚 𝑐𝑎𝑠𝑡,(2) Table 1 Example of provided BF activity values. Register index [h] 1 2 3 4 5 6 7 Value 1.0 0.5 0.0 0.1 1.0 0.4 0.0 Fig. 4. Signal indicating whether the BF is tapping or stopped. 𝑠𝑡𝑜𝑝 = 𝑀−2 ⋃ 𝑚=0 ( 𝑡2𝑚+1, 𝑡2𝑚+2] = 𝑀−2 ⋃ 𝑚=0 𝑚 𝑠𝑡𝑜𝑝,(3) where we assume that  𝑡2𝑀−1 =𝑡𝑒𝑛𝑑 . At a given time 𝑡∈, we solve a different thermal model depending on whether the BF is stopped or tapping, i.e., if 𝑡belongs to 𝑚 𝑠𝑡𝑜𝑝 or 𝑚 𝑐𝑎𝑠𝑡 for some 𝑚. We denote the corresponding problems as 𝑚 𝑠𝑡𝑜𝑝 and 𝑚 𝑐𝑎𝑠𝑡, respectively. These problems are solved sequentially, meaning that the Problem 𝑚 𝑐𝑎𝑠𝑡 has as initial condition the final temperature value computed in 𝑚−1 𝑠𝑡𝑜𝑝 , for 𝑚= 1,…𝑀− 1. For 𝑚= 0, which corresponds to the first cast, an initial state for the temperature in the trough is supplied as initial condition, which is discussed in Section 3.5. The same applies for 𝑚 𝑠𝑡𝑜𝑝, which uses as initial value the final temperature in 𝑚 𝑐𝑎𝑠𝑡, for 𝑚= 0,1,…, 𝑀 − 2. To generate the splitting among 𝑐𝑎𝑠𝑡 and 𝑠𝑡𝑜𝑝 following (2) and (3), the time values  𝑡𝑐are required, with 𝑐= 0,…,2𝑀− 1. To find them, we use operation data of a BF, supplied by the steelmaking company collaborating in the research project PID2019-105615RB-I00. The data are composed of values ranging from 0to 1, which are equal to the fraction of each hour that the BF was discharging material through the taphole. Assuming that there is only one cast start or stop during each hour, a function which is 1for 𝑡∈𝑐𝑎𝑠𝑡 and 0if 𝑡∈𝑠𝑡𝑜𝑝 is constructed, similar to the BF activity curve displayed in Fig. 2. The function determines the values  𝑡𝑐. For instance, for the set of values gathered in Table 1, we reconstruct the function shown in Fig. 4. In particular, for this simple example, we find that there are two casts. The first one starts and ends at  𝑡0= 0 h and  𝑡1= 1.5h, respectively. The second cast starts at  𝑡2= 3.9h and ends at  𝑡3= 5.4h. 3.2. Computational domain The computational domain is selected as the cross-section depicted in Fig. 5. We consider a fluid subdomain, 𝛺𝐹, which comprises the area occupied by the slag and hot metal pool, and a solid subdomain, 𝛺𝑆, which corresponds to the solid refractory materials. The steel casing is not included in the computational domain, as it incorporates several fins and ribs, which are not compatible with a 2D thermal model. Nevertheless, their effect is considered by appropriately enhancing the heat losses through the boundaries. During the unit cycle stops, the main trough holds the remaining liquids from preceding casts. This allows to preserve heat in the trough, International Journal of Thermal Sciences 173 (2022) 107349 5 P. Barral et al. Fig. 5. Computational domains and their boundaries. Fig. 6. Free surface positions in the fluids. keeping hot metal and slag in a liquid state, which helps separation in subsequent casts (see e.g [21] or [23]). The trough is drained only to perform emergency reparations or at the end of each campaign cycle. The upper surface of the slag and the interface separating hot metal and slag are not fully steady during the campaign cycle, as the liquid level slightly increases during each cast. The level changes are small and very difficult to measure accurately; in the tapping steady state its value is estimated to grow up to a 15% with respect to the level during the stops. Provided that the scope of this work is to assess the influence of the stops between consecutive casts in the heat transfer towards the solid refractory layers, these liquid level variations can be neglected. Hence, we assume that both the slag and hot metal upper surfaces remain at a constant height during the campaign cycle. This value is calculated assuming that the fluids are at rest during the unit cycle stops. Consequently, the approximate heights depend only on the trough design, specifically on the iron and slag runner height, as well as on the fluid density values, as detailed in [12,17]. Using this approach, considering the coordinate axes and origin depicted in Fig. 5, we find that the slag upper surface is located at ℎ𝑠𝑙𝑎𝑔 = 1.5m and the hot metal surface is at ℎℎ𝑚 = 1.15 m, as shown in Fig. 6. The bottom of the fluids is denoted as ℎ𝐵and the corresponding value for the selected crosssection is ℎ𝐵= 1.04 m. Note that the runner has some inclination to facilitate the sliding of the hot liquids, so ℎ𝐵depends on the coordinate 𝑥3of the considered cross-section. During each cast, we follow the strategy proposed in [1] for a single tapping, where it was assumed that a stationary temperature is reached swiftly in the fluids. To model this phenomenon, the temperature in 𝛺𝐹is assumed to follow a time-dependent temperature profile (see Section 3.4), which depends on the initial temperature. Thus, the domain for the Problem 𝑚 𝑐𝑎𝑠𝑡 is 𝛺𝑐𝑎𝑠𝑡 =𝛺𝑆for all 𝑚, as depicted in Fig. 5(a). In the problem corresponding to unit cycle stops, 𝑚 𝑠𝑡𝑜𝑝, the aim is to assess the cooling of the main trough. The complete cross-section is considered, including the fluids that are retained in the runner during the stops. Hence, for 0≤𝑚≤𝑀− 2, the domain of Problem 𝑚 𝑠𝑡𝑜𝑝 is 𝛺𝑠𝑡𝑜𝑝 =𝛺𝐹∪𝛺𝑆∪ 𝛤𝐷, as displayed in Fig. 5(b). The boundary 𝛤𝐷is defined as 𝛤𝐷=𝛺𝐹∩𝛺𝑆. For the sake of simplicity, we denote by 𝛺(𝑡) Table 2 Properties of the materials. 𝑘[W/(m K)]𝑐𝑝[J∕(kg K)] 𝜌[kg∕m3] Working lining 2.6 1212 2500 Safety lining 1.8 1172 2900 Insulation lining 𝑘𝑖𝑛(𝑇) 1050 800 Cover refractory 2.5 1296 2750 Hot metal 16.5 850 7015 Slag 9.7 807 2600 the computational domain at time 𝑡∈: 𝛺(𝑡) = {𝛺𝑠𝑡𝑜𝑝,if 𝑡∈𝑠𝑡𝑜𝑝, 𝛺𝑐𝑎𝑠𝑡,if 𝑡∈𝑐𝑎𝑠𝑡.(4) 3.3. Model equations To model the evolution of the temperature in the trough crosssection, both for the cast and stop problems, we use the transient heat equation (see [24]) 𝜌𝑐𝑝 𝜕𝑇 𝜕𝑡 −div(𝑘(𝑇)∇𝑇)=0,(5) for 𝑡∈and 𝐱∈𝛺(𝑡), with 𝐱= (𝑥1, 𝑥2). The thermal conductivity, the specific heat, and the density are denoted as 𝑘,𝑐𝑝and 𝜌, respectively, and take the values gathered in Table 2, where 𝑘𝑖𝑛 is the following temperature-dependent function (see [1]): 𝑘𝑖𝑛(𝑇) = ⎧ ⎪ ⎨ ⎪ ⎩ 0.13,if 𝑇≤293, 1 98 000 (9𝑇+ 10 103),if 293 < 𝑇 < 1273, 0.22,if 𝑇≥1273. (6) The values in Table 2 have been supplied by the company, in some cases taken from the data sheets elaborated by the material manufacturer. The constant values are given at the expected operation temperature. Only in the case of the insulation lining, temperaturedependent data are available, so the 𝑘𝑖𝑛 function is obtained by fitting the values. On the boundary 𝛤𝑅, which corresponds to the radiation cavity (see Fig. 5), an integral equation is solved to model the radiation heat exchange. Under the assumption that the radiation cavity constitutes an opaque, diffuse and grey surface, whereas the air within is a transparent medium to radiation, the radiosity 𝑅on 𝛤𝑅satisfies, for 𝑡∈, that (see e.g. [25]): (− (1 − 𝜀𝑅(𝐱)))(𝑅(𝑡))(𝐱) = 𝜀𝑅(𝐱)𝜎𝑇 4(𝐱, 𝑡),on 𝛤𝑅,(7) where is the identity operator and 𝜀𝑅is the emissivity, considered constant with value 𝜀𝑅= 0.7. Moreover, 𝜎= 5.67𝑒−8 W/(m2K4) is International Journal of Thermal Sciences 173 (2022) 107349 6 P. Barral et al. Table 3 Heat transfer coefficients. ℎ[W/(m2K)] ℎ𝑅35 ℎ17 ℎ28 ℎ34 the Stefan–Boltzmann constant and denotes the following integral operator (𝑅(𝑡))(𝐱) = ∫𝛤𝑅 𝜔(𝐱,𝐲)𝑅(𝐲, 𝑡)𝑑𝑠𝐲,(8) where the integral kernel 𝜔, known as view factor, is such that (see [26]) 𝜔(𝐱,𝐲) = 𝐧(𝐱)⋅(𝐲−𝐱)𝐧(𝐲)⋅(𝐱−𝐲) 2|𝐱−𝐲|3𝛩(𝐱,𝐲),with 𝐱≠𝐲,(9) 𝐧being the outward-pointing unit normal vector. The visibility or blocking factor 𝛩accounts for obstacles that obstruct the view between points and is defined as 𝛩(𝐱,𝐲) = {0,if 𝐱𝐲 ∩𝛺≠∅, 1,if 𝐱𝐲 ∩𝛺= ∅,(10) with 𝐱𝐲 denoting the segment that connects the points 𝐱and 𝐲. To evaluate the radiative heat flux, 𝑞𝑟𝑎𝑑 , which appears as a term in the boundary condition corresponding to 𝛤𝑅of the heat Eq. (5), the integral operator is introduced, which yields 𝑞𝑟𝑎𝑑 (see [1]): 𝑞𝑟𝑎𝑑 (𝐱, 𝑡) = (𝑅(𝑡))(𝐱) ∶= (−)(𝑅(𝑡))(𝐱).(11) For a more detailed description of the nonlocal radiation model, we refer to [1,27]. A theoretical analysis of the steady coupled conduction and nonlocal radiation can be seen in [28] and the references therein. The transient case is addressed in [29]. Note that in general applications, the hypothesis that are assumed to describe thermal radiation using (7) are far from realistic. This is the case of environments with combustion (see e.g. [30]), where specific methods have to be used to solve the radiative transfer equation (see [25,31]). 3.4. Boundary conditions The boundaries of the domains are decomposed as displayed in Fig. 5 for 𝛺𝑐𝑎𝑠𝑡 and 𝛺𝑠𝑡𝑜𝑝. We start by briefly summarizing the boundary conditions used for the cast problems, which are the same already presented with detail in [1]. Subsequently, the modifications introduced for the stop problems are described. 3.4.1. Cast problems In the following discussion, we detail the boundary conditions considered for each 𝑚 𝑐𝑎𝑠𝑡. The boundary of its domain is independent of 𝑚and is given by 𝛤𝑐𝑎𝑠𝑡 = 3 ⋃ 𝑛𝑐=1 𝛤𝑛𝑐 𝐶∪𝛤𝐷∪ 𝛤𝑅,(12) as depicted in Fig. 5. •On the external boundaries 𝛤𝑛𝑐 𝐶,1≤𝑛𝑐≤3, mixed convection– radiation boundary conditions are considered (see [24]): −𝑘(𝑇)𝜕𝑇 𝜕𝐧=ℎ𝑛𝑐(𝑇−𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐) + 𝜎𝜀𝑛𝑐(𝑇4−𝑇4 𝑒𝑥𝑡,𝑟,𝑛𝑐),(13) for all 𝑡∈𝑐𝑎𝑠𝑡. The values for the heat transfer coefficients, ℎ𝑛𝑐, are displayed in Table 3. The convection and radiation ambient temperatures as well as the emissivity are considered equal (𝑇𝑒𝑥𝑡,𝑟,𝑛𝑐=𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐= 293 K and 𝜀𝑛𝑐= 0.7). Fig. 7. Temperature profile 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡 for a typical initial temperature during the trough campaign cycle. •On 𝛤𝐷=𝛺𝐹∩𝛺𝑆, which is the interface between the fluids and solids, we set a Dirichlet condition, using a known temperature: 𝑇𝑚 𝑐𝑎𝑠𝑡(𝐱, 𝑡) = 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡(𝐱, 𝑡),(14) for all 𝐱∈𝛤𝐷and 𝑡∈𝑚 𝑐𝑎𝑠𝑡. In (14), the temperature 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡 is the profile developed in [1], suitably modifying the initial temperature considered therein so that the temperature at the fluids rises from the temperature at the initial time of the 𝑚th cast problem to a stationary temperature. The stationary temperature satisfies the following expression: 𝑇𝐹 ,𝑐𝑎𝑠𝑡(𝑥2) = −1.52𝑒−9𝑥63.46 2+ 1773,(15) where ℎ𝐵≤𝑥2≤ℎ𝑠𝑙𝑎𝑔 , with ℎ𝐵= 1.04 m, as detailed in Section 3.2. The temperature profile (15) is a fitting of the solution computed for the steady state in [20] at 𝑥3= 8 and 𝑥1= 1.5 m. Note that, as described in [1], any variations of the stationary temperature along the horizontal coordinate are neglected. Specifically, at unit cycle 𝑚, with 1≤𝑚≤𝑀− 1, and time 𝑡∈𝑚 𝑐𝑎𝑠𝑡, we set the slag and hot metal temperature to 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡(𝐱, 𝑡) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑇𝑚−1 𝑠𝑡𝑜𝑝 (𝐱, 𝑡2𝑚) − 3𝜆2(𝑡− 𝑡2𝑚)2 𝛥𝑇 𝑚(𝐱) −2𝜆3(𝑡− 𝑡2𝑚)3 (𝛥𝑇 𝑚(𝐱))2,if 𝑡≤𝑡𝑚 𝑠𝑡𝑒𝑎𝑑𝑦(𝐱), 𝑇𝐹 ,𝑐𝑎𝑠𝑡(𝑥2),otherwise, (16) where  𝑡2𝑚denotes the time value corresponding to the end of the preceding stop, and 𝑇𝑚−1 𝑠𝑡𝑜𝑝 (𝐱, 𝑡2𝑚)is the computed temperature for 𝑚−1 𝑠𝑡𝑜𝑝 at its end time. Moreover, 𝛥𝑇 𝑚(𝐱)stands for the difference among the steady temperature and the initial temperature at point 𝐱, i.e.: 𝛥𝑇 𝑚(𝐱) = 𝑇𝑚−1 𝑠𝑡𝑜𝑝 (𝐱, 𝑡2𝑚) − 𝑇𝐹 ,𝑐𝑎𝑠𝑡(𝑥2).(17) If 𝑚= 0, corresponding to the first unit cycle of the trough campaign, in (16) we consider 𝑇𝑚−1 𝑠𝑡𝑜𝑝 =𝑇0, where 𝑇0is the initial temperature value, discussed in Section 3.5. Furthermore, the time needed to reach the stationary temperature is 𝑡𝑚 𝑠𝑡𝑒𝑎𝑑𝑦(𝐱) =  𝑡2𝑚+ 𝑇𝐹 ,𝑐𝑎𝑠𝑡(𝑥2) − 𝑇𝑚−1 𝑠𝑡𝑜𝑝 (𝐱, 𝑡2𝑚) 𝜆.(18) The parameter 𝜆takes the value 𝜆= 3∕2 [K/s] in (16) and (18). This choice guarantees that in the most extreme scenario (𝑇𝑚−1 𝑠𝑡𝑜𝑝 ≈ 293 K), the stationary temperature is reached in approximately 15 min. For a regular casting regime, the stops are not long enough for the fluids to cool down significantly, and the steady state is reached much faster, as depicted in Fig. 7 using a typical value for the initial temperature 𝑇𝑚−1 𝑠𝑡𝑜𝑝 (𝐱, 𝑡2𝑚), extracted from the numerical results discussed in Section 4. International Journal of Thermal Sciences 173 (2022) 107349 7 P. Barral et al. Fig. 8. Temperature contours of the initial temperature 𝑇0. •On  𝛤𝑅=𝛤𝑅∩𝛺𝑆, a convection contribution (see [24]) due to the forced air flow below the cover is added to the radiation heat flux 𝑞𝑟𝑎𝑑 , obtained by solving (7): −𝑘𝜕𝑇 𝜕𝐧=ℎ𝑅(𝑇−𝑇𝑒𝑥𝑡,𝑐,𝑅) + 𝑞𝑟𝑎𝑑 ,(19) for all 𝑡∈𝑐𝑎𝑠𝑡 and where the temperature of the surroundings is 𝑇𝑒𝑥𝑡,𝑐,𝑅 = 293 K. The heat transfer coefficient, ℎ𝑅, is equal to the value indicated in Table 3. Note that even though the radiation contribution is only applied on the part of the radiation cavity that belongs to the boundary of 𝛺𝑐𝑎𝑠𝑡, the integral Eq. (7) is solved on the complete cavity 𝛤𝑅, as the radiosity on the slag upper surface, 𝛤𝐷 𝑅, is unknown even if the temperature satisfies (16). Therefore, for 𝑚 𝑐𝑎𝑠𝑡,(7) is solved for all 𝑡∈𝑚 𝑐𝑎𝑠𝑡 assuming that 𝑇|𝛤𝐷 𝑅=𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡. 3.4.2. Stop problems For 𝑚 𝑠𝑡𝑜𝑝, their domain, 𝛺𝑠𝑡𝑜𝑝, and therefore their boundary, 𝛤𝑠𝑡𝑜𝑝, are also independent of 𝑚. In this case, we have that 𝛤𝑠𝑡𝑜𝑝 = 3 ⋃ 𝑛𝑐=1 𝛤𝑛𝑐 𝐶∪𝛤𝑅,(20) as displayed in Fig. 5. On the external part of the trough, the boundary conditions are the same for the cast and the stop problems. Thus, for 𝑚 𝑠𝑡𝑜𝑝, on 𝛤𝑛𝑐 𝐶, with 1≤𝑛𝑐≤3, we set (13) for all 𝑡∈𝑠𝑡𝑜𝑝. However, note that the radiation cavity boundary is slightly different for the stop problems when compared to the cast problems, as depicted in Fig. 5, where the boundary becomes 𝛤𝑅instead of  𝛤𝑅, as during the stops the temperature on the fluids is not known and the heat equation is also solved in 𝛺𝐹. 3.5. Initial condition The reparations between campaign cycles last for around 20 days, as displayed in Fig. 2. During this period, the cover is removed, and the refractory layers cool down as reparations are performed. Also, the working lining is heated after being rebuilt to prevent damage from the initial thermal shock. The temperature state after this process, which corresponds to the initial temperature at the start of a campaign cycle, is unknown. Nevertheless, an approximated temperature field must be supplied as initial condition at 𝑡= 𝑡0for the proposed model. As it will be detailed in Section 5, the measurements obtained with three thermocouples located in the considered cross-section of the trough are available. They show that there is a significant amount of residual heat in the trough at the start of the campaign cycles, with the bottom being at a higher temperature than the laterals. To obtain the required initial condition, a possibility would be to consider an averaged homogeneous value using the three thermocouple measurements in the two campaign cycles for which we have data, equal to 322 K. However, it is not very advisable to capture a 2D temperature field in the whole crosssection only using the measurements in these three points, given that the temperature in the safety lining could be substantially larger than in the insulation lining, where the thermocouples are placed. Hence, to find a more realistic initial temperature we propose the following procedure: (1) A steady state problem corresponding to 𝑚 𝑐𝑎𝑠𝑡 is solved to obtain the steady state temperature in the solids. This involves solving the heat Eq. (5) in 𝛺𝑐𝑎𝑠𝑡, assuming a steady state as well as 𝑇𝑚 𝑐𝑎𝑠𝑡 =𝑇𝐹 ,𝑐𝑎𝑠𝑡 on 𝛤𝐷, and satisfying boundary conditions (13) and (19). (2) A Problem 𝑚 𝑠𝑡𝑜𝑝 using the steady state temperature computed in the previous step as initial condition is solved, until reaching a time value 𝑡= 20 days. Following these steps, we obtain the temperature field depicted in Fig. 8. The temperature distribution in the cross-section shows that the residual heat is concentrated in the area close to the bottom of the trough, as predicted by the thermocouple measurements. However, it is observed that the temperatures within the safety lining are substantially higher than those in the insulation lining, where the thermocouples are located. The obtained temperature is denoted as 𝑇0and is used as initial condition for the first cast, i.e., the Problem 0 𝑐𝑎𝑠𝑡. Moreover, as we are solving a sequence of problems alternating between casts and stops, each of these needs an initial condition, which corresponds to the computed solution at the end of the previous problem. 3.6. Strong formulation A summary of the complete formulation of the two problems that are sequentially solved is provided below, including the PDE, the radiation integral equation and the corresponding boundary conditions that were discussed previously: Problem 𝑚 𝑐𝑎𝑠𝑡. Find 𝑇𝑚 𝑐𝑎𝑠𝑡(𝐱, 𝑡)and 𝑅𝑚 𝑐𝑎𝑠𝑡(𝐱, 𝑡)such that: 𝜌𝑐𝑝 𝜕𝑇 𝑚 𝑐𝑎𝑠𝑡 𝜕𝑡 −div(𝑘(𝑇𝑚 𝑐𝑎𝑠𝑡)∇𝑇𝑚 𝑐𝑎𝑠𝑡)=0,in 𝛺𝑐𝑎𝑠𝑡 ×𝑚 𝑐𝑎𝑠𝑡,(21) 𝑇𝑚 𝑐𝑎𝑠𝑡( 𝑡2𝑚) = {𝑇0,if 𝑚= 0, 𝑇𝑚−1 𝑠𝑡𝑜𝑝 ( 𝑡2𝑚),if 𝑚≥1,in 𝛺𝑐𝑎𝑠𝑡,(22) 𝑇𝑚 𝑐𝑎𝑠𝑡 =𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡,on 𝛤𝐷×𝑚 𝑐𝑎𝑠𝑡,(23) −𝑘(𝑇𝑚 𝑐𝑎𝑠𝑡)𝜕𝑇 𝑚 𝑐𝑎𝑠𝑡 𝜕𝐧=ℎ𝑛𝑐(𝑇𝑚 𝑐𝑎𝑠𝑡 −𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐) +𝜎𝜀𝑛𝑐((𝑇𝑚 𝑐𝑎𝑠𝑡)4−𝑇4 𝑒𝑥𝑡,𝑟,𝑛𝑐),on 𝛤𝑛𝑐 𝐶×𝑚 𝑐𝑎𝑠𝑡,(24) −𝑘𝜕𝑇 𝑚 𝑐𝑎𝑠𝑡 𝜕𝐧=ℎ𝑅(𝑇𝑚 𝑐𝑎𝑠𝑡 −𝑇𝑒𝑥𝑡,𝑐,𝑅) + 𝑞𝑟𝑎𝑑 ,on  𝛤𝑅×𝑚 𝑐𝑎𝑠𝑡,(25) 𝑞𝑟𝑎𝑑 =(𝑅𝑚 𝑐𝑎𝑠𝑡),on 𝛤𝑅×𝑚 𝑐𝑎𝑠𝑡,(26) (− (1 − 𝜀𝑅))(𝑅𝑚 𝑐𝑎𝑠𝑡) = 𝜀𝑅𝜎(𝑇𝑚 𝑐𝑎𝑠𝑡)4,on 𝛤𝑅×𝑚 𝑐𝑎𝑠𝑡,(27) where 1≤𝑛𝑐≤3and 0≤𝑚≤𝑀− 1. In (22), at the first cast of the main trough campaign cycle, the initial condition, 𝑇0, is that already described in Section 3.5. In (23), the temperature 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡 is given by (16). Problem 𝑚 𝑠𝑡𝑜𝑝. Find 𝑇𝑚 𝑠𝑡𝑜𝑝(𝐱, 𝑡)and 𝑅𝑚 𝑠𝑡𝑜𝑝(𝐱, 𝑡)such that: 𝜌𝑐𝑝 𝜕𝑇 𝑚 𝑠𝑡𝑜𝑝 𝜕𝑡 −div(𝑘(𝑇𝑚 𝑠𝑡𝑜𝑝)∇𝑇𝑚 𝑠𝑡𝑜𝑝)=0,in 𝛺𝑠𝑡𝑜𝑝 ×𝑚 𝑠𝑡𝑜𝑝,(28) 𝑇𝑚 𝑠𝑡𝑜𝑝( 𝑡2𝑚+1) = 𝑇∗,𝑚 𝑐𝑎𝑠𝑡( 𝑡2𝑚+1),in 𝛺𝑠𝑡𝑜𝑝,(29) −𝑘(𝑇𝑚 𝑠𝑡𝑜𝑝) 𝜕𝑇 𝑚 𝑠𝑡𝑜𝑝 𝜕𝐧=ℎ𝑛𝑐(𝑇𝑚 𝑠𝑡𝑜𝑝 −𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐) +𝜎𝜀𝑛𝑐((𝑇𝑚 𝑐𝑎𝑠𝑡)4−𝑇4 𝑒𝑥𝑡,𝑟,𝑛𝑐),on 𝛤𝑛𝑐 𝐶×𝑚 𝑐𝑎𝑠𝑡,(30) International Journal of Thermal Sciences 173 (2022) 107349 8 P. Barral et al. −𝑘 𝜕𝑇 𝑚 𝑠𝑡𝑜𝑝 𝜕𝐧=ℎ𝑅(𝑇𝑚 𝑠𝑡𝑜𝑝 −𝑇𝑒𝑥𝑡,𝑐,𝑅) + 𝑞𝑟𝑎𝑑 ,on 𝛤𝑅×𝑚 𝑠𝑡𝑜𝑝,(31) 𝑞𝑟𝑎𝑑 =(𝑅𝑚 𝑠𝑡𝑜𝑝),on 𝛤𝑅×𝑚 𝑠𝑡𝑜𝑝,(32) (− (1 − 𝜀𝑅))(𝑅𝑚 𝑠𝑡𝑜𝑝) = 𝜀𝑅𝜎(𝑇𝑚 𝑠𝑡𝑜𝑝)4,on 𝛤𝑅×𝑚 𝑠𝑡𝑜𝑝,(33) where 1≤𝑛𝑐≤3and 0≤𝑚≤𝑀− 2. Moreover, with 𝑇∗,𝑚 𝑐𝑎𝑠𝑡 we denote an extension of the computed temperature for the cast problem using (16): 𝑇∗,𝑚 𝑐𝑎𝑠𝑡(𝐱, 𝑡2𝑚+1) = {𝑇𝑚 𝑐𝑎𝑠𝑡(𝐱, 𝑡2𝑚+1),if 𝐱∈𝛺𝑐𝑎𝑠𝑡, 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡(𝐱, 𝑡2𝑚+1),otherwise.(34) 3.7. Weak formulation The weak form of both 𝑚 𝑐𝑎𝑠𝑡 and 𝑚 𝑠𝑡𝑜𝑝, denoted as 𝑚 𝑐𝑎𝑠𝑡 and 𝑚 𝑠𝑡𝑜𝑝 respectively, is presented below. The weak formulation of the problems is the form that is solved numerically using a Finite Element Method, and is obtained analogously to that proposed in [1] for a tapping problem. Problem 𝑚 𝑐𝑎𝑠𝑡. For each 𝑡∈𝑚 𝑐𝑎𝑠𝑡, find 𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡), defined in 𝛺𝑐𝑎𝑠𝑡, and 𝑅𝑚 𝑐𝑎𝑠𝑡(𝑡), defined on 𝛤𝑅, such that 𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡) = 𝑇𝑚 𝐷,𝑐𝑎𝑠𝑡(𝑡)on 𝛤𝐷and ∫𝛺𝑐𝑎𝑠𝑡 𝜌𝑐𝑝 𝜕𝑇 𝑚 𝑐𝑎𝑠𝑡(𝑡) 𝜕𝑡 𝑣 𝑑𝐴 +∫𝛺𝑐𝑎𝑠𝑡 𝑘(𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡))∇𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡)⋅∇𝑣 𝑑𝐴 +∫ 𝛤𝑅 ℎ𝑅𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡)𝑣 𝑑𝑠 + 3 ∑ 𝑛𝑐=1 ∫𝛤𝑛𝑐 𝐶(ℎ𝑛𝑐𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡) + 𝜎𝜀𝑛𝑐(𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡))4)𝑣 𝑑𝑠 (35) = 3 ∑ 𝑛𝑐=1 ∫𝛤𝑛𝑐 𝐶(ℎ𝑛𝑐𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐+𝜎𝜀𝑛𝑐𝑇4 𝑒𝑥𝑡,𝑟,𝑛𝑐)𝑣 𝑑𝑠 +∫ 𝛤𝑅 (ℎ𝑅𝑇𝑒𝑥𝑡,𝑅 −(𝑅𝑚 𝑐𝑎𝑠𝑡(𝑡)))𝑣 𝑑𝑠, for all test functions 𝑣defined in 𝛺𝑐𝑎𝑠𝑡 being sufficiently regular, and satisfying: (− (1 − 𝜀𝑅))(𝑅𝑚 𝑐𝑎𝑠𝑡(𝑡)) = 𝜀𝑅𝜎(𝑇𝑚 𝑐𝑎𝑠𝑡(𝑡))4,on 𝛤𝑅,(36) for 𝑡∈𝑚 𝑐𝑎𝑠𝑡. Furthermore, 𝑇𝑚 𝑐𝑎𝑠𝑡( 𝑡2𝑚) = 𝑇𝑚−1 𝑠𝑡𝑜𝑝 ( 𝑡2𝑚)for 1≤𝑚≤𝑀− 1. If 𝑚= 0, then 𝑇0 𝑐𝑎𝑠𝑡(0) = 𝑇0. Problem 𝑚 𝑠𝑡𝑜𝑝. For each 𝑡∈𝑚 𝑠𝑡𝑜𝑝, find 𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡), defined in 𝛺𝑠𝑡𝑜𝑝, and 𝑅𝑚 𝑠𝑡𝑜𝑝(𝑡), defined on 𝛤𝑅, such that ∫𝛺𝑠𝑡𝑜𝑝 𝜌𝑐𝑝 𝜕𝑇 𝑚 𝑠𝑡𝑜𝑝(𝑡) 𝜕𝑡 𝑣 𝑑𝐴 +∫𝛺𝑠𝑡𝑜𝑝 𝑘(𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡))∇𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡)⋅∇𝑣 𝑑𝐴 +∫𝛤𝑅 ℎ𝑅𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡)𝑣 𝑑𝑠 + 3 ∑ 𝑛𝑐=1 ∫𝛤𝑛𝑐 𝐶(ℎ𝑛𝑐𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡) + 𝜎𝜀𝑛𝑐(𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡))4)𝑣 𝑑𝑠 (37) = 3 ∑ 𝑛𝑐=1 ∫𝛤𝑛𝑐 𝐶(ℎ𝑛𝑐𝑇𝑒𝑥𝑡,𝑐,𝑛𝑐+𝜎𝜀𝑛𝑐𝑇4 𝑒𝑥𝑡,𝑟,𝑛𝑐)𝑣 𝑑𝑠 +∫𝛤𝑅 (ℎ𝑅𝑇𝑒𝑥𝑡,𝑅 −(𝑅𝑚 𝑠𝑡𝑜𝑝(𝑡)))𝑣 𝑑𝑠, for all test functions 𝑣defined in 𝛺𝑠𝑡𝑜𝑝 being sufficiently regular, and satisfying: (− (1 − 𝜀𝑅))(𝑅𝑚 𝑠𝑡𝑜𝑝(𝑡)) = 𝜀𝑅𝜎(𝑇𝑚 𝑠𝑡𝑜𝑝(𝑡))4,on 𝛤𝑅,(38) for 𝑡∈𝑚 𝑠𝑡𝑜𝑝 and 0≤𝑚≤𝑀− 2. In addition, 𝑇𝑚 𝑠𝑡𝑜𝑝( 𝑡2𝑚+1) = 𝑇∗,𝑚 𝑐𝑎𝑠𝑡( 𝑡2𝑚+1) (see (34)). 4. Numerical results In this section, the numerical results concerning the heat transfer in the BF main trough cross-section are presented. The model is Fig. 9. Mesh of the BF main trough cross-section. solved using the open-source finite element computing platform FEniCS (see [32,33]), which allows to solve PDEs from a high-level framework. The different components of FEniCS enable to automatically perform the steps to generate efficient code to discretize the problem using a finite element method, using user-specified high level abstractions similar to the problem weak formulation as starting point. A triangulation 𝜏ℎof the domain 𝛺𝑠𝑡𝑜𝑝 is constructed such that 𝛺𝑠𝑡𝑜𝑝 = ∪𝐾∈𝜏ℎ𝐾, shown in Fig. 9. The triangulation is conformal to the different materials in the BF main trough. Thus, the triangulation used for the domain 𝛺𝑐𝑎𝑠𝑡 is a subset of 𝜏ℎand is identical to that presented in [1], where a mesh sensitivity study was carried out. It is formed by a total 65 440 triangles whose restrictions to 𝛤𝑅are a total of 400 edges. From the weak formulations detailed in Section 3.7, a method of lines leads to the spatial semi-discretization, whereas the radiosity integral equations are discretized using a collocation method, also following the procedure described in [1]. The H211b adaptive step size controller (see [34]) is used to adjust the step size using the local error estimate provided by the ESDIRK 3/2a RK scheme, which simply becomes the difference between the last two stages in virtue of its properties (see [35]). Note that the aim is to solve a problem with very large time scales (in the order of months) with mostly gentle temperature changes, but also involving short periods of sharp temperature change in time, especially at the beginning of each cast. Therefore, step size adaptivity as well as the use of high-order time discretization schemes are essential to control the committed error while reducing the computational cost. 4.1. Algorithm To couple the radiosity integral equation with the heat equation, we use the segregated approach described in [1], which involves a fixed point algorithm for each RK stage. It consists of three embedded loops; an inner one to compute the radiosity and temperature at a given stage and time step, another for the time stepping and an outer loop to advance throughout the sequence of cast and stop problems. It follows the flow chart depicted in Fig. 10. The fixed point algorithm control is based on the relative incremental residual for both the temperature and radiosity, which must fall below a given threshold. In addition, the step size control rejects and recomputes the current time step if the change ratio is below a given value of 0.9. If the end time of the corresponding problem is reached, then the step size is adjusted accordingly and, if that step is finally accepted, the algorithm advances to the subsequent problem. The same parameters for the H211b step size controller to those reported in [1] are used in this work, with the only difference that a maximum step size 𝛥𝑡 = 300 s is set during the cast problems to improve the stability of the algorithm. International Journal of Thermal Sciences 173 (2022) 107349 9 P. Barral et al. Fig. 10. Flow diagram of the algorithm. Table 4 Start date, end date and total number of casts (𝑀) for each solved campaign cycle. Start date End date Number of casts Campaign cycle 1 28/05/2018 28/07/2018 347 Campaign cycle 2 21/08/2018 18/10/2018 325 4.2. Simulation of two campaign cycles The problem described in Section 3.6, composed of a sequence of cast and stop problems, is solved. This sequence is determined by data spanning two different campaign cycles to assess whether the BF is being cast, as described in Section 3.1. A summary of the main characteristics of the considered campaign cycles is displayed in Table 4. The two solved cycles are consecutive, with the downtime between them removed. Table 5 Maximum step size, number of accepted steps, number of rejected steps and acceptance rate using the H211b controller for the considered campaign cycles. Max. step size (s) Accepted steps Rejected steps Acceptance rate Campaign 1 1.89𝑒+04 21 509 408 9.81𝑒−01 Campaign 2 9.05𝑒+03 21 171 411 9.81𝑒−01 In Fig. 11, the step size value is plotted against the time value for the two considered campaign cycles, from day 10 to day 20. The step size rises quickly from the start size at each cast or stop problem, which is 𝛥𝑡0= 1 s. Frequently, this first step is rejected and recomputed with a smaller value. Furthermore, the presence of long stops during the campaign cycles is evident, as they allow for much larger step sizes without increasing the error beyond the selected tolerance value. In Table 5, the maximum step size, the number of accepted and rejected steps, and the acceptance rate are gathered for the two campaign cycles that are solved. The acceptance rate is above 0.98 in both cases, while the maximum step size in each campaign is determined by the length of the longest stop during the BF activity recorded during the campaign, which lasts for longer than two days, as shown in Fig. 11. In Fig. 12(a), the computed temperature field 𝑇0 𝑐𝑎𝑠𝑡,ℎ, solution of Problem 0 𝑐𝑎𝑠𝑡, after the first cast of the BF in Campaign cycle 1is displayed, i.e. at time  𝑡1. Even though during each 𝑚 𝑐𝑎𝑠𝑡 only the domain 𝛺𝑐𝑎𝑠𝑡 =𝛺𝑆is solved, the fluids are also included in the shown contour plots as their temperature is assumed to be known and equal to (16), as described in Section 3.4. High temperatures are observed only in the zone that is close to the liquids, as most of the solid is still unaffected, showing similar temperature values to those of the initial condition (see Fig. 8). In contrast, the parts of the working lining that are affected by the radiation from the slag surface show significant heating, as well as the refractory cover, reaching temperatures around 1100 K. To showcase the effect of the short stops in the runner activity, the temperature after the subsequent stop following the first tapping, 𝑇0 𝑠𝑡𝑜𝑝,ℎ(⋅, 𝑡2)is depicted in Fig. 12(b). Approximately half an hour has passed since the end of the preceding cast, and most of the cooling takes place in the fluids. On the other hand, the cover and the working lining still remain at a high temperature. In Fig. 12(c), the temperature contours after approximately 10 days of casts and stops are shown. Specifically, the contours correspond to 𝑇76 𝑐𝑎𝑠𝑡,ℎ(⋅, 𝑡153). The temperature field is substantially different to that displayed in Figs. 12(a) and 12(b), as the bulk of the solids displays a much higher temperature. Towards the end of the campaign, the temperature is similar, as depicted in Fig. 12(f) for 𝑡= 57.84 days. However, the temperature field is far from a steady state, as demonstrated by the temperature contours shown in Fig. 12(d). This temperature snapshot corresponds to a time value after a long stop, which lasted for around 2.2 days. The cover and the upper part of the working lining and fluids has substantially cooled down. Nonetheless, due to the large thermal inertia of the trough, the bulk of the solids remains above 1000 K. Fig. 11. Time step values computed by the H211b controller during the two campaign cycles. International Journal of Thermal Sciences 173 (2022) 107349 16 P. Barral et al. Fig. 21. Temperature over time at the corrected thermocouple locations. Despite the simplifying assumptions that were made, the computed temperatures show qualitative agreement with the measurements. Nevertheless, as there is a substantial degree of uncertainty concerning the position of the thermocouples, a technique to find a corrected position inside a feasible region was used. In particular, the hybrid GRSA algorithm was employed to find the position for each thermocouple, minimizing the corresponding cost functional. This position is calculated as the one that yields the closest time-dependent numerical temperature profile to the measurement. The computed placements are similar for the two campaigns and three analysed thermocouples. At these corrected positions, much better levels of agreement were achieved among the measurements and the computed temperatures. Provided that the model does not account for wear or material deterioration, the discrepancies among the computed temperature values and the measurements, which grow towards the end of the campaigns, suggest severe damage and loss of properties in the working lining. The results also indicate that most of the wear is located in the lateral sections of the runner, as the differences are much larger for the lateral thermocouples. In future work, further studies should be conducted to experimentally validate the model. More detailed experimental measurements of the material properties should be carried out to address their temperature dependency. The model could also be improved by introducing the phase change that takes place in the fluids during the stops, especially during those that are prolonged. Using better thermocouple monitoring, both concerning the number of devices and the reliability of their placement, an inverse heat transfer model could be developed, which would allow to use the computed deviation of the numerical solution from the experimental measurements to assess the wear profile in the working lining. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was partially supported by ERDF and Xunta de Galicia funds under the ED431C 2017/60 grant, by the Ministerio de Ciencia, Innovación y Universidades through the Plan Nacional de I+D+i (MTM2015-68275-R) and the grant BES-2016-077228, and by the Agencia Estatal de Investigación through project [PID2019-105615RB- I00/ AEI/10.13039/501100011033]. References [1] P. Barral, L.J. Pérez-Pérez, P. Quintela, Transient Thermal Response with Nonlocal Radiation of a Blast Furnace Main Trough, Zenodo, 2021, http://dx. doi.org/10.5281/zenodo.4923518. [2] B.-Y. Guo, D. Maldonado, P. Zulli, A.-B. Yu, CFD modelling of liquid metal flow and heat transfer in blast furnace hearth, ISIJ Int. 48 (12) (2008) 1676–1685. [3] W. Cheng, E. Huang, S. Du, Numerical analysis on transient thermal flow of the blast furnace hearth in tapping process through CFD, Int. Commun. Heat Mass 57 (2014) 13–21. [4] Y. Kaymak, H. Bartusch, T. Hauck, J. Mernitz, H. Rausch, R. Lin, Multiphysics model of the hearth lining state, Steel Res. Int. 91 (11) (2020) 2000055. [5] P. Stevenson, Q. He, Slug flow in a blast furnace taphole, Chem. Eng. Process. 44 (10) (2005) 1094–1097. [6] L. Shao, H. Saxén, A simulation study of two-liquid flow in the taphole of the blast furnace, ISIJ Int. 53 (6) (2013) 988–994. [7] S. Kuang, Z. Li, A. Yu, Review on modeling and simulation of blast furnace, Steel Res. Int. 89 (1) (2018) 1700071. [8] L. Qiu, Y. Feng, Z. Chen, Y. Li, X. Zhang, Numerical simulation and optimization of the melting process for the regenerative aluminum melting furnace, Appl. Therm. Eng. 145 (2018) 315–327. [9] L. Qiu, Y. Li, Y. Feng, Z. Chen, X. Zhang, Three-dimensional fluid-solid coupling heat transfer simulation based on the multireference frame for a side-blown aluminum annealing furnace, Eng. Appl. Comput. Fluid Mech. 13 (1) (2019) 1036–1048. [10] A. Kumar, S.A. Khan, S. Biswas, A. Pal, Strategic steps towards longer and reliable blast furnace trough campaign — Tata Steel experience, Ironmak. Steelmak. 37 (2010) 15–20. [11] H. Kim, B. Ozturk, Slag-metal separation in the blast furnace trough, ISIJ Int. 38 (5) (1998) 430–439. [12] M.J. Luomala, T.T. Paananen, M.J. Köykkä, J. Fabritius, T. Matti, H. Nevala, J.J. Härkki, Modelling of fluid flows in the blast furnace trough, Steel Res. Int. 72 (4) (2001) 130–135. [13] Q. He, G. Evans, P. Zulli, F. Tanzil, B. Lee, Flow characteristics in a blast furnace trough, ISIJ Int. 42 (8) (2002) 844–851. [14] R.V.P. Rezende, A.F.C. Silva, C.R. Maliska, The blast furnace trough two-phase flow and its influence in the refractory lining wear: mathematical modeling and numerical simulation, in: Proc. of the 19th International Congress of Mechanical Engineering, 2007. [15] M. Kou, S. Yao, S. Wu, H. Zhou, J. Xu, Effects of blast furnace main trough geometry on the Slag-Metal separation based on numerical simulation, Steel Res. Int. 90 (2) (2019) 1800383. [16] M.J. Monteiro de Oliveira, G.F.R. Rodrigues, I.A. Silva, J.J.M. Peixoto, C.A.d. Silva, Modeling of two-phase flow in blast furnace trough, Steel Res. Int. (2020) 2000485. [17] L. Wang, C.-N. Pan, W.-T. Cheng, Numerical analysis on flow behavior of Molten Iron and Slag in main trough of blast furnace during tapping process, Adv. Numer. Anal. 2017 (2017). [18] Y. Ge, M. Li, H. Wei, D. Liang, X. Wang, Y. Yu, Numerical analysis on velocity and temperature of the fluid in a blast furnace main trough, Processes 8 (2) (2020) 249. [19] S. Vázquez-Fernández, A.G.-L. Pieiga, C. Lausín-Gónzalez, P. Quintela, Mathematical modelling and numerical simulation of the heat transfer in a trough of a blast furnace, Int. J. Therm. Sci. 137 (2019) 365–374. International Journal of Thermal Sciences 173 (2022) 107349 17 P. Barral et al. [20] P. Barral, B. Nicolás, L.J. Pérez-Pérez, P. Quintela, Numerical simulation of wearrelated problems in a blast furnace runner, in: Recent Advances in Differential Equations and Applications, in: SEMA SIMAI Springer Series 18, Springer, 2019, pp. 229–244. [21] W. Davenport, I. Cameron, M. Sukhram, K. Lefebvre, Blast Furnace Ironmaking: Analysis, Control and Optimization, Elsevier, 2019. [22] M. Geerdes, R. Chaigneau, I. Kurunov, Modern Blast Furnace Ironmaking: An Introduction, IOS Press, 2015. [23] P. Geyer, Z. Halifa, Blast furnace tapping practice at ArcelorMittal South Africa, Vanderbijlpark works, in: Proc. of the Southern African Institute of Mining and Metallurgy Furnace Tapping Conf., 2014, pp. 97–112. [24] Y.A. Çengel, A.J. Ghajar, Heat and Mass Transfer: Fundamentals & Applications, McGraw Hill Education, 2015. [25] M.F. Modest, Radiative Heat Transfer, Academic Press, 2013. [26] T. Tiihonen, Stefan–Boltzmann radiation on non-convex surfaces, Math. Methods Appl. Sci. 20 (1) (1997) 47–57. [27] A. Bermúdez, D. Gómez, M.C. Muñiz, R. Vázquez, A thermo-electrical problem with a nonlocal radiation boundary condition, Math. Comput. Modelling 53 (1–2) (2011) 63–80. [28] P.-E. Druet, Weak solutions to a stationary heat equation with nonlocal radiation boundary condition and right-hand side in 𝐿𝑝(𝑝 > 1), Math. Methods Appl. Sci. 32 (2) (2009) 135–166. [29] P.-E. Druet, Existence of weak solutions to the time-dependent MHD equations coupled to the heat equation with nonlocal radiation boundary conditions, Nonlinear Anal. RWA 10 (5) (2009) 2914–2936. [30] Q. Gao, Y. Pang, Q. Sun, D. Liu, Z. Zhang, Numerical analysis of the heat transfer of radiant tubes and the slab heating characteristics in an industrial heat treatment furnace with pulse combustion, Int. J. Therm. Sci. 161 (2021) 106757. [31] C. Balaji, Essentials of Radiation Heat Transfer, Springer, 2021. [32] A. Logg, K.-A. Mardal, G. Wells, Automated Solution of Differential Equations By the Finite Element Method: The FEniCS Book, Vol. 84, Springer Science & Business Media, 2012. [33] M.S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M.E. Rognes, G.N. Wells, The FEniCS project version 1.5, Arch. Numer. Softw. 3 (100) (2015) 9–23. [34] G. Söderlind, Digital filters in adaptive time-stepping, ACM Trans. Math. Software 29 (2003) 1–26. [35] A. Kværnø, Singly diagonally implicit Runge–Kutta methods with an explicit first stage, BIT 44 (3) (2004) 489–502. [36] K.F.C. Yiu, Y. Liu, K.L. Teo, A hybrid descent method for global optimization, J. Global Optim. 28 (2) (2004) 229–238. [37] L. Ingber, Simulated annealing: Practice versus theory, Math. Comput. Modelling 18 (11) (1993) 29–57. [38] Y.-J. Wang, J.-S. Zhang, An efficient algorithm for large scale global optimization of continuous functions, J. Comput. Appl. Math. 206 (2) (2007) 1015–1026. [39] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.