Numerical simulation of radiative heat transfer in indoor environments on programmable graphics hardware

International Journal of Thermal Sciences 96 (2015) 345e354

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Numerical simulation of radiative heat transfer in indoor environments on programmable graphics hardware € sler b, Clemens Felsmann b Stephan Kramer a, *, Ralf Gritzki b, Alf Perschk b, Markus Ro a b

€ t Go €ttingen, Lotzestraße 16-18, 37073 Go €ttingen, Germany Institut für Numerische und Angewandte Mathematik der Universita €t Dresden, D-01062 Dresden, Germany Institut für Energietechnik, Technische Universita

a r t i c l e i n f o

a b s t r a c t

Article history: Received 1 April 2014 Received in revised form 16 February 2015 Accepted 17 February 2015 Available online 17 March 2015

The efficient use of energy for heating and cooling of indoor environments requires an accurate prediction and analysis of radiative heat transfer. Therefore it is necessary to use modern computer methods as otherwise the computational costs may become higher than for convective heat transfer, which is known to be a huge computational problem. A key problem in calculations of radiative heat transfer is the problem of mutual visibility arising in the determination of view factors. This is the same challenge that global illumination has to cope with which is one of the fundamental topics in computational graphics. The visibility problem is efficiently solved by modern graphics hardware. Therefore an OpenGLbased algorithm is developed to quickly and accurately calculate view factors for arbitrary, complex geometries. Theoretical and implementation details of the applied methods are given. We demonstrate the advantages of the developed computational method by a virtual test room for a tubular radiator, a heating of a warehouse by ceramic infrared heaters and the heat transfer in a car cabin. © 2015 Elsevier Masson SAS. All rights reserved.

Keywords: Obstructed view factors Radiative heat transfer Parallelization Hemicube OpenGL Radiator Indoor environment

1. Introduction Radiation is a very efficient way of heat transfer, as it is known from many natural and technical processes. It is of paramount importance for heating or cooling of indoor spaces. Accurate prediction of energy consumption in indoor environments requires to solve the combined problem of convective airflow and radiative heat transfer (RHT) on a similar spatial scale. Only this allows to design heating, ventilation and air conditioning (HVAC) systems which offer a guaranteed level of thermal comfort while minimizing the energy consumption. For heating systems for flats and offices low operating temperatures are often preferred to minimize energy losses. Modern isolation techniques reduce the demands on heating power as well, relaxing the requirements on operating temperatures even further. Especially modern, environment-friendly and cost-efficient ways of heating like geothermal heat pumps and passive cooling are restricted in the available temperature differences. To realize the necessary amount of radiant power large-surface heating devices

* Corresponding author. Present address: Max-Planck-Institut €ttingen, Germany. Biophysikalische Chemie, Am Fabberg 11, 37077 Go E-mail address: [email protected] (S. Kramer). http://dx.doi.org/10.1016/j.ijthermalsci.2015.02.008 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.

für

have to be employed, cf. [30]. Another way is to efficiently burn gas in overhead luminous radiant heaters where very high heat flux densities arise which require special attention in a simulation. In both cases a prediction of operative temperatures and other thermal comfort criteria is only possible by detailed computation of the RHT with explicit spatial resolution, i.e. by solving the integral equation. It is common practice to model the indoor airflow as lowturbulent hydrodynamic flow [27,33] and to solve it with the necessary high resolution in space and time. However, with these new requirements RHT may become the computational bottleneck, a role previously played by convective heat transfer. For both phenomena, convection and radiation, it is known that the quality and expense of their simulation can be substantially improved by using nowadays programmable graphics cards as cheap and highly efficient numerical coprocessors [6,42]. In the Lambertian approximation RHT is described by a Fredholm integral equation of the second kind. The resulting matrices are rather well-conditioned such that iterative algorithms converge quickly. Furthermore, using modern compression techniques like adaptive cross approximation [4,24,37,38] or wavelet-based compression [39] the view factor matrix can be stored in an efficient way such that only O ð1Þ matrix entries are necessary per degree of freedom (DoF) instead of the usual O ðNDoF Þ, where NDoF is the number of degrees of freedom. This usually is the number of

346

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

