Computer Simulations of Photocatalytic Reactors

Abstract

Photocatalysis has been considered future technology for green energy conversion and environmental purification, including carbon dioxide reduction, water splitting, air/water treatment, and antimicrobial purposes. Although various photocatalysts with high activity and stability have already been found, the commercialization of photocatalytic processes seems to be slow; it is thought that the difficulty in scaling up photocatalytic processes might be responsible. Research on the design of photocatalytic reactors using computer simulations has been recently intensive. The computer simulations involve various methods of hydrodynamics, radiation, and mass transport analysis, including the Monte Carlo method, the approximation approach–P1 model, and computational fluid dynamics as a complex simulation tool. This review presents all of these models, which might be efficiently used for the scaling-up of photocatalytic reactors. The challenging aspects and perspectives of computer simulation are also addressed for the future development of applied photocatalysis.

1. Introduction

Semiconductor photocatalysis, especially based on TiO₂, has been known for over 40 years as a promising strategy to achieve aims such as environmental remediation (e.g., water treatment and air purification) and renewable energy (e.g., hydrogen production and carbon dioxide conversion to methane and light hydrocarbons) [1,2,3]. There have been many studies on photocatalysis, which can be divided into three main fields: development of new photocatalysts, fundamental studies on photocatalysis, and design of photocatalytic reactors [4,5,6,7,8,9].

Despite the constantly growing number of research studies in photocatalysis and some successful applications [10,11,12,13], a gap remains between research and commercialization. One of the main problems is the scaling-up of photocatalytic reactors [14]. Although some semi-pilot photoreactors exhibit similar or even better performance than laboratory ones, e.g., air-sparged hydrocyclone (ASH), due to irradiation of a thin and highly aerated liquid layer [8,15,16,17,18], the broad application of photocatalytic reactions needs a more systematic approach. This process requires the development of mathematical models such as the radiation absorption–scattering model, the radiation emission model, the kinetic model, and the fluid-dynamic model [19,20,21]. These models are the subdomains of a closely linked, complex numerical system based on differential equations [20]. The main point of the modeling procedure is the solution of the radiation transport equation (RTE). This equation depends on the specific geometry of the photocatalytic reactor and the type of light source [22,23]. The RTE solution requires the application of a proper mathematical model, such as the statistical Monte Carlo method [24,25,26,27] or approximation approach–P1 model [28,29,30]. Another important procedure is the discrete ordination method (DOM), which is recognized as an accurate and flexible procedure [19,22,31,32,33,34,35,36,37,38]. As a consequence of the increase in computing power, more research on the modeling of photocatalytic reactors also involves the concept of computational fluid dynamics (CFD). In general, the CFD method is widely used in the development of catalytic reactors and analysis of the systems involving fluid flows and heat transfer to understand the interactions of the catalyst with the surrounding reactive flow field in catalytic reactors (gas–solid, liquid–solid, and gas–liquid–solid) [39,40]. In the case of photocatalytic reactors, the CFD-aided modeling procedure also introduces the simulation of fluid flows (single- and multi-phase) and solves the RTE [32,37,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57]. In the case of slurry photoreactor systems, there is the concept of local volumetric rate of energy absorption (LVREA), defined as the energy required by photons absorbed per time and volume inside the photoreactor. The photocatalytic reaction rate is the function of the LVREA [23,41]. The solution of the LVREA equation requires the analytical solution of RTE [58]. The commonly applied methods of RTE solution have already been mentioned above.

This review analyzes the present state of knowledge on computing methods in the design of photocatalytic reactors. The continuous development of computer-aided methods seems to be a necessary step to increase the number of successful large-scale applications of technologies based on photocatalytic processes. The aim of this work is to indicate the currently considered numerical methods and corresponding software solutions dedicated to the design of photocatalytic reactors. Notably, the most important advantages of computer-aided design are the substantial reduction in time and costs (less experimental trials) and the ability to study systems in which experiments are difficult or impossible to carry out.

2. Numerical Approaches Supporting Scaling-Up of Photocatalytic Reactors

2.1. Scaling-Up Methodology of Photocatalytic Reactors

