Improved discrete ordinate method for accurate simulation radiation transport using solar and LED light sources

Accepted Manuscript Improved Discrete Ordinate Method for accurate simulation radiation transport using solar and LED light sources José Moreno, Cintia Casado, Javier Marugán PII: DOI: Reference:

S0009-2509(19)30404-X https://doi.org/10.1016/j.ces.2019.04.034 CES 14945

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

4 February 2019 29 March 2019 19 April 2019

Please cite this article as: J. Moreno, C. Casado, J. Marugán, Improved Discrete Ordinate Method for accurate simulation radiation transport using solar and LED light sources, Chemical Engineering Science (2019), doi: https:// doi.org/10.1016/j.ces.2019.04.034

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Improved Discrete Ordinate Method for accurate simulation radiation transport using solar and LED light sources

José Moreno, Cintia Casado, Javier Marugán* Department of Chemical and Environmental Technology, ESCET, Universidad Rey Juan Carlos, C/Tulipán s/n, 28933 Móstoles (Madrid), Spain * Corresponding author: Tel. +34 91 664 7466; E-mail: [email protected].

1

Abstract This work describes adaptive quadrature, a new feature designed to improve the Discrete Ordinate Method (DOM) in parallel and cone-shaped radiation sources. OpenFOAM was chosen as the simulation framework to implement the new feature. It already included a Discrete Ordinate Method (fvDOM), but it focuses on thermal radiation and is limited to isotropic emission, diffuse phenomena, and non-scattering media. The model was completed to cover volumetric and superficial absorption, isotropic and anisotropic scattering (with userdefined phase functions), diffuse and specular reflection, diffuse and parallel transmission, and three types of superficial emission sources, i.e., isotropic, cone-shaped, and parallel. Additionally, the model is prepared to work in grey mode or with wavelength bands. Multiple regions with different optical properties are also allowed. Most of the model has been validated by individual simulations of every feature in simple geometries that permit an analytical solution, with errors between 0% and 6.08%. These simulations were also verified in a comparison with an established version of the standard DOM, with differences between model implementations below 2.5%. Some advantages of the developed adaptive quadrature are also analysed using simulation results for different radiative sources and angular discretization. The main conclusion is that adaptive quadrature better defines view angle and light direction of emission sources compared to the established DOM, improving significantly the accuracy of the simulation of non-isotropic sources.

Keywords: OpenFOAM; radiation; photoactivated; quadrature; simulation; CFD. 1. Introduction Numerous methods exist to model radiation, most of which constitute simplifications that are unable to simulate a particular group of phenomena (e.g., radiosity methods are limited to

2

diffuse light and unable to model scattering, while the harmonic spheres method is almost entirely limited to volumetric phenomena). In this field, simpler is faster, and often, simulations can achieve a reliable level of accuracy in a short time if the correct simplification is chosen. However, in some applications, complete methods must be chosen, incurring a much higher computational cost and time. The Discrete Ordinate Method (DOM) and Photon Monte Carlo (PMC) are the most commonly utilized methods when no simplifications can be achieved. Both methods are able to simulate the behaviour of light on surfaces (absorption, emission, reflection, and transmission) and media in between (absorption, emission, transmission, and scattering). Moreover, PMC is more accurate in specular phenomena, while DOM better models diffuse phenomena. DOM is a commonly implemented method in multiphysics simulation frameworks (e.g., ANSYS Fluent, OpenFOAM, COMSOL, etc.), while PMC is usually founded as an independent software because its task queue differs from the iteration process of those frameworks. Our work is focused on photoactivated processes, a field in which most simulations require complex methods, such as DOM or PMC, because the variety of occurring phenomena hinders simpler methods. DOM is the most commonly employed method to solve radiation in photoactivated processes [1], [2]. The Discrete Ordinate Method (DOM) is based on discretization, not only of space and time, but also of wavelength and solid angles. It solves grey phenomena (those in which no changes in radiation wavelength occur) in discrete solid angles representing their central direction (discrete ordinates). To solve media with wavelength-varying properties, the wavelength continuum is divided in discrete wavelength bands, each of which is solved as grey, with properties based on the central wavelength.

3

Although fully functional for applications in other fields, the previous DOM model in the OpenFOAM repository lacks some important features required to solve photoactivated processes, including: a) Boundary conditions are: i) transparent diffuse wall, which can neither be semitransparent nor exhibit parallel transmission; or ii) opaque grey diffuse wall, including an emissivity coefficient, which splits incident radiation into an absorbed and diffuse reflected one, but cannot be specular reflected. This aspect hinders the simulation of reflectors commonly used in photoactivated processes. b) Surface emission is only diffuse- and temperature-dependent. Perfectly conceived for thermal radiation emission, it does not allow simulation of light sources commonly employed in photoactivated processes (e.g., mercury lamps, LED, sunlight, etc.). c) There is no scattering in the radiative transfer equation. Consequently, it is not possible to simulate the behaviour of slurry systems, in which the volume between surfaces is a participating media comprised of a suspension of particles. For these reasons, the aim of this work is the development of a new radiation model integrated in the OpenFOAM core that is able to consider all of the possible phenomena present in photoactivated processes. fvDOM has been widely used in research works, mainly in the field of thermal radiation [3]– [6]. To the best of the authors’ knowledge, however, there is only one extant work [7] about a DOM-based model developed to simulate a photoactivated process, i.e., an algae reactor. In that work, the pixelation effect was included to adapt to non-structured meshes, and scattering was based on Mie’s theory. However, information about boundary conditions or light sources was not included, and the code is not available in OpenFOAM repositories. Consequently, this model was not available to be used as a basis for the current work.

