# A Radiative Transfer Modeling Methodology in Gas-Liquid Multiphase Flow Simulations

Authors: Gautham Krishnamoorthy; Rydell Klosterman; Dylan Shallbetter
Journal: Journal of Engineering
 (2014)

---

## Abstract

A methodology for performing radiative transfer calculations in computational fluid dynamic simulations of gas-liquid multiphase flows is presented. By considering an externally irradiated bubble column photoreactor as our model system, the bubble scattering coefficients were determined through add-on functions by employing as inputs the bubble volume fractions, number densities, and the fractional contribution of each bubble size to the bubble volume from four different multiphase modeling options. The scattering coefficient profiles resulting from the models were significantly different from one another and aligned closely with their predicted gas-phase volume fraction distributions. The impacts of the multiphase modeling option, initial bubble diameter, and gas flow rates on the radiation distribution patterns within the reactor were also examined. An increase in air inlet velocities resulted in an increase in the fraction of larger sized bubbles and their contribution to the scattering coefficient. However, the initial bubble sizes were found to have the strongest impact on the radiation field.


---

## Body

## 1. Introduction

### 1.1. Challenges in Modeling Multiphase Radiative Transfer

Modeling radiative transfer in multiphase flows is important in several applications such as solid fuel combustors [1, 2], externally irradiated gasifiers [3], photocatalytic reactors [4, 5], and photobioreactors (PBRs) [6, 7]. While the procedure for coupling radiative transfer with the hydrodynamics has been well established in dilute multiphase flows (local dispersed phase volume fractions less than 10%) such as pulverized fuel combustors [8], the effect of radiative transfer is often neglected or grossly simplified in computational fluid dynamic (CFD) simulations where all the phases are present in significant fractions such as bubbling bed and circulating fluidized bed gasifiers [9, 10]. This simplification often takes the form of an “optically thin” radiation exchange between the phases to approximate the radiative source term in the phase energy equations. In the optically thin approximation, a radiation temperature of the phases is computed and is employed in conjunction with the phase thermodynamic temperature, an empirical radiative heat transfer coefficient to compute the radiative source term, and consequently the temperature change in the phases resulting from radiative heat exchange. Since a rigorous solution to the radiative transfer equation (RTE) is not carried out in this approach, the optically thin approximation cannot predict the radiative fluxes at different surfaces or the distributions of incident radiation within a reactor which is an important variable of interest in PBRs and photocatalytic systems.The reasons for this lack of rigorous coupling between radiative transfer and hydrodynamics in nondilute multiphase flows may be attributed to three issues: computational cost, complexity, and model incompatibilities. The spatial, temporal, directional, and wavelength dependencies resulting from the complex interactions of light with different multiphase media can add a significant computational cost to the simulations. Adequate resolution of the flow features in a multiphase reactor typically requires time-step sizes on the order of10-3 to 10-2 seconds for a simulation spanning several seconds/minutes. Furthermore, the phases within a reactor may exhibit strongly directionally dependent scattering characteristics such as air bubbles in a bubble column [11, 12], coal and ash particles in a combustor [8], and the algae in a PBR [13]. Similarly, the fluid media (such as H2O (l), CO (g), CO2(g), H2O (g)) may have wavelength dependent absorption coefficients.Complexity in the simulation of radiative transfer in nondilute multiphase flows is introduced by the need to incorporate additional equation and terms within current RTE solution frameworks of CFD codes. These additional equations and terms may arise from volume averaging procedures that are employed to derive a multiphase/vector version of the RTE [14–16] and the need to account for dependent scattering effects [17, 18] in extremely dense flows.In terms of model incompatibilities, the RTE requires as inputs the absorption and scattering coefficients of the phases which are determined from the absorption and scattering efficiencies, respectively, and the surface area density of the phases [9]. However, multiphase models that are most commonly used in multiphase CFD simulations such as the two-fluid model (TFM) and population balance (PB) and volume of fluid (VOF) models provide as outputs the spatial and temporal variations in the volume fractions and the number densities. Consequently, the surface area density information must first be extracted from these variables to determine the absorption and scattering coefficients.


### 1.2. Scope of the Present Study

Clearly, overcoming all of the above mentioned difficulties associated with a rigorous radiation-multiphase hydrodynamic coupling will require significant increases in computational power and significant modifications to the current frameworks for solving the RTE. However, before embarking on such a huge undertaking, it is important to get a feel for the accuracy gained at each level of modeling enhancement. It is important, for instance, to do the following.(1)
To assess the sensitivities of the radiative transfer predictions to the prediction variations among the multiphase models. Multiphase models can vary considerably in their accuracies and computational cost, from the fast and less accurate mixture model (MM) to the more expensive and accurate PB based models that account for bubble growth and coagulation. While high fidelities in the multiphase modeling [19], turbulence modeling [19, 20], drag laws [21], and bubble breakup and coalescence closures [22, 23] are deemed necessary to accurately predict the gas hold-up and phase velocities, it is important to examine their impact on the radiation field predictions within the reactor.(2)
To assess the sensitivities of the radiative transfer predictions to the radiative properties of the participating medium (bubbles) in the application of interest. The absorption and scattering efficiencies as well as the scattering phase function are the parameters that determine the phase radiative properties. Several state-of-the-art open-source codes [24] are available to obtain the absorption and scattering efficiencies of the dispersed media, for instance, Bohren and Huffman’s Mie subroutines [25] for calculating scattering and absorption by spheres as well as T-matrix based methods [26]. However, these high fidelity models are associated with higher computational costs and it is important to ascertain the importance of determining these properties accurately in our application of interest. For instance, in our previous study of oxy-coal combustion [27], employing Mie theory within a CFD framework to determine the coal and ash particle scattering efficiencies resulted in radiative transfer predictions identical to those obtained from employing a constant scattering efficiency since the radiation was dominated by the gas-phase.Therefore, the goal and scope of this paper are to demonstrate for the first time a coupling methodology for obtaining a solution to the RTE when employed in conjunction with the dense-discrete phase model (DDPM) and MM, TFM, and PB based multiphase models. By taking an externally irradiated bubble column photoreactor as our model system, its multiphase hydrodynamics was simulated employing the four multiphase models. The second goal of this paper is to assess the sensitivities of the radiative transfer predictions to the variations in the phase predictions among the models and scattering properties of the phases. The effects of gas flow rates and bubble sizes on the average light distribution within the reactor are also investigated.The novel aspect of this study is the development of user-defined functions (UDFs) to compute the surface area density and the scattering coefficients of the dispersed phases from the volume fraction, and the fractional contribution of each bubble size to the bubble volume, that were output from these multiphase models in the commercial CFD code ANSYS FLUENT (version 12) [9]. We emphasize that the goal of this study is not validation, but initial results towards the demonstration of a modeling methodology to couple radiative transfer and other multiphase models in a CFD framework. First, the reactor configuration and flow conditions were carefully chosen to result in multiphase flow regimes in which all four multiphase flow models explored in this study can be reasonably applied. Second, experimental data that enable simultaneous validation of multiphase models (gas hold-up, bubble sizes) and radiative transfer variables (directional variations in the radiative intensity) in such reactor configurations are currently lacking. Previous studies have considered radiative measurements in well-controlled, steady-state experiments [28] or have carried out 1D radiative transfer calculations in a fixed ensemble of dispersed phase particles [29].Although only a bubble column photoreactor is investigated in this study with interactions between the light, water, and air bubbles, a natural extension of this study can include simulations of bubble column PBRs with the inclusion of algae absorption and scattering. Bubble columns have shown good mass transfer rates, mixing, low shear stress, and potential for achieving good scalability in addition to low capital cost, lack of moving parts, and high surface area to volume ratios [30–32]. Efficient mixing and mass transfer are enabled by bubbling gas through a sparger with perforated plates employed in large bubble columns to facilitate bubble breakup and coalescence. While the bubble column hydrodynamics are affected by bubble size distributions and the interfacial bubble surface areas impact mass transfer rates, the spatial variation of light within a PBR is one of the most important limiting factors governing algal growth [33]. At optimal intensities, the microalgae growth rate is directly proportional to the light intensity, up till saturation levels. Any further increases in light intensity cause inhibition of cell growth, a phenomenon known as photoinhibition [34]. Through CFD simulations of an internally radiated air-lift PBR, Wheaton and Krishnamoorthy [7] showed that the radiation distribution within the reactor was sensitive to the bubble size and the air flow rate. Therefore, an accurate determination of the radiative transfer within a bubble column PBR is important.


## 1.1. Challenges in Modeling Multiphase Radiative Transfer

