Numerical simulation of fluid flow and heat transfer in

Applied Thermal Engineering 117 (2017) 263–274

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

Numerical simulation of fluid flow and heat transfer in an industrial continuous furnace S. Defaee Rad a,⇑, A. Ashrafizadeh a, M. Nickaeen b a b

Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran Richard D. Berlin Center for Cell Analysis and Modeling, University of Connecticut Health Center, Farmington, CT 06030, USA

h i g h l i g h t s A 3D, transient numerical study of an industrial continuous furnace is proposed. Geometry of furnace load is simplified by a thermally equivalent model. Load movement is avoided by implementing unsteady boundary conditions. Fairly accurate results at a reasonably low computational cost are obtained.

a r t i c l e

i n f o

Article history: Received 30 August 2016 Revised 18 January 2017 Accepted 8 February 2017 Available online 10 February 2017 Keywords: Welding electrode Continuous furnace 3D numerical simulation Load’s temperature history Conjugate heat transfer

a b s t r a c t We present a three-dimensional computational study of the transient heat transfer and turbulent fluid flow inside an arc-welding electrode continuous furnace. The model is implemented in FLUENT, a finite volume commercial code. Large difference in the length scales of the electrodes and the furnace, movement of the electrodes, and existence of various modes of heat transfer are the major factors influencing the accuracy and efficiency of the simulation. Two modeling strategies are used to overcome these difficulties. First, the electrode geometry and material composition are simplified using a thermally equivalent model. Second, the electrode movement inside the furnace is avoided by implementing timedependent boundary conditions applied to a fixed domain. A fairly close agreement is obtained in comparing the load’s temperature history against the experimental data with an absolute relative difference below 2.7%. The space between the electrode trays in the furnace is very limited. The analysis shows that the input air does not circulate effectively between the trays. This lowers the thermal efficiency of the furnace and leads to uneven treatment of the electrodes. The present study provides guidelines to improve future furnace designs. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Arc welding electrodes play an important role in the construction of mechanical structures in various industries. Before the electrode is processed into its final state, it is usually heated in a furnace. The electrode is moved along the length of the furnace and, depending on the application, is thermally treated by various heat transfer mechanisms. This heat treatment is necessary to assure the quality of the welding and to ensure that all mechanical, metallurgical, environmental and economical measures are satisfied. The present paper aims at the performance evaluation of an industrial convective furnace, used to cure arc-welding electrodes. ⇑ Corresponding author. E-mail addresses: [email protected] (S. Defaee Rad), ashrafizadeh@kntu. ac.ir (A. Ashrafizadeh), [email protected] (M. Nickaeen). http://dx.doi.org/10.1016/j.applthermaleng.2017.02.031 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.

Five fans are employed to direct hot air, processed in a shell and tube heat exchanger, into the furnace. Similar to many other applications of continuous furnaces [1–3], the load’s temperature history has considerable importance in this case as well. However, because of the transient, multidimensional and turbulent flow conditions as well as existence of simultaneously different modes of heat transfer, it is hard to carry out an accurate, quick and inexpensive analysis. It has been shown that the preceding goals are achievable if Computational Fluid Dynamics (CFD) is employed [4,5]. Computational Fluid Dynamics provides an efficient tool to obtain the transient flow and temperature fields, to inspect the effects of thermo-physical and geometrical parameters on the curing process and to evaluate the success of various optimization methods. Furthermore, commercial CFD packages, such as FLUENT, CFX, STAR-CD or PHOENICS have increased the computational competency significantly and have

264

S. Defaee Rad et al. / Applied Thermal Engineering 117 (2017) 263–274

Nomenclature c1e, c2e, cl constants in the k–e model cp specific heat capacity at constant pressure h convection heat transfer coefficient j thermal conductivity k turbulent kinetic energy n outward unit normal vector p pressure t time T Temperature ui velocity vector component (i = 1, 2, 3) v velocity magnitude xi Cartesian coordinate component (i = 1, 2, 3) e dissipation rate of turbulent kinetic energy q density