4

The developed software is based on OpenFOAM class architecture, which defines special types of variables, such as fields and templates, to create models to be easily integrated with other OpenFOAM modules. All of the phenomena that can be present in a photoactivated process (and in other process related with radiation transfer) are included (see Figure 1), containing the volumetric absorption of the medium, anisotropic scattering, semi-transparent walls with absorption, parallel and diffuse transmission, and specular and diffuse reflection.

Figure 1. Radiative phenomena included in the enhanced model.

5

Individual tests were carried out for each phenomenon, being compared to equivalent simulations in the well-known CFD software ANSYS Fluent (as a reference of the established DOM), as well as with an analytical solution, when possible. Then, the adaptive quadrature effect was studied to prove its performance. Parallel and cone emissions were simulated in different combinations of direction, view angle, and discretization, in a mesh imported from ANSYS Fluent to avoid differences attributable to spatial discretization. The enhanced model may be utilized in any photoactivated process, such as photolytic or photocatalytic processes. It also allows design, scale, and optimization of photoreactors containing any kind of light source, the presence of reflector panels, semi-transparent walls, or multiple regions with different optical properties.

2. Numerical Model 2.1. Implementation of the Discrete Ordinate Method (DOM) The radiative transfer equation (RTE) (Eq. 1) describes the conservation of radiative intensity in a direction of space. Solving the RTE involves solving all directions, which can be achieved by a previous integration through simplifications (Rosseland [8], model P-n [9], [10], view factors [11], method of moments [12], [13]) or by a subsequent numerical integration, as does the DOM or Monte Carlo models [14], [15]. It consists of separately solving a finite number of directions that are representative of all nearby addresses (discrete ordinates) distributed following a distribution map called a quadrature (see Figure 2). Next, the integral is calculated as the sum of all of the discrete ordinates.

6

(Eq. 1)

where

is the intensity of photons with wavelength , travelling in direction

volumetric absorption coefficient;

is the volumetric scattering coefficient;

Boltzmann constant; T is the temperature; and

is the is the Stefan-

is the phase function, which

describes the distribution of scattered radiation.

Figure 2. Graphical representation of a 15 x 15 quadrature and its elements.

The quadrature scheme used for directional discretization exerts a major influence on the precision and shape of the radiation calculations. Specifically, a low discretization (few divisions in θand φ angles) leads to undesired transformation of ordinates from cones to rays, resulting in a loss of homogeneity in the emission. Moreover, the difference of solid angle modulus (size) between ordinates causes preferential directions, unintentionally curving the emission. Consequently, the choice of a suitable directional discretization could be even more

7

crucial than the choice of the mesh. A correct quadrature orientation, for the same degree of discretization, increases precision, avoiding the saw tooth phenomena (cuts on the outer surface emitted following the divisions between ordinates [16], [17]) and/or mistakes in parallel and cone emission direction. The same quadrature scheme used in ANSYS Fluent [18] was implemented in the model, being the number of divisions in θ angle (traversing the sphere from the z-axis to the xyplane), and a number of divisions in φ angle (from the x-axis to the y-axis) determined by the user. Other studies [19] have used alternative quadrature, such as Gauss, which offers the advantage that all solid angles have the same modulus, but in turn greatly limits the user’s choice of discretization level. Since the chosen quadrature scheme does not generate ordinates of the same solid angle modulus, when transmission between directions exists (by reflection, scattering, and/or transmission), the solid angle module must be considered. To reduce the number of calculations, and since the quadrature is constant between iterations, a fit according to the solid angle is introduced into the scattering matrix (Eq. 2), which is solved only once prior to the beginning of resolution iterations. (Eq. 2) where

,

,

, and

are, as described above, radiant intensity, absorption coefficient,

scattering coefficient, and solid angle of the direction being solved, respectively; and refers to radiation scattered from other directions to the current one, but it is calculated and stored explicitly, using values of the last iteration. This equation is integrated for the volume of each cell, obtaining an algebraic version that varies according to the numerical schemes used. It is then resolved as a matrix linear system in the domain.

8

Once the RTE is solved for all discrete ordinates, the incident radiation in a given point is the sum of intensities in all directions:

(Eq. 3)

where

is the incident radiation in a given wavelength; and

is the discrete ordinate area.

Net radiant flux in boundaries is calculated as a summation of the input of every discrete ordinate. Its intensity normal to boundary surface is:

(Eq. 4) where

is the solid angle of a discrete ordinate (module and direction); and

is the normal

vector of the face where flux is being calculated. Incident radiation and fluxes will constitute the output of the simulation of every combination of quadrature and band. Flux may be incoming, emitted, or net (sum of both). As the third flux may be easily calculated from the other ones, only net and incoming fluxes are given as results. Net was chosen due to its possible use in heat transfer cases, and incoming flux is necessary for diffuse phenomena calculation in the code. Global model results, given as files to the user, are fluxes and incident radiation in each wavelength bandwidth, and the target in photoactivated simulations, i.e., absorbed energy, which is calculated in each band as:

(Eq. 5) where

is the absorption volumetric coefficient (m-1) for a specific wavelength. Besides the

