Study on particle and photonic flux distributions in a magnetically stirred photocatalytic reactor

Abstract. When dispersed in water, nanoscale photocatalysts can aggregate into microscale secondary particles due to their high surface energy. The interaction between the incident photons and the aggregated particles is expected to be significantly changed due to their comparable length scales. Hydrodynamics of the particles after aggregation and the resultant photonic flux distributions in slurry photocatalytic reactors are therefore essential. The magnetically stirred photocatalytic reactor is compact and simple for fabrication, and therefore can be readily employed in lab-scale photocatalytic tests. However, studies toward its optimization are rare. In our study, the photocatalyst distribution was simulated using a Eulerian–Lagrangian approach, while the evolution of free liquid surface led by magnetic stirring was modeled by volume of fluid method. Subsequently, based on the photocatalyst distribution, the photonic flux distribution was obtained through a mean free path-based Monte Carlo method. Outcomes suggest that at a stirring speed of 900 rpm, the 10  μm particles can be well suspended, and a moderate free liquid surface will also be present. Moreover, under sufficient stirring, larger catalyst particles are more likely to be densely distributed in the outer region, which contributes to an increase in the overall photonic absorption even higher than the one with evenly distributed photocatalysts.


Introduction
Photocatalytic technology has been regarded as a prominent approach for hydrogen generation and water treatment owing to its mild reaction conditions and environmental friendliness. 1,2 It is now well established that photocatalytic reactions occur by reactions of photogenerated charges at the photocatalyst surface either directly with adsorbed molecules or by creating reactive intermediates that desorb from the photocatalyst and react with solution species. In-depth study has been made on photocatalyst materials to enhance their photochemical conversion efficiency. 3,4 On the other hand, understanding of reactor-level hydrodynamics and the corresponding photonic flux distribution might provide promising approaches to further increase the applicability of photocatalysis, and thus has received increasing attention.
Generally, photocatalytic reactors are either in slurry (in which photocatalyst particles are dispersed in the aqueous solution) or in immobilized form (with the photocatalyst attached to a surface). As the slurry photocatalytic reactor possesses high photocatalyst surface areas and mass-transfer rates, it is preferred in both lab-scale and scaled-up applications. To guarantee photocatalyst particles are well-suspended and illuminated, hydrodynamics and radiation distribution (both internal and external) are specially focused. On this premise, reactors have been developed for different applications that use different light sources as well as different geometries, including thin films, 5 moving beds, 6 fluidized beds, 7 annular recirculating, 8 top-irradiated cylindrical, 9 side-irradiated irregular, 1 and the scaled-up parabolic trough reactor. 10 Experimental measurement can be involved in the design and optimization of these slurry reactors. However, in some cases, experimental measurement to determine the particle and photonic flux distributions is rather demanding; modeling will then be needed to quantitatively evaluate these reactors.
From the perspective of hydrodynamics, due to interphase interactions, photocatalyst particles will generally migrate and distribute according to the flow regimes in reactors. Semiconductor photocatalyst particles in the micrometer size regime are typically dense and thus disperse poorly, resulting in poor photocatalytic performance; a qualified reactor therefore must generate flow fields that result in good suspension of photocatalyst particles of this size. The fast-emerging computational fluid dynamics (CFD) technique has been employed to study the catalyst-liquid flow. Usually, simulations were based on the time-efficient Eulerian-Eulerian approach with the solid phase considered to be continuous. [11][12][13] Another methodology for the CFD simulation of liquid/solid system is the Eulerian-Lagrangian (E-L) approach, which describes the discrete particle's behavior. 14 In the E-L approach, the interphase forces can be exerted on every single particle, and necessary collisions can be included. Although the E-L approach can be time-consuming, it might be helpful in cases where the discrete phase concentrations are low and the flow domain is small. Besides catalyst-liquid two-phase flow, there might exist two immiscible fluids in some kinds of reactors, in which gas/liquid interface needs to be considered, i.e., the top-irradiated cylindrical and side-irradiated irregular reactors. The shape of free liquid surface can be determined through an interface capturing scheme based on volume of fluid (VOF) method. For such reactors with a relatively small size, stirring is usually employed to maintain suspension. Challenges therefore also exist in modeling swirling flow field, which leads to the free liquid surface. 15 For photocatalytic reactions, participation of light is indispensable. Although most of the photocatalyst powders are expected to be in nanoscale during preparation, they will aggregate into microscale secondary particles when dispersed in water due to their high surface energy. Scattering of radiation by submicron or micron particles, instead of nanoparticles, in slurry photocatalytic reactors performs essential roles. 16 In this case, the length scale of the incident photons will be comparable to the particle size. The modeling radiative transfer processes in such scattering media are considerably more complex than in pure absorbing environments, due to the mutable ray direction and the nonuniform photocatalyst concentration. 17 Considering the radiative transfer as discrete photons passing through a continuous multiphase medium and using Monte Carlo method to determine the photon's direction or/and free path is an efficient approach. Two flux, six flux, and mean free path (MFP)-based Monte Carlo method can be classified into that category. Among them, MFP method, which comprehensively describes the photon's direction and free path using a phase function and a concentration-dependent method, is an efficient tool. 18,19 For lab-scale application, both the top-irradiated cylindrical and side-irradiated irregular reactors have been employed. They share similarities in their small scale, high mass transfer rate, and the existence of free liquid surface led by magnetic stirring. However, in the top-irradiated cylindrical reactor, incident radiation usually has to pass the gaseous region and will be partially absorbed by the water spray. In such a top-irradiated reactor, enough negative pressure has to be kept to avoid droplet formation. In a side-irradiated irregular reactor, the catalyst-liquid suspension is directly illuminated; thus, it can be operated under a normal pressure. The aforementioned magnetically stirred photocatalytic reactor (MSPR) is a typical side-irradiated irregular reactor. Despite the existence of the planar light-receiving window, the curved boundary of the reactor as shown in Fig. 1 enables the growth of flow loops without causing unnecessary vortex. Owing to its convenience and economic applicability, MSPR has been employed in numerous lab-scale experimental studies and achieves an acknowledged performance. 1,3,20,21 In our lab, in particular, MSPR has been integrated with an automatic sampling system for long-term hydrogen generation evaluation. 1 However, the optimization work of such photoreactor is rare. It apparently needs parameters such as the reactor volume, optimal photocatalyst loading, and stirring speed. These can be achieved by numerical simulation methods, which can significantly reduce the laborious experimental work, broaden the range of experimental investigation, and provide some valuable theoretical guidance. In this paper, the flow field with free liquid surface, photocatalyst and photonic flux distributions in MSPR will be proposed comprehensively. It is expected that this work will potentially be an initial step for this lab-scale photocatalytic reactors' performance optimization.