Modeling radiative transfer in multiphase flows is important in several applications such as solid fuel combustors [1, 2], externally irradiated gasifiers [3], photocatalytic reactors [4, 5], and photobioreactors (PBRs) [6, 7]. While the procedure for coupling radiative transfer with the hydrodynamics has been well established in dilute multiphase flows (local dispersed phase volume fractions less than 10%) such as pulverized fuel combustors [8], the effect of radiative transfer is often neglected or grossly simplified in computational fluid dynamic (CFD) simulations where all the phases are present in significant fractions such as bubbling bed and circulating fluidized bed gasifiers [9, 10]. This simplification often takes the form of an “optically thin” radiation exchange between the phases to approximate the radiative source term in the phase energy equations. In the optically thin approximation, a radiation temperature of the phases is computed and is employed in conjunction with the phase thermodynamic temperature, an empirical radiative heat transfer coefficient to compute the radiative source term, and consequently the temperature change in the phases resulting from radiative heat exchange. Since a rigorous solution to the radiative transfer equation (RTE) is not carried out in this approach, the optically thin approximation cannot predict the radiative fluxes at different surfaces or the distributions of incident radiation within a reactor which is an important variable of interest in PBRs and photocatalytic systems.The reasons for this lack of rigorous coupling between radiative transfer and hydrodynamics in nondilute multiphase flows may be attributed to three issues: computational cost, complexity, and model incompatibilities. The spatial, temporal, directional, and wavelength dependencies resulting from the complex interactions of light with different multiphase media can add a significant computational cost to the simulations. Adequate resolution of the flow features in a multiphase reactor typically requires time-step sizes on the order of10-3 to 10-2 seconds for a simulation spanning several seconds/minutes. Furthermore, the phases within a reactor may exhibit strongly directionally dependent scattering characteristics such as air bubbles in a bubble column [11, 12], coal and ash particles in a combustor [8], and the algae in a PBR [13]. Similarly, the fluid media (such as H2O (l), CO (g), CO2(g), H2O (g)) may have wavelength dependent absorption coefficients.Complexity in the simulation of radiative transfer in nondilute multiphase flows is introduced by the need to incorporate additional equation and terms within current RTE solution frameworks of CFD codes. These additional equations and terms may arise from volume averaging procedures that are employed to derive a multiphase/vector version of the RTE [14–16] and the need to account for dependent scattering effects [17, 18] in extremely dense flows.In terms of model incompatibilities, the RTE requires as inputs the absorption and scattering coefficients of the phases which are determined from the absorption and scattering efficiencies, respectively, and the surface area density of the phases [9]. However, multiphase models that are most commonly used in multiphase CFD simulations such as the two-fluid model (TFM) and population balance (PB) and volume of fluid (VOF) models provide as outputs the spatial and temporal variations in the volume fractions and the number densities. Consequently, the surface area density information must first be extracted from these variables to determine the absorption and scattering coefficients.


## 1.2. Scope of the Present Study

Clearly, overcoming all of the above mentioned difficulties associated with a rigorous radiation-multiphase hydrodynamic coupling will require significant increases in computational power and significant modifications to the current frameworks for solving the RTE. However, before embarking on such a huge undertaking, it is important to get a feel for the accuracy gained at each level of modeling enhancement. It is important, for instance, to do the following.(1)
To assess the sensitivities of the radiative transfer predictions to the prediction variations among the multiphase models. Multiphase models can vary considerably in their accuracies and computational cost, from the fast and less accurate mixture model (MM) to the more expensive and accurate PB based models that account for bubble growth and coagulation. While high fidelities in the multiphase modeling [19], turbulence modeling [19, 20], drag laws [21], and bubble breakup and coalescence closures [22, 23] are deemed necessary to accurately predict the gas hold-up and phase velocities, it is important to examine their impact on the radiation field predictions within the reactor.(2)
To assess the sensitivities of the radiative transfer predictions to the radiative properties of the participating medium (bubbles) in the application of interest. The absorption and scattering efficiencies as well as the scattering phase function are the parameters that determine the phase radiative properties. Several state-of-the-art open-source codes [24] are available to obtain the absorption and scattering efficiencies of the dispersed media, for instance, Bohren and Huffman’s Mie subroutines [25] for calculating scattering and absorption by spheres as well as T-matrix based methods [26]. However, these high fidelity models are associated with higher computational costs and it is important to ascertain the importance of determining these properties accurately in our application of interest. For instance, in our previous study of oxy-coal combustion [27], employing Mie theory within a CFD framework to determine the coal and ash particle scattering efficiencies resulted in radiative transfer predictions identical to those obtained from employing a constant scattering efficiency since the radiation was dominated by the gas-phase.Therefore, the goal and scope of this paper are to demonstrate for the first time a coupling methodology for obtaining a solution to the RTE when employed in conjunction with the dense-discrete phase model (DDPM) and MM, TFM, and PB based multiphase models. By taking an externally irradiated bubble column photoreactor as our model system, its multiphase hydrodynamics was simulated employing the four multiphase models. The second goal of this paper is to assess the sensitivities of the radiative transfer predictions to the variations in the phase predictions among the models and scattering properties of the phases. The effects of gas flow rates and bubble sizes on the average light distribution within the reactor are also investigated.The novel aspect of this study is the development of user-defined functions (UDFs) to compute the surface area density and the scattering coefficients of the dispersed phases from the volume fraction, and the fractional contribution of each bubble size to the bubble volume, that were output from these multiphase models in the commercial CFD code ANSYS FLUENT (version 12) [9]. We emphasize that the goal of this study is not validation, but initial results towards the demonstration of a modeling methodology to couple radiative transfer and other multiphase models in a CFD framework. First, the reactor configuration and flow conditions were carefully chosen to result in multiphase flow regimes in which all four multiphase flow models explored in this study can be reasonably applied. Second, experimental data that enable simultaneous validation of multiphase models (gas hold-up, bubble sizes) and radiative transfer variables (directional variations in the radiative intensity) in such reactor configurations are currently lacking. Previous studies have considered radiative measurements in well-controlled, steady-state experiments [28] or have carried out 1D radiative transfer calculations in a fixed ensemble of dispersed phase particles [29].Although only a bubble column photoreactor is investigated in this study with interactions between the light, water, and air bubbles, a natural extension of this study can include simulations of bubble column PBRs with the inclusion of algae absorption and scattering. Bubble columns have shown good mass transfer rates, mixing, low shear stress, and potential for achieving good scalability in addition to low capital cost, lack of moving parts, and high surface area to volume ratios [30–32]. Efficient mixing and mass transfer are enabled by bubbling gas through a sparger with perforated plates employed in large bubble columns to facilitate bubble breakup and coalescence. While the bubble column hydrodynamics are affected by bubble size distributions and the interfacial bubble surface areas impact mass transfer rates, the spatial variation of light within a PBR is one of the most important limiting factors governing algal growth [33]. At optimal intensities, the microalgae growth rate is directly proportional to the light intensity, up till saturation levels. Any further increases in light intensity cause inhibition of cell growth, a phenomenon known as photoinhibition [34]. Through CFD simulations of an internally radiated air-lift PBR, Wheaton and Krishnamoorthy [7] showed that the radiation distribution within the reactor was sensitive to the bubble size and the air flow rate. Therefore, an accurate determination of the radiative transfer within a bubble column PBR is important.


## 2. Method

### 2.1. Reactor Geometry and Boundary Conditions

Sánchez Mirón et al. [30] suggest that the diameter of the bubble column in a PBR should not exceed 0.2 m to ensure adequate light distribution within the reactor. Therefore, the selection of the bubble column for investigation was determined by this criterion as well as the availability of an experimental map of the flow pattern within the bubble column to ensure that the hydrodynamics was being represented adequately for the intended purposes of this study. The rectangular bubble column PBR experimentally investigated by Deen [35] was simulated. Figure 1 shows the geometric dimensions of the reactor geometry along with the boundary conditions. The square bubble column had dimensions of 0.15 m × 0.15 m × 1 m with the water level at 0.45 m. The inlet of the bubble column is at the center of the bottom plate. The experimental inlet consisted of a distributor plate with 49 holes with diameter of 1 mm at a square pitch of 6.25 mm with the plate expected to yield bubbles with a mean diameter of 4 mm [35]. Through CFD simulations of this bubble column, Hansen [36] determined that spatially resolving each hole in the distributor plate was not necessary and the inlet could be modeled as a perforated plate. Subsequently, Hansen [36] carried out several measurements of phase velocities and volume fractions and CFD simulations employing an air inlet velocity of 7.84 cm/s. Therefore, this velocity was adopted in this study as well. The outlet was located at the air-water interface and modeled as a pressure outlet with an air backflow volume fraction specified to be zero and the no-slip boundary condition was applied at the walls for both phases. The inlet was resolved employing 144 cells and the reactor resolved using 46,080 cells which corresponds exactly with the resolutions employed by Deen [35] and Hansen [36] in their study. They both obtained an acceptable agreement between the experimental results and numerical predictions at this resolution [35, 36].Figure 1
Geometry and boundary conditions for the bubble column PBR.In order to represent this bubble column as a PBR, a strip (of width 5 cm) at the center of one of the four side walls of the square bubble column was simulated as a fluorescent radiation source in this study (dark strip in Figure1). The radiation source was modeled as a semitransparent surface and the visible spectrum was divided into three wavelength regions from 400 nm to 700 nm of equal width. The radiation distribution of the fluorescent radiator in each wavelength region was specified as per our previously computed values for a fluorescent lamp [7]. These are shown in Table 1. These correspond to 25% in the wavelength region between 400 nm and 500 nm and 20.65% between 600 nm and 700 nm [37]. The absorption coefficient of the glass was set at 80 m−1having a refractive index of 1.474 in each of the wavelength regions.Table 1
Beam irradiation flux from the radiator and radiative properties of water.
Wavelength (microns)
Water absorption coefficient (m−1)
Water refractive index
Beam irradiation flux (W/m2)
0.4-0.5
0.01
1.34
105.24
0.5-0.6
0.1
1.334
228.29
0.6-0.7
0.4
1.332
86.926


