Use of CFD Codes for Calculation of Radiation Heat Transfer: From Validation to Application Use of CFD Codes for Calculation of Radiation Heat Transfer: From Validation to Application

In complex geometries, computational fluid dynamic (CFD) codes are commonly used to predict the heat and fluid transfer. To justify their use for the applications with dominant radiation heat transfer conditions, the implemented models need to be first appropriately validated on simple benchmark examples where the analytical solutions exist. The practical application discussed in this chapter considers the thermal radiation inside the vac- uum vessel of the fusion reactor. Two representative benchmark examples are used to obtain the analytical solution and assess the accuracy of the real case simulations perfo- rmed by CFD codes. The analytical solutions use the view factor method to calculate the net radiation heat flux on each radiating surface. Several numerical methods are available in the CFD codes to solve the thermal radiation problems. The discrete transfer method (DTM) is considered as one of the most efficient for solving the radiation fluxes between the surfaces in the case of radiatively non-participating fluid. Discussion includes description of fundamentals of analytical and numerical thermal radiation methods, validation of radiative heat exchange in simple enclosure problems, estimation of numerical errors and application to the practical case.


Introduction
There are many applications where thermal radiation plays an important role. To name just a few of them, atmospheric physics, astrophysics, astronautics, remote sensing, nuclear engineering and many other applications may be added on the list. The engineering example in this study considers the thermal radiation exchange in the vacuum vessel of the fusion reactor.
In all cases where the thermal radiation is significant, the proper choice of the thermal radiation model affects the accuracy of the solution and also the amount of time required for the computation. Accurate solutions are computationally very expensive; therefore, a special care must be taken to select an appropriate calculation method that corresponds to the physics of the radiative heat transfer under consideration. For example, different methods can be applied in the situations where the medium between the enclosing solid surfaces behaves as optically dense or it is practically transparent to the heat transferred by radiation. In the first situation the radiative heat can be scattered, absorbed or re-emitted in the medium and at the solid surfaces. The radiation models used for participating medium are usually based on the solution of general thermal radiation equation [1]. In the second situation the medium is either not present (vacuum) or it is transparent to the wavelengths at which the thermal radiation occurs. In this case the radiative heat is transferred solely between the boundary surfaces and depends on the surface properties and geometrical orientation that each surface has to the others. This type of thermal radiation is relevant for most of the engineering applications. For simple geometries analytical solutions based on the view factor evaluation exist [2,3]. In the real cases with complex surface geometries, the exact analytical expressions are not available and approximate numerical methods must be used, such as Monte Carlo [4] or discrete transfer method (DTM) [5]. The DTM method has some advantage due to its computational efficiency, easy application to complex geometries and implementation into the CFD codes. The method can be very accurate, but it has a shortcoming, as it is not easy to assess the accuracy of the net radiative heat flux on the individual surface in the complex enclosure. Even if the total energy balance in the enclosure is conserved, it is not self-evident that the radiation heat flows on individual surfaces are calculated correctly. To assess and improve the accuracy of the CFD simulation by the DTM method, a validation against the simple example with the exact analytical solution that resembles the real case is highly recommended.
The following sections give a description of analytical and numerical methods. Analytical method is based on the view factor calculation, reciprocity and summation relations in an enclosure. Description of numerical method considers DTM method as implemented in the CFD code ANSYS CFX [6]. The accuracy of numerical results is validated on the two analytical examples, cylinder and closed ring. The geometry and dimension of the latter approximately resemble the realistic case of the heat radiation inside the vacuum vessel of fusion reactor, the results of which are presented in the last section. The chapter ends with the summary of main conclusions.