volumetric equations defined previously, a general boundary condition defined has an equation for incident ordinates in the face (Eq. 6) and another one for leaving ordinates (Eq. 7):

9

(Eq. 6) (Eq. 7) where t is the transmission coefficient of the boundary, as r is its reflection coefficient;

is

again the intensity of photons with wavelength , travelling in direction ; and the remaining terms are fully described below. 2.1.1. Definitions at boundaries Optical properties of the material have been chosen to define the boundaries. The equations that define the emission sources are described below: •

To distribute radiation in an isotropic emission (Figure 1A), the following formula is used: (Eq. 8)

where •

is the emission flux declared by the user in boundaries.

To calculate intensity in a cone-shaped (Figure 1C) and parallel emission (Figure 1E): (Eq. 9) (Eq. 10) where

is the emission flux;

is the beam direction;

is the cone opening; and

is

the angle between the cone central direction and , i.e., the face normal vector. •

A power cosine law to better define LEDs, as an alternative to cone-shaped sources (Figure 1C), has been implemented:

(Eq. 11)

10

(Eq. 12)

where q is the emission flux;

is the beam direction;

model exponential factor; and

is the face normal vector.

is the cone opening; g is the

Domain walls can be either interior surfaces between regions with different optical properties or exterior walls. They can also transmit radiation or reflect it. Equations defining these properties are described below. To calculate transmitted radiation between regions, the following equations have been used: •

Transmission (Figure 1f) is composed of parallel and diffuse transmission: (Eq. 13) where

is the diffuse fraction in transmission;

in the same direction, but in a neighbouring region; and

is the radiant intensity is the incident

radiative flux coming from the neighbouring region. In the current version, parallel transmission is not refracted, and thus flux will be transmitted between the same ordinates in different regions. A similar equation was used in reflection (Figure 1B and 1D): •

It is composed of specular and diffuse reflection: (Eq. 14)

(Eq. 15)

where

is the reflection diffuse fraction;

incident radiative flux in the boundary;

is the face normal vector;

is the

is the direction estimated as specular

11

reflect; and

is the intensity in direction .

Figure 3. Analytical direction of specular reflection.

In specular reflection, the new direction must be calculated (Figure 3) and paired with an existing discrete ordinate. If the reflected direction is paired with a single discrete ordinate (named in the model as “first order reflection”), i.e., the nearest one, an error would occur in either the flux or the energy in the domain. The enhanced model chooses to preserve energy, as it is considered a more important variable in photoactivated processes than flux, but it can be easily changed to preserve flux or added as a third option. If it is paired with a couple of discrete ordinates (named as “second order reflection”), intensity is splitted between them to preserve both flux and energy, but there is a larger error in light direction. These errors decrease with angular discretization and are also dependent on quadrature orientation.

12

2.1.2. Volumetric definitions Regarding volumetric processes (Figure 1G and 1H), absorption and out-scattering are included in the RTE (Eq. 2) through the definition of absorption and scattering coefficients for a specific wavelength. Although it is not included in this code version, the model is ready to use non-homogeneous distributions of absorption and scattering coefficients as a future improvement of the model or for user development. This option would simulate multiphase simulations (e.g., bubbles, solid-fluid) or inhomogeneous solutions (e.g., reactive media). The inScatter term of the radiative transfer equation (Eq. 2) is calculated through the phase function

chosen by the user, which defines the distribution of scattering along the

spatial direction. (Eq. 16) where

is the ordinate solid angle; and

is its radiative intensity in each wavelength.

To calculate the phase function, in the case of isotropic scattering, the following is used: (Eq. 17)

where N is the number of discrete ordinates. Additional phase functions may be easily added to the model, inheriting the abstract mother class “phaseFuncModel”. As an example, the Henyey-Greenstein phase function was implemented. Its anisotropy coefficient defines the shape of scattering, in a range between -1 (full back-scattering) and 1 (full forward-scattering). (Eq. 18)

13

where

is the angle between discrete ordinates; and

is the anisotrpy coefficient for

wavelength . To ensure energy conservation, the phase function results of all given light paths are later normalized. In wide-band simulations, there are two additional options to define the phase function. The “Coeff” phase function allows the user to set a single-phase function model for all bands, but different coefficients in each band, while the “Band” phase function uses a different model in each band, declared as in grey simulations, with its own coefficients. The inScatter term is also calculated and stored in each wavelength band. 2.2. Enhancements in DOM The main change in the model is the treatment of quadrature. As previously described, quadrature is the scheme used to discretize spatial directions. It has a determinant importance in model precision, not only in quantitative terms (level of discretization in each spherical coordinate), but also qualitative terms (relative orientation between mesh and quadrature). In the enhanced model, each simulation will use a number of quadrature equal to the number of different sources (two different isotropic sources are a single source, as are two coneshaped sources with the same direction and opening; however, isotropic and cone-shaped sources use different quadrature schemes, or two cone-shaped sources with different directions or openings). Each quadrature adapts to its kind of light source, obtaining higher precision without increasing angular discretization. Quadrature rotates to place a discrete ordinate in the direction of parallel sources (Figure 4B), instead of shifting the source direction to the nearest ordinate, as a standard DOM does. In cone-shaped sources, the “north pole” of the quadrature fits the cone central direction, and the

14

nearest circle of latitude is moved to the cone limit. Figure 4C describes the cone-shaped adaption.