### 2.2. Multiphase Hydrodynamics Modeling

CFD simulations of the bubble column investigated in this study (without the radiating surfaces) have previously been carried out employing unsteady solvers to capture the dynamic characteristics of the bubble plume [35, 36, 38, 39]. However, since the purpose of this study is to demonstrate a multiphase radiative transfer modeling methodology, the time-averaged bubble plume behavior was determined employing different models and the sensitivities of the volume averaged incident radiation within the reactor to the operating conditions and the multiphase modeling methodology were examined. The various multiphase models considered in this study are described next.


#### 2.2.1. The Two-Fluid Model (TFM)

The TFM is an Euler-Euler approach, where momentum and continuity equations are solved for each phase (water and air) in an Eulerian reference frame with a single pressure field being shared by both. In the absence of mass transfer from one phase to another as undertaken in this study, the continuity equation for each phaseq reduces to
(1)∂∂tαqρq+∇·αqρqv→q=0,
where αq and ρq are the phase volume fractions and densities, respectively. Constant densities corresponding to ambient conditions (298 K) were employed for water and air in this study.The momentum equation for each phaseq again in the absence of mass exchange between the phases reduces to
(2)∂∂tαqρqv→q+∇·αqρqv→qv→q=-αq∇p+∇·τ-q+αqρqg→+∑p=12R→pq+F→lift,q.
In (2), τ is the stress-strain tensor for the qth phase, Flift is a lift force, Rpq is an interaction force between the phases, and P is the common pressure shared among both the phases.The interaction forceRpq was computed as
(3)∑p=12R→pq=∑p=12Kpqv→p-v→q.
In (3)  Kpq is the fluid-fluid exchange coefficient defined as
(4)Kpq=αpαqρpfτp.
Here f includes a drag coefficient that is a function of the relative Reynolds number between the phases and τp is the “particulate relaxation time.” The lift force has previously been established as an important interfacial force to capture the behavior of the bubble plume [36]. In this study, the Schiller-Naumann drag law was employed to compute “f” in (4) and a lift coefficient of 0.5 was employed to compute the lift force in (2) [36]. By assuming that the turbulence transfer among the phases would play a strong role in determining the bubble dynamics, the k-ε per-phase turbulence model was employed in the simulations where equations for the turbulence kinetic energy and turbulence dissipation rate were solved for each phase [9]. Although the TFM is generally employed in transient simulations [36], it has also been employed in steady-state simulations of coal gasification in fluidized beds, for instance [40].


#### 2.2.2. The Mixture Model (MM)

The MM is a simplified multiphase model that can be employed in scenarios where the phases move at different velocities assuming local equilibrium over short spatial length scales [9]; that is, the dispersed phase is assumed to reach its terminal velocities over short spatial lengths. The MM solves a continuity and momentum equation for the multiphase mixture, a volume fraction equation for the dispersed phase and employs algebraic expressions for the relative velocities between the phases. The continuity equation for the MM is written as
(5)∂∂tρm+∇·ρmv→m=0,
where vm is the mass-averaged velocity and ρm is the mixture density:
(6)∂∂tρmv→m+∇·ρmv→mv→m=-∇p+∇·μm∇v→m+∇v→mT+ρmg→+∇·∑k=12αkρkv→dr,kv→dr,k.
The drift velocities vdr in (6) are computed from the relative velocities between the two phases. The relative velocity in turn was computed from an algebraic expression that was a function of the drag coefficient f, the particulate relaxation time τp, and an acceleration term (a). Inherent in the usage of this algebraic expression is the assumption of an “algebraic slip assumption” where local equilibrium between the phases is assumed to be reached over short spatial length scales. Again, the drag function “f” was determined from the formulation of Schiller-Naumann drag function and employed in the algebraic expression for the relative velocity. The acceleration is determined by solving the equation
(7)g→-v→m·∇v→m-∂v→m∂t=a→.
The volume fraction equation for the dispersed secondary phase (denoted by subscript s) can be written as
(8)∂∂tαsρs+∇·αsρsv→m=-∇·αsρsv→dr,s.


#### 2.2.3. Dense-Discrete Phase Model (DDPM)

The DDPM in ANSYS FLUENT follows an Euler-Lagrange approach where the fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking the bubbles through the calculated flow field. The DDPM tracks the motion of the dispersed phase in a Lagrangian reference frame but simulates dispersed phase interactions through the kinetic theory of granular flow. Initial studies have shown it to be promising with the ability to employ coarse grids resulting in significant computational savings compared to the TFM approach [41]. In the DDPM, the conservation equations of the continuous phase are also modified by taking into account the volume fraction of the dispersed phase. Therefore in the absence of mass exchange between the phases, the continuity equation for the Eulerian phase when DDPM is enabled takes the same form as the continuity equation for the TFM (cf. (1)). When only a single Eulerian phase is present such as the bubble column studied in this paper, the momentum equation gets transformed as
(9)∂∂tαpρpv→p+∇·αpρpv→pv→p=-αp∇p+∇·αpμp∇v→p+∇v→pT+αpρpg→+KDPMv→DPM-v→p+SDPM,explicit.
Volume fractions and the discrete phase velocities are obtained from the Lagrangian field solution.


#### 2.2.4. The Two-Fluid Model with Population Balance (TFM + PB)

The population balance (PB) model describes the evolution of the bubble size distribution as a result of growth, aggregation, and breakage processes. By introducing the concept of a number density functionf(x,v,t) balance equations are solved describing its evolution [9]. The general governing equation for population balance can be written as
(10)∂∂tfx,ν,t+∇·v→bx,ν,tfx,ν,t=Sx,ν,t.
In (10), S(x,v,t) accounts for the bubble birth rate, death rate, and breakup. The drag force is computed in terms of the local mean bubble density. The bubble sizes in the simulation were divided into six discrete bins (or classes) with bubble diameters ranging between 1 mm and 10 mm with a ratio exponent of 2. In the PB solution method employed in this study, all the bubbles were assumed to move at the ensemble averaged gas-phase velocity vb obtained from the solution of the TFM circumventing the need to solve continuity and momentum equations for each of the individual discrete bins. The Luo-model [9] was employed for the aggregation kernel and the Hagesather model formulation [9] was employed for the breakage with its frequency determined by the Luo-model [9].


### 2.3. Radiation Modeling