Analytical solution method for thermal radiation in an enclosure
The theory of radiation exchange between the surfaces described in this section is based on two assumptions. First, the surfaces form an enclosure and second, they are separated by a medium that does not participate in thermal radiation. Radiatively non-participating media has no effect on the radiation transfer between the surfaces. There is no scattering, emission or absorption in the medium. Such conditions occur in vacuum, and also in monatomic and most diatomic gases at low and moderate temperatures, before the ionization and dissipation occurs. In fact, in many engineering applications medium does not affect the radiation heat transfer [2].
The radiation heat exchange between the two differently oriented black surfaces with final dimensions can be expressed as [2,7]: where σ ¼ 5:67•10 À8 Wm À2 K À4 is the Stefan-Boltzman constant and F ij is the view factor. The net radiation heat exchange between the surface A i at absolute temperature T i and an enclosure of N black surfaces at absolute temperature T j may be described by the energy balance on the opaque surface A i ( Figure 1): Absorption depends on irradiation, which depends on emission from other surfaces including those far away from the observed surface. To calculate the total radiative energy balance the entire enclosure has to be considered. Thus, all radiation contributions are accounted for. An open enclosure is in practice closed by introducing artificial surfaces. For example, opening can be considered as a surface with zero reflectivity or as a radiation source when presenting environmental radiation. Enclosure is usually composed of complex geometries which may complicate the calculations. To deal with the complex geometry, the enclosure can be approximated by several simple surfaces that are assumed to be isothermal, as indicated in Figure 1.
Radiation exchange between the surfaces in addition to their radiative properties and temperatures strongly depends on the surface geometries, orientations and distances between them.  This leads to the development of a geometric function known as view factor. The view factor F ij is defined as the fraction of the radiation leaving the surface A i that is intercepted by the surface A j and is calculated a general expression [3]: where S is the distance between the surfaces A i and A j , θ i and θ j are the angles between the surface normal vectors n i , n j and S. (Figure 2).
Two very useful view factor relations are valid for the enclosure. The first one is the reciprocity . The latter relation follows the conservation requirement that all radiation leaving the surface A i must be intercepted by some other surface A j . As shown later, these two relations are not useful only for the analytical calculations but are also important for the assessment of the accuracy of numerical methods. Reciprocity relation can be used to check the accuracy of individual view factor and the summation rule can be used for validation of the energy conservation. The view factor F ii deserves a special consideration. If the surface is convex, no radiation leaving A i will strike itself so F ii ¼ 0. If the surface is concave, it sees itself and part of radiation leaving A i will be intercepted by itself so F ii 6 ¼ 0.
To calculate the radiation exchange in the enclosure of N surfaces a total of N 2 view factors are required. The view factors on N surfaces can be written in a matrix form:   view factors may be obtained from the reciprocity relations. When the enclosure in the model includes M surfaces that cannot see itself F ii ¼ 0 ð ), the total number of the view factors that need to be calculated directly amounts: In some cases, the use of symmetry can additionally reduce the number of directly calculated view factors.
Even for simple geometries the analytical calculation of view factors, Eq. (3) is not easy. The calculation of view factor between the two finite surfaces requires solving of the double area integral, or fourth-order integration. Such integrals are difficult to evaluate analytically except for very simple geometries. For practical use, the view factors can be generated from already known solutions for the frequently used geometries that are collected in the form of tables and charts. The most complete set of solutions is given in a catalog of Howell [3]. For more complicated geometries, the view factors need to be calculated by numerical integration that can be computationally expensive, depending on the complexity of the geometry.
The energy balance on the selected surface A i in an enclosure of N black surfaces is calculated using Eq. (3), where the emission j i, emission and absorption are expressed as: The net radiation heat from the surface A i at absolute temperature T i to the enclosure of N black surfaces at absolute temperature T j can be expressed as [2]: where F ij is the view factor between the surface A i and one of the enclosing surfaces A j . Radiation properties such as emissivity (ε i ¼ 1 for the black body) and temperatures are set, surface areas are easy to calculate, whereas the calculation of the view factors F ij can be a very difficult task for complex geometries.