Problem Description
Irregular spherical MSPR, with a planar light-receiving window for incident radiation, is often illuminated by a Xenon lamp, which generates collimated beams. The complete system can be seen in Fig. 1(a), while the MSPR under operation, with a catalyst loading of 0.5 g∕L TiO 2 (Degussa P25) is shown in Fig. 1 By magnetic stirring, photocatalyst particles are dispersed in the aqueous solution. When this reactor is employed for hydrogen generation evaluation, void volume initially filled with nitrogen gas must be maintained above the slurry flow to collect the generated gases. Due to the magnetic stirring at the bottom of the reactor, apparent free liquid surface between gas and liquid phases will possibly exist. Figure 1(c) shows the free liquid surface flow under various magnetic stirrers' rotational speeds [simplified as rotational speed(s), N, in further discussion]. To consider such an evolution of free liquid surface, VOF method is employed. On the premise that a convincing flow field is obtained, through the E-L approach, the influence of key parameters, including stirrer's rotational speeds and particle sizes, on photocatalyst distribution will be investigated. Fundamentally, sliding mesh method is applied to describe the magnetic stirrer's rotation in CFD software. Moreover, based on the simulated nonuniform photocatalyst distribution, the home-made MFP-based Monte Carlo program will incorporate the concentration distribution outcomes to solve the photonic flux distribution. 19 3 Computational Model

Modeling of Free Liquid Surface Flow
VOF method, which is based on the Eulerian-Eulerian approach, is well suited to determine the aforementioned free liquid surface in MSPR. 22 In this approach, gas and liquid phases share the same velocity and turbulence fields within the whole computational domain, which can be determined by solving a set of governing transport equations with the volume-weighted mixture ∂ ∂t ðα q ρ q Þ þ ∇ðα q ρ q uÞ ¼ 0: (1)