The finite volume (FV) radiation model in ANSYS FLUENT (which is a conservative variant of the discrete ordinates method) was employed in this study to model the radiative transfer. In this study, the number of divisions in the polar (Nθ) and azimuthal angles (Nφ) for the angular intensity vectors were both set to 4, corresponding to a total of 128 directions (8NθNφ) in which the RTE was being solved. Any further refinement in the angular discretization did not impact our results.If “I” represents the directional intensity, k the absorption coefficient (which is due to water alone), σ the scattering coefficient (due to the secondary phase, i.e., the bubbles alone), Ib the black body emissive power (with the water temperature set at 298 K), and Ф the scattering phase function (assumed to be forward scattering for the bubbles), then the differential equation governing the FV method can be written for each wavelength region (λ) and medium refractive index “n” as [42]
(11)∇·Iλr→,s^s^=-(kλ+σλ)Iλ(r→,s^)+ηkλn2Ibλ(r→,s^)+σλ4π∫04πIλr→,s^Φs^,s^′dΩ′.
In (11), “η” corresponds to the fraction of the total blackbody emission in the wavelength region under consideration [42].In this multiphase study, the absorption and scattering coefficients were multiplied by their corresponding phase volume fractions (αp,i) and employed in a modified version of the RTE for each wavelength region (λ) as
(12)∇·Iλr→,s^s^=-∑i=12αp,i(kλ,i+σλ,i)Iλ(r→,s^)+∑i=12ηiαp,ikλ,in2Ibλ,i(r→,s^)+∑i=12αp,iσλ,i4π∫04πIλr→,s^Φs^,s^′dΩ′.
Therefore, in this two-phase study “effective” absorption and scattering coefficients of the medium were computed through a weighted averaging of the individual phase radiative properties by their corresponding volume fractions. The bubble columns simulated in this study were assumed to be isothermal with thermal equilibrium prevailing between the two phases. The radiative transfer variables reported in this study, the incident radiative flux (qin), and the incident radiation (G) were obtained by integrating the angular intensities over all directions as
(13)qin,λr=∑λ∫4πIλr,s^s^dΩ,(14)Gr=∑λ∫4πIr,s^dΩ.
The incident radiation is a measure of the amount of radiation received by a control volume from all directions (e.g., this might represent the dosage received by algae in a photobioreactor). The volume averaged incident in the reactor was employed in this study to garner insights into the effects of different operating conditions such as air flow rate and bubble size on the light distribution. The scattering coefficient of the bubbles in (12) was computed in the DDPM model as
(15)σs=LimitV→0⁡∑N1-εsn1-fsnAsnV.
In (15), the summation is over N particles within the control volume V, εsn is the bubble emissivity (set to zero), Asn is the projected area of the nth particle, and fsn is the scattering factor associated with the nth particle. Near-field approximations have shown that, in the limitation case of spheres much larger than the wavelength of radiation (such as the bubbles examined in this study) that are embedded in an absorbing medium, the extinction efficiency approaches unity as the effects of diffraction can be neglected [43–45]. Therefore fsn in (15) was set to 0 to ensure that the light incident upon a bubble was scattered with a scattering efficiency of 1.While stand-alone radiation codes for accurately characterizing the scattering properties of spherical particles and bubbles are available [24], in multiphysics simulation choices have to be made to assess first the fidelity required in radiative transfer modeling and second the relative importance of gas and dispersed phase radiative property resolution. Since the computational cost of a wavelength dependent radiative transfer calculation in a transient multiphase simulation is already high, these assessments can help optimize the calculation times. For instance, radiative transfer calculations may be carried out once in several time-steps in the calculations if it was established that the radiative transfer is not significantly impacted by small changes to the phase volume fractions that occur over the hydrodynamic time-scales. The second issue of the relative importance of gas and dispersed phase radiative properties will largely be determined by the reactor configurations and the flow conditions within the reactor. For instance, through coal combustion simulations, Krishnamoorthy et al. [2] showed that the fidelity of the gas-phase radiative property models had little impact on the radiation distributions in a lab-scale reactor whereas differences between gray and nongray model predictions were observed in a full-scale boiler [46]. Similarly, through steam-gasification studies on directly irradiated reactor, von Zedtwitz et al. [3] determined that the radiation absorption by particles was three orders of magnitude higher than that absorbed by the gas-phase.The TFM-PB model provides the spatial variations in the probability distribution function (pdf) of the different bubble classes (bin). The scattering coefficient of each bin in the TFM-PB model calculations was calculated from the bin pdf, projected surface area (As), and the scattering efficiency (assumed to be unity) as
(16)σP,bin=volumefraction∑allclassesbubbleclasspdf×1×Ap,binVp.
The total scattering coefficient (σ) at any location for use in the RTE (cf. (12)) was then computed by a summation of the bin scattering coefficients at that location.The TFM and the mixture models do not account for bubble coalescence and breakup. Therefore, a single “class” is employed in the calculations. However, both the models provide volume fraction information as output. Therefore, in the TFM and mixture model calculations, the number density was first determined from the volume fraction information using(17)Numberdensity=BubblevolumefractionVp,
where Vp is the volume of the bubble assuming a spherical bubble shape and the bubble diameter was assumed to be 4 mm. The scattering coefficient was then determined as
(18)σs=numberdensity×As.
Equations (16)–(18) were implemented as user-defined functions (UDFs) and employed in the simulations. Flow charts depicting the exchange of variables between ANSYS FLUENT and the UDFs are shown in Figure 2. The boundary condition associated with the semitransparent radiator surface was specified as
(19)Qradiator=Qirrad.
In (19)  Qirrad is the beam irradiation flux corresponding to each wavelength region as specified in Table 1. The other three side faces of the bubble columns were set as semitransparent surfaces at a temperature of 298 K.Flow diagram depicting the exchange of variables between the CFD calculations and the user-defined functions (UDFs) for computing the scattering coefficients.
(a)
(b)


## 2.1. Reactor Geometry and Boundary Conditions

Sánchez Mirón et al. [30] suggest that the diameter of the bubble column in a PBR should not exceed 0.2 m to ensure adequate light distribution within the reactor. Therefore, the selection of the bubble column for investigation was determined by this criterion as well as the availability of an experimental map of the flow pattern within the bubble column to ensure that the hydrodynamics was being represented adequately for the intended purposes of this study. The rectangular bubble column PBR experimentally investigated by Deen [35] was simulated. Figure 1 shows the geometric dimensions of the reactor geometry along with the boundary conditions. The square bubble column had dimensions of 0.15 m × 0.15 m × 1 m with the water level at 0.45 m. The inlet of the bubble column is at the center of the bottom plate. The experimental inlet consisted of a distributor plate with 49 holes with diameter of 1 mm at a square pitch of 6.25 mm with the plate expected to yield bubbles with a mean diameter of 4 mm [35]. Through CFD simulations of this bubble column, Hansen [36] determined that spatially resolving each hole in the distributor plate was not necessary and the inlet could be modeled as a perforated plate. Subsequently, Hansen [36] carried out several measurements of phase velocities and volume fractions and CFD simulations employing an air inlet velocity of 7.84 cm/s. Therefore, this velocity was adopted in this study as well. The outlet was located at the air-water interface and modeled as a pressure outlet with an air backflow volume fraction specified to be zero and the no-slip boundary condition was applied at the walls for both phases. The inlet was resolved employing 144 cells and the reactor resolved using 46,080 cells which corresponds exactly with the resolutions employed by Deen [35] and Hansen [36] in their study. They both obtained an acceptable agreement between the experimental results and numerical predictions at this resolution [35, 36].Figure 1
Geometry and boundary conditions for the bubble column PBR.In order to represent this bubble column as a PBR, a strip (of width 5 cm) at the center of one of the four side walls of the square bubble column was simulated as a fluorescent radiation source in this study (dark strip in Figure1). The radiation source was modeled as a semitransparent surface and the visible spectrum was divided into three wavelength regions from 400 nm to 700 nm of equal width. The radiation distribution of the fluorescent radiator in each wavelength region was specified as per our previously computed values for a fluorescent lamp [7]. These are shown in Table 1. These correspond to 25% in the wavelength region between 400 nm and 500 nm and 20.65% between 600 nm and 700 nm [37]. The absorption coefficient of the glass was set at 80 m−1having a refractive index of 1.474 in each of the wavelength regions.Table 1
Beam irradiation flux from the radiator and radiative properties of water.
Wavelength (microns)
Water absorption coefficient (m−1)
Water refractive index
Beam irradiation flux (W/m2)
0.4-0.5
0.01
1.34
105.24
0.5-0.6
0.1
1.334
228.29
0.6-0.7
0.4
1.332
86.926


## 2.2. Multiphase Hydrodynamics Modeling

CFD simulations of the bubble column investigated in this study (without the radiating surfaces) have previously been carried out employing unsteady solvers to capture the dynamic characteristics of the bubble plume [35, 36, 38, 39]. However, since the purpose of this study is to demonstrate a multiphase radiative transfer modeling methodology, the time-averaged bubble plume behavior was determined employing different models and the sensitivities of the volume averaged incident radiation within the reactor to the operating conditions and the multiphase modeling methodology were examined. The various multiphase models considered in this study are described next.


### 2.2.1. The Two-Fluid Model (TFM)

The TFM is an Euler-Euler approach, where momentum and continuity equations are solved for each phase (water and air) in an Eulerian reference frame with a single pressure field being shared by both. In the absence of mass transfer from one phase to another as undertaken in this study, the continuity equation for each phaseq reduces to
(1)∂∂tαqρq+∇·αqρqv→q=0,
where αq and ρq are the phase volume fractions and densities, respectively. Constant densities corresponding to ambient conditions (298 K) were employed for water and air in this study.The momentum equation for each phaseq again in the absence of mass exchange between the phases reduces to
(2)∂∂tαqρqv→q+∇·αqρqv→qv→q=-αq∇p+∇·τ-q+αqρqg→+∑p=12R→pq+F→lift,q.
In (2), τ is the stress-strain tensor for the qth phase, Flift is a lift force, Rpq is an interaction force between the phases, and P is the common pressure shared among both the phases.The interaction forceRpq was computed as
(3)∑p=12R→pq=∑p=12Kpqv→p-v→q.
In (3)  Kpq is the fluid-fluid exchange coefficient defined as
(4)Kpq=αpαqρpfτp.
Here f includes a drag coefficient that is a function of the relative Reynolds number between the phases and τp is the “particulate relaxation time.” The lift force has previously been established as an important interfacial force to capture the behavior of the bubble plume [36]. In this study, the Schiller-Naumann drag law was employed to compute “f” in (4) and a lift coefficient of 0.5 was employed to compute the lift force in (2) [36]. By assuming that the turbulence transfer among the phases would play a strong role in determining the bubble dynamics, the k-ε per-phase turbulence model was employed in the simulations where equations for the turbulence kinetic energy and turbulence dissipation rate were solved for each phase [9]. Although the TFM is generally employed in transient simulations [36], it has also been employed in steady-state simulations of coal gasification in fluidized beds, for instance [40].