surface elements subject to RHT. Based on these compression techniques efficient solvers for the associated system of linear algebraic equations can be devised [2,3]. Therefore, the issue in RHT computations is not so much solving but rather setting up the discretized problem, i.e. to determine the coefficients of the linear algebraic equations. In the standard setting of isothermal diffuse, gray surfaces the coefficients in the equations are the view factors which represent the fraction of radiant energy emitted by a given surface that is intercepted by another one [32]. In the case of unobstructed surfaces view factors can be computed easily from closed formulas [21,22]. However, in realistic setups occlusion effects are not an exception but rather the rule. Thus for each pair of surface elements one has to numerically determine their effective mutual visibility. It is this visibility problem which makes the determination of accurate area-to-area (A2A) view factors in the presence of occlusion effects so costly. Besides the visibility problem the accuracy of a radiative transfer calculation also depends on the quality of the chosen mesh, the validity of the assumption that the radiating surfaces are diffuse and that the enclosed medium does not influence the radiative transport, which is not always the case. Whether surfaces can be considered as diffuse depends on the spectral range of the radiation which is of interest in an application and the accuracy with which the material parameters are known. For enclosed media which interact only weakly with the radiation the setting of diffuse, gray surfaces enclosing a non-participating medium is a sufficiently good approximation. Provided one accepts this all-Lambertian setting the discussion of the influence of the imprecise knowledge of the radiative properties and modeling or discretization errors can be based on solid ground. For discretization errors one can use the standard convergence theory as exemplified by Voigt et al. [44]. Furthermore, Strang's Lemma can be used for an a priori analysis of the influence of errors on the view factors [36]. Modeling errors, in particular errors in the reflection coefficients or omitting the interaction with participating media, can largely be answered by goal-oriented error estimation. For a nice review see the seminal work Becker and Rannacher [5]. Their influence on measurable quantities which are given as functionals of the radiosity can be monitored by the techniques developed for quantifying the influence of stabilizing terms for convection dominated flows [8]. There were similar developments in computer graphics [1,25]. Especially the error analysis for radiosity by Lischinski et al. follows the same lines as goal-oriented error estimation but in a discrete setting and specialized for global illumination. Radiative heat transfer and global illumination, i.e. indirect lighting, are described by the same equation [18]. Hence, the techniques developed in computer graphics (CG) to achieve realtime, photo-realistic image synthesis [10] should be beneficial for building performance simulations (BPS) as well. Yet, despite all mathematical similarities the notion of “accuracy” is different. In the RHT part of BPS the radiant flux between surfaces has to be resolved with great accuracy but not the discontinuities of the visibility function itself, i.e. the local variation of the radiant flux density. In CG we rather face a shape-matching problem: photo realism requires an accurate localization of shadows, especially soft ones. Jagged boundaries between lit and unlit areas are more perceptible than a shadow's shade of black. In CG A2A view factors often are calculated only approximately by using the P2A view factors and a suitable adaptive mesh refinement strategy [20]. Using the graphics processing units (GPU) of dedicated graphics cards [12,28,29] for the view factor computations led to a dramatic improvement of execution times. Although the visibility test

proposed by Coombe et al. [12] was fast, it traded speed for accuracy making it less useful for RHT applications [42]. The main advantage of GPUs is that the fundamental CG algorithms are implemented in silico, i.e. in hardware. Especially depth testing, a precursor to any occlusion query, is done with full hardware acceleration. For validation purposes the problem of the A2A view factor between two arbitrary, unobstructed polygons has been solved in closed form [41] by converting the double area integral (2AI) into two nested contour integrals (2CI) using Stokes' theorem. Another, entirely CPU-based attempt based on the 2AI is by Dobrowolski de Carvalho Augusto et al. [9]. For geometries of particular importance special solutions exist, for instance in reactor design [7,15]. Combined with local mesh refinement to convert partial occlusion into an all-or-nothing property [16], Stokes' theorem can be made applicable to obstructed view factors as well. From a mathematical point of view the 2CI can only be applied for unobstructed pairs of surfaces while the 2AI is formally exact. The hemicube (HC) algorithm [11] provides a practical solution for resolving the visibility issue. Although the restrictions on the 2CI can be relaxed by local mesh refinement it yet remains only approximately correct. The main shortcoming of approaches to obstructed view factors not utilizing the GPU is the need to re-implement many of the algorithms for visibility testing which otherwise would be readily available in hardware-accelerated form. Hence, less computational resources can be spent on the primary aim of evaluating the double integral. The hemicube algorithm [11] in its OpenGL-based form [28] is the key to overcome the occlusion problem while increasing performance. OpenGL provides a set of tailor-made functions for rendering the lighting of a three-dimensional (3D) scene with full hardware acceleration. On nowadays GPUs the resolution of the HC can be chosen almost arbitrarily fine such that errors due to aliasing effects [31] can be neglected. An accurate way of computing obstructed A2A view factors is given by the combination of numerical integration with the HC method, i.e. locating the view ports of the HCs at the quadrature points. Attempts to control accuracy by adaptive quadrature are due to Gershbein [17] and Walton [45]. For the latter there is an opensource implementation available by Walton and Pye [46]. It is the aim of this paper to port some of the high-performance solutions of CG for view factors to the field of RHT and to combine them with task-parallel computing techniques to take advantage of modern trends towards non-uniform multi-core architectures. This enables us to predict thermal comfort with high accuracy at the level of detail comparable to the spatial resolution of the finite element problem for the indoor airflow and to accurately assess the true contribution of RHT to the overall thermal performance of heating devices. Similar to our approach OpenGL has been successfully employed for the accurate computation of direct solar radiation on building surfaces of complex shape at high speeds [23]. Compared to previous implementations [16,35] we gain a striking speedup of three orders of magnitude. The paper is organized as follows. In section 2 we briefly repeat the properties of view factors and recall how to compute P2A view factors from the OpenGL-based HC. Section 3 is devoted to the discussion of the details of the parallelization of the overall computation and the resulting software design. The influence of the method's parameters and the efficiency of the parallelization is discussed in section 4 using three real-life test cases. The conclusions are drawn in section 5.

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

2. Basics of radiative heat transfer

qðrÞ  εðrÞ

2.1. Flux balance

JðrÞ ¼ εðrÞEðrÞ þ 9ðrÞHðrÞ;

(1)

where E(r) ¼ sT4(r) is the blackbody emissive power, s is the StefaneBoltzmann radiation constant, H(r) is the irradiation and 9H is its reflected part. The total irradiation at r is given by

Z

cosðqr Þcosðqr0 Þ

HðrÞ ¼

pjr  r0 j2

~ AðrÞ

Jðr0 ÞdAðr0 Þ

(2)

0

that is, the radiosities from all other visible points r . The angles qr 0 and qr0 are those between the surface normals at r, r and the line 0 connecting r and r . The quantity

dFdA/dA0 ¼

~ AðrÞ

Z

The central quantity governing radiative heat transfer between gray, diffuse surfaces is radiosity J, the total radiation energy leaving a surface element per unit time and unit area. In the case of a gray emitter the local emission ε(r) and absorption coefficient a(r) at a point r on the surface1 A of an enclosure do not depend on the wave length of the thermal radiation and Lambert's cosine law has to hold which states that radiation is emitted in all directions with equal intensity. The local reflection coefficient is given by 9ðrÞ ¼ 1  εðrÞ. Therefore, the total radiant heat flux leaving a surface at position r is

