Iraqi Journal of Chemical and Petroleum Engineering Vol.16 No.4 (December 2015) 67- 78 ISSN: 1997-4884 Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method Tahseen Ali Al-Hattab * , Ali Abd Al-Muhsen Al-Assady ** and Wisam Ahmed Abd Al- Wahid *** * University of Babylon, ** University of Baghdad, *** University of Kufa Abstract The temperature distributions are to be evaluated for the furnace of Al-Mussaib power plant. Monte Carlo simulation procedure is used to evaluate the radiation heat transfer inside the furnace, where the radiative transfer is the most important process occurring there. Weighted sum of gray-gases model is used to evaluate the radiative properties of the non gray gas in the enclosure. The energy balance equations are applied for each gas, and surface zones, and by solving these equations, both the temperature, and the heat flux are found. Good degree of accuracy has been obtained, when comparing the results obtained by the simulation with the data of the designing company, and the data obtained by the zonal method. In the above work, a code of Monte Carlo method is built, and used to overcome the mathematical analysis. Key Words: temperature distributions, degree of accuracy Introduction The representation of the heat transfer processes, inside the industrial furnaces, have taken a very big attention in the literatures, where in the past two simplified cases are used, the well stirred enclosure which assume perfectly mixed combustion products, and the long furnace model which assume plug flow of combustion products with no radial temperature gradient. In order to predict the accurate temperature distribution in the enclosure, the zone method have been developed. Because that the zone method can deal only with furnaces of simple shapes, and that it is suffering from a lack to evaluate the total exchange area between zones can not see each other, Monte Carlo method introduced as a solution for this problem. Monte Carlo method was first reported in 1964 [1], and the text book by Siegel, and Howell [2] presented the fundamental concepts of the method, and some examples analyzing simple systems. All the later literatures using Monte Carlo method are based on the same concepts with very much modification to the method to represent the conditions of the problem. There are many studies which studying the radiative transfer inside furnaces, but a few studies selected and presented here, which show the application of the method on furnaces, where Vercammen, and Iraqi Journal of Chemical and Petroleum Engineering University of Baghdad College of Engineering Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 68 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net Froment [3], introduce an improvement to the zone method, where the concepts of the method still the same, but the total exchange areas among the zones are obtained by Monte Carlo method. Steward, and Guruz [4] apply Monte Carlo method to evaluate the total exchange area among zones of an industrial furnace, with a non gray medium contains a particle radiation, and by neglecting the scatter. Gupta, Wall, and True love [5] apply Monte Carlo method to a cylindrical furnace, both isotropic, and anisotropic scattering mediums are examined. Taniguchi, Kudo, and Sasaki [6], apply Monte Carlo method to study two kinds of furnaces, with two, and three dimensional radiative analysis, both luminous, and non luminous flames are studied in the furnaces. System Description The furnace under consideration is the boiler of Al-Mussaib power plant, which is a Babcock and Wilcox radiant boiler. The combustion chamber is approximately 20m in height, 13m in width, and 11m in length. Steady State Energy Equations After dividing the system into (154) surface zones, and (116) volume zones, each zone will be assumed small enough to consider as isothermal. Then, the steady state energy equation is applied to each zone, to evaluate the temperature, and heat flux distributions, after solving them simultaneously: infhinroutfcgoutr QQQQQQ ,,,,  … (1) cinrdsoutr QQQQ  ,, … (2) where, equation (1) is applied to volume zones, and equation (2) is applied to surface zones. The term Qr refers to radiative heat transfer, Qc is convective heat transfer, Qf is advective heat transfer, Qh is heat generation, and Qd is the heat flux at the surface elements. The terms of equations (1), and (2) can be obtained from: VTkQ ggoutr 4 , 4  … (3) ATQ ssoutr 4 ,  … (4) soutrgoutrinr QSXQGXQ ,,,   … (5)   sgc TThAQ  ... (6) upgppginf TcvmQ ,,  … (7) gpgoutf TcmQ  , … (8) The total exchange area, in the present work, will be obtained by Monte Carlo method. Monte Carlo Method Monte Carlo method is a statistical procedure, which simulates the radiation events (absorption, and emission), in terms of random numbers distributed in the range (0-1). In this method, and for each zone, a finite number of energy bundles will be defined, every bundle emitted from the zone will be governed by the azimuthal angle υ, and the polar angle θ, which are given by 1 2 R  … (9)  zonessurfaceforR 2 cos  …(10)  zonesgasforR 3 21cos  … (11) Tahseen Ali Al-Hattab, Ali Abd Al-Muhsen Al-Assady and Wisam Ahmed Abd Al-Wahid -Available online at: www.iasj.net IJCPE Vol.16 No.4 (Dec. 2015) 69 The probable mean free path is given by:   n RKl ln1 … (12) If the bundle traveled the mean free length (1) and it was less than the maximum path (L), the bundle will not be intercepted by wall, then the bundle is said to be absorbed by the gas, and tally there, if not, another condition applied, which is 5 R … (13) If the above condition is applied, then the bundle is said to be absorbed by the surface, and the bundle is tallied there, and if the condition is not satisfied, then the bundle is said to be reflected, and the point of reflection will be treated as emission point, previous procedure will be repeated for the new emission point. This procedure is repeated for each bundle till all the history of bundles is tallied. This procedure can be represented by figure (l). The total exchange area can be evaluated for each zone, by dividing the number of bundles reached to each zone, by the total number of bundles emitted. Evaluation of the Radiative Properties of the Non Gray Gas The gas fills the boiler will be assumed as non gray gas, and the weighted sum of gray gases model will be used to evaluate the total emissivity, and total absorptivity for each zone according to that zone temperature. Total emissivity for the weighted sum of gray gases model is evaluated from the following expression   kips l i i eTa     1 0 , … (14) The weighting factor for i=0, is evaluated from:    l i i aa 1 ,0, 1  … (15) A convenient representation of the temperature dependency of the weighting factors is a polynomial of order J-l given as follows: 1 1 ,,,    j j j jii Tba  … (16) where bε,i,j can be evaluated from the data of appendix (A). To evaluate the total absorptivity, the irradiation temperature of surface surrounding the gas is also introduced. Hence:   kips s l i i eTTa     1, 0 , … (17) Where 1 1 1 1 ,,,,              j j j k k k skjii TTCa  … (18) again the value of kji C ,,, can be found from appendix (A) . The weighting factors for i=0 is given as:    l i i aa 1 ,0, 1  … (19) Three Gray-Gases Model It has seen in equation (12), that the mean free path is a function of the absorption coefficients, which has a variable value from zone to zone, according to that zone's temperature. In the present work three gray-gases model will be used to represent the non gray-gas included in the system, which Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 70 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net means that the gas is consisting of three gray gases, and clear part. For each gray gas, a different value of aborption coefficient will be used, then four values of total exchange area will be kept in the calculations, Hence, the values of total exchange area will be expressed as:     N n niiinii XYTaXY 0 )( … (20) where, an is the weighting factor at the gray gas i, which is obtained in (16). Results and Discussions It is important now to apply the radiation heat transfer to evaluate the temperature of each gas zone. To do this, the energy balance equations for each gas zone should be solved, which contains terms of convection, advection, and combustion rate, as well as the terms of radiative transfer terms. After solving the energy balance equations, the temperature distribution is determined. Figures (2), and (3) show the temperature distribution at (70%) load of the furnace for the first layer beside the front wall (K=l), and the core layer (K=2), respectively. The figures show that, in general, the second row (I=2) have a general values of temperature less than the first row, this is because that the first row (I=1) is a combustion region, where as the second row is not. The temperatures increase in the third row due to heat liberation by combustion, which will has temperature values higher than, such of zones. The temperature values will decrease when the flue gases going up until reaching its lower value in the furnace just before leaving the furnace. The reason behind the temperature values decreasing when the flue gases going up, is due to the loosing of heat both by radiation and convection to the walls of the furnace. Because of the strong effect of the advection from bottom rows to the upper rows, it can be seen that for each row the middle zone temperature is lower than other zones in the same row, where the zones at the sides are the higher temperature zones. It is also can be seen that the right side of the furnace is slightly higher than the left side, due to the effect of the neck part, which represent the top four gas zones on the left sides, where the flue gases are leaving the furnace. This part of the furnace is colder than other parts, which means it absorbs heat by radiation from the left side of the furnace, more than the right side due to the neighboorhood of that zones to this side. One can also notice that, the temperature values of the core layer (K=2), are higher than the values of the first layer (K=l), this is because of the effect of the front wall, which absorb heat by convection, where the core layer has no convection heat transfer to wall, only on the sides. To check the validity of the previous results, first, the outlet temperature of the flue gases will be checked with the test values of the designing company. The comparisons show that, the outlet temperature of the flue gases at (70%) load is (1090C 0 ), where the rate of the temperature leaving the furnace shown in figures (2) and (3), is (1092C 0 ), then there is a percentage of error of (0.18%), which is a good percentage, and the calculation is said to be good. Another comparison will be made is with Al-Habbubi [7], where the results are obtained by using the zonal method as shown in figures (4) and (5). In general the distribution, and the variation of the distribution, has agreement. The difference is at the local values of the layers, where in the Tahseen Ali Al-Hattab, Ali Abd Al-Muhsen Al-Assady and Wisam Ahmed Abd Al-Wahid -Available online at: www.iasj.net IJCPE Vol.16 No.4 (Dec. 2015) 71 first layer (K=1), the temperature values are higher than, the first layer (K=l), obtained by Monte Carlo method, where as for the opposite is evaluated at the core layer (K=2). The results show that Monte Carlo method has the higher temperatures than of the zonal method. The reason behind that is due to the assumption used, where at [7], half of the furnace is used in the calculations, but in the present work, the whole furnace is used, which make a difference in energy distribution of layers. This in general made the difference between the layers (K=l, and K=2) at Monte Carlo method layer and at the zonal method. Then, the variation in the furnace sizing use made this difference of calculations. Fig. 1, Flow chart of Monte Carlo procedure Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 72 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net Fig. 2, Temperature distribution for gas zones in (k) at layer beside the wall and for 70% load. Tahseen Ali Al-Hattab, Ali Abd Al-Muhsen Al-Assady and Wisam Ahmed Abd Al-Wahid -Available online at: www.iasj.net IJCPE Vol.16 No.4 (Dec. 2015) 73 Fig. 3, Temperature distribution for gas zones in (k) at core layer and for 70% load. Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 74 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net Fig. 4, Temperature distribution (K) for 70% load. Temperature in figure is for layer beside the wall volume zone from [7] Tahseen Ali Al-Hattab, Ali Abd Al-Muhsen Al-Assady and Wisam Ahmed Abd Al-Wahid -Available online at: www.iasj.net IJCPE Vol.16 No.4 (Dec. 2015) 75 Fig. 5, Temperature distribution (K) for 70% load. Temperature in figure is for core layer volume zones from [7]. Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 76 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net Conclusions A three-dimensional Monte Carlo procedure is used to predict the temperature and the heat flux distributions inside the furnace of Al- Mussaib power plant with verious loads by using the available data at the field. From the previous discussions it can be concluded that: 1- The results of normal working conditions give a good percentage of error when comparing this results with the testing data for the focused furnace, which means that Monte Carlo method gives a good results for radiation transfer analysis, and will give a good accurate results for temperature and heat flux distribution, only when applying a good assumptions for other physical processes occur in the furnace. 2- The weighted sum of gray gases model gives a good results to represent the emissivity, and the absorbtivity coefficients of non-gray gas. 3- The postulated flow pattern is a very important character, which effects mainly on the final data. 4- Complex geometries can be studied to evaluate it's total exchange areas for zones, by using Monte Carlo method without much difficulties. 5- Monte Carlo method is suffering from very high computational time, but it still an attractive method. List of Symbols Symbol Description Units a Weighted factor A Area m 2 b Coefficient of the emissivity polynomial c Coefficient of the absorptivity polynomial Cp Specific heat coefficient J/Kg.K D Diameter m GG The total radiative interchange area between two volume zones m 2 GS The total radiative interchange area between any volume zones and any surface zones m 2 GX The total radiative interchange area between any volume zones and any zones. m 2 h Convection heat transfer coefficient W/ m 2 .K I Iteration index J Iteration index K Gas absorption coefficient m -1 L Beam length m m· Mass flow rate Kg/m 2 N Number of bundles Q Heat transfer W R Random number SG The total exchange area between any surface zone and any volume zone m 2 SS The total exchange area between any two surface zones m 2 SX The total exchange area between any surface zone and any zone m 2 T Temperature K V Volume m 3 Tahseen Ali Al-Hattab, Ali Abd Al-Muhsen Al-Assady and Wisam Ahmed Abd Al-Wahid -Available online at: www.iasj.net IJCPE Vol.16 No.4 (Dec. 2015) 77 Greek Symbols Symbol Description Units α Absorbtivity γ Azimuthal angle radius θ Polar angle radius ε Emissivity τ Optical depth atm.m σ Stephan- Boltzmann constant. (5.669x 10~ 8 ) w/m 2 K 4 Subscript Symbol Description A Surface b Black d Flux to wall f Advection g Gas h Generation i General gas number j General zone number in Inward out, g Outward gas zone out, s Outward surface zone r Radiation s Surface up Upstream v Volume α Absorption γ Angle ε Emissivity References 1- Howell, J.R., and Perlmutter, M., "Monte Carlo Solution of Thermal Transfer Through Radiant Media between Gray Walls", Transactions of the ASME J. of Heat Transfer, vol. 86, NO.l, PP. 116-122, 1964. 2- Siegel, Robert, and Howell, John R., "Thermal Radiation Heat Transfer", 2 nd edition, McGraw- Hill, 1981. 3- Vercammen, H.A.J., and Froment, G.F., "Improved Zone Method Using Monte Carlo Techniques for the Simulation of Radiation In Industrial Furnaces", Int. J. Heat Mass Transfer, vol. 23, PP. 329- 337, 1980. 4- Steward, F.R., and Guruz, H.K., "Mathematical Simulation of an Industrial Boiler by the Zone Method of Analysis", Part 1 : Heat Transfer in steady confined Flame, ch. 3, PP. 46-71, 1974. 5- Gupta, R.P., wall, T.F., and Truelove, J.S., "Radiative Scatter by Fly Ash In puIverized-Coal-Fired Furnaces: Application of The Monte Carlo Method to Anisotropic Scatter", Int. J. Heat Mass Transfer, vol. 26, No. 11, PP. 1649-1660, 1983. 6- Kudo, Ki, Taniguchi, Hi, Kuroda, A., and Sasaki, T., "Development of General Purpose Computer Code for two/three-Dimension Heat Transfer Analysis", Proceeding of the seventh int. conf. held in U.S.A., vol. VII, part 1, 1991. 7- Al-Habbubi, Mohanned Mohammed, "Radiative Heat Transfer in Boilers Using the Zonal Method", M.Sc. thesis, University of Babylon, 2002. Appendix The values of the coefficient for emissivity of mixtures are given in table (A-1) and the values of the coefficient for the absorptivity of mixtures are given in table (A-2). Simulation the Radiation Zone of Al-Mussaib Power Plant by Using Monte Carlo Method 78 IJCPE Vol.16 No.4 (Dec. 2015) -Available online at: www.iasj.net Table (A-1), coefficients for emissivity Mixture, 1/  cw PP i i K 1 1,, 10  i b 4 2,, 10  i b 7 3,, 10  i b 11 4,, 10  i b 1 0.4303 5.150 -2.303 0.9779 -1.494 2 7.055 0.7749 3.399 -2.297 3.770 3 178.1 1.907 -1.824 0.5608 -0.5122 1 T P atm, 0.001≤ PS≤10.0 atm-m, 600≤ T ≤ 2400 K Table (A-2), Coefficients for absorptivity of mixtures Kia c ,, 1/  cw PP K i j 1 2 3 4 1 1 0.55657 E-00 -0.62824 E-03 0.31876 E-06 -0.52922 E-10 2 1 0.16676 E-01 0.15769 E-03 -0.10937 E-06 0.19588 E-10 3 1 0.28689 E-01 0.20697 E-03 -0.17473 E-06 0.37238 E-10 1 2 0.32964 E-03 0.27744 E-06 -0.26105 E-09 0.37807 E-13 2 2 0.50910 E-03 -0.76773 E-06 0.40784 E-09 -0.69622 E-13 3 2 0.24221 E-03 -0.55686 E-06 0.34884 E-09 -0.67887 E-13 1 3 -0.53441 E-06 0.33753 E-09 -0.10348 E-12 0.26027 E-16 2 3 0.37620 E-07 0.18729 E-09 -0.15889 E-12 0.30781 E-16 3 3 -0.19492 E-06 0.36102 E-09 -0.21480 E-12 0.41305 E-16 1 4 0.12381 E-09 -0.90223 E-13 0.38675 E-16 -0.99306 E-20 2 4 -0.32510 E-10 -0.26171 E-13 0.29848 E-16 -0.58387 E-20 3 4 0.41721 E-10 -0.73000 E-13 0.43100 E-16 -0.83182 E-20 1 T P atm, 1.0 c P atm, 0.001≤ ( wc PP  ) s ≤ 10.0 atm-m, 600 ≤ T, s T ≤ 2400 K