### 2.2.2. The Mixture Model (MM)

The MM is a simplified multiphase model that can be employed in scenarios where the phases move at different velocities assuming local equilibrium over short spatial length scales [9]; that is, the dispersed phase is assumed to reach its terminal velocities over short spatial lengths. The MM solves a continuity and momentum equation for the multiphase mixture, a volume fraction equation for the dispersed phase and employs algebraic expressions for the relative velocities between the phases. The continuity equation for the MM is written as
(5)∂∂tρm+∇·ρmv→m=0,
where vm is the mass-averaged velocity and ρm is the mixture density:
(6)∂∂tρmv→m+∇·ρmv→mv→m=-∇p+∇·μm∇v→m+∇v→mT+ρmg→+∇·∑k=12αkρkv→dr,kv→dr,k.
The drift velocities vdr in (6) are computed from the relative velocities between the two phases. The relative velocity in turn was computed from an algebraic expression that was a function of the drag coefficient f, the particulate relaxation time τp, and an acceleration term (a). Inherent in the usage of this algebraic expression is the assumption of an “algebraic slip assumption” where local equilibrium between the phases is assumed to be reached over short spatial length scales. Again, the drag function “f” was determined from the formulation of Schiller-Naumann drag function and employed in the algebraic expression for the relative velocity. The acceleration is determined by solving the equation
(7)g→-v→m·∇v→m-∂v→m∂t=a→.
The volume fraction equation for the dispersed secondary phase (denoted by subscript s) can be written as
(8)∂∂tαsρs+∇·αsρsv→m=-∇·αsρsv→dr,s.


### 2.2.3. Dense-Discrete Phase Model (DDPM)

The DDPM in ANSYS FLUENT follows an Euler-Lagrange approach where the fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking the bubbles through the calculated flow field. The DDPM tracks the motion of the dispersed phase in a Lagrangian reference frame but simulates dispersed phase interactions through the kinetic theory of granular flow. Initial studies have shown it to be promising with the ability to employ coarse grids resulting in significant computational savings compared to the TFM approach [41]. In the DDPM, the conservation equations of the continuous phase are also modified by taking into account the volume fraction of the dispersed phase. Therefore in the absence of mass exchange between the phases, the continuity equation for the Eulerian phase when DDPM is enabled takes the same form as the continuity equation for the TFM (cf. (1)). When only a single Eulerian phase is present such as the bubble column studied in this paper, the momentum equation gets transformed as
(9)∂∂tαpρpv→p+∇·αpρpv→pv→p=-αp∇p+∇·αpμp∇v→p+∇v→pT+αpρpg→+KDPMv→DPM-v→p+SDPM,explicit.
Volume fractions and the discrete phase velocities are obtained from the Lagrangian field solution.


### 2.2.4. The Two-Fluid Model with Population Balance (TFM + PB)

The population balance (PB) model describes the evolution of the bubble size distribution as a result of growth, aggregation, and breakage processes. By introducing the concept of a number density functionf(x,v,t) balance equations are solved describing its evolution [9]. The general governing equation for population balance can be written as
(10)∂∂tfx,ν,t+∇·v→bx,ν,tfx,ν,t=Sx,ν,t.
In (10), S(x,v,t) accounts for the bubble birth rate, death rate, and breakup. The drag force is computed in terms of the local mean bubble density. The bubble sizes in the simulation were divided into six discrete bins (or classes) with bubble diameters ranging between 1 mm and 10 mm with a ratio exponent of 2. In the PB solution method employed in this study, all the bubbles were assumed to move at the ensemble averaged gas-phase velocity vb obtained from the solution of the TFM circumventing the need to solve continuity and momentum equations for each of the individual discrete bins. The Luo-model [9] was employed for the aggregation kernel and the Hagesather model formulation [9] was employed for the breakage with its frequency determined by the Luo-model [9].


## 2.2.1. The Two-Fluid Model (TFM)

The TFM is an Euler-Euler approach, where momentum and continuity equations are solved for each phase (water and air) in an Eulerian reference frame with a single pressure field being shared by both. In the absence of mass transfer from one phase to another as undertaken in this study, the continuity equation for each phaseq reduces to
(1)∂∂tαqρq+∇·αqρqv→q=0,
where αq and ρq are the phase volume fractions and densities, respectively. Constant densities corresponding to ambient conditions (298 K) were employed for water and air in this study.The momentum equation for each phaseq again in the absence of mass exchange between the phases reduces to
(2)∂∂tαqρqv→q+∇·αqρqv→qv→q=-αq∇p+∇·τ-q+αqρqg→+∑p=12R→pq+F→lift,q.
In (2), τ is the stress-strain tensor for the qth phase, Flift is a lift force, Rpq is an interaction force between the phases, and P is the common pressure shared among both the phases.The interaction forceRpq was computed as
(3)∑p=12R→pq=∑p=12Kpqv→p-v→q.
In (3)  Kpq is the fluid-fluid exchange coefficient defined as
(4)Kpq=αpαqρpfτp.
Here f includes a drag coefficient that is a function of the relative Reynolds number between the phases and τp is the “particulate relaxation time.” The lift force has previously been established as an important interfacial force to capture the behavior of the bubble plume [36]. In this study, the Schiller-Naumann drag law was employed to compute “f” in (4) and a lift coefficient of 0.5 was employed to compute the lift force in (2) [36]. By assuming that the turbulence transfer among the phases would play a strong role in determining the bubble dynamics, the k-ε per-phase turbulence model was employed in the simulations where equations for the turbulence kinetic energy and turbulence dissipation rate were solved for each phase [9]. Although the TFM is generally employed in transient simulations [36], it has also been employed in steady-state simulations of coal gasification in fluidized beds, for instance [40].


## 2.2.2. The Mixture Model (MM)

The MM is a simplified multiphase model that can be employed in scenarios where the phases move at different velocities assuming local equilibrium over short spatial length scales [9]; that is, the dispersed phase is assumed to reach its terminal velocities over short spatial lengths. The MM solves a continuity and momentum equation for the multiphase mixture, a volume fraction equation for the dispersed phase and employs algebraic expressions for the relative velocities between the phases. The continuity equation for the MM is written as
(5)∂∂tρm+∇·ρmv→m=0,
where vm is the mass-averaged velocity and ρm is the mixture density:
(6)∂∂tρmv→m+∇·ρmv→mv→m=-∇p+∇·μm∇v→m+∇v→mT+ρmg→+∇·∑k=12αkρkv→dr,kv→dr,k.
The drift velocities vdr in (6) are computed from the relative velocities between the two phases. The relative velocity in turn was computed from an algebraic expression that was a function of the drag coefficient f, the particulate relaxation time τp, and an acceleration term (a). Inherent in the usage of this algebraic expression is the assumption of an “algebraic slip assumption” where local equilibrium between the phases is assumed to be reached over short spatial length scales. Again, the drag function “f” was determined from the formulation of Schiller-Naumann drag function and employed in the algebraic expression for the relative velocity. The acceleration is determined by solving the equation
(7)g→-v→m·∇v→m-∂v→m∂t=a→.
The volume fraction equation for the dispersed secondary phase (denoted by subscript s) can be written as
(8)∂∂tαsρs+∇·αsρsv→m=-∇·αsρsv→dr,s.


## 2.2.3. Dense-Discrete Phase Model (DDPM)

The DDPM in ANSYS FLUENT follows an Euler-Lagrange approach where the fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking the bubbles through the calculated flow field. The DDPM tracks the motion of the dispersed phase in a Lagrangian reference frame but simulates dispersed phase interactions through the kinetic theory of granular flow. Initial studies have shown it to be promising with the ability to employ coarse grids resulting in significant computational savings compared to the TFM approach [41]. In the DDPM, the conservation equations of the continuous phase are also modified by taking into account the volume fraction of the dispersed phase. Therefore in the absence of mass exchange between the phases, the continuity equation for the Eulerian phase when DDPM is enabled takes the same form as the continuity equation for the TFM (cf. (1)). When only a single Eulerian phase is present such as the bubble column studied in this paper, the momentum equation gets transformed as
(9)∂∂tαpρpv→p+∇·αpρpv→pv→p=-αp∇p+∇·αpμp∇v→p+∇v→pT+αpρpg→+KDPMv→DPM-v→p+SDPM,explicit.
Volume fractions and the discrete phase velocities are obtained from the Lagrangian field solution.