Conservation equation of momentum:
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 4 6 Here, α is the volume fraction of relevant phase and Γ is the source term, which represents the gas-liquid interaction. Subscript q represents particle-laden liquid phase (pl) or gas phase (g), which satisfies The momentum equation in VOF, shown above, is dependent on the volume fractions of gas and particle-laden phases through their density ρ and viscosity μ.
When the magnetic stirrer's rotational speed reaches 600 rpm, the rotational Reynold's number will be of the magnitude of 10 6 , exceeding the critical point (10 5 ) for turbulence. Various turbulence models have been extensively compared in modeling the swirl flows. Reynold's stress models present to be more precise than k-ε and shear stress models and much more time-efficient than large eddy simulation. 22 Therefore, Reynold's stress model was employed in our simulation.

Modeling of Particle Distribution
Photocatalyst-liquid flow simulated using an E-L approach shares a similar continuity equation with the gas-liquid flow, but a different momentum equation.
Here, α 0 expresses the volume fraction of relevant phase in particle-laden phase, and the subscript s represents liquid phase or particle phase, which satisfies P α s ¼ α p þ α l ¼ α pl . S is the source term, which describes the particle-liquid interaction. In our case, the suspension is quite dilute to enable light penetration, and thus, the frequency of collision may be very small accordingly. Also, the magnetic stirring provides high shear rate. It is thus reasonable to consider that the flow field will dominate the particle distribution. Correspondingly, the following assumptions are made for simplicity of modeling: (1) the particles in the reactor are considered to be spherical with equal diameters, (2) collisions between particles can be neglected due to the low particle loading in our simulated system, (3) particle rotation was also neglected considering that usually it is the particle-particle collision that generates torque and causes particles to rotate. Such assumptions have also been proposed to be reasonable in literature modeling photocatalyst (particle)-liquid two-phase flow. 14,23 The particle's collision and reflection with the magnetic stirrer and the reactor boundary were described by a linear spring-dashpot. 24 Despite the situation that particles contact the stirrer or reactor wall, the transitional motion of the i'th particle is therefore determined by the sum of forces exerting on it.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 8 3 The momentum sink S is a volumetric summation of the forces between the liquid and particles. F p is the force exerting on the fluid in a CFD mesh cell and V cell is the cell volume. In this work, buoyant force F B , drag force F D , pressure gradient force F Pg , and virtual mass force F Vm have been considered. 14 Thus, S can be obtained by For the drag force F D , the free stream model is employed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 6 8 8 ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 2 8 where C D is the drag coefficient. The pressure gradient force here refers to acceleration pressure gradient in fluid, while the virtual mass force relates to the force required to accelerate the surrounding fluid. The force is also called the apparent mass force, because it is equivalent to adding a mass to a particle. Pressure gradient force: 25 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 5 3 2 Virtual mass force: 26 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 4 7 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 4 4 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 4 0 3 where C Vm is the virtual mass force coefficient (dimensionless) and A c is a dimensionless number. When the simulated system is of high particle Reynold's number, particle-fluid interaction forces such as Saffman lift force and Brownian force are negligible. 14,27 Due to the high shear rate and large particle size in our simulated system, when impacts of these forces are considered, particle concentration distribution exhibits very limited difference (<1%). Thus, they are neglected in our simulation.