The traditional approach of the scaling-up procedure is based mostly on the empirical methodology, i.e., it starts from the laboratory set-ups, and then gradually, the photoreactor dimensions increase until it reaches the destination size specified for the particular practical application. This approach is, unfortunately, time-consuming and requires significant investments. Therefore, there is a consensus that the photoreactor design based on the mathematical models and the increased computing power of hardware might be useful in supporting the conventional scaling-up methodology [41,58,59,60,61]. Another approach for semi- and pilot-plan photoreactors has been applied, i.e., the use of other types of reactors as photoreactors after the addition of irradiation sources, e.g., an immersion of lamp(s) inside the reactor or an application of outside irradiation (both artificial and natural) [14]. For example, ASH, an efficient gas–liquid contactor used successfully for aeration, SO₂ removal, flotation, and stripping [17,62,63], has shown high efficiency for oxidative decomposition of organic compounds when equipped with an immersed mercury lamp for both photolysis (UV and UV/H₂O₂) and photocatalytic (UV/TiO₂) reactions [8,15,16]. Although this approach might result in the discovery of an efficient photocatalytic reactor, it is random, and more studies are needed for detailed characterization and further scaling-up of such photoreactors.

Satuf et al. proposed a methodology for scaling-up slurry photocatalytic reactors (the example of 4-chlorophenol photodegradation) [58]. A similar conception was formulated by the same research group for the photocatalytic reactor with TiO₂ immobilized on the reactor walls for the removal of air pollutants [61]. Figure 1 shows the flow chart illustrating this scaling-up methodology, which consists of two main parts. The first (upper part of the scheme) relates to the derivation of an intrinsic kinetic expression that is independent of the present state of reactor configuration, radiation source, and experimental conditions, and corresponds to photocatalytic degradation of model compound in a laboratory scale reactor. The intrinsic kinetics is based on a detailed photocatalytic reaction scheme and introduces the factors resulting from irradiation conditions, the amount of photocatalyst, and initial compound concentration [58]. Furthermore, it would be stated as the rule that a carefully considered laboratory photoreactor could render a reaction kinetics expression without time dependence and with validity for the whole reactor, as long as it is possible by the provided method to safely extrapolate the results to more complex conditions. The estimation of kinetics parameters should be performed by introducing a nonlinear optimization program to adjust model predictions to the experimental data [61].

In the radiation model, the key point is to determine the effect of radiation absorption on the photocatalytic reaction rate. The simultaneous existence of absorption and scattering is important for radiation modeling (the RTE method). The mass balances of the photodegraded compound and the intermediate reaction products are considered to show the theoretical evolution of the involved species. They should include the data from the kinetic and radiation models [58]. In the second step (lower part of Figure 1), changing the process scale requires the usage of the same kinetic model for the reaction that was reported in the laboratory experiments. It is important to consider that the local volumetric rate of photon absorption (LVRPA) must be calculated from a different radiation balance, estimated for the larger photoreactor. Therefore, RTE should be applied to the new reactor to predict the LVRPA in a reaction space. The resulting expression should be included in the mass balance and, by the conduction of simulation procedure, to form the basis for calculating the reactor output variables. Subsequently, the validation of obtained results from the modeling procedure is necessary by comparing the experiments conducted in a larger (e.g., pilot) scale photoreactor [61].

2.2. Hydrodynamics Modeling

The complex (integrated) photoreactor model should also consider the role of hydrodynamic conditions expressed by the momentum balance. For photocatalytic reactors, two types of conditions should be considered: (i) free-flow through the reaction system and (ii) the flow through the porous system [21,23,37]. The Navier–Stokes equations are used for the description of momentum balance at steady-state conditions. For the constant density and viscosity of the fluid, the continuity equation and the balance equation for momentum can be formulated [23,33]. The computational fluid dynamics method (CFD) is increasingly used to estimate the hydrodynamics and mass transfer to predict the velocity and concentration fields in photoreactors. Just as it underlies the CFD method, the Navier–Stokes equations as a set of differential equations describe the motion of viscous fluids [39]. These equations should be solved first before procuring the solution to the RTE [37].

There are different hydrodynamic models depending on the flow character (laminar/turbulent). If the fluid is assumed to be Newtonian, incompressible, non-reactive, with constant physical properties, and under laminar steady-state flow, it is possible to express the following equations [47]:

The continuity equation:

$$\nabla \times \left(\mathrm{U}\right)=0$$

(1)

The equation of momentum conservation:

$$\nabla \times \left(\mathsf{\rho }\mathrm{UU}\right)=-\nabla \mathrm{P}-\nabla \times \mathsf{\tau }$$

(2)

With the stress tensor:

$$\mathsf{\tau }=\mathsf{\mu }\left(\nabla \mathrm{U}+\nabla {\mathrm{U}}^{\mathrm{T}}\right)-\frac{2}{3}\mathsf{\mu }\nabla \times \mathrm{UI}$$

(3)

The equation of species conservation:

$$\nabla \times \left({\mathsf{\rho }\mathrm{Um}}_{\mathrm{i}}\right)=-\nabla \times {\mathrm{J}}_{\mathrm{i}}, \mathrm{i}=1, 2, \dots , \mathrm{N}-1$$

(4)

The diffusive flux of species J_i can be estimated using Fick’s first law of diffusion:

$${\mathrm{J}}_{\mathrm{i}}=-{\mathrm{D}}_{\mathrm{m}}\nabla \left({\mathsf{\rho }\mathrm{m}}_{\mathrm{i}}\right)$$

(5)

In the above-mentioned equations, ρ is density, U is velocity, P is pressure, τ is viscous stress tensor, μ is molecular viscosity, I is unit tensor, m_i is the mass fraction of species i, N is the total number of species, and D_m is the molecular diffusivity of species i in the mixture [47]. This model has been successfully applied in the CFD simulations of photocatalytic reactors [43,45,46,49,50,53,64,65,66,67,68]. This model is preferred for gas-phase photocatalytic reactions.

Considering the case of turbulent flow, under the same assumptions that were stated for laminar flow, the Reynolds-averaged Navier–Stokes (RANS) turbulence modeling approach requires to solve the following time-average equations [33]:

Mass conservation equation:

$$\nabla \times \left(\overline{\mathrm{U}}\right)=0$$

(6)

Momentum conservation equation:

$$\nabla \times \left(\mathsf{\rho }\overline{\mathrm{UU}}+\mathsf{\rho }\overline{\mathrm{uu}}\right)=-\nabla \overline{\mathrm{P}}-\nabla \times \overline{\mathsf{\tau }}$$

(7)

Species conservation equation:

$$\nabla \times \left(\mathsf{\rho }\overline{\mathrm{U}}\overline{{\mathrm{m}}_{\mathrm{i}}}+\mathsf{\rho }\overline{{\mathrm{um}}_{\mathrm{i}}^{\prime }}\right)=-\nabla \times \overline{{\mathrm{J}}_{\mathrm{i}}} \mathrm{i}=1, 2, \dots , \mathrm{N}-1$$

(8)

The overbar indicates a time-average value, and u and m′_i are fluctuating flow velocity and mass fraction of species i. The specification of the apparent stress gradient ($\mathsf{\rho }\overline{\mathrm{UU}})$ is related to the adapted turbulence model [33]. Four hydrodynamic models are most often considered for modeling photoreactors: (a) standard k–ε model (S k–ε) [36,38,45,69,70], (b) the realizable k–ε model (R k–ε) [45,69], (c) the Reynolds stress model (RSM) [45,69], and (d) the low Reynolds number k–ε model (AKN—developed by Abe, Kondoh, and Nagano [21,45,71]). Duran et al. performed the CFD for a single-phase flow mass transfer prediction in annular reactors using different hydrodynamic models. Figure 2 shows fluid velocity magnitude differences between U- and L-type annular reactors [45].

For single-phase photoreactor systems, classical Navier–Stokes equations are efficient to calculate the hydrodynamic properties, but in the case of multi-phase flows, two main models are considered [72]: (a) Eulerian–Eulerian (E–E) and (b) Eulerian–Lagrangian (E–L) models. In the E–E model, the multiple phases are described as interpenetrating continua for which equations representing the conservation mass and momentum are solved. This model should be preferred due to its ability to simulate large-scale systems without excessive computational requirements. For the E–L approach, the trajectories of dispersed phase particles are simulated by solving an equation motion for each dispersed phase particle. The motion of the continuous phase is simulated by a conventional Eulerian approach [72]. Pareek et al. performed the simulation with an integrated model for a pilot-scale annular bubble column photocatalytic reactor using the CFD method [42]. They applied the conventional Eulerian model to simulate three-phase flow in the photoreactor. The model was verified using experimental data for photodegradation of Bayer liquor, and the adequacy of the adapted model was proven. Figure 3 shows the results of the simulation of this three-phase modeling study.