## 2.2.4. The Two-Fluid Model with Population Balance (TFM + PB)

The population balance (PB) model describes the evolution of the bubble size distribution as a result of growth, aggregation, and breakage processes. By introducing the concept of a number density functionf(x,v,t) balance equations are solved describing its evolution [9]. The general governing equation for population balance can be written as
(10)∂∂tfx,ν,t+∇·v→bx,ν,tfx,ν,t=Sx,ν,t.
In (10), S(x,v,t) accounts for the bubble birth rate, death rate, and breakup. The drag force is computed in terms of the local mean bubble density. The bubble sizes in the simulation were divided into six discrete bins (or classes) with bubble diameters ranging between 1 mm and 10 mm with a ratio exponent of 2. In the PB solution method employed in this study, all the bubbles were assumed to move at the ensemble averaged gas-phase velocity vb obtained from the solution of the TFM circumventing the need to solve continuity and momentum equations for each of the individual discrete bins. The Luo-model [9] was employed for the aggregation kernel and the Hagesather model formulation [9] was employed for the breakage with its frequency determined by the Luo-model [9].


## 2.3. Radiation Modeling

The finite volume (FV) radiation model in ANSYS FLUENT (which is a conservative variant of the discrete ordinates method) was employed in this study to model the radiative transfer. In this study, the number of divisions in the polar (Nθ) and azimuthal angles (Nφ) for the angular intensity vectors were both set to 4, corresponding to a total of 128 directions (8NθNφ) in which the RTE was being solved. Any further refinement in the angular discretization did not impact our results.If “I” represents the directional intensity, k the absorption coefficient (which is due to water alone), σ the scattering coefficient (due to the secondary phase, i.e., the bubbles alone), Ib the black body emissive power (with the water temperature set at 298 K), and Ф the scattering phase function (assumed to be forward scattering for the bubbles), then the differential equation governing the FV method can be written for each wavelength region (λ) and medium refractive index “n” as [42]
(11)∇·Iλr→,s^s^=-(kλ+σλ)Iλ(r→,s^)+ηkλn2Ibλ(r→,s^)+σλ4π∫04πIλr→,s^Φs^,s^′dΩ′.
In (11), “η” corresponds to the fraction of the total blackbody emission in the wavelength region under consideration [42].In this multiphase study, the absorption and scattering coefficients were multiplied by their corresponding phase volume fractions (αp,i) and employed in a modified version of the RTE for each wavelength region (λ) as
(12)∇·Iλr→,s^s^=-∑i=12αp,i(kλ,i+σλ,i)Iλ(r→,s^)+∑i=12ηiαp,ikλ,in2Ibλ,i(r→,s^)+∑i=12αp,iσλ,i4π∫04πIλr→,s^Φs^,s^′dΩ′.
Therefore, in this two-phase study “effective” absorption and scattering coefficients of the medium were computed through a weighted averaging of the individual phase radiative properties by their corresponding volume fractions. The bubble columns simulated in this study were assumed to be isothermal with thermal equilibrium prevailing between the two phases. The radiative transfer variables reported in this study, the incident radiative flux (qin), and the incident radiation (G) were obtained by integrating the angular intensities over all directions as
(13)qin,λr=∑λ∫4πIλr,s^s^dΩ,(14)Gr=∑λ∫4πIr,s^dΩ.
The incident radiation is a measure of the amount of radiation received by a control volume from all directions (e.g., this might represent the dosage received by algae in a photobioreactor). The volume averaged incident in the reactor was employed in this study to garner insights into the effects of different operating conditions such as air flow rate and bubble size on the light distribution. The scattering coefficient of the bubbles in (12) was computed in the DDPM model as
(15)σs=LimitV→0⁡∑N1-εsn1-fsnAsnV.
In (15), the summation is over N particles within the control volume V, εsn is the bubble emissivity (set to zero), Asn is the projected area of the nth particle, and fsn is the scattering factor associated with the nth particle. Near-field approximations have shown that, in the limitation case of spheres much larger than the wavelength of radiation (such as the bubbles examined in this study) that are embedded in an absorbing medium, the extinction efficiency approaches unity as the effects of diffraction can be neglected [43–45]. Therefore fsn in (15) was set to 0 to ensure that the light incident upon a bubble was scattered with a scattering efficiency of 1.While stand-alone radiation codes for accurately characterizing the scattering properties of spherical particles and bubbles are available [24], in multiphysics simulation choices have to be made to assess first the fidelity required in radiative transfer modeling and second the relative importance of gas and dispersed phase radiative property resolution. Since the computational cost of a wavelength dependent radiative transfer calculation in a transient multiphase simulation is already high, these assessments can help optimize the calculation times. For instance, radiative transfer calculations may be carried out once in several time-steps in the calculations if it was established that the radiative transfer is not significantly impacted by small changes to the phase volume fractions that occur over the hydrodynamic time-scales. The second issue of the relative importance of gas and dispersed phase radiative properties will largely be determined by the reactor configurations and the flow conditions within the reactor. For instance, through coal combustion simulations, Krishnamoorthy et al. [2] showed that the fidelity of the gas-phase radiative property models had little impact on the radiation distributions in a lab-scale reactor whereas differences between gray and nongray model predictions were observed in a full-scale boiler [46]. Similarly, through steam-gasification studies on directly irradiated reactor, von Zedtwitz et al. [3] determined that the radiation absorption by particles was three orders of magnitude higher than that absorbed by the gas-phase.The TFM-PB model provides the spatial variations in the probability distribution function (pdf) of the different bubble classes (bin). The scattering coefficient of each bin in the TFM-PB model calculations was calculated from the bin pdf, projected surface area (As), and the scattering efficiency (assumed to be unity) as
(16)σP,bin=volumefraction∑allclassesbubbleclasspdf×1×Ap,binVp.
The total scattering coefficient (σ) at any location for use in the RTE (cf. (12)) was then computed by a summation of the bin scattering coefficients at that location.The TFM and the mixture models do not account for bubble coalescence and breakup. Therefore, a single “class” is employed in the calculations. However, both the models provide volume fraction information as output. Therefore, in the TFM and mixture model calculations, the number density was first determined from the volume fraction information using(17)Numberdensity=BubblevolumefractionVp,
where Vp is the volume of the bubble assuming a spherical bubble shape and the bubble diameter was assumed to be 4 mm. The scattering coefficient was then determined as
(18)σs=numberdensity×As.
Equations (16)–(18) were implemented as user-defined functions (UDFs) and employed in the simulations. Flow charts depicting the exchange of variables between ANSYS FLUENT and the UDFs are shown in Figure 2. The boundary condition associated with the semitransparent radiator surface was specified as
(19)Qradiator=Qirrad.
In (19)  Qirrad is the beam irradiation flux corresponding to each wavelength region as specified in Table 1. The other three side faces of the bubble columns were set as semitransparent surfaces at a temperature of 298 K.Flow diagram depicting the exchange of variables between the CFD calculations and the user-defined functions (UDFs) for computing the scattering coefficients.
(a)
(b)


## 3. Results and Discussion

### 3.1. Gas-Phase Volume Fractions

The time-averaged gas-phase volume fractions predicted by the different multiphase models along an axial midplane are compared in Figure3. The MM and DDPM predictions show a cylindrical column of bubbles whereas the TFM and TFM with PB models capture the spreading of the bubble plume that was observed experimentally [35]. Furthermore, the TFM + PB simulation that accounts for bubble coalescence and breakup shows higher gas volume fractions than the TFM simulation near the inlet and lower gas volume fractions in the upper sections of the column. Therefore at the experimental conditions explored in this study, qualitative differences in the gas volume fractions clearly manifest themselves as a result of the different approximations made in modeling the phases as outlined in Section 2.2.Figure 3
Contours of the bubble volume fraction predicted by the different models along an axial midplane in the reactor.


### 3.2. Bubble Scattering Coefficient