The cylinder example
To demonstrate the use of analytical methods, the radiation heat transfer in a simple model of an open cylinder is calculated. Diameter and height of the cylinder are 50 and 150 mm, respectively. The cylinder is open at the top surface to a large surroundings at T sur ¼ 27 C: The bottom and the side surfaces of the cylinder are approximated as black surfaces and are maintained at T F ¼ 1500 C and the opening at the top is approximated as blackbody at the temperature of the surroundings. The goal of the exercise is to calculate the heat exchanged between the surfaces assuming that the thermal radiation is the only heat transfer mode and the outer backs of surfaces are adiabatic. The sketch of the cylinder with the bottom surface A 1 , the side surface A 2 and the top surface A 3 is shown in Figure 3.
To calculate the radiation exchange between the three surfaces, in general, nine view factors are needed. According to Eq. (5), only N NÀ1 ð Þ 2 À M view factors need to be calculated directly. Taking into account the two flat (A 1 and A 3 on Figure 3) which cannot see itself (M ¼ 2) and one concave surface (A 2 ) only one view factor has to be calculated directly, the others can be derived using the view factor relations as follows: where, The unknown view factor F 13 can be calculated analytically by solving the integral in Eq. (3). The inside-sphere method [2] was found to be very convenient for the considered geometry, detailed derivation of the view factor F 13 is provided in [8]. The easiest way to obtain the view factor for the sought geometry is to generate it from the available database of already calculated view factors for similar geometry configurations. In the catalog of Howell [3], the view factor for the two parallel coaxial disks of unequal radius can be found (the case C-41). In our case the disks (surfaces A 1 and A 3 ) have equal radius, the view factor F 13 then yields: where H equals to h=2r. The remaining view factors obtained using the view factor relations (Eq. (9)) are given in Table 1.
Using the Eq. (8), we can calculate the net radiative heat transfer from each surface. The net radiation loss to the surroundings through the open surface A 3 amounts À1521.6 W. Analytical results for the net radiation heat from the surfaces are collected in Table 2. The heat flow balance in the enclosed cylinder is the sum of all radiation heat flows and should be equal to zero.
The view factor  Use of CFD Codes for Calculation of Radiation Heat Transfer: From Validation to Application http://dx.doi.org/10.5772/intechopen.74420

Numerical simulation methods
Thermal radiation exchange in complex geometry configurations is calculated using the numerical simulation methods. The accuracy of numerical simulations depends on the numerical mesh and on the type of numerical methods used. For accurate results it is important to understand the main characteristics of the method to efficiently reduce the accompanying numerical errors. To validate the accuracy of the selected method a comparison of the results on a simpler geometry that has an exact analytical solution is of great importance.
Different numerical methods are available for solving the complex thermal radiation problems. In general, the radiation transfer through a medium is affected by absorption, emission and scattering and can be described by the generic radiation transport equation: where I is radiation intensity which depends on position ( r ! ) and direction (s ! ), I b is blackbody emission intensity, K a and K S are the absorption and scattering coefficients of the medium, ϕ is the scattering function and s is the path length. Analytical solutions for Eq. (11) exist only for very simple cases. There are several numerical methods used to predict the thermal radiation, based on Eq. (11) [1,9].
One of the most efficient methods for solving the thermal radiation between the surfaces is the discrete transfer method (DTM) developed by Lockwood and Shah [5]. The DTM solves the simplified form of the radiation transport equation (Eq. (11)) along rays (see Figure 4). To determine the direction of the rays, the unitary hemisphere over the element face is discretized using spherical coordinates. The span is divided into angles by the ray number, and rays directions are computed to pass through the center of the angles. In total, the square of ray number is traced from an element of the surface.
The rays are fired in prescribed directions from discrete point P i , which is located at the center of a cell face on the boundary surface. Each ray is traced until it hits another boundary and Q j is the impingement point. In general Q j is not the central point of a boundary cell, but it is assumed that radiation intensity at Q j and at the central point P j are equal. The boundary conditions or initial values of radiation intensity at point Q j for gray, diffuse surface depends on the incident radiation intensity at point P j , which further depends on the radiation intensity of all other rays that reach point P j (see Figure 4). Simplified radiation transport equation (Eq. (11)) is discretized in 3D finite volumes along the rays. The path along the ray is discretized using the sections formed from breaking the path at volume boundaries. The physical quantities in each volume are assumed to be uniform.
Due to the fixed sampling and ray discretization the physical quantities can be found at fixed points. The accuracy of the DTM is controlled by the number of rays and by the mesh density.
For accurate results the control volumes (mesh) must be chosen in a way that the irradiation field is reasonably homogeneous inside them. The major problem of the DTM is the lack of error information. Large errors can be produced when the ray sample misses the sensitive area or object. This error is known as the ray effect.
Unless the surfaces are black (ε i ¼ 1), the radiation intensity of selected surfaces depends on the radiation intensity of all other surfaces (see Figure 4) and solution requires iterative calculation procedure. The detailed description of the numerical solution procedure can be found in [9,10].
In this study the DTM method implemented within the ANSYS CFX [6] computational code is tested. The option that includes "non-participating media" is the most appropriate for the  Table 3. The numerical results of the same case are obtained by the ANSYS CFX code using the S2S thermal radiation model for the nonparticipating media. For the purpose of sensitivity analysis two different ray numbers on the single numerical mesh with 72,000 hexahedral mesh elements.
Numerical results presented in the third column of Table 3 are calculated using the ray number 48 (actual number of rays is 48 2 ¼ 2304) on the mesh with 72,000 mesh elements and 7200 mesh faces. This means that 2304 rays are fired from the center of each face. Hence, the whole number of rays involved in the calculations is 7200 Â 2304 ¼ 16 588 800. The results in the Table 3 show that the ANSYS CFX prediction of the radiative heat transfer from the surface A 1 is the least accurate. Better accuracy can be achieved by increasing the number of rays, which is evident from the results for the ray number 128 in the fourth column of the Table 3. The overall heat flow balance in the enclosure is very well preserved for both numerical simulations and as such cannot provide the information about the numerical accuracy for the individual heat exchanging surface. Only the comparison with the analytical solution can give this type of information.