emerged as the dominant tools to perform quick numerical analysis with reasonable accuracy. A detailed literature study shows that CFD tools have not yet been employed to analyze the transport phenomena inside an arc-welding electrode continuous furnace. However, reviewing similar applications of CFD for the study of continuous furnaces in steel [6], food [7] and ceramic [8] industries provides useful information for the analysis of electrode curing process. Almost all of the mentioned applications are followed by a standard methodology. First, the physical processes are described as mathematical models. Second, the governing equations are discretized and solved to obtain a numerical solution of the problem. Sensitivity analysis may then be conducted to evaluate the impacts of design and performance parameters on the simulation results. Here, use of Computational Fluid Dynamics with focus on commercial software in the numerical simulation of transport phenomena in steel reheat furnaces is briefly reviewed. Prieler et al. [9] performed CFD modeling of an industrial walking hearth furnace to predict transient heating characteristics of the billets. Emadi et al. [10] obtained both convective and radiative heat fluxes using zone method and investigated the effects of design parameters on billet temperature behavior. Han and Chang [11] conducted a 3D unsteady numerical simulation of a reheat furnace with the commercial software FLUENT and analyzed the optimal slab residence time. Han et al. [12] carried out a periodically transient numerical simulation of heat transfer to slabs in a reheat furnace by using FLUENT. Due to the default software limitations a User-Defined Function (UDF) was developed to process the slab movement. Morgado et al. [13] followed a numerical strategy similar to that of Han [12] and compared two different thermochemical composition models for the slab. Dubey and Srinivasan [14] applied an unsteady three dimensional numerical model to analyze temperature field in the billet. They discretized the governing equations by finite volume method and solved the equations using an in-house MATLAB code. Kim et al. [15] presented a comprehensive FLUENT-based simulation to understand the turbulent flow and radiative heat transfer in a walking beam furnace. One of the significant implications of their work was the simplifications made in the modeling of the furnace geometry to avoid the unaffordable number of computational cells. Furthermore, other researchers conducted similar CFD studies using alternative software suits such as STAR-CD [16], Phoenics [17], FURMO [18] and CONCERT [19]. The current paper addresses a numerical approach, performed with FLUENT version 6.3 [20], to examine the non-isothermal flow phenomena inside an arc-welding electrode continuous furnace. Two important modeling strategies are employed here that are

l lt rk, re u

viscosity turbulent viscosity coefficient turbulent Prandtl number for k and e scalar

Subscripts i,j coordinate directions Superscripts – Reynolds-averaged value ’ Reynolds fluctuation value

crucial to the success of the numerical approach and considerably reduce computational requirements of the simulation. To the best knowledge of the authors, these simplifying strategies have not been used by others for modeling arc-welding electrode continuous furnaces. First, we consider time-dependent boundary conditions to avoid the vastly expensive moving mesh methods required for modeling electrode movement inside the furnace. Instead of moving electrodes forward, we solve inside a fixed geometry (fixed segment of the furnace) and model the hot air inflow through time-dependent boundary conditions that move backward with the same velocity as the electrodes. Since FLUENT does not support moving boundary conditions with its default functions, a User Defined Function [21] is developed and linked to the main solver. Second, we replace the electrode geometry (cross-section) and material composition with a computationally inexpensive but thermally equivalent model. The existence of diverse length scales between the electrodes and the furnace renders the good quality grid generation and the efficient boundary condition implementation very challenging. Therefore, we propose several simplified models for the electrodes and systematically study them under equivalent furnace operating conditions. Accuracy of the numerical results is verified by refining the computational grid to ensure that a grid independent solution is obtained, and by changing the time step used for the temporal discretization. In addition, the electrode’s temperature history is compared to the experimental data obtained from a k-type thermocouple measuring the temperature in a real furnace operating condition. A maximum relative temperature difference below 2.7% is obtained which verifies the correctness of the simulation results and the assumptions used in the mathematical model. The rest of the paper is organized as follows: The furnace and electrode configurations are described in Section 2. The mathematical and numerical models as well as the computational details are presented in Section 3. The results are discussed in Section 4 and, finally, Section 5 summarizes important conclusions.

2. Furnace configuration and arrangement of electrodes The continuous electrode curing furnace considered here consists of preheating, heating and cooling sections that are 10, 9.8 and 4.7 m long, respectively. The electrodes pass through the entire furnace in 146 min with a velocity of 0.0028 m/s. Therefore, the electrodes are treated for 59.5, 58.5 and 28 min in each of the aforementioned furnace zones, respectively. Air is forced into the furnace through a number of nozzles, exchanges heat with electrodes and then leaves from the bottom. The first two sections

265