Another point is the aspect of stirring as one of the most important factors that is still underrated in the available computational studies. The performance of especially immobilized photocatalytic reactors is limited by the rate of mass transfer of the compounds to the photocatalyst layer. For example, Duran et al. [21] applied mass transfer promoters that increase the turbulence effect, such as repeated ribs on the internal wall of the reactor and static delta wing mixers placed in the middle of the reactor. A CFD simulation showed that repeated ribs improve the turbulence within the near-wall region, whereas delta wing mixers are responsible for increasing turbulence in the core flow [21].

An interesting example of CFD-aided studies on mass transfer improvement is the work concerning the Taylor vortex photocatalytic reactor [73,74]. The flow instability created in this type of reactor enhances the photocatalytic efficiency, as confirmed in both experimental and computational simulation studies.

2.3. Lamp Emission and Radiation Modeling

It is possible to distinguish three types of illumination in photocatalytic reactors: (i) direct illumination, such as in immobilized photoreactors; (ii) with the lamp immersed in the reaction environment (slurry reactors); and (iii) with reflecting devices (parabolic or elliptical). This determines the necessity to create lamp emission models [72]. Adequately, three main emission models are considered: (a) the line source model, (b) the surface source model, and (c) the volume source model. The schematic presentations of these models are shown in Figure 4 [75].

The significance of lamp emission models is connected to the effects of reflection, refraction, and absorption of radiation prior to entering the reaction space being incorporated as boundary conditions during CFD-aided radiation modeling [72].

As already mentioned in the Introduction, the solution of the radiative (photon) transfer equation (RTE) is important for the modeling of the radiation field. For monochromatic radiation, this equation is defined as [33]:

$$\frac{\mathrm{dL}\left(\mathrm{r},\mathrm{s}\right)}{\mathrm{dz}}+\left(\mathsf{\kappa }+\mathsf{\sigma }\right)\mathrm{L}\left(\mathrm{r},\mathrm{s}\right)={\mathrm{j}}^{\mathrm{e}}\left(\mathrm{r}\right)+\frac{\mathsf{\sigma }}{4\mathsf{\pi }}{\int }_{4\mathsf{\pi }}\mathrm{L}\left(\mathrm{r}, \mathrm{s} \prime \right)\mathrm{p}\left( \mathrm{s} \prime \to \mathrm{s}\right)\mathrm{d} \mathsf{\Omega } \prime ,$$

(9)

where L is the photon radiance, r is the position vector, s is the propagation direction vector, z is the path length, κ is the absorption coefficient, σ is the scattering coefficient, j^e is the emission term, p is the phase function for the inscattering of photons, and Ω′ is the solid angle about the scattering vector s′. For photocatalytic reactors (slurry systems), due to light scattering in the presence of the photocatalyst, it is impossible to find an analytical solution for the RTE; therefore, the proper mathematical model is needed. The RTE solution provides the opportunity to calculate the local volumetric rate of energy absorption (LVREA) and the local volumetric rate of photon absorption (LVRPA). LVREA depends on the photon distribution in the reaction medium and is useful for the description of slurry-mode photoreactors. As defined in the photocatalytic process, photocatalytic reaction might proceed as the consequence of photocatalyst excitation by the light of a specified wavelength. The number of reacted molecules is proportional to the number of absorbed photons, which can be evaluated by LVREA. Therefore, this is a linkage between LVREA and photocatalytic efficiency.

The first common approach to solve RTE is the application of the Monte Carlo method. This is the class of methods involving statistical sampling to model complex chemical and physical phenomena. During the Monte Carlo simulation, RTE is not solved explicitly, but the trajectories and fates of all photons emitted from the lamp are traced. They are both decided with the usage of random numbers. Depending on the underlying hypotheses, various versions of the Monte Carlo method can be distinguished [76]. This method is often used to provide simulations of photoreactors [24,25,26,27]. For example, Akach et al. performed simulations of the light distribution in a solar photocatalytic bubble column reactor [27]. Figure 5 shows the polar plots of LVREA depending on photocatalyst loading. Computational simulations might help adjust optimal photocatalyst loading without performing appropriate experiments.