Figure 4. A) Quadrature scheme (15 x 15 discretization); B) quadrature rotation in parallel sources; C) quadrature adaption for cone-shaped sources.

To the best of the authors’ knowledge, CFD frameworks lack specific models to simulate LEDs. Consequently, most papers utilize cone-shaped sources [20], [21], analytical equations in transparent media and boundaries [22], or more simplified options. Power cosine has been proven as a reliable scheme to model LEDs emission [23], and it has been included in the enhanced model to offer a more accurate source to users. In this model, radiation is distributed along all possible directions from a given face (a hemisphere), but intensity decreases as discrete ordinates move away from the central direction, as a power of cosine: (Eq. 19)

15

where direction;

is the intensity in a given discrete ordinate;

is the intensity in the central

is the angle between them; and g is a parameter of the model. Both

and g

can be obtained from the LED opening and net flux (data usually found in the manufacturer’s technical specifications) if the LED has a single peak emission. In LEDs with a double peak emission, such as batwing or similar directivity scheme, they can be simulated as a combination of power cosine sources [23], but in different directions. They could be easily added using the power cosine class as a model. (Eq. 20)

(Eq. 21)

The OpenFOAM code of the novel enhanced DOM model presented in this work is available online at: https://gitlab.com/Jose_Moreno/openfoam_dom

3. Results 3.1. Standard features validation and verification The model was tested in a cubical geometry that allows an analytical solution. Predictions of all individual features included in the model are validated with the analytical solution and verified with the same simulations in ANSYS Fluent, in which the DOM model is established and has been tested in numerous previous works. Preconditioned biconjugate Gradient (PBiCG) with Gauss-Seidel as a smoother and Diagonal Incomplete Cholesky (DILU) as a preconditioner were chosen to carry out the simulations during the model verification.

16

3.1.1. Geometry and mesh A cube with 1 m side and 8000 hexahedral cells (20 x 20 x 20) was used, as shown in Figure 5, with the same regular cells in both software. For transmission simulations, a cube of the same dimensions and mesh was added beneath the first one (Figure 5B). Angular discretization used for the validation tests was 15 x 15 (15 divisions in θ, and 15 in φ, in each sphere octant). Parallel emission and volumetric absorption validations use the same mesh described in 3.2. and 3.3. (Figure 6A).

Figure 5. Meshes used for isotropic emission validation.

17

3.1.2. Individual phenomena validations and verifications Isometric emission combined with absorption, reflection and parallel transmission have been individually validated through the view factor method, having an analytical solution in a cube or geometries of similar simplicity. The top wall is declared as an isotropic source that emits 100 W m-2. Net fluxes in each face were chosen to estimate errors. This group of tests have also been verified in comparison to ANSYS Fluent simulations. Parallel emission has been geometrically validated, checking the radiation distribution in an opposite face and the value of fluxes in both sides. Volumetric absorption was validated in combination with parallel emission, testing incident radiation in the opposite face. Parallel emission was not verified in this section, as it uses a different approach than the standard DOM and will be studied in 3.3. Volumetric absorption was verified using isotropic emission to avoid differences based on source simulation. Finally, volumetric scattering was verified, combined with isotropic emission, in both software. A global summary of errors and differences between the enhanced DOM and the standard DOM are shown in Table 1. Table 1. Common features validation and/or verification Total error in fluxes

Net flux

(compared to analytical solution)

(error in energy conservation)

Enhanced DOM

Standard DOM

Enhanced DOM

Standard DOM

Isometric emission

0.16

0.08

2·10-5

0.007

Diffuse reflection

12.38

12.88

0.016

0.006

Specular reflection

0.00

0.00

0.007

0.004

18

Parallel transmission

0.26

0.21

4·10-4

5·10-4

Parallel emission

0.00

X

0

X

Volumetric absorption

0.008

X

0.008

X

0.18

9·10-5

3·10-4

Volumetric scattering

Difference between models:

3.1.3. Combined features verification This case combines all the defined features (apart from non-isotropic sources) to verify the behaviour of the model when several features are coupled. As there is no analytical solution, the case was just verified in comparison to ANSYS Fluent equivalent simulation. Figure 6 presents the contour of fluxes in both simulations, with minor differences in the top face and interface. Differences in the top face are due to post-processing, as the higher colour in the rainbow scale is not the same in either software. The face fluxes integrals are presented in Table 2 to study the possible origin of differences at the interface. It seems that the values of both interfaces are close to zero, but there is a small region in 6A where values rise over zero.

19

Figure 6. Verification of the combined case. (A) Standard DOM; (B) Enhanced DOM.

Table 2. Combined features verification

Top Higher sides Interface Lower sides Bottom

Boundary flux integral (W) Enhanced DOM Standard DOM Difference -41.04 -41.00 0.10% 13.06 13.64 4.25% 22.67 23.69 4.31% -84.90 -86.66 2.07% 27.95 27.35 2.15%

Region average incident radiation (W m-2) Enhanced DOM Standard DOM Difference Top region 150.32 155.35 3.24% Bottom region 45.80 43.66 4.67%

3.2. Enhancements in LEDs modelling A new mesh (Figure 7A) was designed in ANSYS Fluent with approximately 64 k cells and a circled surface (~ 15 cm radius) in the top face to be declared as a lamp. The mesh was

20