Modeling of the Photonic Flux Distribution
In heterogeneous photocatalytic reactions, radiation absorption and conversion on catalyst particles are of great significance, directly determining the ultimate reaction rate. Photocatalysts are dispersed in a multiphase slurry state, in which the transmission of light is affected by both liquid and catalysts' scattering and absorption. In the MATLAB ® code compiled by our group, the MFP is chosen to describe the radiation absorption phenomena brought by particles' sizes and concentrations. Generally, in this method, if a local point is of high particle mass concentration, the absorption rate of photons here will be high, and the photon's MFP obtained will be smaller accordingly. It means that the photon will linger in this region for more steps and has a larger probability of being absorbed. The flowchart for our MFP calculation process is shown in Fig. 2. Random position on the light-receiving window can be decided by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 1 2 2 where r m is the light-receiving window radius.
The MFP of the incident is obtained according to the local particle concentration via E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 2 6 7 where ξ i (i ¼ 1 to 6) is a random sample from Uð0;1Þ; I 0 and IðrÞ are the intensities of the incident light at the initial stage and at position r (vector), respectively; β is the light extinction coefficient. It should be well addressed that the extinction is proportional to the local concentration of catalyst particles C (m −3 ) and inversely proportional to the particle diameter d p , which indicates our strategies to correlate the photon transport with particle concentration and size. 28 For simplicity, the incident photon is assumed to be emitted from an artificial chromatic light source with a wavelength of 365 nm, where the assumed photocatalyst has an extinction coefficient of 5.61 × 10 4 cm 2 g −1 , and asymmetry parameter g of 0.5026. 29 Assuming the photon will be absorbed if it encounters the particle, photon absorptance by particles at a local point can be derived in the following way. The average volume that each photocatalyst particle occupies in the solution is V 0 ¼ V∕NP total , where NP total is the particle number in the whole reactor; therefore, the illumination area of the exterior square can be expressed as in Eq. (14). Details for the description of this equation can be found in our previous work. 9 Ultimately, the local probability of absorption (local absorption rate) of the incident photon, ξ, is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 6 7 4 where A is the projected area of photocatalyst particle (cm 2 ) and ρ s is the density of catalystliquid suspension, which can be easily correlated with the local particle mass concentration. Therefore, the local photon absorption rate is a function of the local particle mass concentration. As shown in Fig. 2, if ξ 3 is larger than the absorption rate ξ, the photon will be scattered to the surrounding region. The probability distribution function for the MFP of the photons can be derived, if the light intensity is described in terms of numbers of incident photons. It can be described as the fraction of incident photons with propagation step equaling or exceeding jrj, like E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 5 2 7 PfR ≥ jrjg ¼ The photon's scattering length is determined by random sampling from the MFP as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 4 7 2 Moreover, in a heterogeneous medium, the extinction coefficient in the photon path might change. For discrete regions, light intensity changes can be described by As for photon's scattering direction, the directional distribution of the intensity of the scattering light in three dimensions can be described by the Henyey−Greenstein phase function. 30 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 6 ; 3 2 9 Φðθ; ϕ; gÞ ¼ 1 4π where θ is the polar angle of the scattering direction with respect to the incident direction, ϕ is the azimuthal angle of the scattering direction with respect to the incident direction, and g, ranging from −1 to 1, is the asymmetry parameter for light scattering, which describes portion of forward or backward scattering. The sampling expression for θ can be obtained by the sampling transformation method.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 2 2 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 6 ; 1 7 2 ϕ ¼ 2πξ 6 : For simplification, all boundaries of the photocatalytic reactor and the gas-liquid interphase are considered to be transparent for photons. Photons that enter the surrounding region of magnetic stirrer will also be considered to be out of the effective region of photocatalysts' radiation absorption. It is worth noting that considering the demanding requirement of E-L model, the particles added into the simulated system are less than the real situation. When calculating the radiation distribution, the particle concentration distribution is magnified to an average level of 0.5 g∕L and loaded into MATLAB ® to interpolate each point in the concentration field. Particles' sizes and stirrer's rotational speeds are assumed to be the main parameters leading to various light-absorption phenomena and will thus be both considered.