Validation of the closed ring example
As a second validation case a closed ring geometry with a rectangular cross-section was selected (see Figure 5). The ring geometry tends to resemble the realistic model of the fusion reactor torus as far as possible but is still simple enough to allow the analytical solution that enables proper validation of numerical results. The inner radius (6 m), outer radius (12 m) and height (11 m) of the ring ensure that heat radiation surfaces of the ring approximately match the inner surfaces of the DEMO in-vessel components. In the Figure 9 the ring dimensions (represented by the yellow rectangle) are compared with the dimensions of the DEMO tokamak. The ring geometry contains a small surface (red rectangular surface in Figure 5) representing the upper surface of the  detached divertor cassette that is kept at a lower temperature (323 K) than the remaining ring surfaces (373 K). The imposed surface temperatures are the same as the temperatures in the DEMO reactor model (see the next Section 4). All surfaces of the ring enclosure are approximated as blackbodies. Analytical results for the net radiation heat transfer on different ring surfaces (as marked in Figure 5) are presented in Table 4 and are calculated using the Eq. (8). Analytical values of view factors were obtained from the catalog of known view factor solutions [3] or directly calculated from Eq. (3).  For validation purposes, the numerical results are computed on the three different hexahedral meshes using the ANSYS CFX code. The meshes are presented in Figure 6. The Coarse mesh in Analytical results of the net heat transfer on the ring surfaces are compared with the numerical results on Table 4. Numerical simulations are performed by ANSYS CFX using the discrete transfer model (DTM) for non-participating media (S2S). The negative net heat transfer value on the surface A5, which represents the detached cassette at 323 K, means that it receives the thermal radiation from other components that are at higher temperatures of 373 K. The other surfaces with positive net heat flow are net emitters. The simulation results in Table 4 are obtained for three different mesh densities using the ray number 48 (actual number of rays per surface element is 2304).
The total sum of net heat flows (heat flow balance) is also shown in Table 4. Ideally (e.g. analytical solution) the net emitted heat flow in an enclosure is equal to the net absorbed heat flow leading to the zero heat flow balance. As shown the total heat flow balance is very accurately predicted for all mesh densities. This means that the energy flow inside the domain is well preserved in all cases, but does not provide any information about the accuracy of thermal radiation calculation. The total relative error r tot of the net thermal radiation transfer in a model can be estimated using the equation: where P i is the net heat flow on the i-th surface. The relative error r i on the i-th surface is defined as the ratio between the difference of analytical P i, an and numerical P i, num net heat flow  Table 4 show that the total relative error r tot of the numerical simulation reduces with the mesh refinement. The computational time increases with denser mesh (2nd row in Table 4). The Medium mesh provides the best compromise between the accuracy of results and the computational time.
The results for three different meshes and three different ray numbers are shown in Table 5.
The computational time expectedly increases with the number of rays on each mesh. In the case of a Coarse and Fine mesh the increased ray number improves the accuracy of the results, whereas this does not seem to be the case for the Medium mesh.
The dependence of numerical solution on the number of rays, computational time and mesh resolution is presented in Figure 7. The ray numbersquare root of the real number of rays is shown in the abscissa of Figure 7. Corse, Medium and Fine mesh are marked by blue, green and red color, respectively. Medium mesh is tested in detail by using denser sampling of ray numbers, also including the odd number of rays.  From the Figure 7a, it can be seen that the computing time increases with the number of rays and with the mesh density. In general the curve of the computing time dependence has the same shape for all mesh densities only the curves are shifted upwards with the higher number of mesh cells. Finer mesh requires high computational resources already at the low ray number.
The accuracy of the DTM method is controlled by the number of rays. Figure 7b it shows that the relative error decreases with increasing number of rays and that the refinement of the mesh reduces the numerical error. Further it can be seen that the relative error decreases very slowly and then remains at a certain level when the number of rays increases over the specific number of rays (e.g. 30 for the medium meshgreen curve). In the case of Medium mesh it can be seen that the difference in the accuracy of the even (e.g. 64) an odd (e.g. 65) ray number is relatively high. Obviously the odd number of rays decreases the numerical accuracy.
The number of rays has to be determined in a way that the accuracy is acceptable and that the simulation is not computationally too demanding. The number of rays used in the following simulations is obtained by comparing the Figure 7a and 7b. Fine mesh requires more computational time than the Medium mesh for the same relative error. Fine mesh has less than two times better accuracy than the Medium mesh at the same number of rays (ray number 48) but requires almost six times more computing time. With the ray number 48 it is feasible to achieve a good accuracy with the reasonable computing time. The total relative error for the numerical solution is equal to 0.03% and is approximately evenly distributed over all surfaces.