imported to OpenFOAM using Fluent3DmeshToFoam tool to avoid differences based on geometry or cell distribution. To test cone-shaped emission, two general cases were used, one with the cone central direction along the z-axis (Figure 7B), and a sloped cone with the direction contained in the xz-plane (Figure 7C). The first kind of simulation (cone normal to top face, studied in sections 3.2.1 and 3.2.2) does not rotate the quadrature, as the central direction is already in the z-axis. Both software will have the same quadrature orientation, allowing the individual effect of quadrature adjustment to the cone opening to be observed. The second set of simulations (sloped cone, section 3.2.3) use a constant opening of 36º in a constant angular discretization of 15 x 15, avoiding quadrature adjustment. All of the differences observed in simulations will be produced by quadrature rotation.

Figure 7. A) Mesh used for cone-shaped emission validation; B) emission parallel to the z-axis; C) emission to the yz-plane.

3.2.1. Effect of viewing angle The first set of cone emission cases was designed to study the effect of cut adjustment in cone opening. As the standard DOM uses fixed positions of discrete ordinates in the quadrature, the cone opening will behave as a discrete range of values with a finite set of possible cone

21

openings. This model, however, adjusting the ordinates cut, will behave as a continuous range. Figure 8 shows the contours of incident radiation in the bottom face and emission direction in the z-axis (see Figure 7B). Two different viewing angles were used. The simulation of the 36º opening cone is similar in both software, with minor differences based on the numerical solver (Figure 8A and B). The external contour of Figure 8A has been used as a template in the rest of the figures to observe changes in cone shape. Figure 8C (34º opening cone source in the standard DOM) shows the same shape as 8A and 8B (36º opening). Moreover, both 34º and 36º belong to the same discrete range in the standard DOM. Figure 8D presents a smaller cone, as the enhanced model is able to simulate a continuous range.

Figure 8. Study of cone opening: incident radiation in the bottom face, for a viewing angle of 36º ( A) enhanced DOM, B) standard DOM), and a viewing angle of 34º ( C) enhanced DOM, D) standard DOM).

For additional clarity, the graph below (Figure 9) represents the profiles of incident radiation along a line in the centre of the bottom face, for the cases shown in Figure 8. Enhanced DOM cones (lines) are different in height and width (area must be equal to preserve energy). Standard DOM ones (dots) have the same width (discrete behaviour of opening), but different height, leading to a false higher flow in the 34º cone.

22

Figure 9. Profiles in cone opening: incident radiation at the bottom wall.

3.2.2. Effect of discretization The second set of cone emission simulations was designed to elucidate the effect of angular discretization in cone shape and size. Instead of changing the cone opening, three different angular discretization (12 x 12, 15 x 15, and 17 x 17) were used to simulate the same cone opening (36º). Figure 10 shows the results of incident radiation in the bottom face for these cases. Images in column B are the simulations already shown in Figure 8 (model case, 15 x 15 discretization). Simulations with a discretization of 12 x 12 and 17 x 17 (Figure 10A and C, respectively) for the enhanced DOM resulted in the same shape with minor differences in the profiles. Conversely, these simulations for the standard DOM differ not only in size (radius is adjusted to the number of ordinates inside of the cone), but also in shape.

23

Figure 10. Effect of discretization in cone shape and size. Incident radiation in the bottom face. Enhanced DOM (top); standard DOM (bottom). Angular discretization: A) 12 x 12, B) 15 x 15, C) 17 x 17.

By comparing all simulations in the standard DOM, a tendency can be observed in behaviour founded in the cone centre, with a lower illuminated zone decreasing in size with a higher discretization. This lower illuminated zone disappears for the higher discretization (Figure 10C). This is ascribed to a lack of discrete ordinates in the centre of the cone, that correspond to the “south pole” of the quadrature. The enhanced model also lacks this central ordinate and would exhibit the same low lighter centre if distance from the source increases. However, the developed algorithm that moves the cut between ordinates resets the value of

(angular

distance between rows of the quadrature) in both sides of the cut, placing a higher number of ordinates inside of the cone. Figure 11 plots the incident radiation profiles at the bottom face. Dashed lines were chosen in the enhanced DOM 15 x 15 and 17 x 17 profiles, as they follow the same tendency. This is because, in both simulations, fitting the cut between ordinates leads to the same number of 24

ordinates inside of the cone, and thus the simulations are equal inside of the cone. Of course, differences in discretization still affect radiation scattering or reflection. The black line (enhanced DOM 12 x 12) has almost the same width, but a less accurate profile in the centre. This is due to the higher size of inner ordinates (not only do 12 x 12 simulations always have bigger ordinates than 15 x 15, but the adjustment of the cut also makes them bigger), making it more obvious that the cone is, in fact, a bunch of light beams.

Figure 11. Study of angular discretization in cones: incident radiation at the bottom wall.

Standard DOM simulations in 15 x 15 and 17 x 17 are quite similar to the enhanced DOM ones, showing the adjustment in both discretization as a middle-point between both standard DOM simulations. The standard DOM 12 x 12 case exhibits an opposite behaviour to the enhanced DOM 12 x 12. In this case, ordinates inside of the cone are not enlarged to fit the cone opening. Instead, surrounding ordinates are included in the cone. A bigger region is illuminated with a minor flux of light.

25