are furnished with almost similar arrangements of nozzles along the top and side walls. Fig. 1 depicts the cross section of the preheating (or heating) zone. It should be noted that the left wall has one row of nozzles more than the right wall. The nozzle arrangement of the cooling section has two major differences. First, nozzles are located only on the right wall. Second, they are aligned in 10 rows, while in other sections the arrangement of the side wall rows is tighter. An arc-welding electrode has a central rod covered by a particular material. It has a length of 450 mm with inner and outer diameters equal to 4.95 and 7.2 mm, respectively. A series of electrodes are organized in a tray with the width of 700 mm as shown in Fig. 2 (a). Three columns, each one consisted of twenty-eight electrode trays (numbered from bottom to top), are then organized as a wagon, as shown in Fig. 2(b), and discharged to the furnace for thermal treatment. 3. Numerical procedure 3.1. Turbulent air flow and energy equations A fluid flow with velocity magnitude similar to the studied furnace is normally considered laminar [22]. However, air enters the welding electrode furnace from three intersecting paths and consequently, the flow becomes turbulent. In this paper, Reynolds averaging of the governing equations is employed to model turbulence. In the Reynolds averaging technique, the solution variables in the instantaneous Navier-Stokes equations are decomposed into the mean and fluctuating components. Velocities are decomposed as

i þ u0i ui ¼ u

ð1Þ

Fig. 2. Illustration of the arrangement of electrodes in each (a) tray and (b) wagon.

i and u0i are the instantaneous, Reynolds-averaged and where ui, u fluctuating velocity components, respectively. Similar decomposition is used for other quantities

þ /0 /¼/

ð2Þ

where / indicates a scalar variable such as the pressure and enthalpy. Substituting the above expressions into the instantaneous mass, momentum and energy equations of an incompressible fluid and taking a time average yields the following equations using the Cartesian tensor notation:

i @u ¼0 @xi

ð3Þ

i i @u @u @p @ @u @ 0 0 j i ¼ þl qui uj þu @xi @xj @xj @xj @t @xj

q

@T @T j @ @T i ¼ þu @t @xi qcp @xi @xi

!

@ 0 0 uT @xi i

ð4Þ

ð5Þ

The Reynolds stresses, qu0i u0j , are related to the mean velocity gradients under the commonly used Boussinesq hypothesis [23,24] Fig. 1. Cross section of the preheating and heating zones. Dimensions are in millimeter.

qu0i u0j ¼

i @u j @u 2 @u þ qk þ lt k dij 3 @xj @xi @xk

ð6Þ

266


Fig. 3. The exact electrode geometry and material composition (model A), and the suggested models (models B to E).

Table 1 Thermophysical properties of the electrode. Parameter

Unit

Central rod

Cover

Mixed material

Density Specific heat Conductivity

kg/m3 J/kg K W/m K

7801 473 43

3000 800 10

5269.4 645.4 25.6

This hypothesis, however, introduces the turbulence kinetic energy k and its eddy viscosity lt as new variables. Numerous turbulence models are devised to model the new variables. The standard k–e model, introduced by Launder and Spalding [25], has been used extensively to model turbulent flow and heat transfer in practical engineering applications. It is a semi-empirical model based on phenomenological principles. We use the k–e model in our simulations because of its widespread use and also its simplicity. The standard k–e model consists of two transport equations, namely the equations for turbulence kinetic energy and its dissipation rate e. These are described by the following equations

@ @ @ i kÞ ¼ ðqkÞ þ ðqu @t @xi @xj @ @ @ i eÞ ¼ ðqeÞ þ ðqu @t @xi @xj

j lt @k @u qu0i u0j qe rk @xj @xi

lþ

lt @ e e2 c2e q re @xj k

c1e

lþ

e

k

qu0i u0j

j @u @xi

ð8Þ

The turbulence eddy viscosity is defined as the model constants are as follows [25]

c1e ¼ 1:44;

c2e ¼ 1:92;

cl ¼ 0:09;

ð7Þ

lt ¼ qcl k2 =e and

rk ¼ 1; re ¼ 1:3

3.2. Electrode energy equation The electrode material is modeled as a Fourier’s solid with the energy conservation equation written as

@T j @ @T ¼ @t qcp @xi @xi

ð9Þ