cosðqr Þcosðqr0 Þ pjr  r0 j2

(3)

is the unobstructed, differential point-to-point (P2P) view factor ~ and AðrÞ is the part of the surface visible from r. The radiative heat fluxes have to be in balance with all other heat fluxes supplied to the surface, i.e.

qðrÞ ¼ JðrÞ  HðrÞ:

(4)

From this energy balance three different integral equations can be derived for the radiative heat transfer in an enclosure with a gray interior surface ([32]; Sec. 5.3).

Z JðrÞ  9ðrÞ

Jðr0 ÞdFdA/dA0 dAðr0 Þ ¼ εðrÞEðrÞ:

2.3. Discretization To numerically solve for either the spatial distribution of radiosity J(r), surface temperature T(r) or heat flux q(r) we have to discretize the integral equation(s) with boundary elements. We subdivide the surface A into N patches of finite size. At this point we explicitly include the option for curvilinear patches and do not assume that the surface patches were isothermal. Furthermore, each patch Ai is considered as being composed of one or several surface elements. The latter applies to the case when we integrate over patches which are non-simplicial, e.g. quadrilateral. We approximate the sought solutions by low-order polynomial boundary elements. Depending on the polynomials their support is either one patch (globally discontinuous elements) or a collection of adjacent patches like in standard conforming finite elements. The numerical solutions Jh(r), Th(r), qh(r) are considered as a linear combination of polynomial trial functions, e.g.

Jh ¼

X

Jj fj ðrÞ

(8)

j

with coefficient vectors J ¼ (J1,…,JNDoF) and T, q, respectively. The number of degrees of freedom NDoF may be different from the number of patches N, especially for higher order boundary elements, e.g. quadratic (discontinuous) Lagrange elements as in €der [40], or the Lischinski et al. [26], the wavelet methods by Schro general higher order Galerkin approaches by Zatz [48], Troutman and Max [43] and Pellegrini [34]. Inserting this ansatz, for instance into Eq. (6), multiplying with a test function ji with support Bi(r) and integrating gives

X

Z ji ðrÞfj ðrÞdBi ðrÞ

Jj

j



X

Bi

Z

Z ji ðrÞ

Jj

Z

(5)

(7)

~ AðrÞ

fj ðr0 ÞdFdA/dA0 dAðr0 ÞdBi ðrÞ

(9)

~ AðrÞ

Bi

ji ðrÞqðrÞdBi ðrÞ:

¼

If radiosity is to be computed from given temperatures we have to solve

ðEðrÞ  Eðr0 ÞÞdFdA/dA0 dAðr0 Þ:

¼

j

2.2. Equations to solve

347

 1  1 qðr0 ÞdFdA/dA0 dAðr0 Þ εðr0 Þ

Z 

Bi

From this Galerkin approach the numerically less expensive collocation methods can be recovered by a proper choice of the test functions or quadrature rules for the outer integral over Bi.

~ AðrÞ

If the heat fluxes at the surface are known but the temperatures are not we get

Z JðrÞ 

Jðr0 ÞdFdA/dA0 dAðr0 Þ ¼ qðrÞ:

2.4. View factors The most expensive task is the computation of the double area integral which represents a generalized obstructed A2A view factor

(6)

~ AðrÞ

Finally, eliminating radiosity gives an integral equation which establishes a direct connection between heat fluxes and temperatures

1 If the orientation of surfaces (indicated by their normal vectors) matters we use bold-face, capital letters. Regular capital letters only denote the area of a surface.

Z

Z

Bi

~ AðrÞ

Fij :¼

ji ðrÞdFdA/dA0 fj ðr0 ÞdAðr0 ÞdBi ðrÞ

(10)

between the supports of the test function ji and trial function fj. Using piece-wise constant test and trial functions j, f which have only one patch as support the standard concept of view factors between isothermal patches is recovered. In this case the discrete version of Eq. (6) is

348



S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354



dij  Fij Jj ¼ qj ;

(11)

where we use Einstein's summation convention, dij is the Kronecker symbol and qj is the average heat flux on patch j. In this setting there is one degree of freedom per patch which is usually located in the center of mass of the patch. Standard mesh refinement methods can then be used to restore the validity of the assumption of isothermal surface elements. For each pair of surfaces Ai and Aj of finite size the obstructed A2A view factor is

FA /A~ :¼ i

j

1 Ai

Z dFdA /A~ dAi : i

(12)

j

Ai

~ . To compute F The visible part of Aj is denoted as A ~ j by j Ai /A numerical quadrature we have to efficiently compute the P2A view factors

Z dFdA /A~ ¼ i

j

~ j ðxi Þ A

    cos qi xi ; xj cos qj xi ; xj ~ dAj : 2  pxi  xj 

(13)

~ , respectively. For The points xi and xj are positions on Ai and A j the choice of coordinates and the definition of the angles qi, qj cf. Fig. 1. 2.4.1. Evaluation of the inner integral To evaluate the integral in Eq. (13) we use the HC method [11]. Let p be an arbitrary point on dAi and let A⊥ j be the area of the ~ onto the HC, cf. Fig. 2. Furthermore, let DFxy be a projection of A j P2Px view factor on the top of the HC, DFyz a P2Px view factor on one of the faces parallel to the yz-plane and DFxz a P2Px view factor on one of the faces parallel to the xz-plane. The P2Px view factors are computed from a closed expression [29] in a preprocessing step. Then, the obstructed P2A view factor of Aj w.r.t. p is given by the sum of the P2Px view factors of the individual pixels covered by A⊥ j H dFdA ¼ i /Aj