3.2.3. Effect of direction To investigate the effect of cone direction, a new set of simulations was performed with the cone emitting through the yz-plane. For these cases, results are not represented at the bottom mesh surface, but are represented in the wall where the cone-shaped beam finishes (left wall in Figure 7C). Three cases were compared between both software with a small change in the z-direction coordinate (in a normalized direction vector, the z-coordinate is the cosine with a vertical direction) of the emitted light (z = 0.7, z = 0.8, and z = 0.9). The cone emission aperture was kept constant, as well as angular discretization (36º, 15 x 15). Figure 12 shows the obtained results. The first row in Figure 12 represents the differences in the ordinates filling the cone shape due to adaptive quadrature behaviour. The red lines represent the ordinates created by the enhanced DOM, with a higher number of ordinates in the cone centre as a result of quadrature rotation. The black rectangles represent the ordinates chosen by the standard DOM to reproduce cone shape. As can be seen in the graphs, the rotation of quadrature allows the enhanced DOM to fit the cone shape (second row in Figure 12). In addition, the total number of discrete ordinates to distribute the light is 180 in enhanced DOM and 40-44 in standard DOM (Figure 12A).

26

Figure 12. Study of change in direction (cone emission opening 36º). Each column describes the same case (direction z = 0.7, z = 0.8, z = 0.9). A) Distribution of discrete ordinates in both software; B) incident radiation in the side wall in the enhanced DOM; C) incident radiation in the side wall in the standard DOM.

It can be easily concluded that adaptive quadrature enhances cone source simulation. This fact constitutes an important capability for the simulation of LED-based systems, as their energy efficiency and the possibility of almost monochromatic emission make them an important tool in photoactivated processes [24], [25].

27

3.2.4. Development of a power cosine model Power cosine has already been tested as a more accurate model to simulate LEDs than cone emission. It is not affected by quadrature adjustment, as it lacks a discontinuity in emission, and will keep approximately the same number of discrete ordinates emitting with or without quadrature rotation. Both models, power cosine with quadrature rotation and a model with no quadrature modification, have been implemented and tested to estimate the rotation effect. The results were also compared with cone-shaped source simulations in the enhanced DOM. The main differences between cone-shaped and power cosine results (Figure 13A and B, respectively) are the shape of the emission front (Figure 13, bottom), higher radiation near the emission central direction in the cone shape, and a higher affected region in power cosine. These differences come from emission directivity, represented in Figure 14A, over a typical LED datasheet. Figure 13B, which shows the light distribution of non-rotated power cosine emission, presents an irregular distribution in the central region, compared to Figure 13C. When the quadrature is rotated, discrete ordinates grow from the centre of the emission, with higher precision in regions with higher emission. However, in the x-axis emission, with no quadrature rotation, discrete ordinates distribution matches the lowest precision region (biggest ordinates) with highest emission. Both distributions are shown in Figure 14B.

28

Figure 13. LED emission along the x-axis. Net flux in the opposite face (top); incident radiation in the central slice (bottom). (A) Rotated cone; (B) non-rotated power cosine; (C) rotated power cosine

29

Figure 14. (A) Emission directivity; (B) discrete ordinates distributions in quadrature in x-axis emission of power cosine, up rotated, bottom non rotated.

3.3. Enhancements in parallel emission The final kind of emission source studied is based on parallel emission (e.g., sun beams). The previous mesh was selected to test parallel emission. Again, light was simulated using a normal and a sloped beam (see Figure 15), but in this case, quadrature already rotates, even in simulations along the z-axis. Moreover, the quadrature rotates to fit the beam direction with the first ordinate over the quadrature equator. Ordinates around the quadrature equator have the most squared shape, and thus are the most representative of the directions inside of its solid angle (relative standard deviation of direction depends on ordinate shape, with a circle being the most representative, and a triangle the least).

Figure 15. A) Mesh used for parallel emission; B) emission to the z-axis; C) emission to the yz-plane.

3.3.1. Effect of discretization The first set of simulations, with the beam normal to the emission surface (along the z-axis), was selected to study the effect of discretization in beam direction. ANSYS Fluent chooses the ordinate containing the beam direction, and so the beam is reallocated to the ordinate

30

central direction. As ordinates around the z-axis are triangle shapes with a common vertex in the pole, precision in light direction grows with θ angle discretization (as triangle sizes decrease, and their central point approaches the pole). Figure 16 shows the results obtained when four different discretization values are used: 3 x 3, 10 x 10, 15 x 15, and 20 x 20. Black lines represent the analytic solution. As can be seen, error in standard DOM decreases with angular discretization (from Figure 16A to 16D), with a higher error in the y-axis than the x-axis, as discrete ordinates near the zaxis have a narrow shape, and the first ordinate is located next to the y-axis. The orange colour of the circle in D shows how the software is calculating the beam intensity with the declared direction and not with the simulated one. Consequently, the flux in the emission surface and in the illuminated face is lower than the one declared by the user. However, these errors are avoided with the implemented code for OpenFOAM. Specifically, as it rotates the quadrature to fit a central ordinate in the beam direction (an ordinate near the sphere equator, with similar size in θand φ), it obtains a correct beam direction, even with the smallest discretization (3 x 3). Regardless, a higher discretization will be necessary if there is scattering or reflectors.

31

Figure 16. Emission normal to the surface (z-axis) in enhanced DOM (top) and standard DOM (bottom). Study of the discretization effect: A) 3 x 3, B) 10 x 10, C) 15 x 15, D) 20 x 20.