3.3. CFD models 3.3.1. Electrode model The existing length scales in the furnace are very diverse. Electrodes with a few millimeters diameter are cured inside a multimeter-long furnace. Therefore, using the real geometry of the electrodes in the numerical simulation will be challenging in terms of good quality grid generation and efficient boundary condition implementation. The circular cross section and the two-material composition of the electrodes impose further complications. Therefore, it is beneficial to replace the electrode with a computationally inexpensive but thermally equivalent model. In Fig. 3, four alternative cross sections with different material compositions are suggested to simplify the electrode’s modeling. Model A represents the exact geometry and material composition. Models B and C consider two-material square cross sections, while models D and E assume one homogeneous material for the circular and square cross sections, respectively. These models are evaluated under a numerical benchmark problem to select the most appropriate one to be used in the furnace simulation. The assumptions made for this problem are mostly derived from the actual furnace operating conditions. Particularly, a single electrode is analyzed under transient convection-conduction thermal conditions. Since the length of the electrode is more than 60 times larger than its diameter, it is wise to carry out a two dimensional analysis and neglect the longitudinal dimension. The

267

S. Defaee Rad et al. / Applied Thermal Engineering 117 (2017) 263–274 Table 2 Temperature of inlet air to the furnace. Name of the zone

Preheating

Heat exchanger number Temperature (K)

1 353

Heating 2 388

3 393

Cooling 4 393

– 313

Fig. 4. The maximum temperature difference across the electrode during the heat treatment process.

convection heat transfer is modelled as a convective boundary condition applied on the electrode. This needs the free stream temperature and the convection heat transfer coefficient to be specified. Free stream temperature in each section of the furnace is considered equal to the average temperature of the inlet air in the corresponding furnace section (Table 2), i.e. 363 K. Convection heat transfer coefficient in a convective furnace is commonly estimated using the equation h = 1 + 0.225 V [26], which works with the English system of units. As explained later in Section 3.4, hot air enters the furnace with a velocity of 1 m/s (or 3.28 ft/s). Thus, the convection heat transfer coefficient is readily calculated as 10 W/m2 K (or 1.738 Btu/ h ft2 F). The initial temperature of the electrode is set to 300 K. The electrode’s thermal treatment lasts 3600 s. The thermo-physical properties of the electrode are assumed constant as reported in Table 1. The mixed material column presents volume weighted average of the values for the central rod and the cover. It is applied for models D and E. Fig. 4 presents the maximum crosswise temperature difference in the electrode (model A). Considering the fact that the maximum temperature difference across model A during the transient simulation is about 0.07 K, it can be argued that the temperature across the electrode is almost uniform. In addition, the electrode has a very low

Biot number value, i.e. Bi = 0.001. This further confirms that the solid material can be considered as a lumped material [27]. Fig. 5 shows the temperature differences between the average temperature of model A and other models. The results for model B demonstrate that changing the electrode cross sectional shape to a square has no significant effect on the results. The results for model E demonstrate that the mixed material hypothesis leads to tolerable differences as shown in the figure. Therefore, model E, which has a simple square geometry and a homogeneous material composition, is selected as the final electrode model. Furthermore, the low Biot number of the model E makes it possible to consider each tray of electrodes as a two dimensional flat plate of 4.95 mm thickness, hence having much less trouble in the mesh generation. 3.3.2. Furnace model The movements of the electrodes and the coexistence of largely different length scales are two major issues which make the numerical simulation lengthy and cumbersome, even if state-ofthe-art computers are utilized. The first issue requires moving grid algorithms to account for the movement of electrodes. Maintaining high quality grid and hence an accurate solution would typically

Fig. 5. Average temperature variation between model A and the other models.

268


Fig. 6. The modeled furnace with an electrode column.

require remeshing techniques that are time consuming and problem dependent [28]. The second issue considerably increases the number of computational cells, in the order of multi millions that are required to sufficiently resolve small geometric features. In

the current research, efficient approaches are proposed to avoid such difficulties. Furthermore, because every column of electrodes in each wagon experiences identical operational sequence, only one column plus its distance with the neighbouring columns (0.59 m) is considered and then, periodic conditions are implemented to the front and rear side boundaries. Fig. 6 displays the sketch of the furnace model. To avoid the moving grid, we assume that the wagon is stationary but furnace boundaries proceed backward with the same velocity as the electrodes. The logical procedure which is followed in the beginning of each time step is to calculate the central position of nozzles arranged on the top and side walls. Next, they are compared with the location of each computational cell. If the difference is less than the radius of nozzles (12.5 mm), normal velocity of the corresponding cell is set to the inlet air velocity. It should be taken into account that FLUENT is not able to perform this specific task with its common features. Therefore, for each furnace wall, one set of User Defined Functions is developed. User Defined Function is a programming feature that allows the user to customize FLUENT and can significantly enhance its standard capabilities. It is written in C programming language and could be dynamically linked with the FLUENT solver. As mentioned before, the furnace height and the distance between trays are incomparable. Thus, producing a high quality grid that is efficient in terms of the total number of computational cells is not straightforward. To achieve these important goals, the model is divided into several sub-regions so that each of them can be meshed separately with hexahedral structured cells. Where necessary, the procedure allowed changing the grid size to have more control on the desired fineness of the sub-region. Note that the solver ensures continuity of fluxes passing through the nonconformal mesh interface between the sub-regions [20]. Fig. 7 shows two views of the generated computational mesh of 144,000 cells. The numerical grid is generated using Gambit software [29].