X ðx;yÞ2A⊥j

DFxy þ

X ðy;zÞ2A⊥j

DFyz þ

X

Fig. 2. Choice of coordinate system for the computation of point-to-area view factors and projection of a distant surface Aj onto the hemicube. Summation of the point-topixel view factors belonging to the shaded pixels gives the point-to-area view factor of Aj. The projection of all distant surfaces Aj is done by rendering them with OpenGL on the graphics card. The camera is located at p and points into the z-direction.

DFxz :

(14)

ðx;zÞ2A⊥j

The expensive part in this procedure is the determination of which of the distant surfaces Aj is visible in a given pixel on the HC. However, this is the fundamental operation of rendering 3D scenes in CG and what GPUs can do with full hardware acceleration. The implementation of the HC requires only two OpenGL functions in order to set up the render passes which project the other surface elements on the surfaces of the HC. The camera position is set by

Fig. 1. Geometry of radiative heat transfer between a differential and a finite surface. The shadow B of the occluding object depends on the point of view xi.

where T is a placeholder for GLfloat or GLdouble. The first three function arguments eyex, eyey, eyez define the camera position p which is some point on Ai. The viewing direction is indicated by the coordinates centerx, centery, centerz of the point c which is defined in the global coordinate system of the whole scene (aka world coordinates). To fully specify the orientation of the camera one has to fix the rotation around the viewing direction. This is done by specifying the components upx, upy, upz of the vector u. The points c and u depend on the particular surface of the HC. To specify them we need a local coordinate system spanned by the vectors ex, ey, ez. We need as auxiliary quantities the vectors dj: ¼ vj þ 1mod3  vj, j2f0; 1; 2g which describe the edges of the triangle formed by the first three vertices v0, v1, v2 of Ai. Then, ex ¼ d1 =jd1 j and ez ¼ ex  d2 =jex  d2 j. The y-axis is given by ey ¼ ex  ez. To render the top of the HC we use u ¼ ex and c ¼ p þ ez. For the side walls we have to set u ¼ ez and c ¼ p ± ex and c ¼ p ± ey, respectively. The actual shape of the HC depends on the viewing volume which is set by

where fovy (¼90 for the HC) is the angle of the field of view in the y direction of the coordinate system of the camera (i.e. in the direction of u) and aspect (¼1 for the HC) is the aspect ratio of the viewport, i.e. the visible rectangular section of the near clipping plane. The parameters zNear and zFar indicate the distance of the near and far clipping plane from the center of the camera along the axis defined by c  p. The distance to the far clipping plane can be approximated by the diameter of the bounding sphere of the total scene and zNear has to be small enough such that the closest other surface element Aj is still inside the viewing volume. To identify the projected parts of the surfaces each surface element gets a unique color in RGBA format by converting its index el_idx by bit-shift operations from a 32 bit integer into 8 bit blocks which then are used as red, green, blue and a components.

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

This allows for 232 z 4,109 individual surface elements in a scene. Adding all P2Px view factors of pixels of the same color gives the obstructed P2A view factor. A sample HC rendering is shown in the FBO in Fig. 3 below. 2.4.2. Evaluation of the outer integral The final numerical approximation of FA /A~ is obtained by i j numerical quadrature over Ai, i.e. we compute a weighted sum of P2A view factors

FAh /A~ i j

¼

X a

  H wa dFdA ~ x i;a : /A i

j

(15)

The quadrature weights ua contain the Jacobian of the transformation of Ai to the reference element. For triangular surfaces this cancels the factor 1/Ai. The superscript h denotes the diameter of the circle enclosing Ai we integrate over and is a measure for its size. In the presence of occluding objects visibility between two surfaces becomes highly discontinuous. Hence we use low-order quadrature rules and improve accuracy rather by multiple regular subdivisions of Ai. Each subdivision of a triangle introduces four new triangles by connecting the midpoints of the old triangle's edges. On each subtriangle T3Ai we apply the quadrature rule and sum over all subsurfaces

FAh /A~ ¼ i

j

X X T3Ai

a



ðTÞ H wa dFdA ~ xi;a : /A i

j

(16) ðTÞ

Within a subtriangle T a quadrature point is denoted as xi;a . As quadrature rule we typically use the 3-point Gaub-Legendre quadrature for planar triangles [19]. The number of subdivision steps is denoted as ns, yielding 4ns triangles. Quadrilateral elements are rendered as such but get split into two triangles for integration over Ai.

349

2.5. Quality control of view factors In contrast to previous approaches by Walton [45], Dobrowolski de Carvalho Augusto et al. [9], Francisco et al. [16], Bopche and Sridharan [7], Feng and Han [15], which only monitor the error in the energy transport between individual pairs of surfaces, we rather opt for a more complex error control by looking at the values of the weighted column sums of the A2A view factor matrix

X Ai i

Aj

!

FAi /Aj ¼ 1

cj:

(17)

which can be derived from the reciprocity relation

Ai FAi /Aj ¼ Aj FAj /Ai

(18)

The advantage of monitoring the column sums, Eq. (17), is that it is applicable in those cases where no closed-form solution of the 2AI exists. For geometries as they arise in practice this is virtually always the case. Furthermore, the columns of the view factor matrix describe how a given surface (the column index) contributes to the heat fluxes at all other visible surfaces (shoots energy) whereas the rows represent how a given surface (the row index) gathers energy from the other surfaces. By construction the row sums of HC-based view factor matrices always equate to unity. If this is not the case, one or several surface elements must have had the wrong orientation, i.e. their surface normals pointed into the wrong direction. The column sums are the only solution-independent true measure for radiant flux conservation. 3. Parallelization strategy The main observation concerning parallelization is that the calculation of an individual row of the A2A view factor matrix is