3.3.2. Effect of direction The second set of parallel emission simulations were selected to elucidate the relation between light direction and angular discretization. A combination of direction and discretization was chosen to investigate those cases in which light discretization causes an undesired change in direction or prevents a desired change. On the one hand, the first three columns in Figure 17 show the effect of light discretization in beam direction. Having a constant direction declared by the user, the beam should be pointing to the same location, but as standard DOM chooses the nearest discrete ordinate, this changes with angular discretization. On the other hand, the last column shows how a change in the direction declared by the user may be updated (enhanced) or not (standard) in the program. In both ANSYS 15 x 15 cases, the nearest discrete ordinate is the same, and thus the simulation

32

is almost the same (flux in the boundary, and incident radiation in ordinates are calculated with the user-declared direction, which has changed).

Figure 17. Study of the discretization effect in beam direction in enhanced (up) and standard (bottom) DOM. The first three columns show the study of the discretization effect (A) 12 x 12, B) 15 x 15 and C) 17 x 17) with the same direction (z=0.7), while the last column shows a change in direction (z=0.75) with the same discretization as B) (15 x 15). The black lines are the analytical solution.

Again, adaptive quadrature produces a superior performance in parallel light sources. This kind of emission source area is already very important and mainly used in solar simulations [26]. The solar vector is transformed into a continuous variable and not bounded in angular discretization options.

33

3. Conclusions A novel Discrete Ordinate Method module has been implemented in OpenFOAM able to simulate diffuse and specular phenomena in both surfaces and volumetric media. The adaptation of quadrature to the emission sources makes this model as accurate as ray-tracing models for the simulation of non-isotropic sources, while isotropic sources are still better. Whereas the fixed quadrature scheme employed in standard DOM misleads the direction of parallel emission sources, and leads to irregular cuts in the directional distribution of coneshaped sources by using the nearest available ordinates, the power cosine law emission sources and the quadrature adaptation included in the enhanced DOM model improves significantly the accuracy in the simulation of solar or LED radiation sources, being able to adapt to subtle changes in cone direction or opening. The standard features of the model were validated against an analytical solution and verified with simulations in already-tested DOM models, with no appreciable differences found. The main strengths (diffuse reflection or scattering) and weakness (specular reflection) of the model follow the expected trend. The effect of adaptive quadrature was also carefully investigated by individually changing the parameters involved (i.e., angular discretization, central direction, and opening). Future work will include master/slave behaviour for sources, to use the same quadrature for multiple sources (to compare with other software, or combine compatible sources), and scalar/vector storage options for radiant intensity (to maintain precision in specular reflection in low discretization, and as an alternative way to improve source simulations). Based on the results, the enhanced model can be considered as an effective tool not only for photoactivated process simulations with OpenFOAM, but also for an enormous range of applications in which radiation transfer is involved.

34

Acknowledgements This work was supported by the Spanish State Research Agency (AEI) and the Spanish Ministry of Economy and Competitiveness (MINECO) [WATER4FOOD project, CTQ201454563-C3] and Comunidad de Madrid [FOTOCAOS project, Y2018/EMT-5062].

35

Symbol List

Absorbed energy (W m-3) Reflection diffuse fraction Transmission diffuse fraction Henyey-Greenstein shape coefficient Incident radiation (W m-2) Radiant intensity Scattering phase function Lamp emission flux in a wavelength band (W m-2) Flux in a wavelength band (W m-2) Incoming flux in a wavelength band (W m-2) Reflection factor in boundaries Transmission factor in boundaries

Abbreviations

RTE

Radiative Transfer Equation

DOM

Discrete Ordinate Method

Greek letters

Opening angle of cone-shaped sources Angle between cone central direction and surface normal Angle between ordinates in phase function calculation Volumetric absorption coefficient for a specific wavelength band Wavelength Volumetric scattering coefficient for a specific wavelength band Solid angle modulus (representative area of discrete ordinate in the sphere)

Vector symbols

36

Solid angle of discrete ordinate (direction and modulus) Normal vector of face Discrete ordinate normalized direction vector

37

References [1]

C. A. Martin, G. Camera-Roda, and F. Santarelli, “Effective design of photocatalytic reactors: influence of radiative transfer on their performance,” Catal. Today, vol. 48, no. 1–4, pp. 307–313, Jan. 1999.

[2]

A. E. Cassano, C. A. Martin, R. J. Brandi, and O. M. Alfano, “Photoreactor Analysis and Design: Fundamentals and Applications,” Ind. Eng. Chem. Res., vol. 34, no. 7, pp. 2155–2201, Jul. 1995.

[3]

Y. Paixão De Almeida, P. Laranjeira, C. Lage, L. F. Lopes, and R. Silva, “Validation of LES simulation of a turbulent radiant flame using OpenFOAM,” in 22nd International Congress of Mechanical Engineering (COBEM 2013), 2013.

[4]

K. Cid and S. Vianna, “A COMPARATIVE STUDY BETWEEN THERMAL RADIATION MODELS P-1 AND DISCRETE ORDINATES USING CFD SOFTWARE OPENFOAM,” in Anais do Congresso Brasileiro de Fluidodinâmica Computacional, 2016.

[5]

A. Heimsath, F. Cuevas, A. Hofer, P. Nitz, and W. J. Platzer, “Linear Fresnel Collector Receiver: Heat Loss and Temperatures,” Energy Procedia, vol. 49, pp. 386–397, 2014.