Fig. 7. Computational grid. (a) The whole model and (b) the cross section.


269

Fig. 8. Temperature profiles along line AB of the 12th tray for three computational grids at (a) 30, (b) 75 and (c) 105 min.

3.4. Initial and boundary conditions The column of electrodes enters the furnace with initial temperature of 313 K (ambient temperature). At the beginning, the air temperature inside of the furnace is 353 K. Inlet air is forced through the nozzles at 1 m/s. The inlet air temperature varies from one zone to another according to Table 2. Note that there are two

inlet air temperatures in both preheating and heating zones in Table 2 because the air is supplied by two different heat exchangers. Because hot air flows on both sides of furnace walls, it is reasonable to assume that the walls are well-insulated. According to the aforementioned simulation strategy, periodic boundary condition is implemented on the two planes along the column path. Since

270


the fluid velocity or pressure at the furnace outlet is unknown prior to the simulation, it is considered as an outflow-type boundary condition [20]. The energy equation of electrode trays needs either a known temperature or heat flux value as the boundary condition. However, both of them are unknown prior to the calculation. Using principle of heat flux continuity, the energy equation in the solid and fluid phases are coupled to each other at the boundaries. Therefore, the following relation is enforced at the fluid/solid interfaces:

j

dT dT ¼ j dn Fluid dn Solid

ð10Þ

where n is the unit normal vector at the interface. 3.5. Computational details The finite-volume commercial code FLUENT 6.3 is employed to numerically investigate the transient three dimensional fluid flow and heat transfer in the welding electrode continuous furnace. Standard k–e model is used for efficient evaluation of the flow turbulence. The analysis uses coupled pressure-velocity formulation (fully implicit). Second-order upwinding and fully implicit timestepping schemes are used for spatial and temporal discretizations, respectively. The relaxation factors are 0.8 for both pressure and velocity and 1 for the energy equation. The scaled residuals are kept below 5 105 for the continuity and momentum equations and below 1 106 for the energy equation. The simulations are performed on a four-processor Pentium IV desktop (Windows 7 OS, 3 GHz, 2 GBRAM). A typical simulation takes 320 h (around two weeks) to be completed. We also included the load’s temperature history obtained from a simulation with the first-order upwinding scheme (see Fig. 9). In this simulation, the scaled residuals are kept below 1 103 for the continuity and momentum equations. All other simulation parameters are the same. 4. Results and discussion Although the time step size of the employed algorithm is not restricted by any stability constraints, it should be kept sufficiently small in order to obtain accurate results. Here, we tested the effect of time step size on the mean temperature of the 12th tray along its width (line AB in Fig. 2a) using 10 and 20-s time steps. The com-

