Toni Mikkonen
MATRIX HEAT TRANSFER CALCULATION MODEL FOR BACK-PASS TUBE BANK HEAT EXCHANGERS OF
FLUIDIZED BED STEAM BOILERS
Examiners: Professor, D.Sc. (Tech) Timo Hyppänen D.Sc. (Tech) Jouni Ritvanen
Supervisors: M.Sc. (Tech.) Ari Hottinen M.Sc. (Tech.) Jouni Miettinen
School of Energy Systems
Degree Programme in Energy Technology Toni Mikkonen
Matrix heat transfer calculation model for back-pass tube bank heat exchangers of fluidized bed steam boilers
Master’s Thesis 2020
124 pages, 64 figures, 3 tables, and 3 appendices
Examiners: Professor, D.Sc. (Tech) Timo Hyppänen and D.Sc. (Tech) Jouni Ritvanen Supervisors: M.Sc. (Tech.) Ari Hottinen and M.Sc. (Tech.) Jouni Miettinen
Keywords: matrix method, cell method, heat transfer, tube bank heat exchanger
The aim of this Master’s thesis is to develop a heat transfer calculation model for back-pass tube bank heat exchangers of fluidized bed steam boilers to calculate inlet, outlet, and aver- age temperatures for fluids and tube walls in each tube row. Accurate temperature infor- mation of each tube row and thermal stresses determined using them can be used to optimize material and thickness choices for heat exchanger tubes, thereby improving the cost-effec- tiveness of the manufacturing process.
Conventional calculation methods do not provide temperature profiles from inside the heat exchanger because the heat exchanger is calculated as a single unit. The new calculation model is based on the so-called matrix method. In it, tube banks of the heat exchanger are divided into several parts and the group of equations resulting from heat transfer calculations of different parts is solved by matrix calculation. Development of the calculation model re- quires an understanding of heat transfer mechanisms, heat exchanger configurations, energy balances, and pressure drop calculations.
The matrix equations vary depending on the heat exchanger configuration, so an essential part of the calculation model is the logic which creates a correct kind of matrix equation for each situation. Another key point is the calculation procedure in which heat transfer of tube banks and cavities between them is calculated alternately since the radiation emitted by cav- ities into the tube banks depends on tube temperatures and tube temperatures, in turn, depend on the radiation from cavities. Once the calculation model is complete, it still needs to be validated. In the validation process, results obtained from the calculation model are com- pared with the experimental measurements to find fouling factors for different situations.
Correction needed for superheaters is quite small and it is not dependent on the fuel. Econ- omizers need bigger correction which is dependent on the fuel probably due to the effect of condensation fouling mechanism. Results of the matrix calculation model were also com- pared with the results of the conventional calculation model. The comparison showed that the matrix calculation model provides a small improvement to the calculation accuracy.
TIIVISTELMÄ
Lappeenrannan-Lahden teknillinen yliopisto LUT School of Energy Systems
Energiatekniikan koulutusohjelma Toni Mikkonen
Matriisilämmönsiirtolaskentamalli leijupetihöyrykattiloiden takavedon putkipaketti- lämmönsiirtimille
Diplomityö 2020
124 sivua, 64 kuvaa, 3 taulukkoa ja 3 liitettä
Tarkastajat: Professori, TkT Timo Hyppänen ja TkT Jouni Ritvanen Ohjaajat: DI Ari Hottinen ja DI Jouni Miettinen
Hakusanat: matriisimenetelmä, solumenetelmä, lämmönsiirto, putkipakettilämmönsiirrin Tämän diplomityön tavoitteena on kehittää leijupetihöyrykattiloiden takavedon putkipaket- tilämmönsiirtimille lämmönsiirtolaskentamalli, jolla saadaan laskettua yksittäisten putkiri- vien sisäänmeno-, ulostulo- ja keskilämpötilat fluideille sekä putkiseinille. Tarkkoja putki- rivikohtaisia lämpötilatietoja ja niiden avulla määritettyjä lämpörasituksia voidaan käyttää lämmönsiirrinputkien materiaali- ja paksuusvalintojen optimoinnissa, jolloin valmistuspro- sessin kustannustehokkuutta saadaan parannettua.
Perinteiset laskentamenetelmät eivät tarjoa lämpötilaprofiileja lämmönsiirtimen sisältä, koska lämmönsiirrin lasketaan yhtenä kokonaisuutena. Uusi laskentamalli perustuu niin sa- nottuun matriisimenetelmään. Siinä lämmönsiirtimen putkipaketit jaetaan useaan osaan ja eri osien lämmönsiirtolaskentojen tuloksena muodostuva yhtälöryhmä ratkaistaan matriisi- laskennalla. Laskentamallin kehittäminen edellyttää ymmärrystä lämmönsiirtomekanis- meista, lämmönsiirrinkonfiguraatioista, energiataseista sekä painehäviölaskennoista.
Matriisiyhtälöt vaihtelevat riippuen lämmönsiirrinkonfiguraatiosta, joten olennainen osa las- kentamallia on logiikka, joka luo oikeanlaisen matriisiyhtälön kullekin tilanteelle. Toinen keskeinen asia on laskentamenettely, jossa putkipakettien ja niiden välissä olevien onkaloi- den lämmönsiirtoa lasketaan vuorotellen, koska onkaloiden putkipaketteihin lähettämä sä- teily riippuu putkilämpötiloista ja putkilämpötilat puolestaan onkaloista tulevasta säteilystä.
Kun laskentamalli on valmis, se täytyy vielä validoida. Validointiprosessissa laskentamal- lista saatuja tuloksia verrataan kokeellisiin mittaustuloksiin, jotta löydetään likaantumisker- toimet erilaisille tilanteille. Tulistimille tarvittava korjaus on melko pieni eikä se ole riippu- vainen polttoaineesta. Ekonomaiserit tarvitsevat isomman korjauksen, joka on riippuvainen polttoaineesta luultavasti johtuen kondensaatiolikaantumismekanismin vaikutuksesta. Mat- riisilaskentamallin tuloksia verrattiin myös perinteisen laskentamallin tuloksiin. Vertailu osoitti, että matriisilaskentamalli tarjoaa pienen parannuksen laskentatarkkuuteen.
I would like to thank Sumitomo SHI FW Energia Oy for this opportunity to learn a lot about steam boiler technology and software development. Especially, I would like to thank my thesis supervisors Ari Hottinen and Jouni Miettinen for their great advice and mentoring during the project. Also, I would like to thank Rita Marjeta, Ari Kettunen, and Ossi Sippu for the cooperation and support they provided whenever needed.
I would like to thank my thesis examiners Timo Hyppänen and Jouni Ritvanen who gave me guidance and constructive feedback during the Master’s thesis project. Lastly, I would like to thank my family and friends who have supported me during my studies.
Varkaus, December 2020 Toni Mikkonen
TABLE OF CONTENTS
1 INTRODUCTION ... 9
2 PERFORMANCE CALCULATIONS OF HEAT EXCHANGERS ... 10
2.1 Conventional methods for solving heat exchanger performance ... 10
2.1.1 LMTD method ... 10
2.1.2 e-NTU method ... 14
2.2 Matrix method for solving heat exchanger performance ... 16
3 HEAT TRANSFER MECHANISMS ... 19
3.1 Conduction ... 19
3.2 Convection ... 21
3.2.1 Flow over tubes and tube banks ... 25
3.2.2 Flow inside tubes ... 27
3.3 Thermal radiation ... 29
3.3.1 Radiation from surface ... 30
3.3.2 Radiation between surfaces without participating gas ... 34
3.3.3 Radiation between surface and gas ... 37
3.3.4 Radiation between surface and dusty gas ... 42
3.3.5 Radiation between surfaces with participating gas ... 45
4 WORKING PRINCIPLE OF FLUIDIZED BED BOILERS ... 46
4.1 Types of combustion technology ... 46
4.2 Types of working fluid circulation ... 47
4.3 Heat exchangers of fluidized bed boilers back-pass ... 49
4.3.1 Superheaters and reheaters ... 50
4.3.2 Evaporating heat exchangers ... 53
4.3.3 Economizers ... 54
4.3.4 Air preheaters ... 55
5 ENERGY BALANCE OF TUBE BANK HEAT EXCHANGERS ... 56
5.1 The heat from external radiation ... 57
5.1.1 External radiation from cavities ... 58
5.1.2 External radiation from surfaces ... 64
5.2 Heat to the cold fluid ... 64
5.2.1 Outside radiative thermal resistance ... 66
5.2.2 Outside convective thermal resistance ... 68
5.2.3 Tube wall thermal resistance... 73
5.2.4 Foul layers’ thermal resistances ... 74
5.2.5 Inside convective thermal resistance ... 76
5.2.6 Overall heat transfer coefficient ... 79
5.3 Heat to sidewalls, hanger tubes, and losses ... 83
5.4 Heat from the hot fluid ... 85
6 PRESSURE LOSSES OF TUBE BANK HEAT EXCHANGERS ... 86
6.1 Outside pressure losses ... 86
6.2 Inside pressure losses ... 91
6.2.1 Friction pressure losses ... 91
6.2.2 Inlet, outlet, and bend pressure losses ... 95
7 DEVELOPMENT OF THE CALCULATION MODEL ... 97
7.1 Matrix equation formation logic ... 99
7.1.1 Matrix equations for horizontal tube bank heat exchangers ... 99
7.1.2 Matrix equations for vertical tube bank heat exchangers ... 104
7.1.3 Effect of calculation direction on the matrix equations ... 106
7.2 Calculation procedure ... 109
7.2.1 Calculation procedure for horizontal tube bank heat exchangers ... 109
7.2.2 Calculation procedure for vertical tube bank heat exchangers ... 114
7.2.3 Stabilization of the calculation procedure ... 115
8 VALIDATION OF THE CALCULATION MODEL ... 116
8.1 Determination of fouling factors ... 118
8.2 Accuracy of the new model compared to the conventional model ... 121
9 SUMMARY AND CONCLUSIONS ... 123
REFERENCES ... 125
APPENDICES
Appendix 1. Outside matrices Appendix 2. Inside matrices
Appendix 3. Emissivity diagrams for carbon dioxide and sulfur dioxide
𝐴 area m2 𝐴𝑝 specific projection area of the dispersed particles m2/kg
𝑎 transverse pitch ratio -
𝑎𝑓 absorption fraction -
𝑏 longitudinal pitch ratio -
𝐶 heat capacity rate W/K
𝐶𝐷 drag coefficient -
𝑐 diagonal pitch ratio -
𝑐𝑝 specific heat capacity at constant pressure J/kgK
𝐷 depth of the tube bank m
𝑑 diameter m
𝐸 emissive heat flux W/m2
𝑒 effectiveness -
𝐹 view factor -
𝑓 friction factor -
𝑓 fouling factor -
𝑓𝑎𝑙𝑖𝑔 correction factor for aligned tube configuration - 𝑓𝑠𝑡𝑎𝑔 correction factor for staggered tube configuration -
𝑓𝑝 pressure correction factor -
𝑓𝑅𝑁 row number correction factor -
𝑓𝑁𝑅 number of rows correction factor -
𝑓𝑓𝑎 flow angle correction factor -
𝑓𝑡𝑓𝑠 factor typical for the flow situation -
𝑓𝑢 under-relaxation factor -
𝐺 irradiation W/m2
𝐻 height of the tube bank m
ℎ specific enthalpy kJ/kg
ℎ𝑐 convective heat transfer coefficient W/m2K ℎ𝑟 radiative heat transfer coefficient W/m2K
ℎ𝑒𝑥 heat transfer coefficient of the external radiation W/m2K
𝐼𝑒 electric current A
𝐽 radiosity W/m2
𝐾𝑎𝑏𝑠,𝑔 absorption coefficient of the gas phase 1/m 𝐾𝑒𝑚𝑖,𝑔 emission coefficient of the gas phase 1/m 𝐾𝑖 inside correction factor for heat transfer direction - 𝐾𝑜 outside correction factor for heat transfer direction -
𝑘 thermal conductivity W/mK
𝐿𝑐 characteristic length m
𝐿𝑚𝑏 mean beam length m
𝐿𝑝 particle load kg/m3
𝑙 tube length m
𝑁𝑅 number of rows -
𝑁𝑀𝑅 number of main resistance -
𝑝 pressure bar, Pa
𝑄𝑎𝑏𝑠 mean relative absorption efficiency of particles - 𝑄𝑏𝑠𝑐 mean relative backscattering efficiency of particles -
𝑞𝑚 mass flow rate kg/s
𝑅 thermal resistance K/W
𝑅𝑒 electrical resistance 𝛺
𝑅𝑁 row number -
𝑟 radius m
𝑠𝐷 diagonal pitch m
𝑠𝐿 longitudinal pitch m
𝑠𝑇 transverse pitch m
𝑇 temperature ºC, K
𝑡 thickness m
𝑈 voltage V
𝑈 overall heat transfer coefficient W/m2K
𝑉 volume m3
𝑣 kinematic viscosity m2/s
𝑊 width of the tube bank m
𝑤 velocity m/s
𝑤𝑓 low velocity in the free cross section m/s 𝑤𝑛 mean flow velocity in the narrowest cross section m/s
Greek symbols
𝛼 total hemispherical absorptivity -
𝛼𝜆 spectral absorptivity -
𝛼𝜆,𝜃 spectral directional absorptivity -
𝛽 total hemispherical reflectivity -
𝜀 absolute roughness
𝜀 total hemispherical emissivity -
𝜀𝜆 spectral emissivity -
𝜀𝜆,𝜃 spectral directional emissivity - 𝜏 total hemispherical transmissivity -
𝜇 dynamic viscosity Pa⋅s
∆ difference
𝛥𝑝𝑏 bend pressure loss Pa
𝛥𝑝𝑓 friction pressure loss Pa
𝛥𝑝𝑖 tube inlet pressure loss Pa
𝛥𝑝𝑜 tube outlet pressure loss Pa
𝜓 void fraction -
∅𝑒𝑚𝑖 optical thickness for emission -
∅𝑎𝑏𝑠 optical thickness for absorption -
𝜃 flow angle °
𝜌 density kg/m3
𝜙 heat rate W
𝜙′′ heat flux W/m2
𝜎 Stefan-Boltzmann constant 5,67⋅10-8 W/m2K4
Dimensionless parameters
𝑁𝑇𝑈 number of heat transfer units 𝑁𝑢 Nusselt number
𝑃𝑟 Prandtl number 𝑅𝑒 Reynolds number
Subscripts 𝑎 average 𝑎𝑏𝑠 absorbed 𝑎𝑙𝑖𝑔 aligned 𝑏 blackbody
𝑐 cold, cross-sectional 𝑐𝑜𝑛𝑑 conduction
𝑐𝑜𝑛𝑣 convection 𝑐𝑜𝑛𝑡 contact 𝑑𝑦𝑛 dynamic 𝑒𝑓𝑓 effective
𝑒𝑥 external radiation 𝑓 fluid, fouling layer 𝑔 gas
ℎ hot
ℎ𝑡 hanger tubes 𝑖 inside, inlet 𝑙𝑎𝑚 laminar
𝑙𝑚 logarithmic mean 𝑙𝑜𝑠𝑠 losses
𝑚 mean
𝑚𝑎𝑥 maximum 𝑚𝑖𝑛 minimum 𝑚𝑝 mean particle 𝑜 outside, outlet 𝑜𝑣𝑒 overall
𝑝 particles
𝑝𝑠 projected surface 𝑟 ratio
𝑟𝑎𝑑 radiation 𝑟𝑒𝑓 reflected 𝑠 surface 𝑠𝑡𝑎𝑔 staggered 𝑠𝑤 sidewalls 𝑡 tube row wall 𝑡𝑜𝑡 total
𝑡𝑟𝑎 transmitted 𝑡𝑢𝑟𝑏 turbulent 𝑢𝑛𝑐 uncorrected 𝑤 wall
Abbreviations
BFB Bubbling Fluidized Bed CFB Circulating Fluidized Bed 𝐻2𝑂 water vapor
𝐶𝑂2 carbon dioxide 𝑆𝑂2 sulfur dioxide
1 INTRODUCTION
Cost-effectiveness is always desired when heat exchangers are manufactured. In general, when heat exchangers are designed and their steady-state heat transfer performance is inves- tigated, calculations are made so that the overall heat transfer rate and unknown inlet or outlet temperatures are solved for the whole heat exchanger. This usually gives enough ac- curate information about the overall behavior of the heat exchanger, but it does not give accurate information about what is happening in the different parts of the heat exchanger. In particular, it would be beneficial to know the exact local temperatures inside the heat ex- changer and the resulting thermal stresses. This accurate information could be used during the designing and manufacturing process to select optimal construction materials, material thicknesses, and fabrication techniques for the heat exchanger tubes. (Silaipillayarputhur &
Idem 2013, 338)
In this thesis, a heat transfer calculation model is developed for different kinds of tube bank heat exchangers which can be found from back-pass sections of industrial bubbling fluidized bed (BFB) boilers and circulating fluidized bed (CFB) boilers designed by Sumitomo SHI FW Energia Oy. The aim is to obtain intermediate tube and fluid temperatures of the tube rows inside the heat exchangers for tube designing purposes and improve the accuracy of heat transfer performance calculations using a so-called matrix method where tube banks of the heat exchangers are divided into many parts so that individual tube rows are calculated separately. The calculation model is programmed using C# (C sharp) and Fortran languages.
The thesis begins by introducing conventional heat transfer performance calculation meth- ods and the basic idea of the matrix method. Then, different heat transfer mechanisms are covered which are the foundation of the heat transfer calculations. After that, working prin- ciples of fluidized bed boilers and different back-pass heat exchangers are discussed to un- derstand what’s kind of situations the calculation model must be able to handle. In addition, energy balance and pressure losses of the tube bank heat exchangers are considered. In the end, the calculation model is developed and also validated, that is, results obtained from the model are compared to the real measurement data gathered from power plants and also to the results from the earlier conventional model.
2 PERFORMANCE CALCULATIONS OF HEAT EXCHANGERS
In this chapter, conventional methods for heat transfer performance calculations of heat ex- changers are discussed and the idea of the matrix method is introduced.
2.1 Conventional methods for solving heat exchanger performance
When the performance of the heat exchanger is calculated, there are six parameters from which four are known and two are unknown. Those parameters are mass flow rates, inlet temperatures (or specific enthalpies), and outlet temperatures (or specific enthalpies) of both fluids. From now onward specific enthalpies are called just enthalpies in the text. In the following examples, mass flow rates and inlet enthalpies are assumed to be known. En- thalpies are used instead of temperatures because the calculation model must be able to han- dle also situations where one of the fluids evaporates. In those situations, the temperature does not change so enthalpy change is wanted to know. In real cases, there are also pressure changes in both fluid sides but those are neglected still at this point. Traditionally perfor- mance is calculated for the whole heat exchanger, that is, outlet enthalpies (or temperatures) of both fluids and total heat transfer rate are wanted to know. Calculations can be made using either the LMTD (logarithmic mean temperature difference) method or the e-NTU (effec- tiveness-number of heat transfer units) method (Ganapathy 2018, 63-64).
2.1.1 LMTD method
Consider a simple heat exchange situation where parallel flowing fluids are separated by a wall. This so-called parallel-flow heat exchanger is illustrated in figure 1.
Figure 1. Parallel-flow heat exchanger.
Assuming that there are no heat losses to the surroundings, and thus all the heat released from the hot fluid goes to the cold fluid, energy balance is as follows (Incropera et al. 2007, 676-677):
𝜙ℎ = 𝜙𝑐 (2.1)
where 𝜙ℎ is the heat from the hot fluid (W) and 𝜙𝑐 is the heat to the cold fluid (W).
Heat rates can be expressed as follows:
𝜙ℎ = 𝑞𝑚,ℎ𝑐𝑝,ℎ(𝑇ℎ,𝑖 − 𝑇ℎ,𝑜) = 𝐶ℎ(𝑇ℎ,𝑖− 𝑇ℎ,𝑜) = 𝑞𝑚,ℎ(ℎℎ,𝑖− ℎℎ,𝑜) (2.2) 𝜙𝑐 = 𝑞𝑚,𝑐𝑐𝑝,𝑐(𝑇𝑐,𝑜− 𝑇𝑐,𝑖) = 𝐶𝑐(𝑇𝑐,𝑜− 𝑇𝑐,𝑖) = 𝑞𝑚,𝑐(ℎ𝑐,𝑜− ℎ𝑐,𝑖) (2.3) where 𝑞𝑚 is the mass flow rate (kg/s), 𝐶 is the heat capacity rate (W/K) and ℎ is the specific enthalpy (kJ/kg). Subscripts ℎ and 𝑐 are for hot and cold fluids and i and o are for inlet and outlet, respectively.
Notation 𝜙 without subscripts will be used for the heat transfer rates in the following equa- tions because in this simplified case it is the same for both fluids. Specific heat capacities in earlier equations are not constant but a function of temperature and pressure, and also a so- called overall heat transfer coefficient change due to changing fluid properties and due to changing flow conditions. Usually, these changes are not very big in the heat exchangers, so a common practice is to use average values of fluid properties. (Incropera et al. 2007, 677) In the LMTD method, the first step is to guess outlet enthalpies to calculate average en- thalpies and then average thermal properties of both fluids in the heat exchanger. After that average overall heat transfer coefficient for the heat exchanger can be calculated. Overall heat transfer coefficient is defined for the heat transfer systems where many heat transfer mechanisms exist, and it is related to the total thermal resistance of the system (Incropera et al. 2007, 100). Solving the overall heat transfer coefficient is introduced more specifically in chapter 5.
After the overall heat transfer coefficient has been solved, the heat transferred between fluids can be calculated from the following equation which has an analogy to Newton’s law of cooling (Incropera et al. 2007, 678):
𝜙 = 𝑈𝐴𝑠𝛥𝑇𝑙𝑚 (2.4) where 𝑈 is the overall heat transfer coefficient (W/m2K) and 𝛥𝑇𝑙𝑚 is the logarithmic mean temperature difference (K).
Because temperatures of the fluids are not changing linearly in the heat exchanger, which is illustrated in figure 2, the average temperature difference between hot and cold fluid cannot be calculated just by taking the average from the inlet and outlet temperature differences.
Thus, a so-called logarithmic mean temperature difference is used which gives a proper value for the average temperature difference. (Incropera et al. 2007, 676)
Figure 2. Temperature profiles of a parallel-flow heat exchanger. (Incropera et al. 2007, 677)
Logarithmic mean temperature difference for the above parallel flow case is calculated as follows (Incropera et al. 2007, 678):
𝛥𝑇𝑙𝑚 = 𝛥𝑇1−𝛥𝑇2
ln (𝛥𝑇1/𝛥𝑇2)= 𝛥𝑇2−𝛥𝑇1
ln (𝛥𝑇2/𝛥𝑇1) (2.5)
where temperature differences are:
𝛥𝑇1 = 𝑇ℎ,𝑖− 𝑇𝑐,𝑖 (2.6)
𝛥𝑇2 = 𝑇ℎ,𝑜− 𝑇𝑐,𝑜 (2.7)
In the case of the counter-flow heat exchanger, which is illustrated in figure 3, the logarith- mic mean temperature difference is calculated similarly as earlier using equation 2.5. The only difference is that temperature differences 𝛥𝑇1 and 𝛥𝑇2 are calculated using the follow- ing equations because temperature profiles are now different which is shown in figure 4 (Incropera et al. 2007, 679):
𝛥𝑇1 = 𝑇ℎ,𝑖− 𝑇𝑐,𝑜 (2.8)
𝛥𝑇2 = 𝑇ℎ,𝑜− 𝑇𝑐,𝑖 (2.9)
Figure 3. Counter-flow heat exchanger.
Figure 4. Temperature profiles of a counter-flow heat exchanger. (Incropera et al. 2007, 679)
If the configuration of the heat exchanger is more complicated than just a simple parallel or counter-flow heat exchanger, the value of the logarithmic mean temperature difference must be corrected as follows (Incropera et al. 2007, W-37):
𝛥𝑇𝑙𝑚 = 𝐹𝑙𝑚𝛥𝑇𝑙𝑚,𝐶𝐹 (2.10)
where 𝐹𝑙𝑚 is the correction factor for the logarithmic mean temperature difference (-) and 𝛥𝑇𝑙𝑚,𝐶𝐹 is the logarithmic mean temperature difference calculated using counter flow equa- tions.
Equations and graphs for correction factors of different heat exchanger configurations can be found in the literature. There are also some special cases where configuration does not influence the behavior of the heat exchanger, that is, the correction factor 𝐹𝑙𝑚 is 1. This is the case if one of the fluids is evaporating or condensing so the temperature of that fluid does not change. (Incropera et al. 2007, W-39)
After the heat transfer rate is determined using equation 2.4, outlet enthalpies can be calcu- lated from heat rate equations 2.2 and 2.3 as follows:
𝑞𝑚,ℎ(ℎℎ,𝑖 − ℎℎ,𝑜) = 𝜙
=> ℎℎ,𝑜= ℎℎ,𝑖− 𝜙 𝑞𝑚,ℎ
𝑞𝑚,𝑐(ℎ𝑐,𝑜− ℎ𝑐,𝑖) = 𝜙
=> ℎ𝑐,𝑜 = ℎ𝑐,𝑖+ 𝜙 𝑞𝑚,𝑐
If calculated outlet enthalpies differ a lot from the initial guesses, these new outlet enthalpies are used as initial guesses and the calculation procedure is repeated until the outlet enthalpies no longer change significantly.
2.1.2 e-NTU method
If the e-NTU method is used instead of the LMTD method, the first step in a solving process is once again to guess outlet enthalpies and calculate the average overall heat transfer
coefficient. After that, a value for the NTU can be calculated as follows (Incropera et al.
2007, 687-688):
𝑁𝑇𝑈 = 𝑈𝐴𝑠
𝐶𝑚𝑖𝑛 (2.11)
where 𝑁𝑇𝑈 is the number of heat transfer units (-) and 𝐶𝑚𝑖𝑛 is the minimum heat capacity rate (W/K).
𝐶𝑚𝑖𝑛 can be either 𝐶ℎ or 𝐶𝑐 depending on which is smaller. Also, the ratio of heat capacity rates must be determined as follows (Incropera et al. 2007, 688):
𝐶𝑟= 𝐶𝑚𝑖𝑛
𝐶𝑚𝑎𝑥 (2.12)
where 𝐶𝑟 is the heat capacity rate ratio (-) and 𝐶𝑚𝑎𝑥 is the maximum heat capacity rate (W/K).
In special condensation or evaporation cases temperature of the one fluid does not change so the heat capacity rate of that fluid is considered approaching infinity and is then also the maximum heat capacity rate. In this situation, the heat capacity ratio is approaching zero, that is 𝐶𝑟= 0. (Incropera et al. 2007, 688)
Effectiveness of the heat exchanger is a function of NTU and 𝐶𝑟 so it can now be determined.
The effectiveness also depends on the heat exchanger configuration. Correlations of effec- tiveness for different heat exchanger configurations can be found from the literature. For example, the correlation for the earlier discussed parallel flow situation is as follows (In- cropera et al. 2007, 687, 689):
𝑒 =1−exp [−𝑁𝑇𝑈(1+𝐶𝑟)]
1+𝐶𝑟 (2.13)
For the counter flow situation, correlation depends on the value of 𝐶𝑟 as follows:
𝑒 = 1−exp [−𝑁𝑇𝑈(1−𝐶𝑟)]
1−𝐶𝑟⋅exp [−𝑁𝑇𝑈(1−𝐶𝑟)] 𝐶𝑟< 1 (2.14)
𝑒 = 𝑁𝑇𝑈
1+𝑁𝑇𝑈 𝐶𝑟= 1 (2.15)
As the configuration does not affect the behavior of the heat exchanger if there are conden- sation or evaporation, that is 𝐶𝑟 = 0, there is only one single correlation for these situations which is shown below:
𝑒 = 1 − exp (−𝑁𝑇𝑈) (2.16)
After the effectiveness is determined, another expression for the effectiveness is needed to solve the performance. That is as follows (Incropera et al. 2007, 689):
𝑒 = 𝜙
𝜙𝑚𝑎𝑥 = 𝜙
𝐶𝑚𝑖𝑛(𝑇ℎ,𝑖−𝑇𝑐,𝑖) (2.17)
where 𝜙𝑚𝑎𝑥 is the maximum possible heat transfer rate (W).
From equation 2.17 the heat transfer rate can be solved:
𝜙 = 𝑒𝐶𝑚𝑖𝑛(𝑇ℎ,𝑖− 𝑇𝑐,𝑖) (2.18)
After the heat transfer rate is determined, outlet enthalpies can be calculated from equations 2.2 and 2.3 as usual. If calculated outlet enthalpies differ a lot from the initial guesses, these new outlet enthalpies are used as initial guesses and the calculation procedure is repeated until the outlet enthalpies no longer change significantly.
2.2 Matrix method for solving heat exchanger performance
When intermediate enthalpies (and temperatures) inside the heat exchanger are also wanted to know, in addition to outlet enthalpies, the heat exchanger must be divided into many parts.
Consider earlier introduced parallel-flow heat exchanger which is now divided into two dif- ferent parts that is illustrated in figure 5.
Figure 5. Parallel-flow heat exchanger divided into two parts.
The first step of the solving process is to guess outlet enthalpies as usual. Assuming linear enthalpy profiles inside the heat exchanger for both fluids, initial guesses for intermediate enthalpies can also be calculated. After that, average enthalpies, average fluid properties, overall heat transfer coefficients, and heat transfer rates in each part of the heat exchanger can be obtained. In the literature, the e-NTU method is usually used to solve heat rates when the heat exchanger is divided into several parts (VDI 2010, 34-35, Silaipillayarputhur &
Idem 2013).
After heat transfer rates are calculated for each part of the heat exchanger, the equation group of heat rate equations is as follows:
𝑞𝑚,ℎ(ℎℎ,𝑖 − ℎℎ,1) = 𝜙1 (2.19)
𝑞𝑚,𝑐(ℎ𝑐,1− ℎ𝑐,𝑖) = 𝜙1 (2.20)
𝑞𝑚,ℎ(ℎℎ,1− ℎℎ,2) = 𝜙2 (2.21)
𝑞𝑚,𝑐(ℎ𝑐,2− ℎ𝑐,1) = 𝜙2 (2.22)
Above equations can be calculated open and all unknown terms can be moved to the left side and know terms to the right side. Then equations are as follows:
−ℎℎ,1 = 𝜙1
𝑞𝑚,ℎ− ℎℎ,𝑖
ℎ𝑐,1 = 𝜙1
𝑞𝑚,𝑐 + ℎ𝑐,𝑖
ℎℎ,1− ℎℎ,2 = 𝜙2 𝑞𝑚,ℎ
ℎ𝑐,2− ℎ𝑐,1= 𝜙2 𝑞𝑚,𝑐
The group of equations can then be converted to matrix form as follows from which the method takes its name:
[
−1 0 0 0
0 0 1 0
1 −1 0 0
0 0 −1 1
] [
ℎℎ,1 ℎℎ,2 ℎ𝑐,1 ℎ𝑐,2]
= [
𝜙1/𝑞𝑚,ℎ− ℎℎ,𝑖 𝜙1/𝑞𝑚,𝑐 + ℎ𝑐,𝑖
𝜙2/𝑞𝑚,ℎ 𝜙2/𝑞𝑚,𝑐 ]
(2.23)
In the above matrix equation, the square matrix on the left can be called a matrix A, the enthalpy vector can be called a vector h and the known vector on the right-hand side can be called a vector b. The vector h can be solved by inverting the matrix A and multiplying it by the vector b as follows:
𝒉 = 𝑨−1⋅ 𝒃 (2.24)
If calculated enthalpies differ a lot from the initial guesses, calculated enthalpies are used as new initial guesses and the solution is iterated as in the case of conventional methods. The working principle of the calculation model for tube bank heat exchangers will be introduced more accurately in chapter 7.
In the literature, the matrix method is also called a cell method. These methods can be used to solve any kind of heat exchanger configuration. For example, in figure 6 there is a so- called shell-and-tube heat exchanger that is divided into six different parts.
Figure 6. Shell-and-tube heat exchanger.
3 HEAT TRANSFER MECHANISMS
Thermal energy is transferred whenever temperature differences exist. Heat transfer can occur through three basic mechanisms, which are conduction, convection, and thermal ra- diation. (Incropera et al. 2007, 2-3)
3.1 Conduction
Heat transfer by conduction is related to the interactions between the molecules of a sub- stance. Molecules with higher temperatures have stronger random motion and thus higher kinetic energy. Therefore, during intermolecular interactions caused by these random mo- tions, energy is transferred from more energetic molecules to less energetic molecules. In conductor materials, conduction is also associated with the movement of energy-carrying free electrons. (VDI 2010, 18)
The so-called Fourier's law describes the amount of heat transferred by conduction. In the case of one-dimensional heat conduction through the wall in the x-direction, the equation is as follows (Incropera et al. 2007, 4):
𝜙′′ = −𝑘𝑑𝑇
𝑑𝑥 (3.1)
where 𝜙′′ is the heat flux (W/m2), 𝑘 is the thermal conductivity (W/mK) and 𝑑𝑇/𝑑𝑥 is the temperature gradient (K/m).
Under steady-state conditions, the temperature distribution is linear within the plane wall, so the earlier equation can also be represented as follows:
𝜙′′ = −𝑘𝑑𝑇
𝑑𝑥 = −𝑘∆𝑇
∆𝑥 = −𝑘𝑇𝑠2−𝑇𝑠1
𝑡𝑤 = 𝑘𝑇𝑠1−𝑇𝑠2
𝑡𝑤 (3.2)
where 𝑇𝑠1 and 𝑇𝑠2 are the surface temperatures of the plane wall (K) and 𝑡𝑤 is the wall thickness (m).
The heat transfer rate is obtained when the heat flux is multiplied by the heat transfer area as follows (Incropera et al. 2007, 98):
𝜙 = 𝜙′′𝐴𝑐 = 𝑘𝐴𝑐𝑇𝑠1−𝑇𝑠2
𝑡𝑤 = 𝑇𝑠1𝑡𝑤−𝑇𝑠2
𝑘𝐴𝑐
(3.3)
where 𝜙 is the heat transfer rate (W) and 𝐴𝑐 is the cross-sectional area normal to the direction of the heat transfer (m2).
At a general level, thermal conductivity values are higher for solids than for liquids and higher for liquids than for gases. The reason is that in this order, the distance between mol- ecules usually increases, leading to weaker and more rarely occurring interactions between molecules. The value of thermal conductivity also depends on the temperature. Pressure can also affect thermal conductivity in the case of gases if the conditions are close to a vacuum, and in liquids, if the conditions are close to the critical point. (Incropera et al. 2007, 4, 60- 65)
Heat conduction can be seen as analogous to electrical conduction. The electric current can be calculated from Ohm's law as follows (Incropera et al. 2007, 99):
𝐼𝑒 =∆𝑈
𝑅𝑒 (3.4)
where 𝐼𝑒 is the electric current (A), ∆𝑈 is the voltage difference (V) and 𝑅𝑒 is the electrical resistance (𝛺).
When comparing equations 3.3 and 3.4, it can be seen that the term 𝑡𝑤/𝑘𝐴𝑐 in equation 3.3 corresponds to the electrical resistance in equation 3.4. Thus, a so-called thermal resistance for conduction is as follows (Incropera et al. 2007, 99):
𝑅𝑐𝑜𝑛𝑑 = 𝑡𝑤
𝑘𝐴𝑐 (3.5)
where 𝑅𝑐𝑜𝑛𝑑 is the thermal resistance for conduction (K/W).
If the structure through which the heat is conducted consists of several different materials, the temperature may change significantly at the material interfaces. The reason for this tem- perature drop is the fact that there are small gaps at the material interfaces due to the surface roughness. Materials are thus not in perfect contact. A term called thermal contact resistance 𝑅𝑐𝑜𝑛𝑡 is used to describe this effect. (Incropera et al. 2007, 102)
3.2 Convection
The term advection is used when referring to a situation where heat is transferred only by the macroscopic motion of the fluid. In reality, in addition to this macroscopic motion, there are still also random molecular level motions in the substance which affect the heat transfer.
The term convection is used when referring to the combined effect of these two mechanisms.
Regarding the relative strengths of advection and conduction mechanisms in convection heat transfer, conduction has a strong effect near the surface where velocity is low as will be illustrated later, and the effect of advection increases farther from the surface. Convection can be classified as either natural or forced convection. In natural convection, the movement of the fluid is caused by buoyancy forces which are generated by density differences in the fluid. Differences in density, in turn, are due to differences in temperature. In forced con- vection, the movement of the fluid is caused by some external force. For example, the move- ment of the fluid can be caused by a pump. (Incropera et al. 2007, 6-8)
Convection is usually associated with situations where heat is transferred between a solid surface and a moving fluid. Essential terms related to such situations are velocity and thermal boundary layers. Consider a situation in which a fluid flow having constant velocity and temperature profiles in the y-direction begin to flow over an isothermal flat plate in the x- direction. The situation is illustrated in figure 7. (Incropera et al. 2007, 348-350)
Figure 7. Development of velocity and thermal boundary layers on a flat plate. (Incropera et al. 2007, 348- 350)
The velocity of the fluid molecules in contact with the surface can be assumed to be zero and due to shear stresses, the velocity of the adjacent fluid layer in the y-direction also de- creases. When the flow proceeds in the x-direction, more and more fluid layers are affected by earlier fluid layers. A region where there are velocity gradients in the y-direction is called
a velocity or hydrodynamic boundary layer and its thickness increases until it is so thick that the effect of the previous fluid layer on the next one is insignificant. In addition to the ve- locity boundary layer, a thermal boundary layer with temperature gradients also begins to develop if surface and fluid temperatures differ from each other. Fluid molecules in contact with the surface acquire the surface temperature, and then the heat is transferred to the next fluid layers. Like the thickness of the velocity boundary layer, also the thickness of the ther- mal boundary layer increases as the flow proceeds in the x-direction as heat is further trans- ferred to the fluid in the y-direction. (Incropera et al. 2007, 348-350)
The velocity boundary layer can be either laminar or turbulent in nature. In a laminar flow, viscous forces are strong enough to keep the flow regular so fluid molecules move along certain streamlines. In a turbulent flow, viscous forces cannot keep the flow regular so oc- casional vortices are generated. Those vortices mix the flow so that higher speed molecules from the free stream are carried towards the surface and slower speed molecules from the surface are carried towards the free stream. For this reason, the shape of the velocity profile is flatter in turbulent than in the laminar boundary layer which can be seen from figure 8.
The turbulent boundary layer has three different layers, which are a turbulent region, a buffer layer, and a viscous or laminar sublayer. Usually, when the velocity boundary layer is de- veloping on some surface, the boundary layer is first laminar, then comes transition zone where the flow has both laminar and turbulent properties which share varies as a function of time and finally flow turn to be fully turbulent. In figure 8, boundary layer development from laminar to turbulent on a flat plate is illustrated. (Incropera et al. 2007, 359-361)
Figure 8. Turbulent and laminar velocity boundary layers on a flat plate. (Incropera et al. 2007, 360)
A dimensionless parameter called a Reynolds number describes the nature of the flow, that is, whether the boundary layer has time to change from the laminar to turbulent or whether it is constantly laminar. The equation for the Reynolds number is as follows (Incropera et al.
2007, 360):
𝑅𝑒 =𝜌𝑤𝐿𝑐
𝜇 (3.6)
where 𝑅𝑒 is the Reynolds number (-), 𝜌 is the density (kg/m3), 𝑤 is the velocity (m/s), 𝐿𝑐 is the characteristic length (m) and 𝜇 is the dynamic viscosity (Pa⋅s).
Dimension used as a characteristic length in the above equation depends on the case under consideration. For example, in the foregoing case where the fluid flows over the flat plate, a length of the plate is used as a characteristic length. On the other hand, in cases where the fluid flows over or inside a tube, the diameter of the tube is usually used (Incropera et al.
2007, 424, 487). Determination of the flow velocity also differs depending on the flow situ- ation.
Basically, the Reynolds number tells the relative strengths of the inertial and viscous forces acting in the fluid. The boundary layer is called to be laminar when the Reynolds number is small and when the Reynolds number is large enough, the boundary layer is turned to be turbulent. A critical Reynolds number is a term used to describe the border between these two flow types. Because in reality change from laminar to turbulent flow is not sharp but the amount of turbulence grows slowly, any fixed value does not tell the absolute truth. The value range of the critical Reynolds number also depends on the case. For example, if the flow proceeds over the flat plate, the range is from 105 to 3⋅106 and the typical fixed value used in calculations is 5⋅105. (Incropera et al. 2007, 359-361)
The so-called Newton’s law of cooling describes the amount of heat transferred by convec- tion as follows (Incropera et al. 2007, 8):
𝜙′′ = ℎ𝑐(𝑇𝑠− 𝑇𝑓) (3.7)
where ℎ𝑐 is the convective heat transfer coefficient (W/m2K), 𝑇𝑠 is the surface temperature (K) and 𝑇𝑓 is the fluid temperature (K).
The heat transfer rate is obtained when the heat flux is multiplied by the heat transfer area as follows (Incropera et al. 2007, 99):
𝜙 = ℎ𝑐𝐴𝑠(𝑇𝑠− 𝑇𝑓) = 𝑇𝑠−𝑇1 𝑓 ℎ𝑐𝐴𝑠
(3.8)
where 𝐴𝑠 is the surface area (m2).
As in the case of conduction, also in the case of convection, the thermal resistance can be found from the electricity analogy. Thus, resistance is as follows:
𝑅𝑐𝑜𝑛𝑣 = 1
ℎ𝑐𝐴𝑠 (3.9)
where 𝑅𝑐𝑜𝑛𝑣 is the thermal resistance for convection (K/W).
The value of the convection heat transfer coefficient is dependent on the conditions in the boundary layers. Conditions are affected by surface geometry, nature of the flow, and fluid properties like density, viscosity, thermal conductivity, and specific heat capacity. In gen- eral, convection heat transfer coefficients are higher for forced than for free convection and higher for liquids than for gases. (Incropera et al. 2007, 8, 355)
In a region where thicknesses of the boundary layers are increasing, the heat flux and so also the convection heat transfer coefficient are not constant, but they decrease due to decreasing velocity and temperature gradients near the surface. The heat transfer coefficient is bigger for turbulent than for laminar flow because due to turbulent mixing velocity gradients and then also temperature gradients are higher near the surface. The value of the convection heat transfer coefficient can be calculated as follows (Incropera et al. 2007, 350, 361, 371):
ℎ𝑐𝑜𝑛𝑣= 𝑁𝑢⋅𝑘
𝐿𝑐 (3.10)
where 𝑁𝑢 is the Nusselt number (-).
A dimensionless parameter called a Nusselt number describes the dimensionless temperature gradient at the surface and it is a function of the Reynolds number and the Prandtl number (Incropera et al. 2007, 371). Nusselt number correlations for different flow situations have been defined by performing experimental tests (Incropera et al. 2007, 403). The Prandtl
number which is a dimensionless parameter is defined as follows (Incropera et al. 2007, 376):
𝑃𝑟 =𝑐𝑝𝜇
𝑘 (3.11)
where 𝑃𝑟 is the Prandtl number (-) and 𝑐𝑝 is the specific heat capacity at constant pressure (J/kgK).
In the following subchapters 3.2.1 and 3.2.2, relevant flow situations relating to the calcula- tion model development are covered. Those situations are flow over the tubes and tube banks, and flow inside the tubes.
3.2.1 Flow over tubes and tube banks
Consider a situation where the tube is in a cross-flow as in figure 9.
Figure 9. Tube in a cross-flow. (Incropera et al. 2007, 423)
At the front of the tube, there is a so-called forward stagnation point where the fluid velocity is zero and the pressure reaches its maximum. When the flow proceeds from that point to- wards the zenith, pressure decreases, velocity increases, and boundary layers develop as in the case of a flat plate. At the zenith, the pressure reaches its minimum and velocity its max- imum. After that, pressure starts to increase, velocity starts to decrease, and boundary layers continue their development. Because velocity decreases after the zenith, at some point ve- locity gradient becomes zero at the surface and the fluid near the surface starts to flow re- versal direction forming vortices. At this so-called separation point the velocity boundary layer comes off from the surface and a wake is generated behind the tube where the flow is irregular. This is illustrated in figure 10. (Incropera et al. 2007, 423)
Figure 10. Development of velocity boundary layer on a tube surface. (Incropera et al. 2007, 423)
The location of the separation point depends on the Reynolds number. If the Reynolds num- ber is big, the laminar boundary layer becomes turbulent before the separation point which leads to that separation occurs later because in the turbulent boundary layer fluid has more momentum to resist velocity decreasing. If the Reynolds number is small, the boundary layer is laminar all the time before separation. This is illustrated in figure 11. The critical Reynolds number for this kind of case is 2⋅105. (Incropera et al. 2007, 423-424)
Figure 11. Turbulent and laminar velocity boundary layers on a tube surface. (Incropera et al. 2007, 424)
In industrial heat exchangers, several tubes form a tube bank. The tube arrangement of the tube bank can be either staggered or aligned which is illustrated in figure 12.
Figure 12. Aligned and staggered tube bank arrangements. (Modified from source: VDI 2010, 726).
Parameters 𝑠𝑇 and 𝑠𝐿 shown in the earlier figure are used to describe perpendicular distances between centerlines of the adjacent tube rows. They are called transverse and longitudinal pitch, respectively (Incropera et al. 2007, 437). When pitches are divided by a tube outer diameter 𝑑𝑜, a transverse pitch ratio 𝑎 = 𝑠𝑇/𝑑𝑜 and a longitudinal pitch ratio 𝑏 = 𝑠𝐿/𝑑𝑜 are obtained (VDI 2010, 726).
In a tube bank, the efficiency of the heat transfer varies. When the flow proceeds inside the tube bank, the efficiency of the heat transfer increases due to increasing turbulence of the flow. The amount of turbulence increases approximately between the first and fifth tube rows and after that stays constant. Heat transfer enhancement inside the tube bank due to turbu- lence also depends on the tube configuration. Heat transfer is more efficient in staggered arrangements because tubes of the next row are not hidden behind the tubes of the previous row, as is the case of aligned arrangements. Therefore, a larger amount of tube surface is directly exposed to the flow. Aligned tube arrangements, where 𝑠𝐿/𝑠𝑇 < 0,7, that is, adja- cent tube rows in a longitudinal direction are very close to each other, are especially bad from the perspective of heat transfer. (Incropera et al. 2007, 441)
3.2.2 Flow inside tubes
Then, consider a situation where the flow is inside the tubes. When the flow enters the tube, hydrodynamic and thermal boundary layers start to develop as in the case of a flat plate. The hydrodynamic boundary layer is growing until layers from opposite sides of the tube join at the centerline of the tube which is illustrated in figure 13. (Incropera et al. 2007, 486-487)
Figure 13. Development of velocity boundary layer inside a tube. (Incropera et al. 2007, 487)
The beginning part of the tube where the hydrodynamic boundary layer is developing is called a hydrodynamic entrance region. After that region, the flow is called to be hydrody- namically developed laminar or turbulent flow depending on that whether the boundary layer became turbulent before the opposing layers merged. For the flow inside the tubes, the crit- ical Reynolds number range is from 2300 to 10000 and the commonly used fixed value is 2300. The velocity profile of the developed flow is parabolic in the case of laminar flow as in figure 13 but in the case of turbulent flow, it is flatter as was discussed already in chapter 3.2. In a turbulent case, the profile looks similar to what the profile is in the entrance region in figure 13. The entrance region is much shorter in the case of turbulent flow and its end- point cannot be determined so accurate than in the case of laminar flow due to velocity pro- files are so similar in the entrance and developed regions. (Incropera et al. 2007, 486-487) In the case of the thermal boundary layer, similar boundary layer development occurs at the beginning part of the tube than in the case of the hydrodynamic boundary layer. In this case, the beginning part of the tube is called a thermal entrance region. If the surface condition is fixed to constant temperature or heat flux, the flow will reach thermally developed condi- tions at some point. The temperature profile of the developed flow depends on the surface condition. In figure 14 profiles are shown for both constant surface conditions. (Incropera et al. 2007, 492)
Figure 14. Development of thermal boundary layer inside a tube. (Incropera et al. 2007, 492)
3.3 Thermal radiation
Thermal radiation is an energy that all forms of matter release when the temperature of the matter is above 0 Kelvin. This phenomenon called emission is due to oscillations and tran- sitions of the electrons inside the matter. These oscillations occur more often when the in- ternal energy of the matter is higher, so higher temperature leads to a higher amount of emit- ted radiation. In contrast to conduction and convection, the transport of radiation does not require any medium because radiation can be considered to be transported by photons or electromagnetic waves, depending on the theory. Thermal radiation does not include all of the spectrum of the electromagnetic radiation but only wavelengths approximately from 0,1 to 100 micrometers. Only this kind of radiation is generated by the temperature of the matter and affects the temperature of the matter. Thermal radiation can be either a volumetric or a surface phenomenon depending on the medium. Emission from gases or semitransparent solids is a volumetric phenomenon. Emission from most of the solids and liquids, in turn, is a surface phenomenon because inside the matter adjacent molecules absorb most of the ra- diation emitted by interior molecules. Therefore, only radiation from near the surface can reach the space adjacent to emitting solid or liquid. (Incropera et al. 2007, 724-726)
Thermal radiation has two inherent features that make it difficult to solve radiation problems.
Firstly, thermal radiation has a spectral distribution which means that total radiation consists of many wavelengths, and the magnitude of the radiation is different for different wave- lengths. Secondly, thermal radiation has a directional distribution which means that radiation can differ depending on its direction. These features are illustrated in figure 15. (Incropera et al. 2007, 726)
Figure 15. Spectral and directional distributions of thermal radiation. (Incropera et al. 2007, 726)
The radiation heat transfer rate can be expressed using a similar type of equation as the convection heat transfer rate. The equation below describes the radiation heat transfer from a small surface to a much bigger surface surrounding a smaller one (Incropera et al. 2007, 10):
𝜙 = ℎ𝑟𝐴𝑠(𝑇𝑠− 𝑇𝑠𝑢𝑟𝑟𝑜𝑢𝑛𝑑𝑖𝑛𝑔) = (𝑇𝑠−𝑇𝑠𝑢𝑟𝑟𝑜𝑢𝑛𝑑𝑖𝑛𝑔) 1
ℎ𝑟𝐴𝑠
(3.12) where ℎ𝑟 is the radiative heat transfer coefficient (W/m2K).
The equation for resistance can be found as usual and it is as follows:
𝑅𝑟𝑎𝑑 = 1
ℎ𝑟𝐴𝑠 (3.13)
where 𝑅𝑟𝑎𝑑 is the thermal resistance for radiation (K/W).
Expressing the radiation heat transfer similar way than convection heat transfer by using heat transfer coefficient is a good way in practical applications, although the radiation heat transfer coefficient has not a very strong connection to physical reality (Vakkilainen 2016, 133). This is because with the heat transfer coefficient radiation is made proportional to the linear temperature difference. In the following subchapters from 3.3.1 to 3.3.5, different thermal radiation situations are covered, and it will be noticed that the radiation is propor- tional to the difference of two temperatures which have been raised to the fourth power.
(Incropera et al. 2007, 10-11) 3.3.1 Radiation from surface
A so-called ideal radiator, also called a blackbody, needs to be defined when radiation from the surface is under consideration. The blackbody emits radiation perfectly which means that it emits the maximum amount of radiation that is possible for the surface at a certain tem- perature. In the case of blackbody emission, the radiation has no directional distribution.
Thus, the blackbody is a so-called diffuse emitter. The heat flux emitted by the blackbody can be calculated as follows (Incropera et al. 2007, 736-738):
𝐸𝑏 = 𝜎𝑇𝑠4 (3.14)
where 𝐸𝑏 is the emissive heat flux of the blackbody (W/m2) and 𝜎 is the Stefan-Boltzmann constant (5,67⋅10-8 W/m2K4).
Real surfaces do not emit as much radiation as the blackbody. The heat flux emitted by real surfaces is calculated as follows (Incropera et al. 2007, 9-10):
𝐸 = 𝜀𝐸𝑏 = 𝜀𝜎𝑇𝑠4 (3.15)
where 𝐸 is the emissive heat flux (W/m2) and 𝜀 is the total hemispherical emissivity (-).
The value of the emissivity is between 0 and 1 so it describes the relative emissive power of the real surface compared to the black surface at the same temperature (Incropera et al. 2007, 9-10). Emissivity can be different for different wavelengths and directions, so the total hem- ispherical emissivity used in the above equation represents an average value over the entire wavelength and direction ranges (Incropera et al. 2007, 746). The total hemispherical emis- sivity is dependent on the surface temperature (Incropera et al. 2007, 755).
When the total radiation leaving from the surface is considered, also the radiation that hits the surface from its surroundings must be taken into account. Incoming radiation to the surface is called irradiation. If the surface is semitransparent, three options can occur to the irradiation. These three possible options are reflection, absorption, and transmission, as shown in Figure 16. (Incropera et al. 2007, 752-753)
Figure 16. Irradiation to a semitransparent surface. (Modified from source: Incropera et al. 2007, 753)
The radiation balance of the surface in the above situation is as follows:
𝐺 = 𝐺𝑟𝑒𝑓+ 𝐺𝑎𝑏𝑠+ 𝐺𝑡𝑟𝑎 = 𝛼𝐺 + 𝛽𝐺 + 𝜏𝐺 (3.16) where 𝐺 is the total irradiation (W/m2), 𝐺𝑟𝑒𝑓 is the reflected irradiation (W/m2), 𝐺𝑎𝑏𝑠 is the absorbed irradiation (W/m2), 𝐺𝑡𝑟𝑎 is the transmitted irradiation (W/m2), 𝛼 is the total hemi- spherical absorptivity (-), 𝛽 is the total hemispherical reflectivity (-) and 𝜏 is the total hemi- spherical transmissivity (-).
Once again, it must be noted that the shares of the different components strongly depend on the wavelength and the direction of the radiation. Therefore, components in the above equa- tion are averaged values over all spectrum and all directions. Temperature cannot be consid- ered to have a large effect on these portions (Incropera et al. 2007, 755). In most engineering applications, the reflection is usually assumed to be diffuse instead of specular because this assumption is close to the truth for rough not-mirror-like surfaces. The difference between these two reflection types is illustrated in figure 17. (Incropera et al. 2007, 752-757)
Figure 17. Reflection from a rough (on the left) and a mirror-like (on the right) surface. (Incropera et al. 2007, 757)
In the case of a blackbody, the surface absorbs radiation perfectly which means that it com- pletely absorbs all the radiation that hits it, no matter what the wavelength or incoming di- rection is. Thus, there is not any reflection or transmission in the case of a blackbody (In- cropera et al. 2007, 736).
The total radiation that leaves from the surface is the sum of the emission and reflection and it is called radiosity (Incropera et al. 2007, 823):
𝐽 = 𝐸 + 𝛽𝐺 (3.17) where 𝐽 is the total radiosity from the surface (W/m2).
The net radiation heat transfer rate that leaves from the surface is the difference of radiosity and irradiation as follows (Incropera et al. 2007, 823):
𝜙 = 𝐴𝑠(𝐽 − 𝐺) (3.18)
For an opaque surface (no transmission so 𝛼 = 1 − 𝛽), the net power can also be expressed using emission and absorption (Incropera et al. 2007, 823):
𝜙 = 𝐴𝑠(𝐸 − 𝛼𝐺) (3.19)
If the opaque surface can be assumed to be gray (𝛼 = 𝜀) and diffuse, the radiosity can also be expressed as follows:
𝐽 = 𝐸 + 𝛽𝐺 = 𝐸 + (1 − 𝛼)𝐺 = 𝜀𝐸𝑏+ (1 − 𝜀)𝐺 (3.20)
From equation 3.20, 𝐺 can be solved as follows:
𝐺 =𝐽−𝜀𝐸𝑏
(1−𝜀) (3.21)
When 𝐺 is placed to the equation 3.18, the following is obtained:
𝜙 = 𝐸𝑏−𝐽
(1−𝜀)/𝜀𝐴𝑠 (3.22)
Equation 3.22 is useful for solving complex radiation exchange problems between many surfaces because it is analogous to electrical networks. In a numerator, there is a driving potential, and in a denominator, there is a so-called surface radiative resistance (Incropera et al. 2007, 823-824).
When equation 3.22 was derived, two assumptions were made which help to solve radiation problems. Firstly, the surface was assumed to be diffuse which means that emissivity and absorptivity are independent of direction. Secondly, total emissivity and total absorptivity
were assumed to be equal. This is true only in specific circumstances. Spectral directional emissivity and absorptivity are always equal (𝜀𝜆,𝜃 = 𝛼𝜆,𝜃) because they are inherent proper- ties of the surface. Thus, the spectral and directional distributions have no effect on them.
Spectral emissivity and absorptivity are equal (𝜀𝜆 = 𝛼𝜆) if the irradiation that hits the surface or the surface itself is diffuse, that is, independent of direction. Irradiation can be assumed to be diffuse for many surfaces existing in engineering applications. Especially for noncon- ducting materials, this is the truth at a reasonable level. Total hemispherical emissivity and absorptivity are equal if the irradiation to the surface equals the emission from the same temperature blackbody or if the surface is a so-called gray surface, that is, emissivity and absorptivity are not dependent on wavelength distribution. The surface can be assumed to be gray if emission from the surface and irradiation to the surface are between the same wavelength region where spectral emissivity and absorptivity are roughly constant. This sit- uation is illustrated in figure 18. (Incropera et al. 2007, 764-766)
Figure 18. A situation where a gray surface assumption can be applied. (Incropera et al. 2007, 766)
3.3.2 Radiation between surfaces without participating gas
When a radiation exchange is considered between two or more surfaces in an enclosure without gas participating to the radiation exchange, a view factor is an essential term to know. For example, the view factor Fij describes the portion of the radiation originating from the surface i which ended up to the surface j. So, the radiation energy transferred from the surface i to the surface j is as follows (Incropera et al. 2007, 812-820):