Application on the vacuum vessel of the fusion reactor
The results of the validation cases can provide a useful information on the mesh density and the ray number for the DTM method. Based on the closed ring case, a similar mesh density and ray number are applied in the ANSYS CFX input model for the real DEMO tokamak geometry.
The presented practical example focuses on the heat loading of the divertor cassette immediately after the shutdown of the DEMO fusion reactor. The demonstration fusion reactor DEMO is planned to be the last major step before the commercial fusion power plant (Figure 8a). During the reactor operation the divertor (Figure 8b) is subjected to a high incident heat flux of removed plasma particles with the values above 10 MW=m 2 . Such heat loads may eventually cause severe damaging and consequently the need for regular replacement of the divertor cassettes that is envisaged at a 2-year cycle.
The component of interest is one of the 54 divertor cassettes that is being replaced during the regular maintenance. The cassette under replacement (Figure 8c) is unplugged from the cooling pipes, while the remaining cassettes and the blanket remain actively cooled. Because of the lack of internal cooling the detached cassette is heated up due to the decay heat in activated structure materials. Thermal radiation from the colder surfaces of the surrounding divertor cassettes and blankets represents the only cooling mode of the detached cassette.
The assumed conditions during the maintenance (see Figure 9) are defined as follows [11]: • Blankets are actively cooled at T ¼ 373 K.
• All the remaining divertor cassettes are actively cooled at T ¼ 373 K: • Vacuum vessel is set at room temperature T ¼ 300 K.
Boundary conditions on the detached divertor cassette (see Figure 10) are set for the purpose of this study: • At the plasma facing surface the temperature is either fixed at 323 K or the passively cooled surface is assumed.
• Temperature of the side surfaces facing actively cooled components is set to 373 K.
• At surface facing the vacuum vessel the temperature is set to 300 K.
• The decay heat just after the shutdown is adopted as 25 kW=m 3 [11].  Taking into account the detached divertor cassette the toroidal symmetry of the tokamak is not preserved, therefore the full 360 geometry of the vessel has to be modeled. Since the heat load on the detached divertor cassette needs to be evaluated, it is sufficient to consider only the invessel components (blanket and divertor) and vacuum vessel without the manway ports. The ANSYS CFX model of DEMO vacuum vessel with one detached divertor cassette (out of 54) is shown in Figure 9. Components included in the model are additionally simplified: all gaps between the actively cooled divertor cassettes, blanket and vacuum vessel are neglected as they have a negligible effect on the solution. The gap around the detached divertor cassette has been considered with the appropriate boundary conditions. The simulation model is set up as a closed cavity radiation problem with surfaces approximated as gray bodies (ε i ¼ 0:25) [11].
In addition to the thermal radiation between the system surfaces also the heat conduction in the detached cassette and internal heat generation due to the decay heat in activated structure materials has to be considered.
Three cases have been simulated by applying different boundary conditions on the detached cassette. The simulation results are presented in Table 6. The most realistic case is the Case 1, where the internal heat generation is considered inside the cassette body and its top surface is modeled as a passive boundary condition. The second case (Case 2) is the most similar to the conditions of the closed ring validation case, therefore it does not include the internal heat generation and the top surface temperature of the detached cassette is set to 323 K (50 C). In the Case 3, besides the internal heat generation, an additional surface cooling of the detached cassette is modeled by imposing its top surface temperature to the value 323 K. Comparing to the Case 1, the Case 3 aims to evaluate the effect of external surface cooling on the reduction of maximum temperature inside the cassette body. The results presented in Table 6 are calculated by the ANSYS-CFX using the surface-to-surface thermal radiation model with the ray number 48. The mesh shown on Figure 9 has 698,898 mesh cells.
The Table 6 includes the heat flow balance inside the tokamak enclosure that shows perfect energy conservation in the tokamak (the total sum of emitted and received heat flows is zero in all cases). The overall error represents only 10 À6 % or less of all exchanged radiation heat inside the tokamak. Based on the results obtained by comparing the numerical and analytical solutions of the closed ring example, we are aware that the accurately calculated heat flow balance does not mean that the heat flows on an individual surface are equally precisely predicted. Taking into account the results of the closed ring example, we may assume that the estimated error for an individual surface is of the same order as the error in the benchmark model, considering that the sufficiently dense mesh is used. Also, it is assumed that the emissivity does not affect the error, since the emissivity is not a geometry dependent parameter. Relative error for the individual net heat flow is thus estimated at 0.03%.
The calculated heat flows in the Cases 2 and 3 are exactly the same, as the cases have the same radiative boundary conditions. The negative heat flow on the detached divertor cassette in cases 2 and 3 means that the detached cassette (323 K) receives the radiation heat from other components at a higher temperature (373 K), which is also in accordance with the results of the closed ring example. In the Case 1, the heat flow on the detached cassette is positive, which means its top surface is hotter than the surrounding surfaces. In this case, the detached cassette is passively cooled by the thermal radiation. However, the temperature at the surface is still rather high (420 K). Total sum of heat flows 0.000 0.000 0.000 Table 6. Numerical results and estimated heat flows for three different cases.
Use of CFD Codes for Calculation of Radiation Heat Transfer: From Validation to Application http://dx.doi.org/10.5772/intechopen.74420 In addition to the net heat flows, Table 6 also shows the maximum temperature inside the detached cassette body and on its top surface. In the Cases 1 and 3, it can be seen that the maximum temperature increases due to the internal heat generation and amounts to 423 and 395 K for the Case 1 and Case 3, respectively.
The temperature distribution inside the detached cassette for the Cases 1 and 3 is presented in Figure 10. Due to the external cooling (imposed temperature of 323 K on the top surface), the region of maximum temperature in the Case 3 is lower and displaced toward the inside of the cassette body (Figure 10a). The temperature peak for the Case 1 is significantly higher (423 K) and located closer to the upper surface (see Figure 10b).

Conclusions
To justify the use CFD codes for the applications with dominant radiation heat transfer, the implemented models need to be first appropriately validated on simple analytical examples.
The practical application under consideration is the heat load on the detached divertor cassette inside the DEMO vacuum vessel after the reactor shutdown. The detached cassette is subjected to the internal heating and is cooled solely by thermal radiation from the surrounding in-vessel components. The thermal analysis is performed with the CFD code using the discrete transfer method (DTM) for the thermal radiation modeling. method applied within the ANSYS CFX code gives accurate predictions of the thermal radiation in the complex geometry configurations.