parison revealed that the relative difference between the two solutions is always less than 0.5%. As the former needs much more computational resources, the latter is chosen for temporal discretization. The accuracy of the numerical simulation is strongly affected by the grid size. Generally, finer computational cells result in a more accurate solution. However, the significant increase in both computer memory requirement and calculation time puts strict limitations on the grid fineness especially in 3D simulations. Therefore, numerical simulations are often tested using several grid sizes to study grid independency. In this paper, three grids with 60,000, 80,000 and 144,000 control volumes are considered. The total number of control volumes is mainly controlled by changing the resolution of control volumes located in trays and the spaces between them. The temperature profile along line AB of the 12th tray is shown in Fig. 8 at (a) 10 min into the preheating process, (b) 75 and (c) 105 min into the heating process. The maximum variation between the temperatures of the finest grid and the two other grids are below 1% and 0.5%, respectively. Since the grid with 144,000 cells leads to a more accurate solution and is affordable with the available computational resources, it is used in the remainder of this work. It should be noted that the temperature profiles in Fig. 8 are not symmetric due to the fact that the nozzle arrangements at the side walls are not completely identical. Because no analytical or similar numerical solution exists for welding electrode continuous furnaces, experimental data is used to validate the simulation results. A data logger assembled with a K-type thermocouple measures the 12th tray temperature at z = 10 cm along line AB (see Fig. 2a). The numerical result is compared against the experimental data in Fig. 9. The temperature curves fairly follow each other during the whole process with an absolute relative error below 2.7% for the second-order upwind scheme. We also included the temperature curve obtained by a first-order upwind scheme. The two numerical simulation curves are very similar (the maximum temperature variation is 2.4 °C), however the maximum absolute relative error is around 3% for the first-order scheme. Please note that the first-order upwind simulation took approximately 75 h, less than a quarter of the simulation time for the second-order scheme. This suggests that the firstorder upwind scheme with its robustness and lower computational cost could be used in the first trials of such industrial simulations. The temperature curves in Fig. 9 are of bell-shaped type which is common in such furnaces [30]. At the furnace entrance, the air temperature is 313 K, i.e. the temperature of the air surrounding the furnace. Temperature increases sharply at the start of

Fig. 9. Comparison of simulated and measured temperatures at z = 10 cm along line AB of the 12th tray. Vertical dashed lines mark different zones.


Fig. 10. Air path lines inside the furnace at (a) 75 min (heating zone) and (b) 135 min (cooling zone).

preheating until approximately one third of the zone length. Then, as the furnace is influenced by the second heat exchanger, the air temperature experiences another rapid increase, coming to 380 K at the end of the preheating zone. The heating rate starts to decline in the heating zone and this condition is maintained until a temperature of around 390 K is reached. In the cooling zone, the temperature decreases approximately linearly and at last, the electrodes leave the furnace at about 376 K. Fig. 10 presents path lines, i.e. the trajectory of air particles from the nozzles to the furnace outlet, at t = 75 (heating zone) and 135

271

(cooling zone) minutes. The path lines are colored according to the wall from which they are originated. Because of the short distance between the trays, only a small proportion of the air penetrates to the inter-tray areas of the wagon. The majority of the airflow passes around the wagon and is directed to the furnace outlet, see Fig. 10a. Fig. 10b shows that the air simply enters the cooling zone from the right wall (there is no nozzle on the left wall in the cooling zone) and after passing around the wagon exits the furnace. Heat flux distributions of the electrode trays are demonstrated at t = 75 and 135 min in Fig. 11. By convention the negative sign means that heat is transferred to the surface. Note that the highest and lowest trays have the maximum heat flux values. Due to the fact that the whole top wall and a part of the side wall nozzles expel hot air nearly perpendicular to the 28th (highest) tray, it is subjected to a high heat flux rate. However, mass conservation requires the velocity of air to increase at the furnace bottom. This accounts for the increase in convection heat transfer coefficient and hence heat flux around the 1st (lowest) tray. In comparison with these trays, all the other trays show a rather similar behavior and their heat flux declines considerably as we move from the edges toward the center. This is mainly associated with the decreased air velocity magnitude between the trays (as explained in the discussion of Fig. 10). Note that in the cooling zone (Fig. 11b), heat transfer occurs in the opposite direction, i.e. from the load to the furnace, and hence the heat flux is positive. Moreover, due to the asymmetric nozzle arrangement in the cooling zone, a one-sided heat transfer is observed for almost all of the trays in this zone. Fig. 12a–c compares the temperature of the 1st, 2nd, 12th, 27th and 28th trays along line AB at t = 10 (preheating zone), 75 (heating zone) and 146 (end of the cooling zone) minutes. In the preheating zone, there is no major temperature difference between trays. However, as time advances, the differences become more distinct. Moreover, because of the increased heat transfer near the lowest and highest trays, they always have higher temperatures compared to the other trays in the preheating and heating zones. Since similar heat transfer condition governs the temperature variation of the other trays (2nd, 12th, 27th), they have similar temperature profiles with values close to one another. Leaving the heating zone, the 28th tray has the highest temperature, however it gets a high volume of cold air flow coming from five rows of nozzles above it. Thus, it is effectively cooled down. In fact, it has the lowest average temperature at the end of the cooling zone as shown in Fig. 12c. The 12th and 27th trays on the other hand have the lowest temperatures entering the cooling zone. The available cold air flow passing around them is sufficient to lower the temperature of their right sides. The 1st and 2nd trays also have relatively high temperatures when they enter the cooling zone. The first tray gets proportionally more cold air flow passing around it while the air flow around the 2nd is restricted due to the short distance between the trays. Therefore, the 2nd tray has a higher temperature as it leaves the furnace. Finally, the air velocity increases around the left side of lower trays due to the mass conservation. The increased convection heat transfer coefficient on the left side combined with the continuous flow of cold air coming from the right wall creates the local maximum in the temperature profile of the 1st and 2nd trays in Fig. 12c.