Fig. 3. Pipelining and parallelization of the hemicube rendering and the summation of the P2Px view factors. Details see text.

350

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

completely independent of all other rows. Furthermore, the individual P2A view factors needed to compute the A2A view factors by numerical quadrature are also independent of each other. Hence, the overall problem of computing the A2A view factor matrix is highly parallel on several levels. In addition to the approach shown here, it could easily be distributed over several compute nodes. The rendering of the HC on the GPU and the summation of the P2Px view factors on the CPU is parallelized according to the producer-consumer model [47]. The implementation is based on the multithreading module and the QSemaphore example of the QT library. The GPU produces rendered hemicubes, i.e. large arrays of P2Px view factors and the CPU consumes them to form the P2A view factors which ultimately get condensed into the A2A view factors of one row. In stage (0) the hemicube is rendered into a framebuffer object (FBO). Stage (1) of the pipeline consists of copying the content of another FBO into a pixel buffer object (PBO). Stage (2) takes care of transferring data from another PBO from the GPU to the CPU-side of the PCI-Express bus by means of a direct-memory access transfer. Parallel to these buffer operations the contents of a CPU-side copy of a PBO are summed up to get the P2A view factors which is stage (3). All stages can be executed in parallel because each stage of the pipeline has its own set of FBO, PBO and CPU-side PBO-mirror. The tiles in the latter indicate how different parts of the image of the hemicube are distributed among several threads for computing the P2A from the P2Px view factors. A graphical representation of the algorithm is given in Fig. 3. The FBOs display a sample image of a rendered hemicube.

4. Results There are two basic ways to warm up an indoor space, either by convection or by radiation. From the point of view of thermal comfort the advantages of radiation are that heat can be supplied to the receiver without any air movement and without any remarkable delay. This is often made use of by modern heating devices. In practice both phenomena, convection and radiation, occur together. Therefore it is necessary to accurately predict the individual contribution of radiation and convection to the heat supply. Moreover in most cases only the exact prediction of both phenomena, radiation and convection, leads to useful simulation results.

Fig. 5. Details of the finite element surface mesh for the flow field around a radiator.

4.1. Test cabin compliant with EN 442 As a first example the performance of a well known and widely spread radiator is tested in detail. The selected tubular radiator has two or three rows, twelve segments and a height of 0.6 m. Its thermal output is assessed according to DIN EN 442, part 2 [13] in a virtual test cabin based on a coupled building- and flow simulation, see Fig. 4. The test cabin itself has a square base area with an edge length of 4 m and a height of 3 m. All walls of the cabin are cooled, except the wall directly behind the radiator. For more details related to the test set up see Ref. [13]. The calculation procedure includes the fluid flow and the three dimensional heat conduction within the metal parts of the radiator, the indoor air flow and heat conduction of the test cabin and the detailed RHT between the radiator and the test cabin by the procedure given in the previous sections. Therefore both the test cabin, including the outer radiator surfaces and the radiator interior zone are geometrically modeled, meshed and simulated using different bidirectional coupled simulation tools. Figs. 5 and 6 show some details of the different grid sizes used for the involved simulation tools. To ensure the correct radiant heat balance all geometrical details of the radiator have to be modeled, which leads to a mesh for the RHT problem of a similar resolution as for the flow problem, cf. Fig. 6.

Fig. 4. Surface temperature distribution and flow field around a radiator with three rows and twelve segments in a test cabin according to EN 442, cf. last row in Table 1.

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

351

investigated so far agree well with the manufacturer's data, see Table 1. 4.2. Warehouse

Fig. 6. Details of the surface mesh for radiative heat transfer from a radiator.

For more details concerning the simulation tools and the coupling algorithm see Ref. [27]. During the simulation the energy balance between the radiator and the test cabin walls is calculated. Under well-defined boundary conditions the thermal output is generated and documented. In addition to the data given by manufacturers (which are used as a reference) we can specify the part of heat transfer due to radiation. Because of the accuracy and flexibility the simulation tool is applicable to any type of radiator and any type of heating/cooling problem. The results presented here and for all radiator types

Large indoor spaces need special care in heating to use energy in an efficient way. Here radiation also plays an important role. One way is to work with ceramic infrared heaters. They directly burn natural gas in a special ceramic device such that the heating areas can be kept small because of the high operating temperatures. This allows for very compact designs. The prediction of thermal comfort in a large space heated by ceramic infrared radiators is only possible by a precise radiation calculation. Fig. 7 shows a sketch of a warehouse heated by radiation based on ceramic infrared heaters. This example shows a typical case, where small faces are crucial for the quality of the results. The extreme differences in size of the heating devices (with surface areas of 0,25 m2) and the whole warehouse (base area: 1500 m2) in combination with the high wall temperatures of the heater require a very accurate calculation of the view factors. The part of RHT for those infrared heaters lies between 70 and 80 percent of the overall heat supply. Hence, small inaccuracies in the corresponding view factors would lead to significant errors of the radiation balance. 4.3. Car cabin Another example where radiation has to be calculated in detail is the simulation of the indoor environment and the thermal comfort in a car cabin, see Fig. 8. Fig. 8 shows the inner surface

Fig. 7. Surface mesh of the warehouse. Some of the elements have been removed for view on the inside. The details of the geometry of the ceramic infrared heaters are shown on the lower left.

352

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

Fig. 8. Inner surface temperatures of the simulated cabin.