Figure4 shows the distribution of the bubble scattering coefficient from the various models that were computed employing the UDFs. First, the results show that the spatial variations in scattering coefficient are analogous to the gas-phase volume fraction predictions. The MM and DDPM models show a high central core scattering coefficient that is reflective of the high bubble volume fractions at the center noticeable in Figure 3. The scattering coefficients predicted by the TFM + PB model are slightly lower than those predicted by the TFM in the upper section of the column and may be attributed to the reduced volume fraction of the bubbles observed in Figure 3. The PB modeling approach was carried out by employing the discrete method to solve the PB equations by dividing the bubble sizes into six classes [9]. The scattering coefficient distribution associated with each bin in the 3D column at steady state is shown in Figure 5. The different bin sizes are also indicated in Figure 5 with the minimum bin size set at 0.001 m. The scattering coefficient in the TFM + PB modeling approach is a function of the gas-volume fraction, bubble size as well as the probability of that bubble size occurring at that spatial location. Consequently, the highest value of the scattering coefficient is observed close to the inlet and corresponds to an injection bubble diameter of 4 mm before any coalescence or breakup occurs. However, some contribution to the scattering coefficient is also observed from larger sized bubbles as a result of coalescence whereas the scattering coefficient associated with the smallest bubble diameter is small indicating a lack of significant breakage at these experimental gas flow velocities. The sensitivities of the extent of coalescence and breakage to the gas inlet velocities are explored in the next section. Therefore, the observed variations between the TFM and TFM + PB model scattering coefficient predictions observed in Figure 4 may be attributed to these coalescence and breakage phenomena. Table 2 summarizes the results from radiative transfer calculations in the bubble column at inlet air velocities of 7.84 cm/s and inlet bubble diameters of 4 mm from employing different multiphase models. Results from no bubble participation (through attenuation) in the radiative transfer are also provided as a reference. The results of interest in the radiative transfer analysis are the incident radiation (cf. (14)) which is a measure of the average dosage received by the algae and the surface incident radiative flux which would be measured by a surface sensor located outside the reactor. In the absence of any light attenuation by the bubbles, the wall opposite to the radiator receives more radiation than the side walls (the two vertical walls in Figure 1 that are adjacent and perpendicular to the wall containing the light source) due to the view angle. However, for the bubble column and operational conditions examined in this study the variations in the radiative transfer predictions among the multiphase models were not significant.Table 2
Summary of radiative transfer calculations at inlet air velocities of 7.84 cm/s and inlet bubble diameters of 4 mm.
Model
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
No bubble participation
202
62
49
DDPM
205
45
54
MM
205
49
53
TFM
207
44
53
TFM + PB
206
50
52Figure 4
Contours of the scattering coefficient (in m−1) predicted by the different models along an axial midplane in the reactor.Figure 5
Scattering coefficients (in m−1) of the different sized bubbles in the TFM + population balance simulations.


### 3.3. Effect of Inlet Gas Velocity

Previous studies have shown that as the gas inlet velocity is increased the bubble size becomes larger and the Sauter mean diameter increases within the turbulent churn regimes encountered in a cylindrical bubble column of 16.2 cm diameter and gas inlet velocities of 8 cm/s to 30 cm/s [22]. In order to investigate the impact of the changes to the flow rate and the bubble size on the radiative transfer, two additional simulations employing the TFM + PB modeling approach and an initial bubble diameter of 4 mm were carried out for gas inlet velocities of 3.92 cm/s and 15.68 cm/s which correspond to half and twice our original inlet velocity values, respectively. The volume fraction and scattering coefficients at the different inlet velocities are compared in Figure 6. The volume fraction of the bubbles near the inlet increases with the increase in inlet velocities and the bubble plume is also observed to spread more with an increase in velocity. Figure 6 shows the contribution to the scattering coefficient by the 10 mm diameter bubbles and the 4 mm diameter bubbles at the different inlet velocities. While the 4 mm diameter bubbles contribute significantly to the scattering coefficient, the contribution from the larger sized bubbles to the scattering coefficient is not insignificant at higher inlet gas velocities. Figure 7 shows the volume averaged bubble class pdf for the different bubble sizes in the bubble column at different air inlet velocities. With an increase in velocity, there is an increase in the fraction of the large sized bubbles in agreement with the results reported in Chen et al. [22]. At the highest velocities explored, the largest sized bubbles are found almost in the same fraction as the 4 mm bubbles. However, since the scattering coefficient is directly proportional to the bubble class fraction and inversely proportional to the volume, the contribution of the larger bubbles to the scattering coefficient is still lower than that of the 4 mm sized bubbles as previously noted (cf. Figure 6). Table 3 summarizes the sensitivities of the radiative transfer predictions to these flow rate variations. As illustrated in Figure 6, the attenuation of radiation is still dominated by the 4 mm bubbles across all flow rates. Consequently, the impacts of the flow rate variations on the volume averaged incident radiation, in this reactor configuration, at these operating conditions, are minimal. However, upon extending the radiative transfer modeling methodology proposed in this study to a stirred-tank reactor configuration we were able to discern a 10% variation in the volume averaged incident radiation as a function of air flow rates, the results of which will be published in a later study.Table 3
Sensitivity of radiative transfer predictions to the inlet gas velocities in the simulations carried out with the TFM + population balance approach and initial bubble diameter of 4 mm.
Model [air inlet velocity]
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
TFM + population balance (3.92 cm/s)
206
55
50
TFM + population balance (7.84 cm/s)
206
50
52
TFM + population balance (15.68 cm/s)
207
48
50Variations in the scattering coefficients (in m−1) of the different sized bubbles as a function of inlet gas velocities in the TFM + population balance simulations: (a) 4 mm bubbles; (b) 10 mm bubbles.
(a)
(b)Figure 7
Variations in the volume averaged bubble probability density functions as a function of the inlet gas velocity for different bubble sizes (TFM + PB model).


### 3.4. Sensitivity to Initial Bubble Diameter

Our calculations so far have assumed an initial bubble size of 4 mm. The sensitivities of the radiative transfer predictions to the initial bubble size were investigated next to take into account variations in the air sparger make and design. Simulations for bubble diameters of 10 mm and 1 mm were carried out employing the TFM approach. Air volumetric flow rate of1.1×10-4 m3/s was maintained for all three bubble diameters. Figure 8 shows the volume fraction and scattering coefficient predictions from the different simulations. The smaller sized bubbles travel higher within the reactor (when the flow rate is maintained constant) and contribute to an increase in gas volume fractions at higher elevations. Furthermore, as noted previously, since the scattering coefficient is inversely proportional to the bubble diameter it increases with a decrease in bubble diameter. The results from employing different initial bubble diameters in the simulations are summarized in Table 4. A decrease in initial bubble diameter contributes to a decrease in the radiative flux incident upon the wall opposite to the radiator (due to stronger attenuation of light), an increase in volume averaged incident radiation resulting from increased scattering and redistribution of light due to the increase in bubble surface area (since the air volumetric flow rate is fixed). The radiative transfer predictions therefore seem to be more strongly impacted by the initial bubble size than the air flow rates or the multiphase model for the flow regimes encountered within this reactor.Table 4
Sensitivity of radiative transfer predictions to the initial bubble diameter.
Model (bubble diameter)
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
TFM (1 mm)
215
33
55
TFM (4 mm)
207
44
53
TFM (10 mm)
205
56
50Variations in the bubble volume fractions (a) and scattering coefficients (in m−1) (b) as a function of the bubble diameter in the TFM simulations.
(a)
(b)


## 3.1. Gas-Phase Volume Fractions

The time-averaged gas-phase volume fractions predicted by the different multiphase models along an axial midplane are compared in Figure3. The MM and DDPM predictions show a cylindrical column of bubbles whereas the TFM and TFM with PB models capture the spreading of the bubble plume that was observed experimentally [35]. Furthermore, the TFM + PB simulation that accounts for bubble coalescence and breakup shows higher gas volume fractions than the TFM simulation near the inlet and lower gas volume fractions in the upper sections of the column. Therefore at the experimental conditions explored in this study, qualitative differences in the gas volume fractions clearly manifest themselves as a result of the different approximations made in modeling the phases as outlined in Section 2.2.Figure 3
Contours of the bubble volume fraction predicted by the different models along an axial midplane in the reactor.


## 3.2. Bubble Scattering Coefficient

Figure4 shows the distribution of the bubble scattering coefficient from the various models that were computed employing the UDFs. First, the results show that the spatial variations in scattering coefficient are analogous to the gas-phase volume fraction predictions. The MM and DDPM models show a high central core scattering coefficient that is reflective of the high bubble volume fractions at the center noticeable in Figure 3. The scattering coefficients predicted by the TFM + PB model are slightly lower than those predicted by the TFM in the upper section of the column and may be attributed to the reduced volume fraction of the bubbles observed in Figure 3. The PB modeling approach was carried out by employing the discrete method to solve the PB equations by dividing the bubble sizes into six classes [9]. The scattering coefficient distribution associated with each bin in the 3D column at steady state is shown in Figure 5. The different bin sizes are also indicated in Figure 5 with the minimum bin size set at 0.001 m. The scattering coefficient in the TFM + PB modeling approach is a function of the gas-volume fraction, bubble size as well as the probability of that bubble size occurring at that spatial location. Consequently, the highest value of the scattering coefficient is observed close to the inlet and corresponds to an injection bubble diameter of 4 mm before any coalescence or breakup occurs. However, some contribution to the scattering coefficient is also observed from larger sized bubbles as a result of coalescence whereas the scattering coefficient associated with the smallest bubble diameter is small indicating a lack of significant breakage at these experimental gas flow velocities. The sensitivities of the extent of coalescence and breakage to the gas inlet velocities are explored in the next section. Therefore, the observed variations between the TFM and TFM + PB model scattering coefficient predictions observed in Figure 4 may be attributed to these coalescence and breakage phenomena. Table 2 summarizes the results from radiative transfer calculations in the bubble column at inlet air velocities of 7.84 cm/s and inlet bubble diameters of 4 mm from employing different multiphase models. Results from no bubble participation (through attenuation) in the radiative transfer are also provided as a reference. The results of interest in the radiative transfer analysis are the incident radiation (cf. (14)) which is a measure of the average dosage received by the algae and the surface incident radiative flux which would be measured by a surface sensor located outside the reactor. In the absence of any light attenuation by the bubbles, the wall opposite to the radiator receives more radiation than the side walls (the two vertical walls in Figure 1 that are adjacent and perpendicular to the wall containing the light source) due to the view angle. However, for the bubble column and operational conditions examined in this study the variations in the radiative transfer predictions among the multiphase models were not significant.Table 2
Summary of radiative transfer calculations at inlet air velocities of 7.84 cm/s and inlet bubble diameters of 4 mm.
Model
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
No bubble participation
202
62
49
DDPM
205
45
54
MM
205
49
53
TFM
207
44
53
TFM + PB
206
50
52Figure 4
Contours of the scattering coefficient (in m−1) predicted by the different models along an axial midplane in the reactor.Figure 5
Scattering coefficients (in m−1) of the different sized bubbles in the TFM + population balance simulations.