5. Conclusions A three-dimensional numerical simulation of turbulent fluid flow and heat transfer inside an arc-welding electrode continuous furnace has been carried out using FLUENT software. To the best knowledge of the authors, there is no previous numerical

272


Fig. 11. Heat flux distributions shown from two different angles at (a) 75 min (heating zone) and (b) 135 min (cooling zone). Heat flux is shown in W/m2.

simulation for this particular thermal problem despite its practical importance. Large number of electrodes, disparity between the size of electrodes and the furnace dimensions, transient and threedimensional nature of the analysis as well as the presence of diverse modes of heat transfer contributed to the complexity of the simulation. In this paper, two main simplifying assumptions were applied. First of all, the geometry and material composition of the electrode were replaced with a computationally inexpensive and thermally equivalent model. Second, the electrode movement was modelled using transient boundary conditions applied to a fixed computational domain. This was implemented by linking our own developed User Defined Functions to the main solver. The computational domain was decomposed into several subregions and was meshed separately for a better control on the local refinements. After performing grid-independency tests, a computational grid consisting of 144,000 structured cells was selected. The conservation equations of mass, momentum and energy were discretized using finite volume method. Numerical simulation was performed on a Pentium IV, 3 GHz processors with a 2 GB access

memory, on which a complete solution takes about 320 h. Some of the most important results are summarized here. 1. There is a fairly close agreement between the numerical and experimental temperature history of the load. The maximum relative temperature difference is always less than 2.7% for the second-order upwinding scheme. This validates the approach employed in the simulation. Moreover, our comparison shows that the first-order upwinding scheme could be used in the first trials of similar industrial simulations to obtain an acceptable preliminary solution with a relatively lower computational cost. 2. Because electrode trays are arranged at very short distances from each other, only a small proportion of the input air penetrates between them. The remaining air flows around the trays and leaves the furnace before it could effectively exchange heat with the electrodes. This leads to a significant decline in the heat transfer efficiency. Therefore, it is necessary to re-design the arrangement of trays and improve the thermal performance of the furnace.


273

Fig. 12. Temperature curves along line AB of five trays at (a) 10 (preheating zone), (b) 75 (heating zone) and (c) 146 (end of the cooling zone) minutes.

3. The maximum heat flux occurs at the lowest and the highest trays. Consequently, these trays experience the highest temperatures in the preheating and heating zones. Intermediate trays in the column have considerably smaller heat fluxes and are cured under relatively identical conditions.

4. The proposed numerical approach and the underlying assumptions deliver fairly accurate results at a reasonably low computational cost. Most of the strategies employed here are general and can be adopted in similar industrial applications of continuous furnaces to ensure a near optimal performance is achieved or to go even further and effectively redesign future furnaces.

274