temperatures of the modeled car cabin. In this case the relation between the occupied space and the control volume is close to one. To calculate the influence of the very differing surface temperatures onto the occupants body parts it is useful to model the occupant in detail, see Fig. 9. In such cabin simulations the thermal comfort always has to be evaluated in detail because of considerable differences between head and feet of the driver. The contribution of RHT is very sensitive to changes of the current situation inside and outside the car. The situation shown in Fig. 10 is a medium winter day with temperatures around 0  C. The air heating has to compensate the heat losses in order to ensure a defined operative temperature at the sensor point, which is in the instrument board close to the driver. Fig. 11 displays some stream traces of the calculated flow field inside the car cabin, driven by the air inlets. The streamtraces are colored according to the air temperature. In addition Fig. 12 shows some vertical contour plots of the operative temperature which is used for evaluating the local level of thermal comfort. This operative temperature is a combination of the field values of the air temperature and the radiation temperature. 4.4. Benchmarking against Dudka's RRV The performance is tested on a Macbook Pro Retina powered by a Core i7-3820QM and an NVidia GeForce GT 650 m running under Mac OSX 10.9.5. As comparison we used the RRV package by Kamil Dudka [14]. It is a serial implementation of the radiosity

Fig. 9. Different surface meshes of the head of the manikin for CFD simulation and radiative heat transfer.

Fig. 10. Inner surface temperatures and flow field with inlet devices of the simulated car cabin.

method and computes the view factors using the HC method implemented in OpenGL. We chose this particular implementation because it is easy to install and comes with a set of reasonable test scenes. The results of our tests are shown in Table 2. The first column lists the numbers of pixels per dimension Npx for the HC used, the second gives the CPU wall time t1000 per 1000 triangles of RRV, the third the bandwidth utilization of the PCI express (PCIe) bus by our implementation, the fourth column contains the CPU wall time t1000 per 1000 triangles of our implementation. The speedup due to our parallelization is given in the last column. We ran our implementation with a 1-point quadrature rule for the outer integral, such that both programs render one HC per surface element. As test geometry we used the unrefined room2.xml with 3936 triangles from RRV's collection of sample scenes. For unknown reasons RRV was particularly slow for Npx ¼ 1024. Using both implementations as given, ours is 6e23 times faster. How2 pixels per HC, ever, while our implementation needs only 3Npx 2 RRV needs 9Npx putting more pressure on the PCIe bus and wasting two thirds of the pixels copied from the GPU to CPU. Both implementations use glReadPixels to copy pixel data. This function is known to limit the maximum transfer rate and therefore our 1.65 GB/s can be considered as exhausting the PCIe bus, although the theoretical transfer for PCIe 2.0 should be 4 GB/s. However, this value is for a GPU which does not do any rendering in between.

Fig. 11. Stream traces of the simulated flow field in the cabin.

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

353

Fig. 12. Operative temperature in the car cabin e combination of air temperature and radiative temperature in vertical planes.

4.5. Accuracy of numerical radiative heat transfer The geometries in practical cases are often very complex, such that Eq. (17) cannot be fulfilled exactly. In such situations it is necessary to check where the elements are located, which generate a deviation in column sum. The influence of the dimension Npx of the hemicube and the refinement level ns of the quadrature rule is shown in Fig. 13. The final criterion is to fulfill the energy balance, which is demonstrated by the warehouse with infrared heaters. The dependence of the quality of the energy balance on the method parameters is given in Fig. 14. 5. Conclusion In order to predict or analyze radiative heat transfer in indoor environments in an efficient and detailed way it is necessary to apply modern computer methods and technologies. It is shown that view factors can be calculated by use of an OpenGL-based procedure. One promising path, from basics of radiative heat transfer up to computational details, is given in the paper. Compared to previous implementations accuracy is improved and

