Copyright is owned by the Author of the thesis. Permission is given for a copy to be downloaded by an individual for the purpose of research and private study only. The thesis may not be reproduced elsewhere without the pennission of the Author. Mathematical Model of the Forced Cooling of Anodes used in the Aluminium Industry A thesis presented in partial fulfilment of the requirements for the degree of Master of Science Christopher Charles Palliser May 1994 in Mathematics at Massey University ABSTRACT The aluminium industry consumes large amounts of electrodes, especially anodes, to operate the smelters. These anodes must be baked at high temperatures in order to give them certain mechanical and electrical properties, after which they are cooled. Baking is done in large furnaces made up of pits inside which the anodes are placed in layers and surrounded by packing coke. The furnaces are of two types - open and closed. In a closed furnace, the pits are lined with refractory bricks inside which flues run vertically and large covers are used to close over parts of the furnace. This thesis presents a mathematical model of part of a forced cooling section of a closed furnace, where air is being sucked or blown through the flues by fans, so that the anodes cool more rapidly. Both one- and two-dimensional models are developed in order to calculate the transient temperature distribution in the anodes, packing coke and side flue wall. For the two-dimensional model, the transient temperature and pressure distributions of the air in the side wall flues and fire shafts are also calculated. After exploring an analytical method for the one-dimensional case, numerical techniques are used thereafter. Given initial block and air temperatures, the two-dimensional model allows calculation of the appropriate temperature and pressure distributions for various mass flows of air in the side wall flues and fire shafts. The results show that for a sufficiently high mass flow, the anodes can be cooled enough so that they can be safely removed from the pits after three fire cycles (the length of time the anodes are exposed to forced cooling). ACKNOWLEDGEMENTS This thesis is dedicated to my wife, Alataua. Without her love, patience, encouragement, hard work and prayers over the last five years, none of this would have been possible. Fa' afetai lava la'u pele. My love and thanks to my parents, Ray and Trixie, who never lost faith in me and have given me so much. Lots of alofa and kisses for my children, Matthew, Carissa, Joshua, Susan and Samuel. My special thanks to my supervisor, Dr R. McKibbin, for his expertise, guidance and patience. I am also grateful to Professor G. C. Wake and other colleagues in the Department of Mathematics for their encouragement and support. I thank New Zealand Aluminium Smelters Ltd for the scholarship. Thanks to Mr R. Braithwaite and his colleagues for their assistance and hospitality. This thesis script has been typed by Mrs G. Tyson, whose skill and care I gratefully acknowledge. Thank you to my sisters, brothers-in-law and families: Ruth and Ian Campbell, Claire and Peter Bruin, Catherine and Peter Kelly and Keneti Faulalo. Thanks Peter Bruin for your expert assistance with the computer. I also thank my friends who have supported and prayed for me, especially Haig and Pam Thomson, Bill and Felicity Kimberley, Tony and Maraia Statham. CONTENTS ABSTRACT ACKNOWLEDGEMENTS NOMENCLATURE 11 iii VI CHAPTERl CHAPTER2 CHAPTER3 CHAPTER4 CHAPTERS CHAPTER6 INTRODUCTION 1.1 Background 1.2 Previous work 1.3 This work - an outline DEVELOPMENT OF THE MODELS 1 2 4 2.1 Introduction 7 2.2 Calculation of cooling areas around the anode 9 2.3 One-dimensional model 13 2.4 Two-dimensional model 15 DERIVATION OF THE HEAT EQUATION IN THE BLOCK 3.1 The heat equation 3.2 The boundary conditions 3.3 Summary THERMAL PROPERTIES OF THE BLOCK 4.1 Introduction 4.2 Carbon anode 4.3 Packing coke 4.4 Flue wall ONE-DIMENSIONAL ANALYTICAL SOLUTION 5.1 Introduction 5.2 Solution with different boundary conditions ONE-DIMENSIONAL NUMERICAL SOLUTION 6.1 Introduction 6.2 Constant thermal conductivities 6.3 Non-constant thermal conductivities 6.4 Results 6.5 Discussion of results IV 17 19 20 22 22 25 26 28 28 34 34 50 52 55 CHAPTER 7 TWO-DIMENSIONAL NUMERICAL SOLUTION 7.1 Introduction 57 7.2 Air flow in the flue 57 7.3 The heat equation in the block 63 7.4 Heat balance at the flue wall 81 7.5 Results 86 7.6 Discussion of results 87 CHAPTERS THE AIR PRESSURE IN A SECTION 8.1 Introduction 94 8.2 Calculation of the pressure 94 8.3 Results 99 8.4 Discussion of results 103 CHAPTER9 SUMMARY AND CONCLUSION 9.1 Introduction 104 9.2 One-dimensional model 104 9.3 Two-dimensional model 105 9.4 Air pressure 105 9.5 Conclusion 106 REFERENCES 107 V NOMENCLATURE Alternative uses are separated by a semi-colon. Dimensions are in square brackets. A,B a, b,c,d,e,f Bi cm, m = 1 to 9 E,e ea, ep, ew estTf F, f I\ Fo,Fo Foa, Fop, Fow G h hc,hr k ka, kp, kw kf defined functions constants lengths associated with a flue and a flue wall [m] Biot number [ - ] constants specific heat capacity at constant pressure of the block and the fluid or air [J/kg K] specific heat capacity at constant pressure of the anode, packing coke and flue wall [J/kg K] total and internal energy per unit mass of the block [J/kg] internal energy per unit mass of the anode, packing coke and flue wall [J/kg] estimated fluid temperature [K] defined functions Fourier [ - ] and 'hatted' Fourier number of the block [1/K] Fourier number of the anode, packing coke and flue wall [ - ] defined function enthalpy per unit mass of the block [J/kg]; heat transfer coefficient [W/m2K] heat transfer coefficient for convection and radiation [W/m2K] thermal conductivity of the block [W /mK]; real constant thermal conductivity of the anode, packing coke and flue wall [W /mK] thermal conductivity of the fluid [W /mK] VI m max. mm. Nu p percentdiff pfl, pfs Pr Q Qw Re s T t Ta, Tp, Tw Tf V V w X x,y,z characteristic length [m] width of the block [m] depth of the block [m] mass flow of the fluid [kg/s] maximum minimum Nusselt number [ - ] pressure in the block; pressure of the fluid or air [N/m2] percentage difference between fluid temperatures pressure of the fluid in a flue and a fire shaft [N/m2] Prandtl number [ - ] rate of heat flow per unit area of the block [J/sm2 ] rate of heat flow per unit area at the flue wall [J/sm2J Reynolds number [ - ] surface of an elemental volume in the block [m2 ] temperature of the block [K] time [s] temperature of the anode, packing coke and flue wall [K] fluid or air temperature [K] elemental volume in the block or a flue [m3] velocity of the fluid [mis] width of a 'large' flue [m] defined function spatial coordinates vii Greek I\ cx,cx cxa, exp, aw ~Q ~Qa, ~Qp, ~Qw ~Qf ~t ~ta, ~tp, ~tw ~a,~xp,~xw ~y µ p pa, pp, pw pf thermal diffusivity [m2/s] and 'hatted' thermal diffusivity of the block [m2 /sK] thermal diffusivity of the anode, packing coke and flue wall [m2/s] defined functions change in the pressure of the fluid; pressure difference across a fan [N/m2 ] heat transferred from the block in time ~t [J] heat transferred from the anode, packing coke and flue wall in time ~t [J] heat gained by the fluid or air in time ~t [J] overall time step length [s] length of time step in the anode, packing coke and flue wall [s] distance between mesh points in the x-direction [m] distance between mesh points in the anode, packing coke and flue wall [m] distance between mesh points in they-direction [m] emissivity of the flue wall [ - ] ; roughness of the flue wall [m] defined functions defined functions real constant; defined function; friction factor [ - ] dynamic viscosity of the fluid [kg/ms] density of the block [kg/m3 ] density of the anode, packing coke and flue wall [kg/m3 ] density of the fluid [kg/m3] viii CT Subscripts 1, i = 1 to Nx J, j = 1 to Ny p q, q = 1 to Nxa r, r = 1 to Nxp s, s = 1 to Nxw w n Superscripts n Stefan-Boltzmann constant [ W/m2K4] integration variable defined function mesh points in the x-direction mesh points in the y-direction constant pressure [N/m2] mesh points in the x-direction in the anode mesh points in the x-direction in the packing coke mesh points in the x-direction in the flue wall flue wall natural number time steps, n = 1 to Nt IX CHAPTERl INTRODUCTION 1.1 Background The aluminium industry consumes large amounts of electrodes, especially anodes, to operate the smelters. These carbon anodes are made of petroleum coke held together by a pitch binder. They must be baked to a given temperature, approximately 1200 °C, following a given temperature profile (no more than 10 -15 °C/hour) in order to end up with the required mechanical and electrical properties. Baking is done in large furnaces made up of pits inside which the unbaked anodes are placed in layers and surrounded by packing coke. These furnaces are of two types, one of which is the Riedhammer (vertical ring, or closed) furnace, a schematic of which is shown in Figure 1.1. In the Riedhammer furnace, the pits are lined with refractory bricks inside which flues run vertically and large covers are used to close over parts of the furnace. A typical Riedhammer furnace consists of two or three fire trains grouped together on a rectangular shaped ring. As shown in Figure 1.1, each fire train comprises about fourteen sections, or sets of pits, and consists of three zones - preheat, fire and cooling zone. Hot combustion gases flow through the flues in the fire and preheat zones, whilst air flows through the flues in the cooling zone. The cooling zone is divided into two parts - natural and forced cooling. In the natural cooling part, the anodes are just left to cool. The forced cooling sections have big fans which either blow or suck air through the flues to increase the rate of cooling of the anodes (see Figure 1.1). There is one fan per section. The rate of cooling is not constrained and may be done as quickly as possible. The fans, fire ramps and exhaust manifold are moved in the fire direction by one section every 32 or 36 hours. This time period is called the fire cycle. 1 ~ Fire direction Cooling fans Exhaust manifold blower cover off Fire ramps 14 Forced Natural Fire Preheat cooling cooling Figure 1.1 A schematic longitudinal view of a typical fire train arrangement in a Riedhammer furnace at NZ Aluminium Smelters Ltd (see Bourgeois et al., 1990). 1.2 Previous work Mathematical modelling of ring furnaces started seriously in 1980 with Furman and Martirena (1980). Since then others have developed different aspects of ring furnace ( open and closed) mathematical modelling. The advent of computational fluid dynamics packages has enabled the development of more elaborate models. Due to the ring furnace's large dimensions and time constant (2½- 3 weeks), experimentation on a real furnace is not only impractical, but also risky, lengthy and costly. Hence the need for mathematical models in order to analyse and predict performance. Furman and Martirena ( 1980) used a three-dimensional finite difference model. The total duration of the baking cycle was simulated. Therefore the time period was long, 300 - 400 hours. In order to save computation time, most of the time steps used were correspondingly long. The first 10 hours were divided in steps of 0.1, 0.2, 0.5, 2.2, 3.0 and 4.0 hours, all further steps were 7 .0 hours long. To ensure the stability of the calculations with these long temporal steps, an implicit Crank-Nicholson scheme was used. The equations were solved using successive relaxation. 960 nodes were used - 4 in the y­ direction (x-direction in this thesis), 20 in the z-direction (y-direction in this thesis) and 12 in the x-direction (z-direction in this thesis). A 'sensitivity' analysis was performed. This involved introducing a significant variation of a 2 given property and observing the corresponding temperature calculations. It was found that the thermal conductivity of the packing coke and the vertical gradient of the gases' temperature along the flues were the parameters which decisively influenced the calculations. It was claimed, therefore, that these were the only parameters which needed to be known accurately. The temperature-dependent thermal properties of the anodes, packing coke, flue walls and gases were not adjusted as the temperature varied. The paper by de Fernandez et al. (1983) used the identical set of nodes as used in the paper of the previous paragraph, so it was a three-dimensional model. However, no mention was made of the solution method. Initially the temperature difference between the top and bottom layer of anodes in a pit was calculated. It was found that this difference was in much closer agreement with the experimental difference, when the thermal properties of the anodes, packing coke and flue walls were adjusted according to the temperature reached at the end of the last time step. The 'negative' image of the heating temperature distribution was used to try and model the cooling temperature distribution. It did not work and this was attributed to two factors, one of which was the occurrence of natural convection. It was suggested that a battery of fans be used on uncovered sections in order to overcome this natural convection. Clearly, fans were not being used at this particular smelter when this paper was written. The transient two-dimensional model presented by Bourgeois et al. (1990) neglected the heat transfer in the longitudinal direction, that is, from the fire shaft to head wall. This helped to simplify the model and keep CPU time down. It was claimed that experimental studies had shown that this longitudinal temperature variation was small compared with the flue wall to anode-centre and vertical variation. Unlike the models from the earlier papers, this one incorporated pressure measurements, namely the draught profile along a fire train. One of the limitations of the model was that it could not determine the temperatures in the forced cooling sections (they were considered disconnected from the main fire train). There was fairly good agreement between calculated and experimental results. Bui et al. (1992) divided a furnace section into four zones - fire shaft, under-lid, pit and under-pit. For gas flow distribution, it was stated that the 3 underlid zone was the most important; whereas for heat transfer to, and therefore presumably from, the anodes, the pit zone was the most important. A three-dimensional model of heat transfer and fluid flow for the under-lid zone of any section of the fire train was developed. For the purposes of validation, the under-lid zone of the first covered cooling section was simulated. The solution procedure was not discussed in detail, but the general purpose computational fluid dynamics PHOENICS code was used as a solver. A larger number of nodes was used, namely 23180. The calculated results followed reasonably well the trend of the measured ones. 1.3 This work - an outline A heat transfer and pressure distribution model of part of a forced cooling section (from now on called section) of the Riedhammer furnace is presented in this thesis. The cover has been removed, the packing coke is still in place and a blowing or sucking fan (blower or sucker) is positioned over the fire shafts - see Figure 1.2. Each uncovered section is a separate entity disconnected from the main fire train. end flue wall I I y foundation Figure 1.2 l head wall fire shaft I direction of air flow Schematic longitudinal view of a flue wall in a forced cooling section 4 The model presented here is used to determine the effect that different mass flows of air have on: (a) the block (anodes, packing coke and side flue wall) temperature, (b) the air temperature in the side wall flues, and ( c) the air pressure in the side wall flues and the fire shafts. Heat is transferred from the anodes into the air in the side wall flues via the packing coke and the side flue wall. The transfer of heat is taken to be by conduction in the anodes and flue wall, and is assumed to be by conduction in the packing coke. The heat is then transferred into the air by convection and radiation from the surface of the flue wall. The section is three-dimensional. Simplified one-dimensional and two­ dimensional models are studied in Chapter 2. This is done by concentrating on the anodes, packing coke, side flue walls and side wall flues part of the section. Simplifications involved in developing the models are justified by a result from Bui et al. (1992) and dimensional considerations. In Chapter 3, the heat conduction equation is derived in the present context, to give a partial differential equation. The consequences of assuming constant thermal conductivities or otherwise is examined. Using the models from Chapter 2, boundary conditions are then added to the partial differential equations to give boundary value problems. The thermal properties of the anodes, packing coke and side flue walls are discussed and calculated in Chapter 4. The one-dimensional heat equation is solved analytically for three different sets of boundary conditions in Chapter 5. The last set are those of the model developed in Chapter 2. In Chapter 6, explicit numerical methods are used to solve the boundary value problem. Using these, the problem is solved for the case where the thermal conductivities are constant within the anodes, packing coke and side flue wall, to give the transient temperature distribution in the block. In this constant 5 thermal conductivities case, the thermal conductivities are the only thermal property that is not being adjusted as the temperature varies with time; whereas in the non-constant case all thermal properties are being adjusted as the temperature varies with time (it is assumed these properties are dependent on temperature). For the non-constant thermal conductivities case, the boundary value problem is set up, but not solved due to the introduction of non-linear terms. For this one-dimensional case, the air temperature in the side wall flues is assumed to be constant. The two-dimensional model is developed in Chapter 7. Only the constant thermal conductivities case is considered. This builds on the work done in the previous chapter on the one-dimensional model. Unlike the one-dimensional model, the temperature of the air in the side wall flues is changing with time and space as heat is transferred into it from the block. This is modelled using an implicit numerical method and combined with the two-dimensional heat equation for the block to give the transient one-dimensional temperature distribution of air along the flues and the transient temperature distribution in the block. As a check on the working, the heat given out by the block and the heat gained by the air in the flues is calculated. Because of the set-up of the model, these should be approximately equal if the calculations are done correctly. The transient one-dimensional pressure distribution in the side wall flues and fire shafts is calculated in Chapter 8. 6 CHAPTER 2 DEVELOPMENT OF THE MODELS 2.1 Introduction Ideally the model should be three-dimensional and provide the temperature distribution that is calculated throughout the section, starting from the fire shafts, moving into the head wall, then the gallery, then the area between the foundation and the bottom of the anodes and finishing in the flues, flue walls, packing coke and anodes. However, this would be too complicated and time consuming, so the model is simplified by: (i) restricting it to the flues, flue walls, packing coke and anodes zone, which is the most important one for anode cooling (see Bui et al., 1992); (ii) making it one- and two-dimensional. Most of the cooling of the anodes occurs through them transferring heat into nearby air. Therefore (ii) of the model simplification is done by finding which sides of the anodes have the most exposure to air. It is then assumed that the anodes transfer most of their heat from these sides. Heat transfer from the remaining sides is then ignored, enabling the development of a two-dimensional model. In order to understand and solve this two­ dimensional model, a one-dimensional model is developed first. 7 end wa flue I!__ Pit containi 20 anodes ng - Packin coke g ~ end flue wall end flue wall \ Head wall I outer side wall flue outer side 9ue wall I / I DDD □ DDDDDDD1DDD □□ \ I □ □ ~ □ □ Outer pit □ llJ- i--- □ □ □ □□□p □□□ q □□□□ DODD □ □ ll / \ □ □ □ inner side inner side Inner pit □ wall flue flue wall 0.820 □ \ I □ m □ □ DD □ DDDDDDDDDDDDD □ l □ ~0.616--l □ m □ □ b.s10 Inner pit □ m □ □ □ □ DDDDDDDDDDDDDDDD □ □ □ □ □ Inner pit □ ff □ □ □ DDDDDDDDDDDDDDDD □ □ □ □ □ Outer pit □ □ □ □ □ DDDDDDDDDDDDDDDD \. outer side flue wall Figure 2.1 Schematic plan view of a section with the packing coke removed from the top to expose the top row of anodes. I, I r;; D.495 In ~ '- - -- - '- - '- (The flues are assumed to be rectangular, whereas in fact the corners of each flue are rounded.) 8 Fire shaft 2.2 Calculation of cooling areas around the anodes The anodes are rectangular blocks. There are twenty anodes stacked in a pit - five rows of four. For the purposes of modelling, these twenty anodes are treated as one big anode. This big anode has six sides, each of which is adjacent to air via packing coke and/or a flue wall: (a) one side is adjacent to the air which is above the anode; (b) one side is adjacent to the air which is flowing between the bottom surface of the anode and the foundation; (c) two sides are adjacent to the air flowing in the end wall flues; (d) two sides are adjacent to the air flowing in the inner/outer side wall flues. (As can be seen from Figure 2.1, if the pit is an inner one, then these two sides are adjacent to the inner side wall flues; whereas if the pit is an outer one, then one of these sides is adjacent to the inner side wall flues and the other side is adjacent to the outer side wall flues.) The areas of these different sides of the anode that are adjacent to air are now calculated. Since most of the anode cooling occurs through heat transfer into nearby air, the area between flues is ignored. A line of symmetry, an adiabatic boundary, is assumed to run through the centre of each flue, see Figure 2.3, so the appropriate areas are calculated on half flue measurements. Since the dimensions of each pit/flue wall are identical, the calculations are done for one pit/flue wall area (for measurements refer to Figures 2.1 - 2.3). (a) The area of the side that is adjacent to the air which is above the anode = (length of anode) x (width of anode) = (0.616 X 5) X 0.810 = 2.5 m2 . (b) The area of the side that is adjacent to the air which is moving between the bottom surface of the anode and the foundation = area as in (i) = 2.5 m2 • 9 (c) The area of the two sides that are adjacent to the air flowing in the end wall flues = (length of one end wall flue + 2 x ½ width of one end wall flue) x (depth of anode) x (number of flues in one end wall) x (number of end walls) = (0.120 + 2 X 0.06) X (1.166 X 4) X (5) X (2) = 11.2 m2 ( d) The area of the two sides that are adjacent to the air flowing in the inner/outer side wall flues = (length of one inner/outer side wall flue + 2 x ½ width of one inner/outer wall flue) x (depth of anode) x (number of flues in one inner/outer side wall) x (number of inner/outer side walls) = (0.178 + 2 X 0.074) X (1.166 X 4) X (16) X (2) = 48.7 m2. Therefore total area of the sides not adjacent to the air in the inner/outer side wall flues = (a)+ (b) + (c) = 16.2 m2 Pit Anode Figure 2.2 Flue wall Schematic cross section of an anode in a pit 10 inner/outer side wall _ lines of symmetry ~~----+---~--\ 0.060 0.178 m DJ m end wall m 0.120m lw 0.065m [!J Figure 2.3 Schematic plan view of a pit comer and flues Comparing this figure with that of (d), namely 48.7 m2 , it is clear that the greatest areas of the anode adjacent to the air are the two sides adjacent to the air flowing in the inner/outer side wall flues. It is therefore assumed that the anode transfers most of its heat into the air flowing in the inner/outer side wall flues. It is further assumed that there is no heat transfer into the air which is: ( a) above the anode; (b) flowing between the bottom surface of the anode and the foundation; ( c) flowing in the end wall flues. This last assumption allows modelling in two dimensions rather than three. This agrees with experimental studies mentioned in Bourgeois et al. (1990). The studies found that longitudinal (from fire shaft to head wall) temperature variation was small compared with the flue wall to anode-centre and vertical variation. The x-direction is taken to be across the pits, and the vertical as the y-direction. (The z-direction is taken to be along the pits - see Figure 1.2.) This is developed further in §2.3 and §2.4. 11 □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ packing coke anode end flue wall I I I I I I I I lines of symmetry :/ I I I I I X end flue wall Figure 2.4 [[] rn m m uJ uJ ~ ~ [] '~ inner/outer side wall flue Schematic plan view of one pit with anode, packing coke, flue walls and flues 12 2.3 One-dimensional model It is assumed that lines of symmetry, which denote adiabatic boundaries, run through the centre line of each pit and flue, as shown in Figure 2.4. All pits are treated alike, that is, there is no difference between outer and inner pits [although it is realized that there is heat loss through the outer side flue walls on either side of the section - see Bui et al. (1992)]. From now on, inner/outer side flue walls are referred to as the flue walls and inner/outer side wall flues as flues. A representative slice of anode, packing coke, flue wall and flue is selected. This slice is adjacent to one flue, as shown in Figure 2.4. As in §2.2, it is assumed that the anode transfers heat into the air in the flue not only along a(= 0.178 m) but also along band c too (both= 0.074 m) see Figure 2.5. This poses a difficulty in regard to the thickness of the flue wall. For example, from the outer edge of the packing coke, d, to the flue wall at a, distance= 0.140 - 0.074 = 0.066 m. But what about the distance from d to the flue wall at b or c? lines of symmetry !-- Fl,re wall ---~; Rue I I I I I I I px I I anode Packing coke Figure 2.5 m om d a J.178m l C ( )10.074m OTJJ I I Schematic plan view of the representative slice showing anode, packing coke, flue wall and flue. 13 In order to overcome this difficulty, each flue is thought of being stretched or elongated, as illustrated in Figure 2.6, so that the sixteen flues in a flue wall become one 'large' flue. This large flue, amongst other things, must have the same cross-sectional area as the sixteen individual flues (for measurements refer to Figures 2.3 and 2.5). Cross-sectional area of sixteen individual flues = ( cross-sectional area of one flue) x 16 = (0.178 X 0.148) X 16 0.421504 m2 Cross-sectional area of 'large' flue (length of 'large' flue) x (width of 'large' flue) = (length of one individual flue x 16 + distance between two individual flues x 15 + distance e + distance f) x W, where W = the width of 'large' flue, = (0.178 X 16 + 0.060 X 15 + 0.030 + 0.030) X W 2 = 3.808 xWm, and this has to equal 0.421504 m2 , ⇒ W = 0.111 m. Therefore the thickness of the flue wall in the case of the 'large' flue is 0.140 - 0 · 111 = 0.085 m (rounded) see Figure 2.6. Jines of symmetry --- --: J : I 'Larae' flue I I ...-("' 0 I I 1 ___ L I I X I I Figure 2.6 Schematic plan view of 'large' flue and surroundings. 14 The block is defined to be the shaded area in Figure 2.6. Lx width of block = (width of anode)+ (width of packing coke)+ (width of flue wall) 0.585 m. Note: In evaluating the Reynolds number and the heat transfer coefficient in Chapters 6, 7 and 8, the flue dimensions of Figure 2.5 are used, that is the 'real' flue dimensions. 2.4 Two-dimensional model The one-dimensional model developed in §2.3 is extended by introducing a vertical or y component as shown in Figure 2.7. axes for sucker axes for blower I line of symmetry I I I Air covering of packing coke I L X I I I packing I anode Ly I coke I I r X 'It Foundation Figure 2.7 flue wall l flue I line of symmetry Direction of air flow for a sucker I I Direction of air flow 1 for a blower I I I I Schematic cross section of the two-dimensional block and flue The same representative slice is taken as for the one-dimensional case, except that here the slice has a y component. In fact the slice extends vertically from the top surface of the anode to the bottom surface of the anode. LY depth of anode = 1.166 X 4 = 4.664 m. 15 The x and y axes are positioned differently according to whether a blower ( L ) or a sucker ( r' ) is used, Blower versus sucker In reality, in the case of a blower operating, the air has to pass through the fire shafts, the gallery and into the area between the bottom of the anode and the foundation prior to it flowing up through the flues. This is not the case if a sucker is operating, since the air has no 'history' like this prior to it flowing down the flues. The air here is drawn from the air above the top surface of the anode (see Figure 1.2). Hence for a blower, the inlet air temperature is determined at the top of the fire shafts, whereas for a sucker it is determined at the top of the flues. But given an inlet air temperature for a blower (air temperature at the top of the fire shafts, which is atmospheric temperature) the two-dimensional model for a blower is unable to determine what the air temperature will be when the air arrives at the bottom of the flues - the inlet air temperature for the model. It is assumed that the inlet air temperature for the model in the case of a blower is atmospheric temperature. Atmospheric temperature is assumed to be 20 °c. There is no such difficulty for a sucker, since the inlet air temperature in reality (atmospheric temperature) is the same as that for the two-dimensional model for a sucker. Therefore for the purpose of the two-dimensional model, it makes no difference to the temperature distribution as to whether a blower or a sucker is operating except that the block temperature distribution is reversed, since y = 0 for a blower corresponds to y = LY for a sucker and y = LY for a blower corresponds toy= 0 for a sucker (see Figure 2.7). 16 CHAPTER3 DERIVATION OF THE HEAT EQUATION IN THE BLOCK 3.1 The heat equation The rate at which heat accumulates in an elemental volume (V) of block material is given by 1t fJf pE dV where p is the density, tis time and E is V the total energy (internal, kinetic and potential) per unit mass of the system. This must be balanced by the rate of heat leaving V. This is given by fJ - Q dS, where Q = dQ/dt is the rate of heat flow (conductive, convective aid radiative) per unit area at the surface (S) of V. First consider 1t fff pEdV. V Since V is closed (mass does not enter or leave V) and stationary then E = e, where e = internal energy per unit mass, and 1t ff f pE dV = ff f a ot (pe) dV V V = fff ae p ot dV V since p is constant with respect to time - see Chapter 4. Introducing the enthalpy h = e + £, where p = pressure (Currie, 1993) and p using the fact that p and p are constant with respect to time gives: d fff p E dV = fff oh dt p ot dV V V fff ah aT (T = temperature) = p oT dt dV V fff oT = p cp dtdV V 17 . c1h h ·t· h . smce aT at constant pressure = cp = t e spec1 1c eat capacity at constant pressure. Now consider ff -Q dS. Since there is no convective or radiative effect, s by Fourier's Law (Rogers and Mayhew, 1992) Q = - k VT (k = thermal conductivity) Therefore ff -QdS=ff k VTdS s s = ff k VT,ndS s = fff V ·(kv'T) dV by Gauss' Divergence V Theorem. Therefore fff c1T ff f V •(kv'T) dV p cp dtdV = V V c1T V ,(kv'T) (3.1) ⇒ p cP dt = If k is constant, then Equation (3 .1) becomes i.e. I ~Tt = an2T I h k/ h 1 d'ff · · . 0 v . w ere a = pep= t erma 1 us1v1ty. If k is not constant and varies with temperature, then in the one-dimensional case Equation (3 .1) becomes 18 _ dk (aT) 2 + k a 2 T - dT ax ax2 Therefore aT (aT) 2 a 2 T " dk at = tx ax + a ax2 where a = dT / pep . In the two-dimensional case Equation (3.1) becomes Therefore pep t~ = fx (k ~J) + fy (k ~~) _ ak aT + k a 2 T + ak aT + k a 2 T - ax ax ax2 ay ay ay2 3.2 The boundary conditions Referring to (a) and (b) on p.11 of §2.2, (a) No heat transfer into the air which is above the anode ⇒ ⇒ { y = 0, adiabatic boundary on y = L y, aT {y = O, - = 0 on ay y = Ly, sucker blower . sucker blower (b) No heat transfer into the air which is flowing between the bottom surface of the anode and the foundation ⇒ ⇒ adiabatic boundary on Y { y=L' y=0, aT {y = Ly, sucker - = 0 on ay y = 0, blower. 19 sucker blower Referring to the assumption made in §2.3, that is, it is assumed that an adiabatic boundary runs through the centre line of the pit e1T ⇒ ax = 0 on x = 0. Suppose the block has some initial temperature, T 1, ⇒ T = T 1 when t = 0. Clearly the temperature at x = Lx (on the flue wall) is decreasing with time because of the air flowing past and heat being transferred from here into the air. Let this decreasing temperature be T w(t), where w = wall and (t) denotes dependence on time 3.3 Summary The boundary value problem for the heat conduction in the block is now specified, both for the one- and two-dimensional models. One dimension e1T dt (x, t) e1T dt (x, t) constantk } non-constant k with aT (0 t) dX ' = 0, t> 0, T(Lx, t) = T w(t), t> 0 and and T(x, 0) =Tl, 0 ::;; x ::;; LX. Two dimensions dT (a2T a2 T) dt(x, y, t) = a ax2 + ay2 constant k t~ (x, Y, t) = & [(~JJ + (~~J] o::;;x:s;Lx, t > 0 o::;;x:s;Lx 0::;; y::;; Ly t>O [a 2T a2 T] + a dX2 + ay2 non-constant k 20 with oT ax (0, y, t) = 0, 0 :=:; y :=:; Ly, t>0, T(Lx, y, t) = Tw(t), 0 :=:; y :=:; Ly, t > 0, T(x, y, 0) = Tl' 0 :=:; x :=:; LX, 0 :=:; y :=:; Ly, and aT aT oy (x, 0, t) = oy (x, Ly, t) = 0, 0 :=:; x :=:; LX, t > 0. 21 CHAPTER 4 THERMAL PROPERTIES OF THE BLOCK 4.1 Introduction The thermal properties of interest are of course those that appear in the heat k a= - and pep equation, namely k, p and cP. Once these are found, 11. dk / a = dT p cp can be calculated. This chapter is included since in the readings for this thesis, values or formulae for the values of k, p and cp were rarely given. The non-referenced values and formulae given here were obtained from Braithwaite (private communication, 1993). In order to simplify expressions and save computing time, linear approximations are used for formulae when approximations are deemed necessary. As most of the block temperatures in the section;:::: 400 K (127 °C), linear approximations seem reasonable. For the linear approximations, the equation p(T) = p(T 0) + (T - T 0)p[T 0, T 1] p(T )-p(T ) is used (Burden and Faires, 1989) where p[T0 , T 1] = 1 T 0 4.2 Carbon Anode Thermal conductivity, ka (W /mK) From Log and Oye ( 1990), 1- 0 ka = ka0 + 0.274-0Zao (-fr - -{Tc;) - 2.8 x 10-\ka0)°-77 (T - T0 ) where T0 = room temperature 293 K say, and ka0 = thermal conductivity of anode at room temperature, in W/mK. ka0 is given by the equation, ka0 = -1.2105 + 0.1744 Le = 4.1959, 0 since Le ( crystallite height) = 31 A for anodes. 22 ka is approximated using the points (400, 4.91) and (970, 6.35) - see Figure 4.1. >, -> - 6 0 :l 5 "C -i:: ~ 4 0 E 0 - 3 ;: 2 ('0 -E 1 ,_ Q) .c: 0 I- 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Ct) s::t- LO (0 I'-- CX) O') Temperature (K) Figure 4.1 Thermal conductivity of the anode [from Log and Oye (1990)] ka = (2.53 X 10-3)T + 3.90 (4.1) d~~) = 0.137-vicao .Jt -2.8 x 10-3 (ka0)0•77 W/mK2 is approximated using the points (400, 5.58 x 10-3) and (970, 5.63 x 10-4) - see Figure 4.2. 0.008 ~ 0.006 .... ·- ('0 > ........ ~ ~ ~ E 0.004 ... "C -.e r:: ;: -o- 0.002 0 (,) 0 0 (') 0 0 0 0 "SI" LO 0 0 (0 0 0 I'-- Temperature (K) Figure 4.2 0 0 co 0 0 O') The rate of change with respect to temperature of the thermal conductivity of the anode 23 d~~) = (- 8.80 X 10-6)T + 9.10 X 10-3 (4.2) The equations for ka and d(ka)/dT are valid for temperatures :S: 973 K (700 °C), as most of the anode temperatures in the section are. If any anode temperature > 973 K, then Equations ( 4.1) and ( 4.2) are still used. Specific heat capacity, cap (J/kgK) c~ = (- 0.13374402 x 10-8)T4 + (0.64604614 x10-5)T3 - (0.11572658 X 10-1)T2 + (0.95663719 X 101)T - 0.13731392 X 104 This is approximated using the points (400, 9.81 x 102 ) and (970, 1.73 x 103) - see Figure 4.3. -~ C) 2000 ~ - 1800 -:, - 1600 .... ctl 1400 Cl) .s:: 1200 (.) 1000 -(.) 800 Cl) 600 c.. Cl) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Cl) ""'" LO co I'- co CJ) Temperature (K) Figure 4.3 Specific heat capacity at constant pressure of the anode cap = (1.3l)T + 4.55 x 102 (4.3) Density. pa (kg/m3) The density of the anode is approximated by 1550. 24 4.3 Packing Coke Thermal conductivity. kp (W/mK) Using data supplied by Braithwaite (private communication, 1993), the average packing coke particle size is found to be 5.7 mm. This is then used in de Fernandez et al. (1983) to approximate kp using the points (473, 4.00 x 10-1) and (1073, 1.16) - see Figure 4.4. >, -> -0 :::; "C - C: ~ 0 E 0 ._ - :ii= res - E ,_ (1) .c: I- 6 5 4 3 2 1 0 400 ---■- -----· 900 Temperature (K) Figure 4.4 -· 1400 Thermal conductivity of the packing coke (interpolated from data supplied by Braithwaite, 1993) kp = (1.27 X 10-3)T - 1.99 X 10-l d~.f) = 1.27 x 10-3 W/mK2 Specific heat capacity. cpp (J/kgK) 4.08 X 107 cpp = 933.0 + (0.916)T - T2 (4.4) (4.5) This is approximated using the points (400, 1.04 x 103) and (1000, 1.81 x 103) - see Figure 4.5. 25 -~ Cl ,:,:: --:, .._, -(13 Q) .c: u -u Q) c.. Cl) 2000 1800 1600 1400 1200 1000 800 600 0 0 (".) 0 0 'Sj- 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ~ ~ ~ ro m o ~ N Temperature (K) Figure 4.5 Specific heat capacity at constant pressure of the packing coke cp = (l.28)T + 5.27 X 102 p Density, pp (kg/m3) The density of the packing coke is approximated by 670. 4.4 Flue Wall Thermal Conductivity, kw (W/mK). (4.6) This is approximated using the points (657, 1.35) and (1265, 1.77) - see Figure 4.6. >, .... > - 6 u :::, 5 "C - C: ~ 4 g E - 3 - 3: 2 (13.._, E 1 ■ ,_ -----------· Q) .c: 0 I- 600 800 1000 1200 1400 Temperature (K) Figure 4.6 Thermal conductivity of the flue wall (interpolated from data supplied by Braithwaite, 1993) 26 kw = (6.91 X 10-4)T + 8.96 X 10-l dt;) = 6.91 x 10-4 W/mK2 Specific heat capacity, cwp (J/kgK) (4.7) (4.8) The specific heat capacity at constant pressure of the flue wall is approximated by 1250. Density, pw (kg/m3) The density of the flue wall is approximated by 2440. 27 CHAPTER 5 ONE-DIMENSIONAL ANALYTICAL SOLUTION 5.1 Introduction k is assumed to be constant, so the boundary value problem is (from §3.3) oT o2T o s x sLx, t > 0 dt = a ax2' with oT t> 0, ax (0, t) = 0, (5.1) T(Lx, t) = T wCt), t> 0, T(x, 0) = T1, 0 s x s LX . Instead of trying to solve (5.1) straightaway, initially boundary value problems with simpler boundary conditions are solved analytically. Using the knowledge gained from solving these simpler problems, (5.1) is then solved. 5.2 Solution with different boundary conditions (a) The ends of the block are maintained at constant temperature, T0, for all time, t > 0, and the block is initially at temperature T 1 for 0 s x s LX. aT a2 T dt = a ax2' with T(0, t) = T(Lx, t) = To, T(x, 0) = T1, Let T(x, t) = 0(x, t) + \jf(x) Homogeneous boundary conditions are obtained for a differential equation involving 0, if \jf = TO. With this choice of \jf, the problem becomes 28 with 0(0, t) = 0(Lx, t) = 0, 0(x, t) = T 1 - T 0, t> 0, (5.2) can be solved using separation of variables. Let 0(x, t) = X(x) G(t) Equation (5.2) becomes G'(t) X(x) = aX"(x)G(t) ⇒ G'(t) X"(x) aG(t) = X(x) = - A say. ⇒ X'(x) + A X(x) = 0 with X(0) = X(Lx ) = 0 and G'(t) + AaG(t) = 0 Solving Equation (5.3), with A< 0, A= 0, A> 0 in tum: (i) (ii) If A = -k2 then X = c 1 cosh kx + c2 sinh kx Now, X(0) = 0 ⇒ c1 = 0 and X(L) = 0 ⇒ c2 = 0 A = - k2 gives a trivial solution. If A = o Now, X(0) = 0 ⇒ c4 = 0 and X(L) = 0 ⇒ c3 = 0 A = 0 gives a trivial solution. 29 (5.2) (5.3) (5.4) (iii) Take 'A, = k2 then X = c 5 eikx + c 6 e-ikx or X = C7 cos ~ X + Cg sin -{i: X Now X(O) = 0 ⇒ c7 = 0 and X(Lx) = 0 ⇒ Cg sin ~ Lx = 0 ⇒ ~Lx = nTC, nEN n2n2 ⇒ A = - 2-, n E N, LX which are the eigenvalues. which are the eigenfunctions; these are orthogonal on [O, LJ. Now, solving Equation (5.4) 0 {B . (nnx) -n21t2at/L2 E sm - e , n L X The most general solution of this form is 00 (nnx) n 2 2 /L 2 0(x, t) = I, Bn sin L e - n at , n=l X where 00 0(x, 0) = T 1 -T0 = I, Bn sin ntx X n=l Using the orthogonality of the eigenfunctions gives 2 rL, . (nnx) T 1 - To n Bn = L Ji (T 1 - T 0) sm L dx = [ 1 - (- 1) ] x o x nn T - T 00 1 (nnx) 2 ? 2 0(x,t) = 1 O L, - [1 - (- 1?] sin - e-n n-at/L, n n Lx n=l Now 0 is known, T is obtained as T = 0 + T0 . 30 (b) A boundary value problem is now solved that is almost the same as (5.1) except that here T(Lx, t) = T0 = a constant, whereas the boundary condition in (5.1) is T(Lx, t) = Tw(t) - varies (decreases) with time. aT a2T 0 :s; x :s; LX, t>0 dt = a ax2' with aT ax (0, t) = 0, t>0, T(Lx, t) = T0 , t> 0, T(x, 0) = T1 , 0 :s; x :s; LX. As in (a), let T(x, t) = 0(x, t) + \jf(x), and \jf = T0 gives homogeneous boundary conditions, so the problem becomes ae a2e 0 :s; x :s; LX , t>0 at = a at2 , with ae ax (0, t) = 0, t> 0, (5.5) 0(Lx, t) = 0, t> 0, 0(x, 0) = Tl - Ta, 0 :s; x :s; LX. Solving (5.5) using separation of variables gives where 00 0(x, 0) = T1 -T0 = L,Cn n=l [ (2n - l)nx] cos 2Lx Using the orthogonality of the eigenfunctions gives C 2fLx(T T) [(2n-l)nx]d 4(T1 -To)(-l)n-l n = Lx Jo 1 - o cos 2Lx x = ; (2n - 1) 31 . 0( ) _:!_(T -T )~ (-lt-l [(2n-1)1txJ-(2n-1)2rc2at14Lx 2 • • X, t - 1 O L..J --- COS ---- e 7t n=l (2n -1) 2Lx Now 0 is known, T is obtained as T = 0 + T 0. (c) See (5.1). Let T(x, t) = 0(x, t) + \jf(t). Homogeneous boundary conditions are obtained for a differential equation involving 0 if \jJ = T wCt). With this choice of \jf, the problem becomes with ae ax (0, t) = 0, t > 0, 0(Lx, t) = 0, t> 0, 0(x, 0) = T1 -Tw(t), o:::;x:s;Lx From (b), ~ [(in - l)nxJ 0(x, t) = Li Dn(t) cos lL n=l x h D (0) - :!_ (T 1 -Tw(O)) (- l)n-1 were n - 7t (2n-1) Therefore Equation (5.6) becomes ~ i.n (t) cos [(ln-l)nx] + ~T (t) Li dt n 2L dt w n=l x = ex f Dn(t)[(2n-1)1tl2 cos [(2n-l)nxl n=l 2Lx 2Lx 32 Using the orthogonality of the eigenfunctions gives ~ D (t) Lx + fLx ~ T (t) cos [(2n - l)rcx] dx dt n 2 Jo dt w 2L X = - a D (t) [(2n - l)rc]2 Lx (5.7) n 2L 2 X d 2L The above integral = dt Tw (t) x (- it-1 (2n - l)rc Hence Equation (5.7) becomes d dt Dn (t) + 11. Dn (t) = f(t) where "l _ [(2n - l)rcJ I\, - a 2L X and f(t) = 1t Tw(t) 4 (- It. (2n - l)rc Therefore i.e. 8 can now be determined and hence T = 8 + T w(t) can be calculated. The difficulty with this analytical solution is knowing T w(t). Numerical methods enable it to be calculated. 33 CHAPTER 6 ONE-DIMENSIONAL NUMERICAL SOLUTION 6.1 Introduction This chapter finds the solution to the one-dimensional boundary value problem with constant thermal conductivities using explicit numerical methods. All the thermal properties, except for the thermal conductivities of the anode, packing coke and flue wall are adjusted according to the temperature reached at the end of the last time step. The non-constant thermal conductivities case is also examined, but not solved due to the introduction of non-linear terms. 6.2 Constant thermal conductivities Recall from §3.3 that the boundary value problem is aT a2 T 0 s x s LX, t>0 dt = a ax2' with aT t> 0, ax (0, t) = 0, (6.1) T(Lx, t) = Tw(t), t> 0, T(x, 0) = T1, 0 s x s LX. aT a2 T dt and ax2 are approximated by using Taylor series expansions. aT 2 T(x, t + ~t) = T(x, t) + ~t dt + O(~t) aT = T(x, t + ~t) - T(x, t) + O(~t) at ~t aT (~x)2 a2T 3 T(x + ~x, t) = T(x, t) + ~x ax + 2! ax2 + O(~x) + ... oT (~x)2 a2T 3 T(x - ~x, t) = T(x, t) - ~x ~ + - 2 - 1- --2 - O(~x) + ... ax . ax 34 a2 T T(x + Lix, t) + T(x - Lix, t) - 2T(x, t) + O(Lix)2 ax2 = (Lix)2 So is approximated by T(x, t + Lit) - T(x, t) = a [T(x + Lix, t) + T(x - Lix, t) - 2T(x, t)] ( 6 _ 2 ) Lit (Lix/ where the error of the left hand side is of order Lit, and the error of the right hand side is of order (Lix)2. Tf is defined as the block temperature in Kelvin at mesh point i, i E z+, at time step n, n E z+. i runs from 1 to Nx where i = 1 is the mesh point on x = 0 and i = Nx is the mesh point on x = Lx. Lix is the distance between mesh points in metres. n runs from 1 to the end of the time period, where n = 1 corresponds to t = 0 and n = Nt corresponds to the end of the time period. Lit is the length between time steps in seconds. Rewriting Equation (6.2) using the subscript/superscript notation gives: T~+I _ T~ I I = Lit a~ [T::_1 + Tf_1 - 2Tf] I (Lix)2 i.e. (6.3) a!1Lit where Fof = - 1 - 2 = the Fourier number of the block at mesh point i at time (Lix) step n and it is dimensionless; af = the thermal diffusivity of the block at mesh point i at time step n, in m2 /s. 35 If Fof > ½, then the higher the value of T7 , the lower will be the resulting value of Tt1 . It can be shown that this may lead to difficulties in regard to the Law of Conservation of Energy (Rogers and Mayhew, 1992). If i = Nx, then the value of Ti+I = TNx+l in Equation (6.3) is not known. On the flue wall, the rate of heat flow per unit area (J/sm2) at time step n = Qw dQwn n aT~x . ' say= dt = - kNx dX, by Founer s Law. But by conservation of heat flow at the flue wall, the rate of heat per unit area flowing out through the wall by conduction = the rate of heat per unit area flowing into the fluid or air by convection and radiation. i.e. (6.4) where hn = the heat transfer coefficient at the flue wall at time step n, in W/m2K, Tfn = the fluid or air temperature at time step n, in K, (Tfn is assumed to be constant for all n, so Tf11 is denoted as Tf) and k~x = the thermal conductivity of the block at mesh point Nx at time step n, in W /mK. aTn hn _fu - =-- (TNn x - Tf) ax - kn Nx (6.5) Now imagine that the side of the flue wall is increased by Llx to Lx + LlX, which corresponds to the mesh point Nx + 1. Suppose that the temperature here varies in such a way that the temperature at i = Nx is always what it should be under the conditions of the actual problem. Using a Taylor series expansion gives 36 i3T~x T~x+l - T~x-1 2 ~x, = ~~~~ + O(Lix) ... 0 2Lix hn = - n (T~x -Tf) from Equation (6.5) kNx T~x+l = where the error of this approximation is (Lix)2. Substituting this into Equation (6.3) gives at i = Nx, where Bf = ~hn = the Biot number of the block at time step n and it is kNx dimensionless. As before, in order to satisfy the Law of Conservation of Energy, it is required that ⇒ Fon < __ 1 __ Nx - 2(1 + Bin) At i = 1, the value of Ti-I= T1_1 = T0 in Equation (6.3) is not known. A Taylor series expansion gives aTn The adiabatic boundary condition is a~ = 0 37 Tn Tn 2 - 0 = 0 where the error of the left hand side is of order (.6x)2. 2.6x Tn _ Tn 2 - 0 · Substituting this into Equation (6.3) gives at i = 1, h F n < 1 w ere o1 _ 2. Determining Fof a~ .6t Fo!l - - 1 - I - (.6x)2 (a) Obtaining af ki is constant within each of the anode, packing coke and flue wall, and equals ka, kp or kw depending on whether the calculations are occurring in the anode, packing coke or flue wall. ki is calculated using the temperature midway between the initial and assumed final average temperature of the block. Suppose this final average block temperature is Tf. For example, suppose the initial block temperature is 873 K and Tf is 373 K. So ki is calculated using 373 + 873 - 373 = 623 K. 2 ka = (2.53 x 10-3) 623 + 3.90, see Equation (4.1) kp = (1.27 x 10-3) 623- 1.99 x 10-1, see Equation (4.4) kw = (6.91 x 10-4) 623 + 8.96 x 10-1 , see Equation (4.7) 38 Pi = pa, pp or pw and ( cp) f = ( cap)f, ( cpp)f or ( cw P)f depending on whether the calculations are occurring in the anode, packing coke or flue wall. For example, suppose the initial and assumed final block temperatures are as before, and the temperature at mesh point i which is in the packing coke at time step n is 685 K. Then kp = (1.27 X 10-3 )623 - 1.99 X 10-l = 5.92x 10-1 W/mK pp = 670 kg/m3 and (cpp)f = (1.28)685 + 5.27 x 102 , see Equation (4.6) = 1.40 X 103 J/kgK Therefore 5.92 X 10-l = ------- m2 /s (670)(1.40 X 103 ) (apf is the thermal diffusivity of the packing coke at mesh point i at time step n.) (b) Obtaining Lix Depending on whether the calculations are occurring in the anode, packing coke or flue wall, Lixa, Lixp or Lixw are used respectively. To obtain Lixa, the number of spatial steps in the anode, Nxa, are specified and then the width of the anode is divided by Nxa-1. That is, A width of the anode DXa = N 1 xa- Similarly for Lixp and Lixw. ( c) Obtaining Lit See later 39 Determining Bin (a) Obtaining kw See (a) on p. 38. (b) Obtaining hn hn has two components - radiation and convection, that is hn = hrn + hen, where ~ = the heat transfer coefficient for radiation at time step n and hen = the heat transfer coefficient for convection at time step n (both in W/m2 K). hen is a function of Tfl , and Tfl is assumed it to be constant for all n. Therefore hen is denoted as he. where and where and ~ = £ cr (T~x - Tf) ( ( T~x)2 - (Tt)2) £ = wall surface emissivity and is dimensionless = 0.97 [Braithwaite (private communication, 1993)] cr = Stefan-Boltzmann constant = 5.67 X 10-8 W/m2 K4 . he = (kf)Nu f, kf = thermal conductivity of the fluid or air, in W/mK, Nu = the Nusselt number (dimensionless) R = the characteristic length 4 x area = . m. penmeter See Figure 2.3 for the dimensions of a flue. Note: Since the air in the flue is being blown or sucked past the end of the block, it is assumed that forced convection is occurring rather than natural or free convection. Buoyancy forces, which are associated with natural convection, are ignored here. 40 where and Nu = 0.023 Re0·8 Pr0.4. Re = Reynolds number (dimensionless) Pr = Prandtl number (dimensionless) (6.6) This equation for Nu assumes that: (i) the flow is turbulent, i.e. Re > 4000 (this is checked in the program); (ii) the flow is fully developed, i.e. the ratio of channel length to f is greater than 20, which it is since length of flue _ S,_ _ 4.664 :=:: f - f - 4 X 0.026344 29; 0.652 (iii) the physical properties are constant, which they are. or Re = (pf)vf (v = the velocity of the air, mis, pf= the density µ mt Re= µ(area) 4m of the fluid or air, kg/m3 , and µ = the dynamic viscosity of the air, kg/ms) (m = the mass flow, kg/s) =----- (6.7) µ(perimeter) Note: There are 5 fire shafts in a section and the centre one is blocked off. So if the total mass flow through all 4 fire shafts is x kg/s, then the mass flow through one flue is x/96 kg/s, since there are 96 flues in a section, and x/4 kg/s per fire shaft. ri1 = (pf)L .6.p where L = the connected load of the drive motor of a fan (1500 W for a blower, 3000 W for a sucker), 41 ~p = the total pressure difference across a fan (200 N/m2 for a blower or sucker), and pf "'" 0.75 kg/m3 - see §8.2, Figure 8.1. Therefore m "'" {5.6 kg/s for a blower 11.3 kg/s for a sucker = cfpµ Pr kf ' where cf P = the specific heat capacity at constant pressure of the fluid or air, in J/kgK. cfp, µ and kf all vary according to the temperature of the air, so he= hc(Tf) Values for kf, Pr and cfp were obtained from Perry et al. (1984). Then (kf) Pr µ = f cp (6.8) can be determined. Linear approximations are used for kf and cf P and a quadratic approximation for Pr. See Figures 6.1 - 6.3 respectively. For the quadratic approximation, the equation p(T) = p(T0 ) + (T-T0)p[T0, T iJ + (T-T0)(T-T 1) p[T0, T 1, T2] is used (Burden and Faires, 1989) where 42 >, -> -0 ::I "C -C: ~ 0 E 0 -3: (1l E ,_ Q) .c: I- 0.07 0.06 ,.,,...,..... ............ 0.05 / 11 /II 0.04 ,,,,.. ,...... ,...... 0.03 .,,;• ,.,,...,..... 0. 0 2 --------------- 300 500 700 900 11 00 Temperature (K} Figure 6.1 Thermal conductivity of air [interpolated from data from Perry et al. (1984)] kf = (5.65 X 10-5) Tf + 1.02 X 10-2 using the points (350, 3.00 x 10-2) and (1000, 6.67 x 10-2). -~ ~ 1180 ::; 1160 - ;· 1140 1 1 0 0 ;•-•-• .... ~ 1120 .c: : 1080 ;•-•-· .., 1060 ~ 1 0 4 0 •-•-•-·-· en 300 500 700 900 1100 1300 Temperature (K} Figure 6.2 Specific heat of air at 1 atmosphere [interpolated from data from Perry et al. (1984)] cfp = (1.38 X 10-1 ) Tf + 9.92 X 102 using the points (400, 1.047 x 103) and (1000, 1.130 x 103). 43 (6.9) (6.10) ~ 0.715f .c 0.71 ............. ~ 0.705 .\ /./ 0.7 • ~ 0.695 \, /. ~ 0.69 ., __ • Cl. 0. 68 5 --------1-----i------+------I 300 500 700 900 11 00 Temperature (K} Figure 6.3 Prandtl number of air at 1 bar [interpolated from data from Perry et al. (1984)] Pr = (1.89 X 10-7) Tf2 - (2.20 X 10-4) Tf + 7.54 X 10-l (6.11) using the points (300, 7.05 x 10-1), (600, 6.90 x 10-1) and (900, 7.09 x 10-1). Determining Lit (a) Obtaining Lita (the length of the time step in the anode) i.e. Foa1:1 < 1. 1-2 ( aaf )(Lita) 1 ----'---- < - (Lixa)2 - 2 (Lixa)2 Lita::::; -- 2aa!1 1 (Foaf is the Fourier number of the anode at mesh point i at time step n) (aa7 is the thermal diffusivity of the anode at mesh point i at time step n) The upper bound on Lita is obtained by calculating the maximum value that aaf can take. n max. aai = ka (pa)(min.(c'lJ )~) 44 Suppose the minimum block temperature equals the assumed final average block temperature which equals Tf (see the discussion earlier in this section). min. (caP)f = (l.31)Tf + 4.55 x 102, see Equation (4.3) (b) Obtaining Lltp (the length of the time step in the packing coke). (Llxp)2 Similarly Lltp ~ -~- and the upper bound on Lltp is obtained by 2apf calculating the maximum value that apf can take. max. ap~ = kp 1 (pp)(min.(cpp)1i) and min. (cpP)f = (l.28)Tf + 5.27 x 102, see Equation (4.6). ( c) Obtaining Lltw (the length of the time step in the flue wall) i.e. Fowf ~ min. (-2 1 , 1 ·n ) (Fowf is the Fourier number of 2(1+B1) the flue wall at mesh point i at 1 =---- 2(1 + Bin) ' (awf )(Lltw) 1 -~---<---- (Llxw)2 - 2(1 +Bin) (Llxw)2 LltW ~ ------ 2awf (1 +Bin) time step n) since Bin > 0 for all n. (awf is the thermal diffusivity of the flue wall at mesh point i at time step n) The upper bound on Lltw is obtained by calculating the maximum values that aw!1 and Bin can take. I 45 n max. awi = aw = kw see §4.4. S. B·n (Lixw)hn . 1 ·n . b mce 1 = kw , the maximum va ue of B1 1s o tained by calculating the maximum value of hn . max. hn = max. hr1 + max. hen max. hrn = ecr[ max. T~x - Tf] [ (max. T~)2 - (Tf)2] = ecr[initial T Nx -Tf] [ (initial T Nxf- (Tf)2 ] = ecr[Tix - Tf] [ (Tix)2 - (Tf)2 ] max.hen = hrl (max.kf)(max.Nu) = f_ = (kf)(Nu) since Tf is assumed to be constant f_ = he ( d) Obtaining the overall Lit The upper bound on the overall Lit = minimum (max. Lita, max. Litp, max. Litw). Temperature on the anode/packing coke boundary The temperature on the anode/packing boundary needs to satisfy the equations for temperature distribution in both the anode and packing coke. Ta~ and Foa ~ are defined as the anode temperature (in Kelvin) and Fourier number respectively at mesh point q, q E z+, at time step n, n E z+. q runs from 1 to Nxa, where q = 1 is the mesh point on x = 0, and q = Nxa is the mesh point on the anode/packing coke boundary. Similarly, Tp~ and Fop~ are defined as the packing coke temperature (in Kelvin) and Fourier number respectively at mesh point r, r E z+, at time step n, n E z+. r runs from 1 to Nxp, where r = 1 is the mesh point on the 46 anode/packing coke boundary, and r = Nxp is the mesh point on the packing coke/flue wall boundary. At q = Nxa. From Equation (6.3) TaN:~ = FoaNxa[TaNxa+l + TaNxa-1 + ( \ - 2JTaNxa] (6.12) FoaNxa Ta~xa+l is an imaginary mesh point outside the anode, which coincides with Tp~ only if Lixa = Lixp. At r = 1. From Equation (6.3) Tpg is an imaginary point outside the packing coke, which coincides with Ta ~xa-l only if Lixa = Lixp. By Fourier's Law, on the anode/packing coke boundary (6.13) Q(t) = dQ/dt = _ k VT = _ k aaT ~ _ k (!(x + Lix, t) - T(x - Lix, t)J, x l 2Lix Q (t) = the rate of heat flow per unit area, in J/sm2. i.e. _ ka r-a~xa+l -Ta~xa-lJ = _ kp r-P~ -Tpg J l 2Lixa l 2Lixp (6.14) Equations (6.12)- (6.14) involve three unknowns: Ta~xa+l' Tpg and Ta~.tia = Tp;i+1• These equations are solved for Tpt1 , the temperature at time step n + 1 on the anode/packing coke boundary, using the fact that TaNxa = Tpf. 47 Fopr {2T n 2[ ka Llxp ]T n Tpn 1 +I = -- P2 + --- aNxa-1 1 1 kp Llxa + [ka Llxp( \ - 2J + (~ - 2J]Tpr} (6.15) kp Llxa FoaNxa Fop1 where Fopr ka Llxp 11 = 1 + n -- . FoaNxa kp Llxa In order to satisfy the Law of Conservation of Energy it is required that the coefficients of Tp2 , Ta~xa-I and Tp7 ~ 0. Since ka, kp, aa~xa and ap1 > 0 for all possible temperatures of interest and Foa~xa :s; ½, Fop1 :s; ½ (the restrictions determined earlier), then the coefficients~ 0. Temperature on the packing coke/flue wall boundary This case is similar to the anode/packing coke boundary one. Therefore where T n+l _ Fowr { 2T n 2[ kp LlXW]T n W1 - -- W2 + --- PNx -1 1 2 kw L1Xp P + [ :~ :;( Fo;N,v - 2) + (Fo~f - 2 J]Twf} 12 = l + Fowr kp LlXW FopNxp kw LlXp (6.16) Tw~ and Fow~ are the flue wall temperature (in Kelvin) and Fourier number respectively at mesh points, s E z+, at time step n, n E z+. s runs from 1 to Nxw where s = 1 is the mesh point on the packing coke/flue wall boundary and s = Nxw is the mesh point on x = Lx. 48 As in the previous case, the coefficients of Tw~, Tp~xp-I and Twf ~ 0, so the Law of Conservation of Energy is satisfied. Summary The numerical solution of (6.1) requires the solution of I Taj+ 1 = Foa~ [ 2Ta~ + (F:af - 2 }•f], 2 = Foa~ [ Ta~+l + Ta~-1 q = 2 to Nxa-1, 3 TaN;a = Tpt1 , see Equation (6.15), 4 Tp~+l = Fop~ [ Tp;+l + Tp;_l r = 2 to Nxp-1, + (-l _ 2JTp~] Fop~ 5 Tp~:~ = Twr1 , see Equation (6.16), 6 T n+ 1 F n [T n T n Ws = OW s Ws+l + Ws-1 + (- 1 _ 2JTw~] Fown s s = 2 to Nxw-1, provided Foan ::; 1 1 to Nxa, q 2 ' q = Fopn ::; 1 1 to Nxp, 2 ' r = p Fown < .!_ s = 1 toNxw-1, s - 2' F n < 1 n > 1. OWNxw - Bin)' 2(1 + 49 6.3 Non-constant thermal conductivities Here the boundary value problem is (see §3.3) with oT t > 0, ax (0, t) = 0, (6.17) T(Lx, t) = T w(t), t> 0, T(x, 0) = T1, 0::; x::; LX. ( aTJ 2 aT a 2 T aT ox has to be approximated as well as at and ox2 . As in §6.2, ox is approximated by using a Taylor series expansion. oT (Lix.)2 o 2 T T(x + Lix, t) = T(x, t) + Lix ~ +-2 , - 2- + O(Lix)3 + ... ox . ox aT (Lix)2 a2 T T(x - Lix, t) = T(x, t) - Lix ox + 2! ox2 - O(Lix)3 + ... ( oT) 2 = [T(x + Lix, t) - T(x - Lix, t)] 2 ax 2Lix Th. . . -" (aT)2 h f d CA )2 IS approx1matlon 10r OX as error O Or er LlX . aT a2 T If at and ox2 are approximated as in §6.2, then the partial differential equation oT " (aTJ 2 o2 T . . d at = a OX + a dX2 IS approximate 50 by T(x, t + Lit) - T(x, t) Lit = & [T(x + Lix, t) - T(x - Lix, t)] 2 2Lix [ T(x + Lix, t) + T(x - Lix, t) - 2T(x, t)] +a 2 , (Lix) where the error of the left hand side is of order Lit, and the error of the right hand side is of order (Lix)2. Therefore, in subscript/superscript notation, where A &~ Llt Fo~ = 1 = the 'hatted' Fourier number of the block at I 4(LlX)2 mesh point i at time step n, in 1/K, and &f = the 'hatted' thermal diffusivity of the block at mesh point i at time step n, in rri2/sK, As in §6.2, Fof ::; ½. n - 2Lixhn ( n ) n and TNx+l = n TNx -Tf + TNx-1 kNx Substituting for T~x+l in Equation (6.18) gives at i = Nx: + Fo,\', [ 2Bi "Tf + 2Tri,_ 1 + ( F:ri, - 2( I+ Bi") JT~,] (6.19) 51 where As in §6.2, The boundary condition at i = 1 gives T ~ = Ti, as before. Substituting this into Equation ( 6.18) gives at i = 1: Tj+ 1 = Foj [TJ-T)j2 +Fof[T)+T)+(F:f - 2}r] = Fof [2T) +(F:f - 2 }r] F n < 1 with o1 _ 2 . Equations (6.18) and (6.19) involve non-linear terms. Because of this, it was decided to leave the non-constant thermal conductivities case at this point. 6.4 Results As input, the program requires an initial block temperature (Ta~, q = 1 to Nxa, Tp;, r = 1 to Nxp, Tw;, s = 1 to Nxw), an air temperature (Tf), the number of mesh points in the anode (Nxa), packing coke (Nxp) and flue wall (Nxw), the mass flow (rh) and finally the length of the time period. The mass flow is the total mass flow through all 4 fire shafts. The program's output is the transient temperature distribution of the block. The initial block temperature is assumed to be equal at all mesh points and is chosen as 600 °c. Tf is chosen to be 20 °c, Nxa = 15, Nxp = 4, Nxw = 4 and the length of the time period = 96 hours. The length of the overall time step (.6.t) is chosen to be 20 seconds. The number of mesh points coupled with the small .6.t is sufficient for convergence, that is, more mesh points and a smaller .6.t do not alter the temperatures significantly. Mesh point 1 corresponds to x = 0 and mesh point 21 to x = Lx . .6.xa, .6.xp and .6.xw are approximately the same - 0.028 to 0.032 m. 52 As shown in § 1. 1, Figure 1. 1, the fire trains being used at the New Zealand Aluminium Smelters Ltd have 3 forced cooling sections. Assuming a fire cycle time of 32 hours, the fans are moved in the fire direction by one section every 32 hours. Therefore each section undergoes 3 x 32 = 96 hours of forced cooling. Figures 6.4 - 6.6 show the temperature profile of the block at 32, 64 and 96 hours for various mass flows. The arrows in the figures indicate the mesh points on the anode/packing coke and packing coke/flue wall boundaries. ~~~~~~~~;g~=;=:.~· 450 .:i:. 400 Q 350 ~ .2 300 0 [ 250 E ~ 200 150 100 50 .J-.-------1--------1------_,_ ____ __; 6 11 Mesh points Figure 6.4 16 21 -■- l.5kg/s --o-- 5 kg/s 10 kg/s -O--l5kg/s Temperature profile of the block at 32 hours (1 fire cycle) 53 =t I 400 -•-•-•-•-•-•-•-•-• 'ff -·-·-·-•-· Q~or "' (]) •-+---+--+--•-•-•--· •, I .2 300 -·-·-♦-♦-♦-♦ "• e ~ ~ ~ 250 ·~ ·-~ ~~ ~ . 150 ~u----, 100 ~. 50 +------+-------;------t------------1 6 11 Mesh points Figure 6.5 16 21 -•- l.5kg/s -a-- 5kg/s 10 kg/s -o--15kg/s Temperature profile of the block at 64 hours (2 fire cycles) 500 450 400 Q 350 I ~ •-•-•-•-•-•-•-•-•-• • 'ff ~ 300 r - -•-•-•-·"' I ~ 250 ~ ~-u--,___,,__ '• y E +--+---+-+-+-+--+-•-• " ~ 200 C=::::§♦-♦ ♦-♦ • ·~ ·---■-. 150 ~ 100 . ~~ ~ 50' : : ' 0 l 6 11 16 21 Mesh points Figure 6.6 -•- l.5kg/s -a-- 5kg/s 10 kg/s -0--15kg/s Temperature profile of the block at 96 hours (3 fire cycles) 54 Figure 6. 7 shows the temperature profile of the block at 96 hours for different thermal conductivities, but a constant mass flow (5 kg/s). Firstly the thermal conductivities of the three materials (ka, kp and kw) are calculated as discussed in §6.2. Then ka is divided by 2, whilst kp and kw are left unaltered. Then kp is divided by 2, whilst ka and kw are left unaltered. Finally kw is divided by 2, whilst ka and kp are left unaltered. 150 100 50 -----;------+-----;------; 6 11 Mesh points Figure 6.7 16 21 -■- ka, kp, kw --a-- ka/2, kp, kw ka, kp/2, kw --<>-- ka, kp, kw/2 Temperature profile of the block at 96 hours with m= 5 kg/s 6.5 Discussion of results Figures 6.4 - 6.6 show that when the mass flows increase, the block cools more quickly. They also show that as time goes on and the temperature difference between the block and the air decreases, then the rate of cooling of the block also decreases. These results confirm that the set-up of the model and the calculations are correct. For this particular set of initial temperatures, it takes a mass flow of 15 kg/sin order for the anodes to reach approximately 200 °Cat 96 hours (3 fire cycles) - see Figure 6.6. This is similar to the mass flows used experimentally. At about this temperature, the anodes can be safely removed from the pits and stacked on the floor. 55 As expected, altering the thermal conductivities effects the cooling rate of the anodes, see Figure 6.7. The greatest effect on anode cooling occurs when kp is altered. Altering ka produces the least effect, whilst altering kw produces an intermediate effect. This result agrees with that of Furman and Martirena (1980). It was found that kp was one of the two parameters which decisively influenced the calculations, see § 1.2. 56 CHAPTER 7 TWO-DIMENSIONAL NUMERICAL SOLUTION 7.1 Introduction Only the constant thermal conductivities case is examined. The main difference, besides the increase in dimension, between this chapter and the previous one, is that here the fluid or air temperature in the flue is not constant. As mentioned in §1.3, the air temperature (and pressure) is changing as heat is transferred from the block out into the air in the flue. 7.2 Air Flow in the Flue An elemental volume (V) of air in the 'large' flue is considered, as shown in Figure 7.1. y air flow I __ ...,., I ➔ d I I I I I I I I - -, \ \ \ \ \ \ Figure 7.1 X Elemental volume of air in the flue V = d W /:iy m3, where d = 0.238 m and W = 0.111 m (see Figure 2.6). Mass of the air in V = pf x d x W x /:iy kg (pf = density of air or fluid in kg/m3 ) Energy of the air in V = mass x specific heat capacity x temperature = ((pf)dW l:iy) X cfp X Tf J 57 By the conservation of mass: The change in the mass of the air in V over time L~.t = the mass of the air entering V during time .M - the mass of the air leaving V during time Lit. i.e. ( (pf)dW Liy )r+ti.t - ( (pf)dW Liy )t = ( (pf)vdW)y Lit - ( (pf)vdW)y+ti.y Lit (v = the velocity of the air in m/s) ⇒ (pf)t+D.t - (pf)t + ((pf)v)y+li.y - ((pf)v)Y = 0 Lit Liy ⇒ a(pf) a - + -((pf)v) = 0 at ay (7.1) Similarly by the conservation of energy: The change in the energy of the air in V over time Lit = the energy or heat gained by the air in V from the block during time Lit + the energy of the air entering V during time Lit - the energy of the air leaving V during time Lit. i.e. ⇒ ⇒ ⇒ But ( (pf)dW Liy (cfp)Tf)t+ti.t - ( (pf)dW Liy (cfp)Tf)t = (Qw)d Liy Lit + ((pf)vdW(cfp)Tf)YLit - ( (pf)vdW(cfp)Tf)y+ti.y Lit a a a a pf at ( cf P Tf) + cf P Tf at (pf)+ (pf)v ay ( cf P Tf) + cf P Tf ay ( (pf)v) 1 . =wCQw) cfPTf gt (pf) + cfPTf :y ((pf)v) = cfPTf (tt (pf) + :y ((pf)v)) = 0, see Equation (7.1). a a 1 . Therefore pf at ( cf P Tf) + (pf)v ay ( cf P Tf) = W ( Qw) 58 and since (pf)v = ii then (7.2) f(Tfd(cfp)aTf + cf aTfJ+ m (Tfd(cfp)aTf + cf aTfJ= _!_(Qw) p dTf dt p dt dW dTf ay p ay w i.e. f aTf m aTf 1 · p at + dW dy = WF ( Qw) d(cfp) where F = Tf dTf + cfp. Th . aTf . I d . h h I . . f h . . e transient term, dt , 1s neg ecte smce t e t erma mertia o t e air 1s very small compared to the thermal inertia of the block, see Thibault et al. (1985). Therefore dTf d(Qw) ay = mF (7.3) From Equation (6.4), Qw = Qw(t) = h(t)[T(Lx, t) - Tf(t)] But for this two dimensional case, Qw = Qw(y, t) = h(y, t) [T(Lx, y, t) - Tf(y, t)] A . . aTf . T I . . . pprox1matmg ay usmg a ay or senes expansion gives aTf = Tf(y + £ly, t) -Tf(y- £ly, t) + O(£l )2 ay 2 £ly Y Equation (7.3) becomes: 59 Tf(y + t:,.y, t) -Tf(y - 1::,.y, t) d h(y, t)[T(Lx, y, t) - Tf(y, t)] = 2 t:,.y mF(y, t) where the error of the left hand side is of order (t:,.y)2. Therefore Tf(y + 1::,.y, t) 2 t:,.y d h(y, t)[T(Lx, y, t) - Tf(y, t)] ( 7 _ 4 ) = . F( ) + Tf(y - 1::,.y, t) m y, t Tf J is defined as the air temperature (in Kelvin) in the flue at mesh point j at time step n, and hj as the heat transfer coefficient (in W/mK2) at the flue wall at mesh point j at time step n, j and n E .z+. j runs from 1 to Ny, where j = 1 is the mesh point in the flue corresponding to y = 0 in the block, and j = Ny is the mesh point in the flue corresponding to y = LY in the block. 1::,.y is the distance between mesh points (in metres) in the y direction and F1} is the value of F at mesh point j in the flue at time step n. Tfj is defined as the block temperature (in Kelvin) at mesh point (i, j), where i and n are defined as in §6.2. j runs from 1 to Ny, and j = 1 corresponds to y = 0 in the block, and j = Ny corresponds toy= LY in the block. Rewriting Equation (7.4) using the subscript/superscript notation gives: = 2 1::,.y d hf (T~x,j - Tff) + Tf~ ri1 pl: J-1 J or Tf ~ = G~ 1 + Tf ~2, J J- J- where Gf_1 = m pl: 1 J- Tff-1) and (7.5) It is assumed that the inlet air temperature at time step n, Tf f , remains constant for all n. For the next mesh point, Tf ~ , Equation (7 .5) cannot be 60 used, since Tf f_2 = Tf g is an imaginary mesh point outside the flue. In this case, it is necessary to find a Taylor series approximation such that the mesh points are inside the flue. The 3 - 4 - 1 approximation outlined as follows is such an approximation. aTf (/1y)2 a2 Tf 3 Tf(y - /1y, t) = Tf(y, t) - /1y a + - 2 - --2 - O(/1y) + ... (7.6) Y ay aTf c211y)2 a2 Tf 3 Tf(y - 2/1y, t) = Tf(y, t) - 2 /1y a + 2 - 2- - O(/1y) + ... (7.7) Y ay 4 x (7.6) minus (7.7) gives: aTf 3 4Tf(y- /1y, t) - Tf(y- 2/1y, t) = 3Tf(y, t) - 2/1y ay + O(/1y) + ... ⇒ 3Tf(y, t) - 4Tf(y - /1y, t) + Tf(y - 2/1y, t) 2 + O(/1y) ... 2/1y Switching to subscript/superscript notation gives: aTf 3Tff - 4 Tff _1 + Tff_2 ay = 211y where the error of the right hand side is of order (11y)2. So in this case Equation (7.3) becomes: 3Tf1:1 - J 4Tff_1 + Tf1:1 2 J- d h1:1(T~ · - Tff) J X,J = 2/1y rh p1:1 J or 3Tff+1 4Tff + Tf1:1 1 d hf+1(T~x,j+l - Tff+1) J- = 2/1y . pn m j+l i.e. Tf n_ 1 on 3 Tf n 1 Tf n J = - 4 j+l + 4 j+l + 4 j-1 (7.8) and for j = 2, Equation (7 .8) becomes Tf ~ = - ¼ G~ + ¾ Tf ~ + ¼ Tf 7 (7.9) 61 But in Equation (7.5), the G j~l are calculated at the previous spatial step. For example, for Tf ~, G~ is calculated and for Tf~Y' G~y-l is calculated. But in Equation (7.9), the GT(= G~) is calculated at the next spatial step, since for Tf ~ , G~ is calculated. Also, note that G~ is used twice, once for calculating Tf 2 in Equation (7 .5) and again for calculating Tf ~ in Equation (7.9). So in order to be consistent, instead of calculating G ~ for Tf ~ in Equation (7.9), G 7 is calculated. Therefore Summary 1 2 3 where and Tf n = constant 1 Tf ~ = - ¼ G~ + ¾ Tf ~ + ¼ Tf ~ 2 11y d h1:1(T~ · -Tf1:1) G~ = J X,J J J m Fl: J d(cfp)i dTf = 0.138, see Equation (6.10). The calculation of Tf at the next time step (7.10) (7.11) (7.12) Relaxation is used and this is outlined as follows. Given the air temperatures along the flue at say time step n, that is given Tf J' j = 1 to Ny, the air temperatures along the flue at the next time step are estimated using Equations (7.11) and (7.12). Call these estimates estTf1J- It is not necessary to estimate Tf f since Tf 7 = constant for all n = atmospheric temperature = 20 °C ( see §2.4). 62 So estTf~ = - ¼ G f + ¾ Tf ~ + ¼ Tf f and estTfj = G }-I + Tfr-2 , j = 3 to Ny. Next the percentage differences between Tf f and est Tf j relative to Tf J are calculated for j = 2 to Ny. i.e. percentdiff1 jTff - estTff I = ~---~ X 100, j = 2 to Ny. Tfl:1 J If the maximum percentage difference for all j 2:: 0.1, then the above process is repeated, replacing the Tq with est Tff ,j = 2 to Ny. New estTf f and hence a new maximum percentage difference is calculated. This process keeps repeating itself until the maximum percentage difference < 0.1. When this inequality is satisfied, then Tf j+ 1 = estTf J , j = 2 to Ny. Once the Tfy+1 are calculated, then this enables the calculations in the block to be done for time step n + 1. After these temperature calculations in the block are completed, then the calculation of the T~ at the next time step begins again and so on. 7 .3 The heat equation in the block Recall from §3.3 that the boundary value problem is OT ( 0 2 T 0 2 T J 0 S XS Lx, at(x, y, t) = a ax2 + ay2 ' 0sysLy, t>0 with aT -(0, y, t) = 0, ax 0sysLy, t > 0, T(Lx, y, t) = T wCt), osy sLy, t > 0, (7.13) T(x, y, 0) = T1, 0sxsL, osysLy, aT aT 0 S XS Lx, t>0. - (x, 0, t) = ay (x, Ly, t) = 0, ay Discretising (7.13) as in Chapter 6, gives in subscript/superscript notation: 63 (7.14) i = 1 to Nx, j = 1 to Ny, n > 1, a,:1.~t where T~. and Fo~. = ~ 2 are the temperature and the Fourier number of the 1,J 1,J (~y) block respectively at mesh point (i, j) at time step n, with dT1n· __ ,_J dx 0, j = 1 to Ny, n > 1, dT~x,j - hf( n n) • = -- TN · - Tf · J = 1 to Ny, n > 1, dx kn x,J J ' Nx,j 1 T· · = T 1,J 1' i = 1 to Nx, j = 1 to Ny, dT:11 __ 1,_ dy 0, i = 1 to Nx, n > 1. Using Figure 7.2 as a guide and Equation (7.14), Ttr1 is now determined for all i and j, with the partial derivatives of the boundary conditions approximated as was done in §6.2. lines of symmetry / clT "" dY=O axes for a sucker y, I I I adiabatic: I clT : dX =0 I I I :Y adiabatic X T(x, y, t) Block axes for a blower .....,_ ____ ~ , x adiabatic I Figure 7.2 ' ' ~ue ' Tf(y, t) ' ' : ¥x= -~(T-Tf) Schematic cross-section of the block and flue showing axes and boundary conditions 64 i=l,j=l aTf · oT!\ The boundary conditions __ ,J = 0 and - 1 -' = 0 give T 3, 1 ax oy = Tf, 2 respectively, therefore provided Fof,1 :::; 1 i = 2 to Nx-1, j = 1 oT:11 1, n n The boundary condition dy = 0 gives Ti,O = Ti,z• Hence provided Tn+l i,l Fo:11 :::; I, Foti{(~:)' (Tf+1,1 + T:'-1,1) + 2Tt2 + [F~f.i - 2( (!)' + I) ]Tf1} (7.16) 1 i = Nx,j = 1 oT:11 - 1 -' = 0 gives TNx,O = TNx,z and as in §6.2 ay -h1:1 ( ) ~ TNx,j - Tff , j = 1 to Ny, Nx,j gives TNx+l,l - 2(Lix)hf ( n n) n T Nx,l - Tf1 + T Nx-1,1 kNx,l 65 [ 1 + --- F0Nx,1 (7.17) provided (&)h~ where Bif = -----"- = the Biot number at mesh point j = 1 at time step n. kN~,1 i = 1, j = 2 to Ny-1 aTf,i . n n ax = 0 gives To,j = T2,j. provided Fo1n· :::; ,J [ 1 + --- Fo1,j 1 i = 2 to Nx-1,j = 2 to Ny-1 provided [ 1 + --­ Fo!1, 1,J 66 (7.19) i = Nx, j = 2 to Ny-1 This is very similar to the case i = Nx, j = 1 studied earlier. Using this earlier case as a guide gives: [ 1 + ---- Fo~x,j provided 1 fu~JS-2-~-!-~_)_2_+_(_!_~_)_2_B_i_f_+_1~). i = 1, j = Ny oTn. oTn a~·1 = 0 and oyNy = 0 give T3,Ny = T~,Ny and TtNy+l = Tf,Ny-1 respectively. Tn+l F n l,Ny = 01,Ny 11y n n { 2 2(11x) T 2,Ny + 2Tr,Ny-l provided 1 Fol,Ny S ------ { (!:Y + 1) i = 2 to Nx-1, j = Ny oTrNy . Tn n oy = 0 gives i,Ny+l = Ti,Ny-1. [ 1 + n Foi,Ny 67 (7.21) (7.22) provided 1 i = Nx, j = Ny This is very similar to the case i = Nx, j = 1 studied earlier. Using this earlier case as a guide gives: [ (( )2 ( )2 1 b..y b..y •ll + n - 2 - + - B1Ny + FoNx,Ny b..X b..x (7.23) provided Determining b..y b..y is obtained by specifying the number of spatial steps required in they direction, Ny, and then the depth of the block is divided by Ny-1. ~ b..y = Ny-1 m Determining Fo~j A very similar method to that discussed in §6.2 is used. As in §6.2, the ki (thermal conductivies) are calculated using the temperature midway between the initial and assumed final average temperature of the block. Suppose this final average block temperature is Tf} , for any j between 2 and Ny, since all the T~1 , j -:f::. 1, are chosen to be equal. Recall from §2.4, that Tf ~ = atmospheric temperature = 20 °C for all n. Determining BiJ Similar to the method in §6.2, but here h = h(y, t), so 68 Lix h~ Bi~ = J = J kn Nx,j Lixw h~ kw where h~ = (hr~ + he~ ) W/m2K J J J hr~ = the heat transfer coefficient for radiation of the flue wall at J mesh point j at time step n, and he~ = the heat transfer coefficient for convection of the flue wall at J mesh point j at time step n. (kff )(Nuf) and he~ = J .e. kff = the thermal conductivity of the fluid or air at mesh point j in the flue at time step n, in W /mK, Nuj = the Nusselt number at mesh point j in the flue at time step n and is dimensionless, .e. = the characteristic length in m, see §6.2. Determining Lit A similar method to that used in §6.2 is followed. (a) Obtaining Lita i.e. (Liy )2 Lita ~ -2-aa_p_. (-'-(--'.0.'-'-y-)_2 _+_1_J . ,J Lixa 69 (Foaf,j is the Fourier number of the anode at mesh point (i, j) at time step n) ( aaf,j is the thermal diffusivity of the anode at mesh point (i, j) at time step n) The upper bound on Lita is obtained by calculating the maximum value that aa.0 1 . can take. I, max. aa!l. = ka , see §6.2 for this calculation. As in §6.2, the 1 ' 1 (pa)(min.(cap )r,} minimum block temperature is taken to be the assumed final average block temperature which equals TfJ for any j between 2 and Ny (see the discussion earlier in this section). (b) Obtaining Litp Similarly ( c) Obtaining Litw F n < . owi,j _ mm. 1 1 2(c;J + 1)' 2((:J +(~)'Bij +1] 1 = ( 2 2 J ' 2 (:w) + (Li;w) Bif+l since Bi J > 0, '7 j and n. i.e. Litw :s; (( )2 ( )2 J n Liy Liy ·n 2awi • -- + -- B1- +1 ,J Lixw Lixw J The upper bound on Litw that satisfies the above inequality is obtained by calculating the maximum value that aw!!. and Bi0 1 . can take. See 1,J §6.2 for the calculation of max. aw!1 1 .. I, Since Lixw h 1 1: Bi~ = 1 kw 70 the maximum value of Bi J is obtained by calculating the maximum value of hf. max. hf = max. hr'] + max. hcJ max. hr~ = £cr[max. T~xw,j - min. Tff] [(max. T~xw)2 - (min. TfT)2 ] = £cr[initial T Nxw,j - min Tf j] [(initial T Nxw)2 - (min. Tf J )2] = £CT [T~xw,j - min. Tf7] [(Tkxw)2 - (min. Tf j )2] Assume that min. Tf j , the minimum fluid temperature at mesh point j at time step n, is 293 K (atmospheric temperature) for all n. n _ (max. kf1 )(max. Nuf) max. hcj - _f, max. kfj = (5.65 x 10-5) max. Tff + 1.02 x 10-2 , see Equation (6.9). Assume that max. TfJ , the maximum fluid temperature at mesh point j at time step n, is 773 K for all n. max. NuJ = 0.023 (max. Ref )°-8 (max. Prf )°-4, 4 ri1 max. Ref = (min. µ J )(perimeter) (min. kf~)(min. Prf ) min. µJ = · · max. (cfP )J see Equation (6.6) see Equation (6.7) see Equation (6.8) max. Pr~= (1.89 x 10-7) (773)2-(2.20 x 10-4) 773 + 7.54 x 10-1 , J see Equation (6.11) min. kff = (5.65 X 10-5) 293 + 1.02 X 10-2 min. Prj = (1.89 x 10-7)(293)2- (2.20 x 10-4) 293 + 7.54 x 10-1 max. (cfp)j = (1.38 x 10-1)773 +9.92 x 102, see Equation (6.10). 71 (d) Obtaining the overall ~t. The upper bound on the overall ~t = minimum (max. ~ta, max. ~tp, max. ~tw). Temperature on the anode/packing coke boundary Ta~,j and Foa ~,j are defined as the carbon anode temperature (in Kelvin) and Fourier number respectively at mesh point ( q, j) at time step n, q is defined as in §6.2. j runs from 1 to Ny- Similarly Tp~j and Fop~j are defined as the packing coke temperature (in Kelvin) and Fourier number respectively at mesh point (r, j), at time step n. r is defined as in §6.2. At q = Nxa. From Equation (7 .16) Tn+l Nxa,l n [ 1 + 2TaNxa,2 + n - FoaNxa,l (7.24) and from Equation (7 .19), for j = 2 to Ny - 1, Ta~xa,j n [ 1 + TaNxa,j-1 + n - FoaNxa,j 2( (!J + I) ]T•~x..i} (7.25) and from Equation (7.22) Ta~-:;,Ny = Foa~x,,Ny{ (::.r (Ta~xa+l,Ny + Ta~x,-1,Ny) 72 n [ 1 + 2TaNxa,Ny-l + --n-- FoaNxa,Ny (7.26) Ta~xa+l,j, j = 1 to Ny, is an imaginary mesh point outside the anode, which coincides with Tptj only if Llxa = Llxp. At r = 1. From Equation (7 .16) (7.27) and from Equation (7.19), for j = 2 to Ny- 1, + [- 1 _ 2((~J 2 + 1J]rpr.} Fopn . Llxp ,J l,J (7.28) and from Equation (7.22) + [ ln - 2((~J 2 + 1J]Tpf,Ny} (7.29) Fopl,Ny Llxp Tpg,j, j = 1 to Ny, is an imaginary mesh point outside the packing coke which coincides with Ta~xa-l ,j only if Llxa = Llxp. 73 As was done in §6.2 , a heat balance is done on the anode/packing coke boundary using Fourier's Law to give: - ka xa+ ,J xa- ,J = - k ,J ,J · = 1 to N (7.30) ( TanN 1 . - TanN 1 . ] (TP2n. - Tpno.] 2L1xa p 2L1xp ' J y' j = 1. Equations (7.24), (7.27) and (7.30) involve three unknowns: Ta~xa+l,l , Tp8, 1 and Ta~:;, 1 = Tpf~1. These equations are solved for Tpt/, using the fact that Ta~xa, 1 = Tpr, 1 and Ta~xa,2 = Tpf,2 . Tpn+l = Fopf.1{2( ~ )\ n + 2 [( ~ r kaL\xp ]Ta" 1,1 A L1 P2,1 L1 k L1 Nxa-1,1 1 xp xp p xa [ ka L\xa ] n + 2 --- + 1 Tp1 2 kp L1xp ' [ ka L\xa( I + kp L1xp FoaNxa,l - {U:J + 1JJ + [ Fo~P,1 - 2l U:J + 1J J ]TPP,1} (7.31) where A1 = 1 + Fopf,1 ka L1xa FoaNxa,l kp L1xp In order to satisfy the Law of Conservation of Energy it is required that the coefficients of Tp~, 1, Ta~xa-l,l' Tpf,2 and Tp?, 1 2:: 0. ka and kp > 0. aa~xa,l and ap1,1 > 0 and hence Foa~xa,l and Fop1,1 > 0. n 1 n FoaNxa 1 ~ ( 2 J and Fop11 , 2 ( L1y) + 1 , L1xa the restrictions determined earlier. 74 Therefore the coefficients :2:: 0. j = 2 to Ny- 1 Similarly Tpn-!-1 l,J [ ka .6.xa ] n [ka .6.xa ] n + --- + 1 Tplj+l + --- + 1 Tpl j-1 kp .6.xp ' kp .6.xp ' [ ka .6.xa ( 1 2(( .6. Y ) 2 1]] + kp .6.xp Foa Nxa,j - .6.xa + [ 1 + n Fop1,j (7.32) where A2 Fopf,j ka .6.xa = 1 + --~ --- FoaNxa,j kp .6.xp j =Ny Similarly Tpn+l _ l,Ny - Fopf,Ny {2(_EL_J 2 T n + 2[(_EL_J 2 ka .6.xp]Tan A A P2,Ny A k A Nxa-1,Ny 3 D.Xp D.Xp p D.Xa [ ka .6.xa ] n + 2 --- + 1 TplNy-1 kp .6.xp ' + [ ka .6.xa ( 1 _ 2]((EL)2 + 1]] kp .6.xp Foa Nxa,Ny .6.xa + [ ln - 2((E'LJ 2 + 1J]lTpf,Ny} (7.33) Fopl,Ny .6.xp + Fopf_Ny ka .6.xa where A3 = 1 ---"-'-'-'-'- k Foa~xa,Ny P .6.xp 75 Temperature on the packing coke/flue wall boundary This case is similar to the anode/packing coke boundary one. Twns. and Fow~ 1 . are defined as the flue wall temperature (in Kelvin) and ,J , Fourier number respectively at mesh point (s, j), at time step n. sis defined as in §6.2. j = 1 Twn+l _ Fowf.1 { 2( Liy ) 2 Twn + 2[( Liy ) 2 kp LixwJT n 1,1 - 0 A 2,1 A k A PNxp-1,1 1 D.XW D.XW W D.Xp [ kp Lixp ] n + 2 --- + 1 Tw1 2 kw Lixw ' kp Lixp 1 2 Liy l [ [( ) 2 ]] + [ kw Llxw FopN,p,I - Llxp + [ 1 + n Fow1,1 (7.34) where 0 1 Fowf,1 kp Lixp 1 + FopNxp,l kw Lixw j = 2 to Ny- 1 where Twn+l _ Fowf,j {2( Liy ) 2 T n 2[( Liy ) 2 kp Lixw]T n 1,1· - ---=- -- W2 · + -- ---- PNx -1 · 0 2 Lixw ,J Lixw kw Lixp P ,J + [ kp Lixp + 1]Twr ·+1 + [ kp Lixp + 1]Twr •-1 kw Lixw ,J kw Lixw ,J I kp Lixp [ 1 + l kw Lixw FopNxp,j - {(::J + 1]] 2(U;J + 1)J]Twf,j} (7.35) [ 1 + ---­ Fow!1. 1,J Fowf j kp Lixp = 1 + ----=--- --- Fopn kw Axw Nxp,j D. 76 j = Ny T n+l = Fowf,Ny {2( L'.iy )2 Twn + 2 [( L'.iy )2 kp L'.ixwJT n WNy 0 A.. 2,Ny A.. k A PNxp-1,Ny 3 L.J.A W L.J.A W W LlXp where [ kp L'.ixp ] n + 2 --- + 1 TwlNy-1 kw L'.ixw ' + ---- ( 1 Fowf,Ny Fowf,Ny kp L'.ixp 03=1 + ---- Fop~xp,Ny kw L'.ixw (7.36) For the same reasons as before, the coefficients of the various mesh points in Equations (7 .34) - (7 .36) ~ 0, so the Law of Conservation of Energy is satisfied. Summfil:):'. The numerical solution of (7 .13) requires the solution of: (1) Tan+l Foaf 1 { 2( /1y )' Ta~ 1 +2Taf 2 1,1 ' L'.ixa ' ' +[-1 -2((~)'+1J}ari}, Foan L'.ixa ' 1,1 (2) Tan+l q,1 = Foal,, {(::J (ra~+I,I + Ta;_,,,) + 2Tal,2 + [ 1, - 2( ( ~ )' + I J ]ra;,1} , q = 2 Foaq,l L'.ixa to Nxa-1, (3) T n+l Tpf,t1 , see Equation (7.31), aNxa I = , 77 (4) Tpn+l r,1 = Fop~1{( /:J.y J 2 (TP~+11 + Tp~-11) + 2Tp~2 , 1:J.xp , , , + [ 1 n - 2((~J 2 + 1J]TP~,I}, r = 2 to Nxp- 1, Fopr,l /:J.xp (5) Twrf1 , see Equation (7.34), , (6) Tw~j 1 = Fow~,1 { ( :w r ( Tw~ + 1,1 + Tw~-1,1) + 2Tw~.2 + [ 1 n - 2((~) 2 + 1J]Tw~,1}, s = 2 to Nxw - 1, Fow5,1 /:J.xw (7) Tw~:~.1 = Fow~xw,I { ( ... ':w r [ Bif Tff + Tw~xw-1,1] + 2TwNxw,2 + [ ~ - 2(( /:J.y ) 2 FowNxw,l /:J.xw + (:JBi( + 1J]T~xw,1}, (8) Ta(,r1 = Foa~;{2(::.)' Tat; + Taf,j+I + Ta~j-l (10) T n+I aNxa,j + [-l -2((~) 2 + 1J]Tar ·}, J = 2 to Ny - 1, Foan • /:J.xa ,J 1,J [ 1 + n Foaq,j - 2( (::J + 1J] Ta~,i }• q = 2 toNxa-1, j = 2 to Ny- 1, = Tpr,j1 , see Equation (7 .32), 78 + [ 1 n. - 2((~J 2 + 1J]TP~,j}, r = 2 to Nxp- 1, Fopr,J ~xp j = 2 to Ny- 1, (12) T n+l - Twn+l E . (7 35) PNxp,j - 1,j , see quat10n . , (14) T n+l WNxw,j [ 1 + n Fows,j -2((~)2 + ,)]Tw~+ s = 2toNxw-1, j = 2 to Ny 1, n n [ 1 + TwNxw,j+l + TwNxw,j-1 + n FowNxw,j j = 2 to Ny-1, (15) Ta/,l:/y = Foa/,Ny{ 2(:.J TatNy + 2Taf.Ny-l + [ In - 2((~) 2 + 1)]TaP,Ny} , Foal,Ny ~xa 79 (16) (17) (18) Tan+l q,Ny = Foa~,Ny{ (:ar (Ta~+I,Ny + Ta:-1,Ny) T n+l aNxa,Ny n [ 1 + 2Taq,Ny-1 + n Foaq,Ny -2((::J + 1)]Ta~,Ny}, q = 2 to Nxa - 1, = Tpf,Uy , see Equation (7.33), Tpn+l r,Ny = Fop;,Ny{( :p J (Tp;+l,Ny + Tp;_I,Ny) n [ 1 + 2TPr,Ny + 1 F0Pr,Ny -{U:J + 1J]Tp;,Ny}, r = 2 to Nxp - 1 , (19) TpNt~,Ny = Twf,t~, see Equation (7.36), (20) (21) Twn+l s,Ny = Fow~Ny{ c,;w )'{T:'i-1,Ny + Tw~-1,Ny) T n+l WNxw,Ny n [ 1 + 2Tw s,Ny-1 + n Fows,Ny -2((:J + 1)]Tw~,N+ s=2 to Nxw-1, n [ 1 + 2TwNxw,Ny-l + --0-- FowNxw,Ny - 2((-1:!,.y ) 2 + (-1:!,.y ) 2 BiN + l:!,.xw l:!,.xw Y I) ]Tw~xw,Ny} , provided Foa~.j :s; 1 q = 1 to Nxa, 80 Fop~,j 1 :s: 2[(:J f + r = 1 to Nxp Fow~,j 1 :s: {(d;J f + s = 1 to Nxw- 1 Fow~xw,j :S: 7 .4 Heat balance at the flue wall Since there is no heat flow across the boundaries x = 0 y = 0 and y = L ' Y' then all the heat flow from the block occurs along the boundary x = Lx , i.e. the flue wall. The heat transferred from the block across this boundary should equal the heat gained by the air in the flue, i.e. LiQ + LiQf = 0, where LiQf = the heat gained by the fluid or air in time Lit and LiQ = the heat transferred from the block in time Lit. / axes for a sucker / / Figure 7.3 plane of symmetry / / / ✓ / / / / / / ✓ Schematic three-dimensional view of the block showing axes Calculation of LiQf Consider the block, where the z direction is as shown in Figure 7.3. Let Q(t) be the heat in the block at time t. Then 81 see §3.1, = 1t ( J~' f ~y J; pe dz dy dx) where dV = dx dy dz and d is as shown in Figure 7.3, dQ rL, rLY a dt = d Jo Jo at (pe) dy dx (pe is independent of z) = d f ~' f ~y (P cp i~) dy dx see §3.1, = d f ~' f ~y V ·(kVT) dy dx see §3.1, = d iL, iLY [j_(k aT) + j_(k aT)] dy dx 0 0 ax ax ay ay dQ iL · - = d y -Qw dy dt o (7.37) . aT aT smce ax = 0 on x = 0, ay = 0 on y = 0, LY and the rate of heat flow per aT unit area of the flue wall = Qw = k ax on x = Lx , see §6.2. Now consider the air. Recall Equation (7.2) that is, ⇒ 82 ( d( cf 0) oTf oTfJ m a 1 . ⇒ pf Tf dTf dt + cfp dt) + dW oy (cfPTf) = W (Qw) ⇒ oTf m a 1 . dt + dWF(pf) oy (cfPTf) = WF(pf) (Qw) d(cfp) where F = Tf dTf + cf P. (7.38) As discussed in §7.2, the transient term o~f is neglected. Therefore Equation (7.38) becomes . d . m dy (cfPTf) = d(Qw). Integrating with respect to y gives: = - dd~ from Equation (7 .37) Calculation of LiQ Consider an elemental volume V of block, as illustrated in Figure 7.4, centred at mesh point (i + ½, j + i), i = 1 to Nx - 1, j = 1 to Ny - 1. Ny-I Nx-1 [ ] Then LiQ = ~ ~ Vi+½,j+½ (pe)~:½,j+½ - (pe)~+½,j+½ J=l 1=! where V. ½. ½ = Lix Liy d m3 for all i andj, see Figure 7.4, l+ 2,J+ 2 and e = internal energy per unit mass in J/kg - see §3.1. 83 , , , , , plane of symmetry , , Aue wall Figure 7.4 Schematic three-dimensional view of the block showing elemental volume V centred at mesh point (i + ½, j + ½) p = pa, pp or pw and e = ea, ep or ew depending on whether the calculations are occurring in the anode, packing coke or flue wall. ea, ep and ew are the internal energies per unit mass of the anode, packing coke and flue wall respectively. oh oe Recall from §3.1, that cpCT, p) = oT = oT at constant pressure. Therefore e(T, p) = J: cpC-t, p)d"'C + f(p) where f(p) is some function of pressure. If cp(O, p) is assumed to be zero, then e(O, p) = f(p). Then assuming that e(O, p) = 0, gives e(T, p) = s: cp("'C, p)d"'C. Now consider cap = (1.3l)"'C + 4.55 x 102 , see Equation (4.3) Therefore ea = s: [(l.3l)"'C + 4.55 x 102 ]d"'C = (1] 1 )T2 + (4.55 x 102 )T J/lcg 84 Similarly ep = (1J8 )T2 + (5.27 x 102 )T J/kg ew = 1250 T J/kg and .6.Q = .6.Qa + .6.Qp + .6.Qw where .6.Qa = the heat transferred from the anode in time Lit in J, .6.Qp = the heat transferred from the packing coke in time Lit in J, and .6.Qw = the heat transferred from the flue in time Lit in J. Consider .6.Qa Ny-I Nx-1 ( J .6.Qa = (Lixa)(Liy)d(pa) L, L, ean++l ·+! - ean 1 . 1 j=l i=l q 2 ,J 2 q+2,J+2 ean+I1 . 1 and ean 1 . 1 can not be easily calculated. They are approximated q+2,J+2 q+2,J+2 using the average of ean+I and ean calculated at mesh points (q, j), (q + 1, j), (q, j + 1) and (q + 1, j + 1) - see Figure 7.5. V (q,j+l).~------,f, ,.- (q+l,j+l) , , , , , ,'J I •,(>(+2, i 2l , , ~ ,_' ___________ _ (q+l,j) 2 ,'( .j) ~+i 2 , 2 , 't:,J_ ,' 2 Figure 7.5 A more detailed view of V in the anode showing mesh points assuming that a blower is operating 85 ean+I 1 . 1 q+-J+- 2' 2 ean 1 . 1 "'" 0.25 [(ea~j 1 - ea~,j) + (ea~!Lj - ea~+I,j) q+-J+- 2' 2 ( n+l n ) ( n+l n )] + eaq,j+I - eaq,j+I + eaq+I,j+I - eaq+l,j+l Therefore LiQa "'" (Lixa)(Liy)d(pa)(0.25){ ( eatt1 - eat1) + ( eaNtt1 - eaNxa,l) + 2 [ Nq~=a2-l((eaqn+,11 n ) ( n+l n )) L.J - eaq,l + eaq,Ny - eaq,Ny Similarly for LiQp and LiQw. 7.5 Results The program requires similar input to the one for the one-dimensional case, but here the air temperature required to be entered is the initial air temperature (for the one-dimensional case, the air temperature was assumed to be constant). As discussed before (see §7.3), the TfJ are chosen to be equal for allj = 2 to Ny, with Tf? = atmospheric temperature = 20 ° C for all n. Also the number of mesh points in they direction (Ny) has to be entered. The program's output is the transient temperature distribution of the block and of the air in the flue. The initial block temperatures (Ta~,j' Tp;,j, Twl) are chosen to be 600 °C. The initial air temperatures (Tf} , j = 2 to Ny) are chosen to be 100 °c. Nxa = 15, Nxp = 4, Nxw = 4 and Ny= 160. Mesh point (1, 1) corresponds to x = 0, y = 0; mesh point (1, 160) corresponds to x = 0, y =LY; mesh point (21, 1) corresponds to x = Lx , y = 0 and mesh point (21, 160) corresponds to x = Lx, y =LY. The total number of mesh points is 3360. The computer used for the calculations could not handle too many more than this number of mesh 86 points. Lixa, Lixp, Lixw and Liy are approximately the same - 0.028 to 0.032 m. The length of the overall time step (Lit) is 20 seconds, as in the one-dimensional case. Either a blowing or sucking fan is operating, since as has been discussed before (see §2.4), it makes no difference to the temperature profile calculated by this model, as to whether a blower or a sucker is used. Figures 7.6 and 7.7 show the temperature profile of the block in the x direction at y = 0 and y = LY respectively at 32 hours for various mass flows. Figure 7 .8 shows the corresponding temperature profile of the air in the flue at 32 hours for the same mass flows. Similarly for Figures 7 .9 - 7 .11, except that time is at 64 hours, and for Figures 7 .12 - 7 .14 except that time is at 96 hours. The arrows in the figures indicate the mesh points on the anode/packing coke and packing coke/flue wall boundaries. Heat balance at the flue wall The percentage difference between L\Q and L\Qf relative to L\Q is calculated at . L\Q- L\Qf the end of each time step. Percentage difference = ---- x 100. L\Q The percentage difference between total L\Q and total L\Qf relative to total L\Q is also calculated at the end of each time step. Total L\Q (in J) is the total heat transferred from the block from n = 1 up to the present time. Similarly total L\Qf (in J) is the total heat transferred into the fluid or air from n = 1 up to the present time. . total L\Q - total L\Qf Total percentage difference = ------- x 100. In all cases, even at total L\Q 96 hours, both the percentage difference and total percentage difference< 0.5%. 7.6 Discussion of results The temperature profile of the block and of the air in the flue The results here are very similar to those of the one-dimensional model. As the mass flows increase, the block cools more quickly, and as time goes on and the temperature difference between the block and the air decreases, then so does the 87 rate of cooling of the block. As before, these results confirm that the set-up of the model and the calculations are correct. For this particular set of initial temperatures, it takes a mass flow of 15 kg/sin order for the anodes to reach approximately 200 ° C at 96 hours. This mass flow compares well with the experimental ones (5.6 kg/s for a blower, 11.3 kg/s for a sucker - see §6.2). As expected, the temperature of the block in the y direction increases as y goes from Oto LY - compare Figures 7.6 and 7.7, Figures 7.9 and 7.10, and Figures 7 .12 and 7 .13. This is because the air temperature in the flue is also increasing as y goes from O to LY (see Figures 7 .8, 7 .11 and 7 .14 ). However this vertical ( or y) temperature variation of the block is quite small, even for the smallest mass flow of 5 kg/s. This is because the air temperature does not increase greatly above the inlet air temperature (20 °C), even as the air flows along the flue and gains heat from the block. See Figures 7 .8, 7 .11 and 7 .14. This is even more so as time goes on and less heat is being transferred from the block, due to the block cooling that has already occurred. This is due to the fact that the air is being sucked or blown through the flue at a rate that enables a rapid replacement ( especially with larger mass flows) of the air in the flue. Hence, a particular quantity of air is not remaining in the flue long enough for it to gain much heat, and therefore increase in temperature. The transient variation of the air temperature is quite small and as time goes on, becomes smaller. See Figures 7.8, 7.11 and 7.14. This is consistent with the assumption made in §7.2, that the transient or t~ term in the partial differential equation for the air temperature can be ignored (since the thermal inertia of the air is so much less than that of the block). Heat balance at the flue wall The small discrepancy, especially between the cumulative heat given out by the block and the cumulative heat transferred into the fluid, indicates that the calculations in the block and air are correctly related and quite accurate. 88 500 ;--•-•--•-•-•~ _ I ~ ~---· 'Ill" -+-+-♦- -■ 450 ~-•-•-•-~~-••• ., ·-·-•-. J\ 400 -♦--~\ ,.... 350 () ; 300 5 o 250