The P1 approximation method is another approach, which assumes that the angular distribution of radiation intensity is almost isotropic within the medium. Accordingly, the radiative energy that reaches any given point inside the reactor comes almost equally from all possible directions. This assumption is important in the view of the scattering phenomenon because it is possible to eliminate the directional behavior (from the radiation field) that originated from its sources [29]. This method has also been applied in photoreactor simulations [28,29,30].

The third important numerical method of RTE solution is the discrete ordination method (DOM). The DOM discretizes the infinite number of directions involved in Equation (9) to a finite number of directions, s_i, customized to the geometry of the system [29]. The presented description is too simplified because an important issue is to establish the proper boundary conditions, directions, and spatial discretization for any change in the reactor geometry. This flexible and accurate method has also been extensively studied for photocatalytic reactors [19,22,31,32,33,34,35,36,37,38].

2.4. Kinetic Modeling

Most kinetic models are based on the Langmuir–Hinshelwood (L–H) rate form, which includes the phenomenon of the adsorption of the reagent (usually pollutant) on the photocatalyst surface [76,77]. In slurry photocatalytic systems, degradation of the compound depends on the operating conditions, such as pH, temperature, irradiance, oxygen concentration, initial concentration of the compound, and photocatalyst loading (W_cat) [72,78,79,80,81,82]:

$$-\mathrm{r}=\mathrm{k}\left(\mathrm{f}\left[\mathrm{P}\right]\right),$$

(10)

$$\mathrm{k}=\mathrm{f}\left(\mathrm{pH}, \mathrm{T}, \mathrm{E}, \left[{\mathrm{O}}_2\right],{ \mathrm{W}}_{\mathrm{cat}}\right),$$

(11)

where f[P] is the function of the compound concentration and is usually of the first-order or L–H form. The following issue is to consider the effect of light intensity on the reaction rate. When the value of LVREA (corresponding to L^a) and its relation to light intensity are considered, then Equation (10) can be modified as follows [72,76]:

$$-\mathrm{r}=\mathrm{k}{({\mathrm{L}}^{\mathrm{a}})}^{\mathrm{m}}\left(\mathrm{f}\left[\mathrm{P}\right]\right).$$

(12)

At large photocatalyst loading, the reaction rate usually decreases due to a higher electron-hole recombination rate. As a result, the volumetric rate of recombination is directly proportional to photocatalyst loading (W_cat). Then, the final form of the kinetic equation can be presented as [42,72]:

$$-\mathrm{r}={\mathrm{k}}_1\left(\mathrm{f}\left[\mathrm{P}\right]\right)\left[{\mathrm{L}}^{\mathrm{a}}-{\mathrm{k}}_2{\mathrm{W}}_{\mathrm{cat}}\right].$$

(13)

3. Complex Photocatalytic Reactor Modeling Systems—The Available Software Tools

The above-described steps of photoreactor modeling might be presented as shown in Figure 6. The CFD method is an important tool to proceed with the complex modeling of photocatalytic reactors. It might combine all physicochemical phenomena involved in the process: fluid dynamics, species transport, radiation transport, and including its coupling with the chemical reaction rate, which has been described above. There are various CFD simulation frameworks available. Some of them are commercial, such as ANSYS^® Fluent (ANSYS Inc., USA) and COMSOL Multiphysics^® (COMSOL Inc., USA), and others are open-source CFD software, e.g., OpenFOAM and FEATFLOW.

A suitable example of complex modeling of photocatalytic reactors by available general-purpose CFD software (ANSYS^® Fluent) is presented in the work of Casado et al. [83]. The simulation relates to the annular photocatalytic reactor. The analyzed model includes the description of hydrodynamics, radiation transfer, mass transport, and chemical reaction rate based on a mechanistic kinetic model. The results of photocatalytic activity, using methanol as a model compound, show agreement between model predictions and experimental data, with errors between 2% and 10%, depending on the photocatalyst loading (Figure 7) [83].