Funding sources This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. References [1] A.G. Fedorov, K.H. Lee, R. Viskanta, Inverse optimal design of the radiant heating in materials processing and manufacturing, J. Mater. Eng. Perform. 7 (1998) 719–726. [2] M.-J. Huang, C.-T. Hsieh, S.-T. Lee, C.-H. Wang, A coupled numerical study of slab temperature and gas temperature in the walking-beam-type slab reheating furnace, Numer. Heat Transf., Part A 54 (2008) 625–646. [3] N. Therdthai, W. Zhou, T. Adamczak, Two-dimensional CFD modelling and simulation of an industrial continuous bread baking oven, J. Food Eng. 60 (2003) 211–217. [4] S.-Y. Wong, W. Zhou, J. Hua, CFD modeling of an industrial continuous breadbaking process involving U-movement, J. Food Eng. 78 (2007) 888–896. [5] H.K. Versteeg, W. Malalasekera, Introduction to Computational Fluid Dynamics, first ed., Longman Group, Harlow, 1995. [6] V.K. Singh, P. Talukdar, P.J. Coelho, Performance evaluation of two heat transfer models of a walking beam type reheat furnace, Heat Transf. Eng. 36 (2015) 91– 101. [7] S.-Y. Wong, W. Zhou, J. Hua, Designing process controller for a continuous bread baking process based on CFD modeling, J. Food Eng. 81 (2007) 523–534. [8] R. Oba, T.S. Possamai, V.P. Nicolau, Thermal analysis of a tunnel kiln used to produce roof tiles, Appl. Therm. Eng. 63 (2014) 59–65. [9] R. Prieler, B. Mayr, M. Demuth, B. Holleis, C. Hochenauer, Prediction of the heating characteristic of billets in a walking hearth type reheating furnace using CFD, Int. J. Heat Mass Transf. 92 (2016) 675–688. [10] A. Emadi, A. Saboonchi, M. Taheri, S. Hassanpour, Heating characteristics of billet in a walking hearth type reheating furnace, Appl. Therm. Eng. 63 (2014) 396–405. [11] S.H. Han, D. Chang, Optimum residence time analysis for a walking beam type reheating furnace, Int. J. Heat Mass Transf. 55 (2012) 4079–4087. [12] S.H. Han, D. Chang, C.Y. Kim, A numerical analysis of slab heating characteristics in a walking beam type reheating furnace, Int. J. Heat Mass Transf. 53 (2010) 3855–3861.

[13] T. Morgado, P.J. Coelho, P. Talukdar, Assessment of uniform temperature assumption in zoning on the numerical simulation of a walking beam reheating furnace, Appl. Therm. Eng. 76 (2015) 496–508. [14] S.K. Dubey, P. Srinivasan, Development of three dimensional transient numerical heat conduction model with growth of oxide scale for steel billet reheat simulation, Int. J. Therm. Sci. 84 (2014) 214–227. [15] J.G. Kim, K.Y. Huh, I.T. Kim, Three-dimensional analysis of the walking-beamtype slab reheating furnace in hot strip mills, Numer. Heat Transf., Part A 38 (2000) 589–609. [16] C.-T. Hsieh, M.-J. Huang, S.-T. Lee, C.-H. Wang, Numerical modeling of a walking-beam-type slab reheating furnace, Numer. Heat Transf., Part A 53 (2008) 966–981. [17] Y. Tang, J. Laine, T. Fabritius, J. Harkki, The modeling of the gas flow and its influence on the scale accumulation in the steel slab pusher-type reheating furnace, ISIJ Int. 43 (2003) 1333–1341. [18] W.A. Fiveland, R.A. Wessel, FURMO: a numerical model for predicting performance of three-dimensional pulverized-fuel fired furnaces, AIAA/ASME Thermodynamics and Heat Transfer Conference, Boston, MA, June 2–4, 86-HT-35, 1986. [19] S.M. Correa, W. Shyy, Computational models and methods for continuous gaseous turbulent combustion, Prog. Energy Combust. Sci. 13 (1987) 249–292. [20] FLUENT 6.3 User’s Guide, FLUENT Inc, Lebanon, NH, 2006. [21] FLUENT 6.3 UDF Manual, FLUENT Inc, Lebanon, NH, 2006. [22] F.M. White, Fluid Mechanics, seventh ed., McGraw-Hill, New York, 2011. [23] J.O. Hinze, Turbulence, second ed., McGraw-Hill, New York, 1975. [24] D.C. Wilcox, Turbulence Modeling for CFD, third ed., DCW Industries, California, 2006. [25] B.E. Launder, D.B. Spalding, Lectures in Mathematical Models of Turbulence, first ed., Academic Press, London, 1972. [26] F. Kreith, D.Y. Goswami, The CRC Handbook of Mechanical Engineering, second ed., CRC Press, New York, 2005. [27] F.P. Incropera, D.P. DeWitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat and Mass Transfer, sixth ed., John Wiley & Sons, New Jersey, 2007. [28] L. Zhang, X. Chang, X. Duan, Z. Zhao, X. He, Applications of dynamic hybrid grid method for three-dimensional moving/deforming boundary problems, Comput. Fluids 62 (2012) 45–63. [29] GAMBIT 2.3 User’s Guide, FLUENT Inc, Lebanon, NH, 2006. [30] P.S. Mirade, J.D. Daudin, F. Ducept, G. Trystram, J. Clément, Characterization and CFD modelling of air temperature and velocity profiles in an industrial biscuit baking tunnel oven, Food Res. Int. 37 (2004) 1031–1039.