Numerical modeling
The physical model is identical with the geometry of the real reactor shown in Fig. 1. As shown in Fig. 3(a), the simulated object is a sphere (R ¼ 40 mm, center of sphere is the coordinate origin) with three cross-sections on its bottom, top, and right side, respectively. The reactor was agitated by a cylindrical magnetic rod, lying at the bottom of the tank, which has a total length of 25 mm. The radius of half spheres at the two ends of the stirrer is 4.5 mm.
To reveal the rotation of the stirrer well, the mesh in the moving fluid region is six times denser than the stationary one in volume. The reactor geometry is meshed with 1.6 × 10 5 unstructured tetrahedral cells, as shown in Fig. 3(b). For mesh-dependency analysis, velocity magnitude along the axial direction was calculated with various grid numbers. As shown in Fig. 4, modeling over 1.6 × 10 5 cells and 3.2 × 10 5 cells shares almost identical results. It is therefore assumed that the 1.6 × 10 5 cells should be sufficient for our simulation.
In total, 2 × 10 6 TiO 2 particles with a uniform diameter of 3, 5, or 10 μm have been randomly added into the liquid domain of the simulated system after the fluid domain has reached a quasi steady state. The initial free liquid surface position is at y ¼ 2 cm, and the total volume of the solution inside is 199.8 cm 3 . Main physical parameters employed in simulations are shown in   The fluid phase is treated as continuity in an Eulerian reference frame, while the particle phase is treated as discrete phase in a Lagrangian frame. The sliding mesh method is applied to simulate the rotation of the magnetic stirrer. The surrounding region around the stirrer is defined as the moving fluid, while other region outside the moving fluid is defined as stationary. These two regions exchange data through defined interfaces, and the mesh of the moving region is reconstructed every time step. The operation ensures the true revealing of the interaction between the stirrer and the irregular reactor wall due to the change of their relative position. No slip condition is employed for the reactor wall and the stirrer, and the zero shear condition is employed for the upper boundary of gas phase.
The continuum model of liquid and gas is solved by the pressure implicit with splitting of operator method, and its momentum equations are solved by the second-order implicit time integration. The explicit time integration method is used to solve the transitional motions of particles. For every time step during simulation, the magnetic stirrer will rotate no more than 0.2 deg.

Model validation
To validate the accuracy of our model for prediction of radiation absorption, we applied the experimental condition reported in literature to simulate the forward photon flux at the outer wall of photoreactor under different TiO 2 (Degussa P25) loadings. 31 As shown in Fig. 5, the experimental photonic flux intensity was compared with our MFP simulated result. It is seen that our simulation is in good agreement with the experimental results.

Flow Field in the Reactor
In the irregular MSPR, the relative position of the stirrer and the wall varies from time to time. Therefore, obtaining an absolute steady state of flow field is almost impossible. Nevertheless, the fluctuation of flow field in our simulation is within an acceptable range. The case with a fluctuation of velocity within 3% of the average value can be assumed to be within a quasi steady state. The flow patterns were captured once the free liquid surface is stabilized and the mean flow is statistically steady. The predicted flow patterns at a rotational speed of 600 rpm, showing the mean velocity vectors projected on XY plane and YZ plane of the coordinate system, respectively, are given in Fig. 6. Figure 7 shows the predicted flow pattern projected on XY plane at N ¼ 900, 1200, and 1500 rpm, respectively. The concave free liquid surface induced by high rotational speeds can be clearly observed, which can be fundamentally validated by the photographs in Fig. 1  (c). As can be seen, the cylindrical stirrer produces a radial flow, which impinges on the reactor walls and deflects toward the liquid surface. This stream then returns to the stirrer and a recirculation loop is formed. Such a trend can be well validated with previous modeling of agitated vessels. 15,32 But the vortices shapes are significantly controlled by the reactor boundary. As shown in Figs. 6(a) and 7, due to the existence of flat wall, which is designed for effective reception of incident light, the vortex at this side is severely deformed and suppressed, which will be quantitatively explained later. In contrast, as shown in Fig. 6(b), two vortices observed on YZ plane are almost symmetrical and well developed. It is very obvious that the existence of the light-receiving window will exert considerable impact on the flow pattern in the reactor.
As shown in Fig. 8, the positions of the free liquid surface nadir were plotted versus their corresponding stirring speeds. The deformation of the free surface with the increase of rotation speeds can be clearly observed. For the following reasons, the concave free liquid surface would not be an ideal choice. In the first place, the photocatalytic reactor design must guarantee that incoming photons are sufficiently utilized and do not escape without having intercepted particles in the reactor. The existence of concave free liquid surface indeed increases the area of the slurry flow's upper boundary and disturbs photon's transport pathway in the middle region. It therefore leads to photon losses. Also, the unsteady motion of free liquid surface can cause irrational results in experiments.
Under a magnetic stirrer's rotational speed of 600 rpm, axial, radial, and tangential velocities of the flow at various heights on XY plane are shown in Fig. 9, respectively. These velocities were all divided by the tip speed of the stirrer U tip . Analysis of axial and radial velocities corresponds well with the recirculation loop flow pattern. In Fig. 9(a) (axial velocity distribution), upward flow on the reaction boundary region and downward flow on the central region can be revealed. It is found that the impingement with the flat wall (light-receiving window) retards the upward flow on this side. In Fig. 9(b) (radial velocity distribution), the outward flow produced by the stirrer on the bottom region and the inward flow on the upper region can be revealed. The flow is accelerated by the stirrer radially and reduces its speed after being injected from the stirrer. In Fig. 9(c) (tangential velocity distribution), force-vortex and free-vortex region can be easily distinguished. From the central region, forced by magnetic stirring, the tangential   velocity will be accelerated and reaches a velocity peak. After the tangential velocity peak, the tangential velocity decreases gradually and that region is called the free-vortex region. Here, the tangential flow near the light-receiving window is accelerated due to the existence of the lightreceiving window. Velocity profiles here are generally similar to the magnetically stirred cylindrical vessel, but the existence of light-receiving wall has non-negligible effects on these velocity components. To achieve a comprehensive understanding of the effect of reactor geometry, further efforts in consideration of reactor's aspect ratio will be of use.