[6]

M. Hettel, E. Daymo, and O. Deutschmann, “3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO,” Comput. Chem. Eng., vol. 109, pp. 166–178, Jan. 2018.

[7]

B. Kong and R. D. Vigil, “Simulation of photosynthetically active radiation distribution in algal photobioreactors using a multidimensional spectral radiation model,” Bioresour. Technol., vol. 158, pp. 141–148, Apr. 2014.

[8]

A. Habibi, B. Merci, and G. J. Heynderickx, “Impact of radiation models in CFD simulations of steam cracking furnaces,” Comput. Chem. Eng., vol. 31, no. 11, pp. 1389–1406, Nov. 2007.

38

[9]

W. Ge, R. Marquez, M. F. Modest, and S. P. Roy, “Implementation of High-Order Spherical Harmonics Methods for Radiative Heat Transfer on OpenFOAM,” J. Heat Transfer, vol. 137, no. 5, p. 052701, May 2015.

[10] D. Radice, E. Abdikamalov, L. Rezzolla, and C. D. Ott, “A new spherical harmonics scheme for multi-dimensional radiation transport I. Static matter configurations,” J. Comput. Phys., vol. 242, pp. 648–669, Jun. 2013. [11] D. J. Drake, “View-factor method for solving time-dependent radiation transport problems involving fixed surfaces with intervening, participating media,” J. Comput. Phys., vol. 87, no. 1, pp. 73–90, Mar. 1990. [12] M. Avila, R. Codina, and J. Principe, “Spatial approximation of the radiation transport equation using a subgrid-scale finite element method,” Comput. Methods Appl. Mech. Eng., vol. 200, no. 5–8, pp. 425–438, Jan. 2011. [13] V. Vikas, C. D. Hauck, Z. J. Wang, and R. O. Fox, “Radiation transport modeling using extended quadrature method of moments,” J. Comput. Phys., vol. 246, pp. 221– 241, Aug. 2013. [14] P. Kuczyński and R. Białecki, “Radiation heat transfer model using Monte Carlo ray tracing method on hierarchical ortho-Cartesian meshes and non-uniform rational basis spline surfaces for description of boundaries,” Arch. Thermodyn., vol. 35, no. 2, pp. 65–92, Jan. 2014. [15] R. Kong, M. Ambrose, and J. Spanier, “Efficient, automated Monte Carlo methods for radiation transport,” J. Comput. Phys., vol. 227, no. 22, pp. 9463–9476, Nov. 2008. [16] J. C. Chai, H. S. Lee, and S. V. Patankar, “Ray effect and false scattering in the Discrete Ordinates Method,” Numer. Heat Transf. Part B Fundam., vol. 24, no. 4, pp. 373–389, Dec. 1993. [17] V. K. Pareek and A. A. Adesina, “Light intensity distribution in a photocatalytic

39

reactor using finite volume,” AIChE J., vol. 50, no. 6, pp. 1273–1288, 2004. [18] I. ANSYS, “ANSYS User and Theory Guide, ANSYS Fluent, Release 14.5,” Cecil Townsh., 2012. [19] J. Marugán, R. van Grieken, A. E. Cassano, and O. M. Alfano, “Intrinsic kinetic modeling with explicit radiation absorption effects of the photocatalytic oxidation of cyanide with TiO2 and silica-supported TiO2 suspensions,” Appl. Catal. B Environ., vol. 85, no. 1–2, pp. 48–60, 2008. [20] C. Casado et al., “Design and validation of a LED-based high intensity photocatalytic reactor for quantifying activity measurements,” Chem. Eng. J., vol. 327, pp. 1043– 1055, Nov. 2017. [21] A. Sergejevs et al., “A calibrated UV-LED based light source for water purification and characterisation of photocatalysis,” Photochem. Photobiol. Sci., vol. 16, no. 11, pp. 1690–1699, Nov. 2017. [22] Z. Wang, J. Liu, Y. Dai, W. Dong, S. Zhang, and J. Chen, “CFD modeling of a UVLED photocatalytic odor abatement process in a continuous reactor,” J. Hazard. Mater., vol. 215, pp. 25–31, 2012. [23] I. Moreno and C.-C. Sun, “Modeling the radiation pattern of LEDs,” Opt. Express, vol. 16, no. 3, p. 1808, Feb. 2008. [24] M. R. Eskandarian, H. Choi, M. Fazli, and M. H. Rasoulifard, “Effect of UV-LED wavelengths on direct photolytic and TiO2 photocatalytic degradation of emerging contaminants in water,” Chem. Eng. J., vol. 300, pp. 414–422, Sep. 2016. [25] A. Uhl, B. W. Sigusch, and K. D. Jandt, “Second generation LEDs for the polymerization of oral biomaterials.,” Dent. Mater., vol. 20, no. 1, pp. 80–7, Jan. 2004. [26] Y. Achdou, “Numerical optimization of a photocell,” Comput. Methods Appl. Mech. Eng., vol. 102, no. 1, pp. 89–106, Jan. 1993.

40

HIGHLIGHTS 

Novel multiregional DOM with anisotropic scattering developed for OpenFOAM



Successful validation vs analytical solution and verification vs standard DOM



Quadrature rotation reduces significantly the directional error in parallel beams



Power cosine or cone shaped adjusted quadrature for improved simulation of LED

41