## 3.3. Effect of Inlet Gas Velocity

Previous studies have shown that as the gas inlet velocity is increased the bubble size becomes larger and the Sauter mean diameter increases within the turbulent churn regimes encountered in a cylindrical bubble column of 16.2 cm diameter and gas inlet velocities of 8 cm/s to 30 cm/s [22]. In order to investigate the impact of the changes to the flow rate and the bubble size on the radiative transfer, two additional simulations employing the TFM + PB modeling approach and an initial bubble diameter of 4 mm were carried out for gas inlet velocities of 3.92 cm/s and 15.68 cm/s which correspond to half and twice our original inlet velocity values, respectively. The volume fraction and scattering coefficients at the different inlet velocities are compared in Figure 6. The volume fraction of the bubbles near the inlet increases with the increase in inlet velocities and the bubble plume is also observed to spread more with an increase in velocity. Figure 6 shows the contribution to the scattering coefficient by the 10 mm diameter bubbles and the 4 mm diameter bubbles at the different inlet velocities. While the 4 mm diameter bubbles contribute significantly to the scattering coefficient, the contribution from the larger sized bubbles to the scattering coefficient is not insignificant at higher inlet gas velocities. Figure 7 shows the volume averaged bubble class pdf for the different bubble sizes in the bubble column at different air inlet velocities. With an increase in velocity, there is an increase in the fraction of the large sized bubbles in agreement with the results reported in Chen et al. [22]. At the highest velocities explored, the largest sized bubbles are found almost in the same fraction as the 4 mm bubbles. However, since the scattering coefficient is directly proportional to the bubble class fraction and inversely proportional to the volume, the contribution of the larger bubbles to the scattering coefficient is still lower than that of the 4 mm sized bubbles as previously noted (cf. Figure 6). Table 3 summarizes the sensitivities of the radiative transfer predictions to these flow rate variations. As illustrated in Figure 6, the attenuation of radiation is still dominated by the 4 mm bubbles across all flow rates. Consequently, the impacts of the flow rate variations on the volume averaged incident radiation, in this reactor configuration, at these operating conditions, are minimal. However, upon extending the radiative transfer modeling methodology proposed in this study to a stirred-tank reactor configuration we were able to discern a 10% variation in the volume averaged incident radiation as a function of air flow rates, the results of which will be published in a later study.Table 3
Sensitivity of radiative transfer predictions to the inlet gas velocities in the simulations carried out with the TFM + population balance approach and initial bubble diameter of 4 mm.
Model [air inlet velocity]
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
TFM + population balance (3.92 cm/s)
206
55
50
TFM + population balance (7.84 cm/s)
206
50
52
TFM + population balance (15.68 cm/s)
207
48
50Variations in the scattering coefficients (in m−1) of the different sized bubbles as a function of inlet gas velocities in the TFM + population balance simulations: (a) 4 mm bubbles; (b) 10 mm bubbles.
(a)
(b)Figure 7
Variations in the volume averaged bubble probability density functions as a function of the inlet gas velocity for different bubble sizes (TFM + PB model).


## 3.4. Sensitivity to Initial Bubble Diameter

Our calculations so far have assumed an initial bubble size of 4 mm. The sensitivities of the radiative transfer predictions to the initial bubble size were investigated next to take into account variations in the air sparger make and design. Simulations for bubble diameters of 10 mm and 1 mm were carried out employing the TFM approach. Air volumetric flow rate of1.1×10-4 m3/s was maintained for all three bubble diameters. Figure 8 shows the volume fraction and scattering coefficient predictions from the different simulations. The smaller sized bubbles travel higher within the reactor (when the flow rate is maintained constant) and contribute to an increase in gas volume fractions at higher elevations. Furthermore, as noted previously, since the scattering coefficient is inversely proportional to the bubble diameter it increases with a decrease in bubble diameter. The results from employing different initial bubble diameters in the simulations are summarized in Table 4. A decrease in initial bubble diameter contributes to a decrease in the radiative flux incident upon the wall opposite to the radiator (due to stronger attenuation of light), an increase in volume averaged incident radiation resulting from increased scattering and redistribution of light due to the increase in bubble surface area (since the air volumetric flow rate is fixed). The radiative transfer predictions therefore seem to be more strongly impacted by the initial bubble size than the air flow rates or the multiphase model for the flow regimes encountered within this reactor.Table 4
Sensitivity of radiative transfer predictions to the initial bubble diameter.
Model (bubble diameter)
Volume averaged incident radiation (W/m2)
Incident radiative flux on wall opposite to the radiator (W/m2)
Incident radiative flux on the side walls (W/m2)
TFM (1 mm)
215
33
55
TFM (4 mm)
207
44
53
TFM (10 mm)
205
56
50Variations in the bubble volume fractions (a) and scattering coefficients (in m−1) (b) as a function of the bubble diameter in the TFM simulations.
(a)
(b)


## 4. Conclusions

Currently, there is a lack of rigorous coupling between radiative transfer and hydrodynamics in simulations of nondilute multiphase flows. This is attributable to the high computational cost, complexities arising from the need to incorporate additional equations and terms within current RTE solution frameworks of CFD codes, and incompatibilities in the output and input variables that are employed by the radiation and multiphase models. This study demonstrates for the first time a methodology for performing wavelength dependent, multiphase radiative transfer calculations in CFD simulations of gas-liquid flows. The multiphase hydrodynamics within a bubble column PBR was simulated employing four different multiphase modeling approaches. By approximating the bubbles as spheres, add-on modules were developed and employed to compute the scattering coefficients from the spatial variations in the number density, volume fraction, and bubble bin fraction information that was output from the different multiphase models. The impacts of the multiphase hydrodynamic modeling strategy, initial bubble diameter, and operational parameters on the radiation distribution patterns within the reactor were also examined. Based on the results from this study the following conclusions may be drawn.(1)
Only the TFM and TFM with population balance multiphase models that solved continuity and momentum equations for the individual phases were able to capture the bubble plume spreading behavior that was observed experimentally. Consequently, the scattering coefficients predicted from the various models were significantly different from one another and were very similar to their gas-phase volume fraction distributions.(2)
The differences in scattering coefficient predictions among the models while having some impact on the wall incident radiative fluxes had a minimal effect on the volume averaged incident radiation within the reactor.(3)
Employing a population balance modeling approach, we were able to discern changes to the bubble size distributions within the reactor that were consistent with previous findings; that is, the fraction of the larger sized bubbles increased with an increase in air velocities. This resulted in a corresponding increase in the scattering coefficient of the larger sized bubbles. However, over the flow rates that were investigated, the initial bubble diameter of 4 mm was the most significant contributor to the scattering coefficient. Consequently, the sensitivity of the volume averaged incident radiation within the reactor to the flow rate variations was minimal.(4)
Investigations of the sensitivities of the radiative transfer predictions to the initial bubble size revealed that smaller sized bubbles travelled higher within the reactor (when the flow rate was maintained constant) and contributed to an increase in gas volume fractions at higher elevations. This resulting increase in the scattering coefficients and bubble surface area with decreasing bubble size translated to a greater attenuation of light and a slight increase in the volume averaged incident radiation within the reactor.Therefore, for the modest sized bubble column (0.01 m3) and flow regimes investigated in this study, the radiative transfer predictions were seen to be more strongly impacted by the initial bubble size than the air flow rates or the multiphase model. The demonstrated methodology for radiation coupling in these multiphase CFD simulations is currently being extended to study the algae growth in a stirred-tank reactor.



---

*Source: 793238-2014-11-06.xml*