Particle Distributions
A general overview of particle concentration distributions can be seen in Fig. 10(a) (d p ¼ 5 μm, N ¼ 900 rpm). To demonstrate the trend of particle concentration distribution in such an irregular reactor, the average-weighted face value of number density NP was taken at different planes along the coordinate axis. Figures 10(b)-10(f) show the contour of the particle concentration distributions on various cross-sectional planes of the reactor vertical to the Y axis. Due to the centrifugal effect of the swirl flow, the particles are more concentrated at the outer region of the reactor, leaving an almost vacant vortex zone around the central domain. Enrichment of particles near the front and back wall of the vessel can also be viewed.
Particle concentration distributions in the reactor were also investigated by considering various particle sizes and rotational speeds, as shown in Figs. 11 and 12. A general trend is that with the increase of rotational speeds and the decrease of the particle sizes, average particle concentrations at different heights in the reactor become increasingly closer to each other, demonstrating a better mixing performance. Unlike cases of annular or cylindrical photocatalytic reactors, radial particle distribution in our case demonstrates anisotropic properties. Due to the centrifugal effect of the stirrer, particles tend to gather in the outer region around the wall. As shown in Fig. 11, particles with larger sizes have a larger concentration gradient along X or Z axis, while smaller particles are  distributed more uniformly. The existence of the light-receiving window increases the concentration gradient along the light-incident direction (X axis), because the tangential flow is accelerated near the window and a better suspension effect is reached. As mentioned, the recirculation brought by the flow loops performs an essential role in mixing the suspension flow. Larger particles are difficult to lift by the upflow along the light-receiving window. Also, the aforementioned axial velocity near the window is smaller, which is not favored to lift large particles. Though there exist those negative effects for particle suspension, at the given maximum particle sizes (d p ¼ 10 μm), particles can suspend well at N ¼ 900 rpm in general.
As shown in Fig. 12, higher rotational speeds lead to lager particle concentration gradient, with particles being more concentrated in the outer region. Along Y axis, rotational speeds of 900 and 1200 rpm generate almost the same particle distributions. In both cases, particle distributions are quite even, indicating that the rotational speed of 900 rpm seems to be sufficient for well suspension. Our results show that photocatalyst particles with diameters of up to 10 μm could be well distributed at a rotational speed around 900 rpm with a mild free liquid surface.