Table 1 Comparison of simulated and measured heat supply (manufacturer's data) of different radiator types for temperature levels 75/65/20 and 55/45/20. All temperatures are in +C. Tubular radiator type, temperatures

Heat supply calculated in W total (conv./rad.)

Heat supply by manufacturers total in W

2 2 3 3

284 542 385 735

280 528 375 720

rows, rows, rows, rows,

55/45/20 75/65/20 55/45/20 75/65/20

(174/110) (333/209) (257/128) (496/239)

to to to to

Fig. 13. Deviation from unity of the weighted column sums in case of the warehouse, cf. Fig. 7. We plot the extrema as function of the hemicube dimension Npx and the level of quadrature subdivision ns.

runtime is considerably reduced. This opens the door for detailed investigations of radiative heat transfer in any geometrical configuration. The advantages of the developed computational method is demonstrated by three different cases: a virtual test room for a tubular radiator, a heating of a ware house by ceramic infrared heaters and the heat transfer in a car cabin. In all three cases it is crucial to get very accurate view factors for the calculation of radiative heat transfer. This is discussed by comparing the simulated heat output of a radiator with the manufacturer's measurements, by the energy balance concerning radiation in a ware house and by operative temperatures in a car cabin.

297 564 407 781

Table 2 Comparison of the overall performance of our implementation and of the publicly available RRV package [14]. Npx: number of pixels on the hemicube per dimension, t1000: CPU wall time per 1000 triangles, PCIe: effectively used bandwidth. Npx

256 512 1024 2048

RRV

Our implementation

t1000/s

PCIe/GB/s

t1000/s

Speedup

5.84 21.3 136 215

0.61 1.35 1.65 1.65

0.96 1.75 5.73 22.9

6.08 12.2 23.7 9.39

Fig. 14. Relative error of the radiant flux balance for the warehouse, cf. Fig. 7, for the 3point Gaub rule.

354

S. Kramer et al. / International Journal of Thermal Sciences 96 (2015) 345e354

GPU-computing is a quickly developing research field which offers new challenges for the solution of scientific problems. Radiative heat transfer is only one area of interest, but because of the similarity to computer graphics it offers a quick access to modern computational technologies. References [1] J. Arvo, K. Torrance, B. Smits, A framework for the analysis of error in global illumination algorithms, in: Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1994, pp. 75e84, http://dx.doi.org/10.1145/192161.192179. [2] M. Bebendorf, Approximation of boundary element matrices, Numer. Math. 86 (2000) 565e589, http://dx.doi.org/10.1007/PL00005410. [3] M. Bebendorf, Hierarchical lu decomposition-based preconditioners for bem, Computing 74 (2005) 225e247, http://dx.doi.org/10.1007/s00607-004-00996. [4] M. Bebendorf, S. Rjasanow, Matrix compression for the radiation heat transfer in exhaust pipes, in: A.M. S€ andig, W. Schiehlen, W. Wendland (Eds.), Multifield Problems, Springer Berlin Heidelberg, 2000, pp. 183e192, http:// dx.doi.org/10.1007/978-3-662-04015-7_20. [5] R. Becker, R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer. 10 (2001) 1e102. URL: http://journals.cambridge.org/article_S0962492901000010. €oder, Sparse matrix solvers on the gpu: [6] J. Bolz, I. Farmer, E. Grinspun, P. Schro conjugate gradients and multigrid, ACM Trans. Graph 22 (2003) 917e924, http://dx.doi.org/10.1145/882262.882364. [7] S.B. Bopche, A. Sridharan, Determination of view factors by contour integral technique, Ann. Nucl. Energy 36 (2009) 1681e1688, http://dx.doi.org/ 10.1016/j.anucene.2009.09.007. URL: http://www.sciencedirect.com/science/ article/pii/S0306454909003016. [8] M. Braack, A. Ern, A posteriori control of modeling errors and discretization errors, Multiscale Model. Simul. 1 (2003) 221e238, http://dx.doi.org/10.1137/ S1540345902410482. [9] L. Dobrowolski de Carvalho Augusto, B. Giacomet, N. Mendes, Numerical method for calculating view factor between two surfaces, in: 10th Conference of the International Building Performance Simulation, Taylor and Francis, 2007. [10] M. Cohen, J. Wallace, Radiosity and Realistic Image Synthesis. Morgan Kaufmann Series in Computer Graphics and Geometric Modeling, Academic Press Professional, 1993. URL: http://books.google.de/books?id¼7JiYl9m3Y6YC. [11] M.F. Cohen, D.P. Greenberg, The hemi-cube: a radiosity approach for complex environments, Comput. Graph. (1985) 3e40. [12] G. Coombe, M.J. Harris, A. Lastra, Radiosity on graphics hardware, in: Proceedings of Graphics Interface 2004, Canadian Human-Computer Communications Society, School of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, 2004, pp. 161e168. URL: http://dl.acm.org/ citation.cfm?id¼1006058.1006078. [13] DIN EN 442e2, Radiators and Convectors e Part 2: Test Methods and Rating, German version ed., Beuth Verlag, 2003 [14] K. Dudka, RRV: Radiosity Renderer and Visualizer, 2007. URL: http://dudka.cz/ rrv. [15] Y. Feng, K. Han, An accurate evaluation of geometric view factors for modelling radiative heat transfer in randomly packed beds of equally sized spheres, Int. J. Heat Mass Transf. 55 (2012) 6374e6383, http://dx.doi.org/ 10.1016/j.ijheatmasstransfer.2012.06.025. URL: http://www.sciencedirect. com/science/article/pii/S001793101200453X. [16] S.C. Francisco, A.M. Raimundo, A.R. Gaspar, A.V.M. Oliveira, D.A. Quintela, Calculation of view factors for complex geometries using Stokes' theorem, J. Build. Perform. Simul. 6 (2013) 1e14. URL: http://www.tandfonline.com/ doi/abs/10.1080/19401493.2013.808266. [17] R. Gershbein, An Adaptive Gauss Method for Computing Irradiance Coefficients of Galerkin Radiosity Systems, Technical Report, Department of Computer Science, Princeton University, 1995. [18] C.M. Goral, K.E. Torrance, D.P. Greenberg, B. Battaile, Modeling the interaction of light between diffuse surfaces, in: Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1984, pp. 213e222, http://dx.doi.org/10.1145/800031.808601. [19] P.C. Hammer, O.J. Marlowe, A.H. Stroud, Numerical integration over simplexes and cones, Math. Tables Other Aids Comput. 10 (1956) 130e137. [20] P. Hanrahan, D. Salzman, L. Aupperle, A rapid hierarchical radiosity algorithm, SIGGRAPH Comput. Graph. 25 (1991) 197e206, http://dx.doi.org/10.1145/ 127719.122740. [21] J. Howell, A Catalog of Radiation Configuration Factors, McGraw-Hill Ryerson, Limited, 1982. http://books.google.de/books?id¼uzSdMAEACAAJ. [22] J.R. Howell, A Catalog of Radiation Heat Transfer Configuration Factors, 2010. URL: http://www.thermalradiation.net/indexCat.html. [23] N.L. Jones, D.P. Greenberg, K.B. Pratt, Fast computer graphics techniques for calculating direct solar radiation on complex building surfaces, J. Build. Perform. Simul. 5 (2012) 300e312, http://dx.doi.org/10.1080/ 19401493.2011.582154.

[24] S. Kurz, O. Rain, S. Rjasanow, The adaptive cross-approximation technique for the 3d boundary-element method, Magnetics IEEE Trans. 38 (2002) 421e424, http://dx.doi.org/10.1109/20.996112. [25] D. Lischinski, B. Smits, D.P. Greenberg, Bounds and error estimates for radiosity, in: Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1994, pp. 67e74, http://dx.doi.org/10.1145/192161.192176. [26] D. Lischinski, F. Tampieri, D.P. Greenberg, Discontinuity meshing for accurate radiosity, IEEE Comput. Graph. Appl. 12 (1992) 25e39, http://dx.doi.org/ 10.1109/38.163622. €sler, Stabilized finite element [27] G. Lube, T. Knopp, G. Rapin, R. Gritzki, M. Ro methods to predict ventilation efficiency and thermal comfort in buildings, Int. J. Numer. Methods Fluids 57 (2008) 1269e1290, http://dx.doi.org/ 10.1002/fld.1790. [28] E. Lum, K.L. Ma, N. Max, Interactive radiosity using mipmapped texture hardware, in: S. Gibson, P. Debevec (Eds.), 13th Eurographics Workshop on Rendering: Pisa, Italy, June 26e28, 2002, Eurographics Association, 2002. [29] E.B. Lum, K.L. Ma, N. Max, Calculating hierarchical radiosity form factors using programmable graphics hardware, J. Graph. GPU Game Tools 10 (2005) 61e71, http://dx.doi.org/10.1080/2151237X.2005.10129210. [30] M. Maivel, J. Kurnitski, Low temperature radiator heating distribution and emission efficiency in residential buildings, Energy Build. 69 (2014) 224e236. URL: http://www.sciencedirect.com/science/article/pii/S0378778813006762. [31] N. Max, Optimal sampling for hemicubes, IEEE Trans. Vis. Comput. Graph. 1 (1995) 60e76, http://dx.doi.org/10.1109/2945.468388. [32] M.F. Modest, Radiative Heat Transfer, second ed., Academic Press, 2003. [33] P. Nielsen, F. Allard, H. Awbi, L. Davidson, A. Sch€ alin, Computational Fluid Dynamics and Ventilation Design, REHVA Guidebook 10, REHVA, 2007. URL: http://www.rehva.eu/publications-and-resources/guidebooks-shop/ guidebook/?tt_products[backPID]¼6&tt_products[product]¼ 37&cHash¼5cbbabb7a0e196769c71b38e730a39bb. [34] M. Pellegrini, A geometric approach to computing higher order form factors, in: Proceedings of the Fifteenth Annual Symposium on Computational Geometry, ACM, New York, NY, USA, 1999, pp. 69e78, http://dx.doi.org/10.1145/ 304893.304914. [35] A. Perschk, Geb€ aude-Anlagen-Simulation unter Berücksichtigung der hygri€udewa €nden, Ph.D. thesis, Technische Universit€ schen Prozesse in den Geba at Dresden, 2000. URL: http://nbn-resolving.de/urn:nbn:de:swb:14993456440484-94152. [36] N.A. Qatanani, A. Daraghmeh, Asymptotic error analysis for the heat radiation boundary integral equation, Eur. J. Math. Sci. 2 (2013) 51e61. [37] S. Rjasanow, Adaptive cross approximation of dense matrices, in: Int. Association Boundary Element Methods Conf., IABEM, 2002, pp. 28e30. [38] S. Rjasanow, O. Steinbach, The Fast Solution of Boundary Integral Equations. Mathematical and Analytical Techniques with Applications to Engineering, Springer, 2007. URL: http://books.google.de/books?id¼RX5CAAAAQBAJ. [39] C. Scheiblich, V. Kolitsas, W. Rucker, Compression of the radiative heat transfer bem matrix of an inductive heating system using a block-oriented wavelet transform, Magnetics IEEE Trans. 45 (2009) 1712e1715, http:// dx.doi.org/10.1109/TMAG.2009.2012793. € der, Wavelet radiosity: wavelet methods for integral equations, in: [40] P. Schro ACM SIGGRAPH '96 Course Notes e Wavelets in Computer Graphics, 1996, pp. 143e165. € der, P. Hanrahan, On the form factor between two polygons, in: [41] P. Schro Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1993, pp. 163e164, http:// dx.doi.org/10.1145/166117.166138. [42] H. Takizawa, N. Yamada, S. Sakai, H. Kobayashi, Radiative heat transfer simulation using programmable graphics hardware, in: Proceedings of the 5th IEEE/ACIS International Conference on Computer and Information Science and 1st IEEE/ACIS International Workshop on Component-based Software Engineering,Software Architecture and Reuse, IEEE Computer Society, Washington, DC, USA, 2006, pp. 29e37, http://dx.doi.org/10.1109/ICISCOMSAR.2006.70. [43] R. Troutman, N.L. Max, Radiosity algorithms using higher order finite element methods, in: Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1993, pp. 209e212, http://dx.doi.org/10.1145/166117.166144. [44] A. Voigt, N. Hanssen, C. Weichmann, The radiosity equation for solving global heat transfer in industrial furnaces, Math. Comput. Model. 39 (2004) 145e150. URL: http://www.sciencedirect.com/science/article/pii/ S0895717704900039. [45] G.N. Walton, Calculation of Obstructed View Factors by Adaptive Integration, Technical Report NISTIR-6925, National Institute of Standards and Technology, Gaithersburg, MD, 2002. [46] G.N. Walton, J. Pye, View3d, 2011. URL: view3d.sourceforge.net. [47] G. Wilson, P. Lu, Parallel Programming Using Cþþ. Scientific and Engineering Computation, MIT Press, 1996. URL: http://books.google.de/books? id¼n3kA3FzBTnEC. [48] H.R. Zatz, Galerkin radiosity: a higher order solution method for global illumination, in: Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, ACM, New York, NY, USA, 1993, pp. 213e220, http://dx.doi.org/10.1145/166117.166145.