In the newest research performed by Moreno-SanSegundo et al., numerical simulations of three photocatalytic reactors were demonstrated: an annular reactor illuminated by a mercury fluorescent lamp, a tubular reactor coupled to a compound parabolic collector illuminated by sunlight, and a tubular photoreactor illuminated by LEDs [84]. The authors proposed a novel discrete ordinate model integrated with open-source OpenFOAM CFD software. The model predictions were successfully validated by comparing simulations conducted by commercial ANSYS^® Fluent software (Figure 8). The proposed simulation method allows for substantial improvements, particularly in the modeling of the LED photoreactor. This is especially important from the point of view of the rapid development of this type of light source-driven photocatalytic reactors.

The application of available CFD software to perform the simulations of photocatalytic reactors is based on an adaptation of the existing simulation modules to conduct simulations of photoreactors. This configuration must be operated by a user familiar with CFD software; therefore, this might generate problems for non-experts in the design of photoreactors. Another approach is to create the complex simulation tool by direct programming of the proposed mathematical photocatalytic reactor models, but it can be difficult to perform for researchers because of the complexity of the project and the extended knowledge of programming languages. On the market, it is possible to find commercial simulation packages dedicated for the chemical industry, such as Aspen HYSYS^® (Aspen Technology Inc., USA), but they do not have implemented modules for photocatalytic reactors.

Therefore, there is huge potential for complex simulation tools dedicated only for photocatalytic reactors. A promising example in this direction is PHOTOREAC, established as an open-access application developed in the graphical user interface of MATLAB^® software [85]. Acosta-Herazo et al. performed the simulations with the PHOTOREAC environment, which includes three configurations of pilot-scale solar photoreactors: a flat plate photoreactor (FPP), a compound parabolic collector photoreactor (CPCP), and the offset multi-tubular photoreactor (OMTP). Each photoreactor system operates in recirculation with a flow-through mode with water passing through an external tank [85]. As a model photocatalyst, Evonik P25 TiO₂ was selected [83], as it is commonly used as a standard due to one of the highest photocatalytic activities among titania photocatalysts [86,87,88,89].

The PHOTOREAC system consists of two main modules: (i) the photon absorption-scattering module (radiation field modeling) and (ii) the kinetic module (estimation of the radiation-independent kinetic parameters). In the case of the first module, a six-flux absorption-scattering model (SFM) was implemented. This is a simplified model, in which the main hypothesis is that scattering only occurs in the six Cartesian directions. Its application allows for short computational time and low mathematical complexity. The authors incorporated the SFM variant coupled with the Henyey–Greenstein scattering phase function (SFM-HG). Figure 9 shows a screenshot of the PHOTOREAC application window with a working photon absorption-scattering module [85,90].

It is thought that this promising open-access software approach, according to the authors’ recommendations [85], might be useful for introduction to photoreactor engineering in situations when a quantitative margin of error is still acceptable or qualitative results are the main aim of the work, and when the parametric space in the study is extensive. The authors are aware of the software limitations; thus, a plan to extend the capabilities of this CFD tool, dedicated for photocatalytic reactors, was proposed, e.g., the extension of simulation possibility to other photocatalysts and user-defined kinetics expressions [85].

4. Concluding Remarks

Heterogeneous photocatalysis, especially under natural solar radiation, seems to be the Earth’s future for energy conversion and environmental purification. Therefore, novel technologies of photocatalyst use must be developed for broad and easy application. Accordingly, the scaling-up of photocatalytic reactors is an important step to commercialize photocatalytic processes. The design of photocatalytic reactors with their particular complexity should be strictly interlinked with the proper computational tools for scaling-up. There are many processes and parameters that must be considered during the scaling-up procedure. Therefore, methods with a properly constructed system of mathematical models (hydrodynamics, radiation, and kinetics) based on the fundamentals of chemical engineering with an efficient algorithm should be constructed to provide reliable predictions of the performance of designed photoreactors. Furthermore, the constantly growing computing power of available hardware is another factor that should favor the development of these methods. Currently, the computational methods used in research on photoreactors are mainly based on general-purpose CFD tools, which are not strictly dedicated for this application. These simulations require much time and specific knowledge to be adopted in given conditions, to introduce the data, to launch proper models, and to perform the simulations. Nevertheless, the available research papers show that these computational tools provide increasingly accurate results. The latest successful and promising example of the PHOTOREAC MATHLAB^®-based application outlines the recommended direction of research in this field, which is to design computational tools taking into account the specific character of photoreactors to make the simulation procedure simpler with an appropriate level of accuracy.