Photonic Flux Distributions
To our knowledge, the catalyst particle distribution significantly affects the effective penetration of the incident photons and therefore the distribution of radiation field. 33 Higher rotational speed also means a higher energy cost. Therefore, our discussion of the photonic flux distribution will be carried out under the circumstances of 600 and 900 rpm, respectively, in which the free liquid surface can be assumed to be mild. For comparison, uniform particle distribution was also considered. As shown in Fig. 13, in such a homogeneous medium with uniform particle distribution, the intensity of the incident radiation along the direction of photon propagation attenuated exponentially. These results are consistent with the literature. 34 Fig. 12 Normalized particle density along (a) X , (b) Y , and (c) Z axis, respectively, for 5 μm photocatalyst particle suspension under various rotation speeds. Fig. 13 Photon absorption rate per unit area (cm 2 ) at different planes in the direction of incoming photons (X axis), under a uniform particle distribution condition (C ¼ 0.5 g∕L).
Due to the strong absorption and scattering properties of the photocatalyst particles and reacting solutions, rapid decay of incident radiation might exist. Compared with the uniform catalyst distribution, the nonuniform distribution with a high concentration near the light-receiving window proves to be better, as shown in Figs. 14 and 15. This can be more clearly seen as we calculated the overall photon absorption rates for three cases, namely, uniform distribution, rotational speeds of 600 and 900 rpm, respectively, as shown in Fig. 16. Apparently, due to the participation of photons in the reaction, the optimal particle distribution is completely different. Radiative transfer in participating reacting media is the result of the interaction of a material multicomponent continuum with an immaterial photon phase. The two phases coexist in a given region in space and interact according to modes that are defined through the constitutive assumptions made for each medium. It is expected that in the regions where the photocatalyst particles are densely distributed, the light absorption is also high, and vice versa. 35 From this perspective, the outcome from our MSPR is considered to be rational.
As for the effects of particles size, the major local absorption difference caused by the particle size should be attributed to the particle concentration distributions. The particle concentration is strongly correlated with the rotational speed. Although particles with a diameter of 10 μm have a poor performance in suspension at N ¼ 600 rpm, its relatively high concentration near the window well compensates. At N ¼ 900 rpm, radiation absorption from suspended particles of 10 μm in size exceeds that of 5 μm particles. Our simulation demonstrates that photocatalytic reactors with smaller particle sizes benefit from the uniform concentration distribution in axial  position, while higher rotational speeds could result in nonuniform particle distribution in the direction of incoming photons, with the outside particle concentration being higher. Apparently, such nonuniform distribution is obviously preferred than the uniform one as for radiation absorption enhancement. Thus, it will be rational to consider the impact on photonic flux distribution also from photocatalyst concentration gradient instead of photocatalyst concentration alone when evaluating the performance of a specific photocatalyst or photocatalytic reactor. 9,19 According to our evaluation, the MSPR can perform well as a good photocatalytic reactor, as it can circulate particles well so that they effectively use incoming photons. However, the suspension state of particle is strongly correlated with the particle size, which implies the complexity in quantitatively characterizing the effects of particle size. As a further discussion, it is assumed that there exists an optimal particle size of photocatalyst particles in specific photocatalytic reactors, which is a result of competing effects of radiation absorption efficiency, specific surface area, and photochemical kinetics. 36 Also, even though smaller particle might benefit from its larger specific surface area and is preferred to enhance radiation absorption by increasing scattering, it is demanding to describe the intensity of the scattered light, which largely depends on the incident radiation wavelength when the particle size is of submicron level. 37 Besides, the photochemical reactions that take place on the catalyst's surface can be size-dependent as well, especially in the case with phase transformation (i.e., hydrogen generation) where the surface energy should be studied. Considering these challenges, however, finding the optimal size of particles and their distribution is believed to be of significance for further design of photocatalysts and photocatalytic reactors.

Conclusions
The aim of our present work is to evaluate the particle and photonic flux distributions in the MSPR of irregular shape by numerical methods. For better description of the discrete distribution of catalyst particles of micrometer sizes in a continuous reacting media, E-L approach is employed, while the photonic flux distribution is solved using MFP-based Monte Carlo method. Our numerical results elucidate the flow pattern, particle distribution, and therefore photonic flux distributions under various magnetic stirrer's rotational speeds and particle sizes. Generally, flow field and particle concentration along with photonic flux are strongly coupled factors in determining the energy efficiency of photocatalytic process. Rational flow field should lead to higher particle concentration in the light-incident direction, while the high particle concentration accompanied with the high radiation could contribute to a high energy conversion efficiency. Specific conclusions to be noted are listed below.
1. Larger rotational speed is preferred to maintain well-dispersed suspension of larger particles without causing a too severe free liquid surface. 2. Under the studied conditions, TiO 2 particles with diameters up to 10 μm could be well distributed under rotational speed of 900 rpm with an acceptable free liquid surface.
3. Smaller particles possess excellent overall light absorption performance, while larger particles are more likely to be densely distributed in the outer region, leading to the enhancement of local light absorption. 4. Overall photon absorption rate for photocatalyst particles distribution at the rotational speed of 900 rpm is even higher than the case for an ideal homogeneous particle distribution. This result can be potentially useful for the design of MSPR in the future.