Int. J. Anal. Appl. (2022), 20:7 Differential Equations Models and Their Applications in Metallurgy Ridha Selmi1,3,4,∗, Muflih Alhazmi1, Abir Sboui2,4, Hanen Louati1, Amel Touati1, Hechmi Hattab5,6 1Department of Mathematics, College of Sciences, Northern Border University, P.O. Box 1321, Arar, 73222, KSA 2Department of Mathematics, Faculty of Sciences and Art (TURAIF), Northern Border University, KSA 3Department of Mathematics, Faculty of Sciences of Gabes, 6072, Gabès, Tunisia 4Laboratory of partial differential equations and applications (LR03ES04), Faculty of sciences of Tunis, University of Tunis El Manar, 1068 Tunis, Tunisia 5FS Sfax, Univercité de Sfax, Route de la Soukra km 4 - Sfax - 3038, Tunisia 6ISIM Gabès, Université de Gabès, Cité Erriadh, Zrig, 6072, Gabès, Tunisia ∗Corresponding author: Ridha.selmi@nbu.edu.sa, ridhaselmiridhaselmi@gmail.com Abstract. This paper focus on the heat recovery from the metallurgical and mining wastes. We propose and study a new and more realistic mathematical model for heat recovery from molten slag. Our model is based on time delay differential equations. In the theoretical part, we prove that a unique solution exists to the mathematical problem. In the numerical part, we establish an algorithm based on explicit fourth order Runge-Kutta method with delay; the new feature is that the delay must be larger enough than the step of integration. Compared to the classical model (without time delay), the numerical test proves that our model is more efficient and industrially more profitable. 1. Introduction In metallurgy, heat recovery from molten slag is nowadays one of the relevant manners to valorize the huge amount of thermal energy provided by such mining waste. This very high temperature liquid Received: Oct. 7, 2021. 2010 Mathematics Subject Classification. 34K05; 65L03; 65L20; 97M50. Key words and phrases. heat transfer; heat recovery from molten slag; delayed convection; mathematical model; time delay-differential equations; Runge-Kutta method; numerical solution. https://doi.org/10.28924/2291-8639-20-2022-7 ISSN: 2291-8639 © 2022 the author(s). https://doi.org/10.28924/2291-8639-20-2022-7 2 Int. J. Anal. Appl. (2022), 20:7 slag is refrigerated then solidified by quenching procedure via fluids (gas or liquid), in general. After that, the heat recuperated is used for many purposes that range from industrially to domestically use. Such huge thermal energy can be used for producing electrical energy, or heating buildings, for example. A review paper about this subject is [13]. The report [11] deals with this matter in Australia. The paper [10] presents a global study of the problem, with focus on South Africa mining and metallurgical industry. In the literature, from mathematical point of view, models dealing with heat recovery from molten slag are formulated using ordinary differential equations that are taking in account only instantaneous properties of the molten slag. Thus, they naturally consider the ordinary time derivative. Such classical models were given by the ordinary differential equation of the form: dT dt (t) = f (t,T (t)), where T is the temperature of the molten slag, t is the time variable and f is a given function that describes the process of energy transfer. From physical point of view, as stated in [8], models above are based on the fact that convection and radiation are mainly the two manners of heat dissipation, when liquid slag is quenched. The heat transfer between the quenching fluid and the very hot slag is governed by the energy equality cV ρ dT dt (t) = Ah(T (t) −Tf ) + Aεσ0(T4(t) −T4f ), (1.1) where T is the unknown temperature of the molten slag, it depends only on the time t. The parameters V and A are respectively the volume and the surface of the molten slag; Tf is the fixed temperature of the incoming quenching fluid; c and ε state for the heat of the alloy and its integrated radiant remittance, ρ is the melt density, h is the heat transfer coefficient of the interface and finally σ0 refers to the blackbody radiation coefficient. However, we think that models above are neither realistic nor so much efficient. In fact, from physical and industrial point of view, there are two intersecting issues (a property and a phenomenon) that should be considered: • Memory effects caused by the viscoelasticity of the molten slag: As for such property, based on physical experiences and chemical analysis the detailed review paper [9] and references therein putted in evidence the viscoelastic response of molten slag, where according to the authors "at high temperatures (which is our case), the time dependent constitutive relations are needed not only for the stress-stress relationship, but also for the heat flux vector"(page 26). Also, they reviewed and discussed "the various existing implicit constitutive models for non-linear viscoelastic materials that can be used to model the rheological characteristics of slag" (page 28). Int. J. Anal. Appl. (2022), 20:7 3 • Time delay phenomena occurring in the interconnections of different parts of a system: As for such phenomena, authors in [6] used a space-averaging technique and the method of characteristics to propose a time-delay system modelling the flow temperatures of a heat exchanger; they believe that "time delay phenomena naturally occur in the interconnections of different parts of a system, as propagation of matter is not instantaneous. In particular, it occurs in tubular heat exchangers, which are very common devices in industry". Tubular heat exchangers are exactly the case of quenching systems used for molten slag. In [12], a dual-phase-lagging model of the micro-scale heat conduction is re-derived analytically from the Boltzmann transport equation. Then, based on such model, a delay-advanced partial differential equations governing the micro-scale heat conduction are established, as a more realistic model. To take in account facts signaled above, we will opt for a model governed by a time-delay equations. Mathematically speaking, this is based on delay differential equations (DDEs), where the derivative of the unknown molten slag temperature T at a time t is expressed by the values of the temperature at previous times that is dT dt (t) = f (t,T (t −τ)), (1.2) where T (t−τ) is the delay term that takes in account the memory effects of the system, and τ is the delay time. In the first chapter of the book [1], it was discussed, in details, how delayed differential equations are a reliable mathematical tool to modelize systems with memory. Based on (1.2), to make the model of heat recovery from molten slag more realistic, we propose the following new governing equation: cV ρ dT dt (t) = [Ah(T (t −τ1) −Tf ) + Aεσ(T4(t −τ2) −T4f )], (1.3) where we performed a correction via different time delays τ1 and τ2, respectively in the convective term and the radiation term. However, to simplify this model, we recall the following: By one hand, as radiation is an electromagnetic wave travelling, heat transfer by radiation occurs at the speed of light. By the other hand, as convection is a phenomenon related to the physical medium, heat transfer by convection takes place at the speed of convective medium. Thus, we are lead to consider that τ2 << τ1, because heat transfer by radiation is faster than heat transfer by convection. Mathematically speaking, we can neglect τ2 compared to τ1. Thus, model (1.3) reads cV ρ dT dt (t) = [Ah(T (t −τ1) −Tf ) + Aεσ(T4(t) −T4f )]. (1.4) Above, we mean by "faster" the way to transfer energy from one medium to another in least amount of time. 4 Int. J. Anal. Appl. (2022), 20:7 From numerical point of view, it is already known, in the literature, that numerical investigation of delayed differential equations is so delicate and some times problematic. Among others, authors in the recent publication [4] discussed several problems that arise while implementing numerical methods to solve delay-differential equations, as overlapping, difficulties with error estimation, discontinuities of solution derivatives and their detection, with focus on the Runge–Kutta methods. Let us explicitly say that, one of serious problems in such investigation is the fact that integration in time variable will be perturbed by some shift effects due to delays. In [2, 7], this anomaly was overcome by interpolation. Here, we will take the step of time integration so small compared to the delay. So, the time shift will not be in the interval of integration. This allows us to establish an explicit delayed fourth order Runge-Kutta algorithm (DFOR-K) valid for our proposed model. To the best of our knowledge, our idea is new and original. The remainder of this paper will be as follows. In section two, we will prove that our model is mathematically well posed; that is a unique solution exists to the corresponding Cauchy problem. Section three deals with the numerical investigation, we will establish an explicit forth order Runge- Kutta algorithmi, after which numerical test will be presented to validate the model. This paper will be achieved by a closely conclusion. 2. Functional setting and theoretical results In this section, we will prove that a solution to our model exists and it is unique. Following mathematical material in [5], we recall the following. • The naturel functional setting used to study retarded differential equations is the set of con- tinuous functions mapping the interval [−r, 0] into Rn, r > 0, denoted by C = C([−r, 0],Rn), and endowed with the usual norm defined by |φ| = sup −r≤θ≤0 |φ(θ)|, for any given function φ ∈ C. • If we consider σ ∈ R, A ≥ 0 and x ∈ C([σ − r,σ + A],Rn), then for any t ∈ [σ − r,σ + A], we let xt ∈ C be defined by xt(θ) = x(t + θ), −r ≤ θ ≤ 0. • If D is a subset of R×C, f : D → Rn is a given function and "." represents the derivative, then we say that the relation (DDE) ẋ(t) = f (t,xt) is a delayed (retarded or advanced) differential equation (also said functional differential equa- tion) on D and will denote this equation by (DDE). Equation (DDE) is a very general type Int. J. Anal. Appl. (2022), 20:7 5 of equation and includes ordinary differential equations ẋ(t) = f (t,x(t)) whenever r = 0 (and thus θ = 0); differential equations in the form ẋ(t) = f (t,x(t),x(t −τ)), 0 ≤ τ ≤ r, for example which is our case, in this paper. • A delayed initial value problem is simply a delayed differential equation supplemented by an initial value φ at σ: (DIV P ) { ẋ(t) = f (t,xt) xσ = φ, where σ ∈R, φ ∈ C are given. • A function x is said to be a solution of a delayed differential equation on [σ−r,σ + A) if there are σ ∈ R and A > 0, such that x ∈ C([σ − r,σ + A),Rn), (t,xt) ∈ D and x(t) satisfies equation (DDE) for t ∈ [σ,σ + A). It follows that x(σ,φ,f ) is a solution to the delayed initial value problem (DIV P ), that is a solution of (DDE) with initial value φ at σ or simply a solution through (σ,φ), if there is an A > 0, such that x(σ,φ,f ) is a solution of equation (DDE) on [σ − r,σ + A) and xσ(σ,φ,f ) = φ. For interested readers, [5] and [1] are complete references about this topic. Using notation above, we introduce the existence result (Theorem 2.3 page 42), in [5]: Theorem 2.1. Suppose Ω is an open subset in R×C and f : Ω → Rn is continuous, and f (t,φ) is Lipschitzian in φ in each compact set in Ω. If (σ,φ) ∈ Ω, then there is a unique solution of equation (DDE) passing through (σ,φ). We use theorem above to deal with existence and uniqueness of solution to our model. (1) Existence of solution: the function T : t → T (t) should belong to C(J), so it should be continuous in time t on some time interval J ⊂ R. This interval will be defined in practice by the duration of the quenching process. Thus, the operator f : J ×C(J) → R, such that f (t,T (t)) = aT (t−τ) + bT4(t) + c is continuous in (t,T ) as a polynomial in the variable T, here we note that there is no explicit dependence of the operator f on the time t and such dependence is implicitly through the temperature T. Above, the coefficients a, b and c are given by the parameters of the model after standard computation. So, at this step a solution to our model (1.4) exists. (2) Uniqueness of solution: Also, as a polynomial, the operator f is Lipschitzian in T in each compact set in J × C(J). Thus, if we supplement our model by the initial value Tσ = Tin, theorem above asserts that our model has a unique solution T (σ,Tin, f ) passing through (σ,φ = Tin). Here, we note that Tin is the temperature of the molten slag before the 6 Int. J. Anal. Appl. (2022), 20:7 beginning of the quenching process and we recall that Tσ is defined by Tσ(θ) = T (σ + θ), −r ≤ θ ≤ 0, r > 0, so that Tσ ∈ C as it is a constant function of time equal to Tin everywhere. The value of σ can be taken equal to zero, without loss of generality. 3. Numerical investigation 3.1. The delayed fourth order Runge-Kutta algorithm. The classical fourth Runge-Kutta method consists in integrating a differential equation of the form dT dt (t) = f (t,T (t)), using the formula: tn+1 = tn + h, where h is the time step and Tn+1 = Tn + 1 6 (K1 + 2K2 + 2K3 + K4) , for n = 1, 2, . . . ,m− 1, (3.1) where K1 = hf (tn,Tn) K2 = hf (tn + h 2 ,Tn + K1 2 ) K3 = hf (tn + h 2 ,Tn + K2 2 ) K4 = hf (tn + h,Tn + K3) (3.2) and m is the number of integration’s steps. Let us consider the following delay differential equation: dT dt (t) = f (t,T (t −τ),T (t)) = 1 cV ρ [Ah(T (t −τ) −Tf ) + Aεσ(T4(t) −T4f )] = αT (t −τ) + βT4(t) + γ, (3.3) where α = (6 ∗ h)/(ρ ∗ c ∗ d), β = (6 ∗ ε ∗ σ)/(ρ ∗ c ∗ d) and γ = (α ∗ Tf ) + (β ∗ (T4f )), and we approximate the molten slag grains by spheres of radius d/2. We claim that the diameter d and thus the surface A can be adjustable, by industrial methods to increase or decrease the transfer interface between molten slag and the fluid. In reference ( [8]), the coefficients α, β and γ are calculated based the physical parameters of model (1.1), to obtain α = −5.208 10−6 β = −4.626 10−3 γ = 38381.568. For sake of comparison of our delayed model to the non delayed model introduced in ( [8]), we are using the same coefficients in our numerical study. Also, as in ( [8]), we take the initial temperature of the slag to be Tin = T (0) = 1773K. Thus, in our following numerical investigation, we will consider Int. J. Anal. Appl. (2022), 20:7 7 the delayed initial value problem (S) { dT dt (t) = −5.208 10−6T (t −τ) − 4.626 10−3T4(t) + 38381.568 Tin = 1773K. If the delay τ is larger enough than the step of integration, the above fourth Runge-Kutta method can be extended to a fourth Runge Kutta algorithm with delay by expressive constants K̃1, K̃2, K̃3 and K̃4 with the delayed form of the f (t,T (t −τ),T (t)). That is K̃1 = hf (T (tn)) = h [ αT (tn −τ) + βT4(tn) + γ ] ; K̃2 = hf (T (tn) + 1 2 K̃1) = h [ α ( T (tn) + 1 2 K̃1 ) (tn−τ) + β(T (tn) + 1 2 K̃1) 4 + γ ] = h [ α ( (T (tn −τ) + 12(K̃1)(tn−τ) ) + β(T (tn) + 1 2 K̃1) 4 + γ ] , (3.4) where (K̃1)(tn−τ) = hf (T (tn −τ)) = h [ αT (tn − 2τ) + βT4(tn −τ) + γ ] ; (3.5) K̃3 = hf (T (tn) + 1 2 K̃2) = h [ α ( T (tn) + 1 2 K̃2 ) (tn−τ) + β ( T (tn) + 1 2 K̃2 )4 + γ ] = h [ α ( (T (tn −τ) + 12(K̃2)(tn−τ) ) + β ( T (tn) + 1 2 K̃2 )4 + γ ] , (3.6) where (K̃2)(tn−τ) = hf ( T (tn −τ) + 12(K̃1)(tn−τ) ) = h [ α ( T (tn − 2τ) + 1 2 (K̃1)(tn−2τ) ) + β ( T (tn −τ) + 1 2 (K̃1)(tn−τ) )4 + γ ] (K̃1)(tn−2τ) = hf (T (tn − 2τ)) = h [ αT (tn − 3τ) + βT4(tn − 2τ) + γ ] (3.7) and K̃4 = hf (T (tn) + K̃3) = h [ α ( T (tn) + K̃3 ) (tn−τ) + β ( T (tn) + K̃3 )4 + γ ] = h [ α ( (T (tn −τ) + (K̃3)(tn−τ) ] + β ( T (tn) + K̃3 )4 + γ ] , (3.8) 8 Int. J. Anal. Appl. (2022), 20:7 where (K̃3)(tn−τ) = hf ( T (tn −τ) + 12(K̃2)(tn−τ) ) = h [ α ( T (tn − 2τ) + 1 2 (K̃2)(tn−2τ) ) + β ( T (tn −τ) + 1 2 (K̃2)(tn−τ) )4 + γ ] (K̃2)(tn−2τ) = hf ( T (tn − 2τ) + 12(K̃1)(tn−2τ) ) = h [ α ( T (tn − 3τ) + 1 2 (K̃1)(tn−3τ) ) + β ( T (tn − 2τ) + 1 2 (K̃1)(tn−2τ) )4 + γ ] (K̃1)(tn−3τ) = hf (T (tn − 3τ)) = h [ αT (tn − 4τ) + βT4(tn − 3τ) + γ ] . (3.9) Then, Tn+1 = Tn + 1 6 ( K̃1 + 2K̃2 + 2K̃3 + K̃4 ) , for n = 1, 2, . . . ,m− 1, (3.10) where tn+1 = tn + h, h is the time step and m is the number of integration’s steps. We will name the algorithm (3.4)-(3.10) the delayed fourth order Runge-Kutta algorithm that we will denote by (DFOR-K). From the above computation, it is clear that since equation (3.3) is with one constant delay τ, the (DFOR-K) depends on the delay of time equals to 4τ via (K̃1)(tn−3τ), in the last line of (3.9). This proves the following theorem. Theorem 3.1. Let τ be the delay in equation (3.3) and h be the integration step in the (DFOR-K) (3.4)-(3.10). If τ ≥ h, then the (DFOR-K) (3.4)-(3.10) is a valid algorithm. Moreover, it depends on delays τ, 2τ, 3τ and 4τ. 3.2. Numerical results. We present two numerical results. The first concerns the comparison be- tween the model with delay, that it is physically more realistic, and the one without delay. The second result consists in establishing the variation of the heat transfer as a function of the exchange surface between the molten slag and the quenching fluid. In Figure 1 and Table 1, we fix the thermal exchange surface of the molten slag with the fluid used for the recovery of the heat to be A = 200 units of the sphere of diameter d modeling the grains of the molten slag and we present the decrease in the temperature of the molten slag as a function of time, with an integration step Timestep = 0.5 s, a delay τ = 4∗ Timestep, where the number of steps of integration is m = 40. So, the simulation is in a time equal to 20 s. Clearly, we can see that, when taking into account the delay, the decrease in temperature is more important, which implies a saving of the time while transferring heat from the molten slag. Also, we see that the behaviour of the temperature as a function of time presents some asymptote that seems to be the same for the two models (without and with delay). This leads to think about the optimal time to spend in the heat recovery process, in industry. Int. J. Anal. Appl. (2022), 20:7 9 Table 1. Table for the delay τ = 4 Pastemps = 4 steps of integration. pastemps=0.5 s, A =200 sphere unities, number of integration’s steps m = 40. Iteration Times (s) Delayed temperature Not delayed temperature 0 0.0000000000000000 1773.0000000000000 1773.0000000000000 1 0.50000000000000000 1496.3627665694594 1487.6130483245572 2 1.0000000000000000 1228.2200943503881 1260.5919326885185 3 1.5000000000000000 1011.7376907722477 1078.3396255965374 4 2.0000000000000000 843.90934123125589 931.22909171791252 5 2.5000000000000000 715.08690804465118 812.09212122248778 6 3.0000000000000000 616.41586854571949 715.41277990356048 7 3.5000000000000000 540.85202877542747 636.85701598717617 8 4.0000000000000000 482.97008411070431 572.97465271184217 9 4.5000000000000000 438.62089061602052 520.99665252349666 10 5.0000000000000000 404.63356920917545 478.68924422998879 11 5.5000000000000000 378.58332275097030 444.24459090912433 12 6.0000000000000000 358.61453884968427 416.19650804817417 13 6.5000000000000000 343.30637504591556 393.35426511243531 14 7.0000000000000000 331.57045335400164 374.74995043192519 15 7.5000000000000000 322.57283358988946 359.59628591699783 16 8.0000000000000000 315.67440783299958 347.25264075992033 17 8.5000000000000000 310.38531550573572 337.19755665628946 18 9.0000000000000000 306.33005421868648 329.00648616513820 19 9.5000000000000000 303.22076278047979 322.33372710777883 20 10.000000000000000 300.83675480171695 316.89774647706361 21 10.500000000000000 299.00883657838762 312.46924901901673 22 11.000000000000000 297.60728872134200 308.86147204479352 23 11.500000000000000 296.53265440274487 305.92228806723062 24 12.000000000000000 295.70867813729103 303.52777670666762 25 12.500000000000000 295.07689267128490 301.57699142269786 26 13.000000000000000 294.59246909462632 299.98769831306481 27 13.500000000000000 294.22103526397939 298.69290600508037 28 14.000000000000000 293.93623652689627 297.63803951277356 29 14.500000000000000 293.71786551843422 296.77863839112746 30 15.000000000000000 293.55042824457450 296.07848181653878 10 Int. J. Anal. Appl. (2022), 20:7 31 15.500000000000000 293.42204466040243 295.50806134363791 32 16.000000000000000 293.32360570627554 295.04333682354223 33 16.500000000000000 293.24812697431037 294.66472295529155 34 17.000000000000000 293.19025313632625 294.35626369636543 35 17.500000000000000 293.14587796555929 294.10495969758477 36 18.000000000000000 293.11185298858754 293.90022039122590 37 18.500000000000000 293.08576409385415 293.73341762396927 38 19.000000000000000 293.06576024566658 293.59752201192765 39 19.500000000000000 293.05042215001640 293.48680668521655 40 20.000000000000000 293.03866155349442 293.39660593217508 Our second numerical test focus on the effect of the variation of the contact surface between the quenching fluid and the molten slag surface on the heat transfer rate from this molten slag. It is clear that increasing this exchange surface has as an effect to increase the speed of the heat recovery from the molten slag, in both cases of modeling with delay (Figure 2) and without delay (Figure 3). Also, it is clear that while passing from surface A = 100 to surface A = 200, temperature decreases faster than while passing from A = 200 to A = 300 an so on, for the same abscise of time, see for example for t = 5s of (Figure 2). This leads to think about the optimally profitable surface that should be considered in industrial protocols. Moreover, the behavior of the temperature as a function of time presents some horizontal asymptote that should define the stopping time of the heat recovery procedure, in practice. Combining the above two statements, as a result, even more heat exchange speed is gained one takes into account the delay and increases the surface of the contact. Remark 3.1. The fact when passing from surface A = 100 to surface A = 200, temperature decreases faster than when passing from A = 200 to A = 300 an so on, may be du to the hypothesis we made in the beginning that is τ2 << τ1. In a forthcoming paper, we will deal with the model (1.3) and discuss the situation based on this model that we think to be a compete one. Int. J. Anal. Appl. (2022), 20:7 11 Figure 1. Decrease in the temperature, τ = 4∗ timestep, timestep= 0.5s and m = 40. Figure 2. Decrease in the temperature for different surface of heat transfer: A=100, 200, 300, & 400; τ = 4∗ timestep, timestep= 0.5s and m = 40. Figure 3. Decrease in the temperature for different surface of heat transfer: A=100, 200, 300, & 400, without delay, timestep= 0.5s, m = 40. 12 Int. J. Anal. Appl. (2022), 20:7 4. Conclusion Heat recovery and re-use of the huge amount of thermal energy offered by the very high temper- ature of the metallurgical and mining waste, in molten slag allow to built a circular economic model and prevent environmental damages. Such industrial procedure should be an efficient process both in the sense of duration of heat recovery and the optimal logistic protocol. In this paper, we proposed a time delayed model that is more realistic then the instantaneous one. Our model is mathematically well posed and numerically solvable using a variant of Runge-Kutta method which is popular. First, we proved that our model is more efficient compared to the classical one, in the sense that the heat recovery is faster; this saves time in industry. Also, we proved that there should exist both a stopping time to the quenching process and an optimal surface of the contact between the molten slag and the quenching fluid, so that the heat recovery procedure is maximally profitable, in industry. Acknowledgement: The authors extend their appreciation to the Deputyship for Research & Innova- tion, Ministry of Education in Saudi Arabia for funding this research work through the project number ”3400_2020_IF ”. Conflicts of Interest: The author(s) declare that there are no conflicts of interest regarding the publication of this paper. References [1] O. Arino, M.L. Hbid, E. Ait Dads, Delay differential equations and applications, Springer, Netherlands, (2006). [2] N. A. Ayunni Sabri, M. bin Mamat, Solving delay differential equations (DDEs) using Nakashima’s 2 stages 4th order Pseudo-Runge-Kutta Method, World Appl. Sci. J. 21 (2013), 181-186. [3] L. Cheng, M. Xu, L. Wang, From Boltzmann transport equation to single-phase-lagging heat conduction, Int. J. Heat Mass Transfer. 51 (2008), 6018–6023. https://doi.org/10.1016/j.ijheatmasstransfer.2008.04.004. [4] A. S. Eremin, A. R. Humphries, A. A. Lobaskin, Some issues with the numerical treatment of delay differential equations, AIP Conf. Proc. 2293 (2020), 100003. https://doi.org/10.1063/5.0027149. [5] J. K. Hale, S. M. Verduyn Lunel, Introduction to functional differential equations, Springer Verlag, Berlin, 1993. [6] S. Hamze, E. Witrant, D. Bresch-Pietri, C. Fauvel, Estimating heat-transport and time-delays in a heat ex- changer, in: 2018 IEEE Conference on Control Technology and Applications (CCTA), IEEE, Copenhagen, 2018: pp. 1514–1519. https://doi.org/10.1109/CCTA.2018.8511359. [7] F. Ismail, R. A. Al-Khasawneh, A. S. Lwin, M. B. Suleiman, Numerical treatment of delay differential equations by Runge–Kutta method using Hermite interpolation, MATEMATIKA: Malaysian J. Ind. Appl. Math. 18 (2002), 79-90 [8] C. Liu, H. Wu, J. Chang, Research on a class of ordinary differential equations and application in metallurgy, in: R. Zhu, Y. Zhang, B. Liu, C. Liu (Eds.), Information Computing and Applications, Springer Berlin Heidelberg, Berlin, Heidelberg, 2010: pp. 391–397. https://doi.org/10.1007/978-3-642-16339-5_52. [9] M. Massoudi, P. Wang, A brief review of viscosity models for slag in coal gasification, DOE/NETL-2012/1533, National Energy Technology Laboratory, Pittsburgh, PA, 2011. [10] E. Matinde, G.S. Simate, S. Ndlovu, Mining and metallurgical wastes: a review of recycling and re-use practices, J. South. Afr. Inst. Min. Metall. 118 (2018), 825-844. https://doi.org/10.17159/2411-9717/2018/v118n8a5. https://doi.org/10.1016/j.ijheatmasstransfer.2008.04.004 https://doi.org/10.1063/5.0027149 https://doi.org/10.1109/CCTA.2018.8511359 https://doi.org/10.1007/978-3-642-16339-5_52 https://doi.org/10.17159/2411-9717/2018/v118n8a5 Int. J. Anal. Appl. (2022), 20:7 13 [11] D. Xie, Y. Pan, R. Flann, B. Washington, S. Sanetsis, J. Donnelley et al., Heat recovery from slag through dry granulation, in: 1st CSRP Annual Conference. Melbourne (Australia), vol. CSIRO Minerals, pp. 29–30, 2007. [12] M. Xu, L. Wang, Dual-phase-lagging heat conduction based on Boltzmann transport equation, Int. J. Heat Mass Transfer 48 (2005), 5616–5624. https://doi.org/10.1016/j.ijheatmasstransfer.2005.05.040. [13] H. Zhang, H. Wang, X. Zhu, Y.-J. Qiu, K. Li, R. Chen and Q. Liao, A review of waste heat recovery technologies towards molten slag in steel industry, Appl. Energy 112 (2013), 956-966. https://doi.org/10.1016/j.apenergy. 2013.02.019. https://doi.org/10.1016/j.ijheatmasstransfer.2005.05.040 https://doi.org/10.1016/j.apenergy.2013.02.019 https://doi.org/10.1016/j.apenergy.2013.02.019 1. Introduction 2. Functional setting and theoretical results 3. Numerical investigation 3.1. The delayed fourth order Runge-Kutta algorithm 3.2. Numerical results 4. Conclusion References