Numerical Simulation of Energy and Mass Transfer in a Magnetic Stirring Photocatalytic Reactor

Abstract

Hydrogen production via photocatalytic water splitting is one of the promising solutions to energy and environmental issues. Understanding the relationship between hydrogen production in suspended photocatalytic reactions and various influencing factors is crucial for expanding the scale of the system. However, the complexity of physical and chemical factors involved in hydrogen production via photocatalytic water splitting makes systematic research of this technology challenging. In recent research, the simulated light source reactor has become a preferred study object due to its strong controllability. This paper presents a comprehensive energy and mass transfer model for the suspended photocatalytic reaction in a magnetically stirred reactor. The mutual impacts between the flow field, radiation field, and reaction field are analyzed. The simulation results show that the rotating speed of the stirring magneton in the reactor has a significant influence on the flow field. The rotation of the stirring magneton generates a vortex in the central axis area of the reactor, with the relationship between the depth of the vortex f(s) and the rotating speed of the magneton s described as f(s) = 0.27e^0.0032s. The distribution of radiation within the reactor is influenced by both the incident radiation intensity and the concentration of the catalyst. The relationship between the penetration depth of radiation g(i) and the incident radiation intensity i is described as g(i) = 10.73ln(i) − 49.59. The relationship between the penetration depth of radiation h(c) and the particle concentration c is given as h(c) = −16.38ln(c) + 15.01. The radiation distribution in the reactor has a substantial impact on hydrogen production, which affects the concentration distribution law of hydrogen. The total amounts of hydrogen generated in the reactor are 1.04 × 10⁻³ mol and 1.35 × 10⁻³ mol when the reaction times are 1.0 s and 2.0 s, respectively. This study serves as a foundation for the future scaling of the system and offers theoretical guidance for the optimization of the photocatalytic reactor design and operating conditions.

1. Introduction

Solar energy is a widely available, pollution-free, and abundant source of renewable energy [1]. However, due to its intermittence and instability, an energy storage mode is necessary to connect solar energy and consumers [2]. Hydrogen energy is a suitable option [3,4]. The conversion of solar energy into hydrogen energy can be achieved through various technologies, including photoelectricity and water electrolysis [5], solar thermochemistry [6], photobiology [7], photoelectrochemical [8], and photocatalysis [9]. Among these, solar photocatalytic hydrogen production is a promising technology for renewable energy hydrogen production due to its low energy consumption and low investment requirements [10].

The photocatalytic reaction system can be categorized into two types based on the state of the catalyst: the suspended reaction system and the immobilized reaction system [11,12,13]. Compared with the immobilized reaction system, the suspended reaction system is widely used due to its higher ratio of the active catalytic surface area to the reaction volume and reduced mass transfer limitation [14]. Tang et al. synthesized g-C₃N₄ dual-functionalized using an integrated strategy of spatially restricted growth [15]. The experimental results showed that the dissociation of excitons could be greatly promoted by the storage-free electron behavior of carbon nitride dual-functionalized with cyanoamino and K⁺. Yang et al. found that Bi doping in ZCS significantly enhanced light absorption and facilitated the separation and transport of photogenerated carriers through experiments [16]. Kargupta et al. prepared a hybrid photocatalyst to improve hydrogen production [17]. The optimum amount of catalyst was determined to be 30 mg through a comparison of photocatalytic performance. At present, most studies on suspended photocatalytic reactions are based on experimental methods, which allow for an immediate evaluation of the reaction system. However, numerical simulation research is also crucial in order to comprehensively analyze the impact of various factors, gain intuitive visual data, and minimize research costs [18].

The photocatalytic process is a multiphase reaction and flow process involving mass and energy balance. Numerical simulation of this process is challenging due to the diversity of reaction substances, light sources, flow conditions, and catalyst forms [19]. Kumar and Bansal simulated the photocatalytic process in an immobilized reactor using computational fluid dynamics for the first time [20]. The relationship between the photocatalytic degradation rate and the initial concentration, Reynolds number, and reaction time was established on the premise that the radiation distribution in the reactor was uniform and constant. Ren et al. conducted a numerical simulation study on the suspended photocatalytic reaction process under the condition of a compound parabolic concentrator [21]. A two-dimensional radiation model was established considering the influence of the concentrator on the radiation distribution inside the reactor. The effects of different catalyst concentrations, particle sizes, and inlet flow rates on the radiation distribution were analyzed. Yuan et al. modeled the single-phase multi-component flow reaction of photocatalytic reduction of carbon dioxide in a double-tube reactor with the catalyst fixed on the inner wall [22]. The simulation results were analyzed to determine the optimal structural parameters of the reactor. Bagheri and Mohseni established a two-dimensional flow reaction model for photocatalytic water treatment in a ring reactor [23]. The model, which was verified via experimental results, was capable of predicting the degradation rate under the assumption that the radiation distribution in the reactor was independent of the fluid motion parameters. Khodadadian et al. modeled the photocatalytic degradation of toluene in an annular reactor with an inner-wall-immobilized catalyst [11]. They emphasized the significant impact of the flow pattern in the reactor on the photocatalytic reaction and recommended its consideration in future studies. Castedo et al. carried out a simulation study of the complete reaction flow process of photocatalytic hydrogen production in an immobilized reactor. They analyzed the impact of light intensity and channel length on hydrogen production rate [24].

In conclusion, the complexity of physical and chemical factors in photocatalytic water-splitting systems has hindered the industrialization of hydrogen production [25]. Therefore, it is imperative to establish an energy and mass transfer model for the photocatalytic reaction process and gain an understanding of the flow field, radiation field, and reaction field in the reactor. The use of a simulated light source in photocatalytic water splitting has become a key tool for investigating the relationship between hydrogen production and various factors due to its strong operability [26]. In this study, the Eulerian–Lagrangian model is applied to simulate the flow field within the reactor. The volume of fluid method is used to track the moving phase interface for gas and liquid phases. Subsequently, the radiation field distribution in the reactor is calculated using the Beer–Lambert Law and the optical properties of the reaction solution, taking into account the distribution of catalyst concentration. Finally, the hydrogen production distribution in the reactor is obtained via combining the flow field and radiation field through the use of the reaction kinetics equation obtained from experiments.

2. Methodology

2.1. Introduction of System

The laboratory setup of the photocatalytic hydrogen production reaction is shown in Figure 1, including a simulated light source, a photocatalytic reactor, and a magnetic stirring apparatus. The simulation light source is the PLS-SXE300 xenon lamp light source. The spectral range of the light source is 300~2500 nm. The instability of the radiation intensity is less than ±0.5%. The spot diameter is 30~63 mm. The beam-divergence angle of the light source is 5~8°. The lamp has a minimum service life of 1000 h. The photocatalytic reactor is a spherical cap-shaped reactor made of Pyrex glass. The incident light window of the reactor is made flat to minimize changes to the optical path caused by the refraction of the circular glass wall. The diameter of the upper cover of the reactor is 33.0 mm, and the diameter of the incident light window is 50.0 mm. The magnetic stirring apparatus has a speed range of 100~1600 rpm with a display accuracy of 1 rpm. The stirring magneton is a cylindrical B620 type made of PTFE material.

2.2. Computational Model

2.2.1. Geometric Model

The geometric model of the numerical simulation is depicted in Figure 2. The parameters include the gas outlet radius (r₁ = 16.5 mm), the incident light window radius (r₂ = 25.0 mm), and the sphere radius (r₃ = 39.5 mm). The total length (l = 20.0 mm) and the thickness (h = 6.0 mm) of the stirring magneton are also shown. The bottom radius of the reactor is equal to the gas outlet radius.

The total volume of the reactor, as calculated via the following volume equation of a spherical cap, is approximately 245.9 mL:

$$V=\frac{\pi }{3}\left(3R-H\right)H^2,$$

(1)

where R is the radius of the sphere, and H is height of the spherical cap.

In this study, the sliding grid model is used due to the movement of the simulated domain caused by the stirring magneton in the reactor. The simulated domain is divided into two cells by two grid interfaces, as shown in Figure 3. Cell zone 2 is the fluid region that contains the stirring magneton. Grid interface 1 is the circular surface parallel to the xz plane, and grid interface 2 is the lateral surface of the cylinder, as shown in Figure 3b,c. During the simulation, the cell regions slide along the grid interfaces in discrete steps. The wall surface of the stirring magneton is set as a moving wall surface, and its motion mode is to rotate synchronously with the adjacent cell zone.

2.2.2. Flow Model

The actual volume of the photocatalytic reaction is 200 mL, which does not occupy the entire reactor, leading to a clear stratification of the gas and liquid phases. Additionally, the volume fraction of catalyst particles is minuscule, as demonstrated by a catalyst concentration of 0.5 g·L⁻¹, which results in a volume fraction of only 1.0 × 10⁻⁴ in the liquid phase. To simulate the flow fields of the continuous phases (i.e., gas and liquid) and discrete phase (i.e., catalyst particles), the Eulerian–Lagrangian model is employed in ANSYS Fluent software. However, the interaction between the catalyst particles and their impact on the liquid phase is not taken into account.

The VOF model is based on the equation that multiphase fluids are not interspersed with each other. The sum of volume fractions of all phases in each control volume is 1.0. This allows the variables and their properties in a given cell to either represent the mixing of phases or simply represent a single phase, depending on the volume fraction of the phases. Based on the volume of fluid method [27], the interface between the gas and liquid phases is determined via solving the continuous equation of the multiphase volume fraction, which is expressed as follows:

$$\frac{\partial {\alpha }_p}{\partial t}+v\cdot \nabla {\alpha }_p=0,$$

(2)

where α_p is divided into the liquid phase volume fraction α_l and the gas phase volume fraction α_g.

The constraint of volume fraction is expressed as follows:

$${\alpha }_l+{\alpha }_g=1.$$

(3)

The momentum equation that applies to all phases is expressed as follows [28]:

$$\frac{\partial }{\partial t}\left(\rho v\right)+\nabla \cdot \left(\rho vv\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla v+\nabla v^T\right)\right]+\rho g+F,$$

(4)

where ∇p is the pressure term, ∇[μ(∇ν + ∇ν^T)] is the viscous stress term, ρg is the gravity term, and F is the other force term.

The density and viscosity of the mixed fluid in the cell can be calculated using the following equations:

$$\rho ={\alpha }_l{\rho }_l+\left(1-{\alpha }_g\right){\rho }_g;$$

(5)

$$\mu ={\alpha }_l{\mu }_l+\left(1-{\alpha }_g\right){\mu }_g.$$

(6)

The energy equation can be expressed as follows:

$$\frac{\partial }{\partial t}\left(\rho E\right)+\nabla \cdot \left[v\left(\rho E+p\right)\right]=\nabla \cdot \left(k_{eff}\nabla T-{\displaystyle \sum _j{h}_j{J}_j}+{\tau }_{eff}\cdot v\right)+S_h,$$

(7)

where the right side of the equation denotes the energy transfer caused by heat conduction, diffusion, and molecular viscous dissipation, respectively; S_h denotes heat transfer due to the chemical reaction; E denotes the average variable of mass, which can be expressed as follows [29]:

$$E=\frac{{\alpha }_l{\rho }_l{E}_l+{\alpha }_g{\rho }_g{E}_g}{{\alpha }_l{\rho }_l+{\alpha }_g{\rho }_g}.$$

(8)

In this study, the standard k-ε model is utilized to describe the turbulent momentum flux inside the reactor, as the flow state inside the reactor is turbulent under the studied rotation speeds (250~1000 rpm). This model, based on the concept of Boussinesq isotropic eddy-viscosity, has a linear stress–strain relationship and is widely used to simulate the turbulence in stirred vessels [30]. The standard k-ε model is used to capture the turbulence phenomenon in this study [31].

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{v}_i\right)=\frac{\partial }{\partial x_i}\left[\left(v+\frac{v_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon +S_k,$$

(9)

$$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho \epsilon v_i\right)=\frac{\partial }{\partial x_i}\left[\left(v+\frac{v_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}+S_{\epsilon },$$

(10)

where G_k represents the portion of turbulent kinetic energy produced by the laminar velocity gradient; G_b represents the portion of turbulent kinetic energy produced by buoyancy; σ_k and σ_b represent turbulence Prandtl numbers, which σ_k = 1.0 and σ_b = 1.3.

The trajectory of discrete phase particles is predicted by integrating the particle force differential equation in Laplace coordinate system. In this study, Cd₀.₅Zn₀.₅S is used as the photocatalyst particle. The force balance equation of particles in x direction, represented in the Cartesian coordinate system, is expressed as follows [32]:

$$\frac{d{v}_s}{dt}=F_D\left(v-v_s\right)+\frac{g_x\left({\rho }_s-\rho \right)}{{\rho }_s}+F_x,$$

(11)

where v is the fluid phase velocity; v_s is the solid particle velocity; ρ is the density of the fluid; ρ_s is the density of solid particles; g_x is the acceleration of gravity in the x direction.

The drag force can be expressed in the following equations:

$$F_D=\frac{18\mu }{{\rho }_s{d}_s^2}\frac{C_D\mathrm{Re}}{24},$$

(12)

$$\mathrm{Re}=\frac{\rho d_s\left|v_s-v\right|}{\mu },$$

(13)

$$C_D=a_1+\frac{a_2}{\mathrm{Re}}+\frac{a_3}{\mathrm{Re}}$$

(14)

where μ is hydrodynamic viscosity; d_s is the particle diameter, which is 10 μm in this study.

The additional mass force F_x includes consideration of the virtual mass force, the pressure gradient force, and the Saffman lift force, which can be expressed as follows:

$$F_{x-mass}=\frac{1}{2}\frac{\rho }{{\rho }_s}\frac{d\left(v-v_s\right)}{dt},$$

(15)

$$F_{x-pressure}=\left(\frac{\rho }{{\rho }_s}\right)v_s\frac{\partial v}{\partial x},$$

(16)

$$F_{x-lift}=\frac{2K{\upsilon }^{0.5}\rho d_{ij}}{{\rho }_s{d}_s{\left(d_{lk}{d}_{kl}\right)}^{0.25}}\left(v-v_s\right),$$

(17)

where d_ij is the fluid deformation rate tensor.

The top of the reactor serves as a gas outlet, while the other wall surfaces are set as non-slip boundary conditions, meaning that the fluid velocity at these walls is equal to zero. This is expressed mathematically as follows:

$$v(x,t)=V(x,t),x\in \Gamma ,$$

(18)

where V is the velocity of the solid wall, and Г is the boundary of the fluid field in contact with the wall.

Considering that the initial state of catalyst particles is precipitation at the bottom of the reactor when there is no stirring magneton, catalyst particles are selected to be injected into the simulated domain from the bottom of the reactor. To simplify the mathematical model, the photocatalyst particles are assumed to be ideal spheres of same size and evenly distributed at the bottom of the reactor. The random trajectory model is used to simulate the random effects of fluid turbulence on photocatalyst particles. The particle diffusion caused by fluid turbulence is calculated by considering the velocity fluctuations of the fluid. The steady Lagrangian discrete phase model is appropriate for solving problems with well-defined boundary conditions but not for describing particle flow suspended in the continuous phase. Thus, this study employs the unsteady particle discrete phase model to address the particle flow within the reactor.

2.2.3. Radiation Model

Since the beam-divergence angle of the xenon lamp light source is less than 8°, it is assumed that the rays of the light source are parallel to the incident light window of the reactor, as shown in Figure 4.

The flow model is simulated by discretizing the governing equation in the spatial region. Similarly, the radiation model established in this study also employs grid discretization to obtain the radiation field within the reactor.

Under the condition of ignoring the scattering effect of photocatalyst particles on radiation, Bill–Lambert Law is used to construct the radiation control equation. The optical parameters of the reaction solution were measured experimentally, and the governing equation of the radiation field distribution can be expressed as [12]:

$$I_{o,i}=I_{i,i}\cdot {10}^{\left[30.69\ln (C_i)-59.548\right]P_i{C}_{ci}},$$

(19)

where I_o,i is the radiation intensity after passing through the ith grid; I_i,i is the radiation intensity before passing through the ith grid; C_ci is catalyst concentration in the ith grid; P_i is the penetration depth of the ray passing through the ith grid along the ray incident direction.

2.2.4. Reaction Model

The mass fraction of the substance can be estimated using the convection–diffusion equation in a chemical reaction. The conservation equation is expressed as follows [33]:

$$\frac{\partial }{\partial t}\left({\rho }_p{Y}_{pi}\right)+\nabla \cdot \left({\rho }_{p}v{Y}_{pi}\right)=-\nabla \cdot J_{pi}+R_{pi},$$

(20)

where Y_pi is the mass fraction of component i in phase p; J_pi is the diffusion flux of component i in phase p; R_pi is the net production rate of component i in phase p.

In turbulent flow, the diffusion flux can be calculated using the following equation:

$$J_{pi}=-\left(\rho D_{pi}+\frac{{\mu }_t}{S{c}_t}\right)\nabla Y_{pi},$$

(21)

where D_pi is the diffusion coefficient of component i in phase p.

The overall reaction equation for the photocatalytic reaction system with the sacrifice agent studied in this study can be expressed as follows:

$$2H_{2}O+N{a}_{2}S+N{a}_{2}S{O}_3\to N{a}_2{S}_2{O}_3+2NaOH+H_2\uparrow .$$

(22)

The initial concentrations of the sacrifice agents are 0.7 mol·L⁻¹ Na₂S and 0.5 mol·L⁻¹ Na₂SO₃. Ignoring the intermediates, the liquid phase consists of five components: Na₂S, Na₂SO₃, NaOH, Na₂S₂O₃, and H₂O. Prior to the reaction, the reaction solution does not occupy the entire reactor; therefore, the gas phase comprises two components: air and hydrogen.

The calculation of the net production rate can be expressed as follows:

$$R_{pi}=M_{pi}{\displaystyle \sum _{i=1}^N{r}_{pi}},$$

(23)

where M_pi is the molecular weight of component i in phase p; r_pi is the rate of production or decomposition of component i in the reaction, which can be expressed as [12]:

$$r_{pi}=k^{*}\cdot I\cdot {\mathrm{e}}^{\raisebox{1ex}{$-E_a$}\!\left/ \!\raisebox{-1ex}{$RT$}\right.}\cdot \frac{\left(K_1\cdot C_1\right)\left(K_2\cdot C_2\right)}{1+K_1\cdot C_1+K_2\cdot C_2},$$

(24)

where k^* is the rate constant of intrinsic photocatalytic hydrogen production; I is the local radiation intensity; E_a is the apparent activation energy; K₁ and K₂ represent the adsorption equilibrium constants of Na₂S and Na₂SO₃ on the catalyst surface; C₁ and C₂ are the concentrations of Na₂S and Na₂SO₃. During the reaction, the reactants are present in the liquid phase, and the products are present in both gas and liquid phases, as shown in Equation (22). Therefore, the user-defined function DEFINE_HET_RXN_RATE is utilized to simulate the rate of the heterogeneous reaction.

2.2.5. Multi-Field Coupling

Numerical simulation of energy and mass transfer of hydrogen production via suspended photocatalytic water splitting containing sacrificial agent involves the flow field, radiation field, and reaction field. The multi-field coupling procedure is illustrated in Figure 5. To increase the efficiency of the simulation, the discrete phase model is added on the basis of the stable liquid flow in this study. The distribution of the discrete phase concentration and the initial incident radiation intensity within the simulated domain plays a crucial role in determining the radiation field distribution. Similarly, the distribution of the liquid phase and radiation field within the simulated domain are decisive factors in shaping the reaction field distribution. The radiation distribution in the reactor is solved iteratively, utilizing the output file of the discrete phase concentration distribution and user-defined incident radiation intensity. Combined with the concentration distribution of reactants in the liquid phase, the flow field and radiation field are coupled by the hydrogen production rate formula of the reaction model. During the iterative calculation, the convergence of the calculation is determined by the residual value of each physical variable. In this study, the convergence criterion for the residual values of energy is set below 10⁻⁶, while the convergence criterion for the residual values of other variables is set below 10⁻³.

2.2.6. Grid Division

The simulated domain of the reactor is discretized into 3.0 × 10⁵ hexahedral cells, as shown in Figure 6. In the grid sensitivity analysis, the velocity along the y-axis is calculated for different grid numbers. The results obtained with 3.0 × 10⁵ cells and 6.0 × 10⁵ cells are found to be nearly identical. Therefore, it is concluded that 3.0 × 10⁵ cells are sufficient for this study.

3. Results and Discussion

3.1. Model Validation

The numerical simulation of stirred reactors with particles has been widely researched. To validate the flow model in this study, the particle concentration distribution in a similar reactor with a rotation speed of 900 rpm is selected for comparison [26]. The conditions include the following: the concentration and diameter of the particle are 0.5 g·L⁻¹ and 10 μm; the liquid phase is pure water; the radius of the spherical cap reactor is 40.0 mm; the material density of the particle is 4800 kg·m⁻³; the volume of the liquid phase is 199.8 mL. The normalized particle concentration in the z-direction is shown in Figure 7. The blue scatter points represent the reference values, while the orange scatter points show the simulation values obtained using the flow model in this study. As seen from Figure 7, both the reference values and the simulation values exhibit symmetric about z = 0. The relative errors between the simulated values and the reference values are less than 8.0% for the same abscissa, demonstrating the accuracy of the flow model in this study. The computational fluid dynamics model for reactive multiphase flow involving radiation in a tubular photocatalytic reactor has been validated through experiments [9,12]. The simulation process in this study is consistent with the references. However, unlike the references, this paper considers the presence of a stirring magneton and employs the Lagrange model to analyze the particle motion state.

3.2. Distribution Characteristics of Flow Field

The addition of the stirring magneton in the reactor serves to impart energy and create a favorable flow state. The rotation of the stirring magnet transfers mechanical energy to the reactant and forms a fully turbulent mixing zone near the stirring magneton. The generated high-speed jet will push the reactant to circulate within the reactor. Therefore, the rotating speed of the stirring magneton has a great influence on the flow within the reactor.

3.2.1. Effect of Rotating Speed on Continuous Phase

The nephograms of gas–liquid phase distribution with different rotating speeds on the z = 0 section of the reactor are shown in Figure 8, where red represents the liquid phase, and blue represents the gas phase. As seen from this figure, as the rotation speed increases, the center of the free liquid surface exhibits a concave shape, and the liquid near the wall surface bulges to form a vortex. The depth of the vortex increases with the increase in rotating speed, which are 0.63 mm, 1.31 mm, 3.06 mm, and 6.86 mm, respectively. The relationship between rotating speed and vortex depth can be observed through a fitting curve, which can be expressed as f(s) = 0.27e^0.0032s. It indicates that the vortex depth increases exponentially with the rotating speed. This is because the rotation of the stirring magneton makes the liquid in the reactor rotate around the central axis, leading to the liquid rushing towards the inner wall surface and rising under the centrifugal force, causing a drop in the liquid level in the central region.

The nephograms of velocity distribution inside the reactor with different rotating speeds are shown in Figure 9, Figure 10, Figure 11 and Figure 12. The figures show that the velocity at the central axis of the reactor is minimum on the same cross section, except for the location of the stirring magneton. This is due to the tangential flow caused by the rotation of the stirring magneton, creating a cylindrical rotary region at the central axis. On the right side of the reactor is the flat light window, leading to an asymmetrical distribution of the internal flow field. Additionally, the asymmetric velocity phenomenon in the reactor gradually decreases with increasing rotating speed. The velocity near the reactor wall is relatively low due to the high viscous force and the presence of a boundary layer. The stirring magneton is located at the bottom of the reactor, leading to a large velocity gradient and intense flow at the bottom. As seen from the velocity nephogram of the y = 0 section, the maximum velocity appears at the end face of the stirring magneton. Moreover, the velocity at the end face increases gradually with the increase in rotating speed. The maximum velocities obtained by simulation for the four rotating speeds are 0.29 m·s⁻¹, 0.62 m·s⁻¹, 0.92 m·s⁻¹, and 1.22 m·s⁻¹, respectively.

The circumferential velocity at the end face of the rotor in the rotary stirring system can be described by the following equation:

$$w=\frac{\pi ln}{60},$$

(25)

where n is the rotating speed, and w is the circumferential velocity.

When the rotating speeds are 250, 500, 750, and 1000 rpm, the theoretical circumferential velocities calculated by Equation (25) are 0.26 m·s⁻¹, 0.52 m·s⁻¹, 0.79 m·s⁻¹ and 1.05 m·s⁻¹, respectively. It is observed that as the rotating speed increases, the difference between the simulated circumferential velocity and the calculated value becomes increasingly larger. This is because Equation (25) is appropriate for the axisymmetric cylindrical reactor. However, for the asymmetric spherical reactor studied in this research, Equation (25) could reflect the influence trend of the rotating speed on the circumferential velocity at the end face of the rotor to some extent, but it is not precise.

3.2.2. Effect of Rotating Speed on the Discrete Phase

The catalyst particles are introduced into the reactor through a surface jet source from the bottom. The trajectory of catalyst particles is shown in Figure 13. Different colors are used to represent the movement trajectory of different catalyst particles in the reactor at the same time. It can be seen that when the particles enter the reactor from the bottom, they follow a spiral motion. At first, the catalyst particles spiral upward along the wall, as illustrated in Figure 13a. Upon reaching a certain height, the axial velocity, due to the effect of gravity, causes the particles to spiral downward near the central axis, as shown in Figure 13b. When reaching the bottom of the reactor, they spiral upward along the wall once again, as depicted in Figure 13c, repeating this cycle.The tangential component velocity leads to a circular motion of the particles as they spiral upward from the bottom, keeping them close to the wall. On the other hand, when the particles spiral downward, the velocity at the central axis is lower, causing the particles to approach the central axis.

Figure 14 displays the concentration distribution of catalyst particles on the z = 0 section of the reactor. Particle concentration can be divided into four regions. Region one is the area with y > 12.0 mm and x < 15.0 mm in the reactor. The particle concentration in this region is low due to the tangential flow, which results in a lower velocity at the central axis, weaker turbulence intensity, and poor particle mixing. Region two is the area with 15.0 mm < x < 25.0 mm in the reactor. The particle concentration in this area is also low. The catalyst particles follow an upward–downward–upward spiral movement in the reactor, being close to the wall during the upward movement and close to the central axis during the downward movement. The particle trajectory is less noticeable near x = 20.0 mm, contributing to the low concentration in this region. Region three is the area with y < 12.0 mm and x <15.0 mm in the reactor. The particle concentration in this region is higher than the average value of 0.5 g·L⁻¹. This is due to the insufficient lift force generated via the rotation of the stirring magneton to keep the particles fully suspended in the reactor. Additionally, as the rotating speed decreases, particle sedimentation becomes increasingly evident. Region four is the area with x > 25.0 mm in the reactor. The reactor’s structure significantly impacts the velocity in this region, leading to changes in particle concentration. As seen from the figure, the area with high particle concentration at the bottom of this region decreases with the increase in rotating speed. This is due to the increase in turbulent intensity at the bottom, which leads to a better mixing of particles as the rotating speed increases.

Generally speaking, an increase in rotating speed leads to an increase in the pressure head generated during rotation, as well as an improvement in the energy supplied to the liquid and an increase in turbulence, which is more favorable for particle dispersion. However, it is important to note that the rotating speed cannot be increased without limits. To prevent the liquid from overflowing the reactor, the equation for calculating the maximum rotating speed has been established as follows:

$$n\le \sqrt{\frac{4gM}{D^2}{\left(\frac{30}{\pi }\right)}^2\frac{x^4}{x^2-\ln \left(x\right)-3/4}},$$

(26)

$$x\approx \frac{D}{0.625l_s},$$

(27)

where n is the rotating speed; g is the gravitational acceleration; M is the penetration depth of the stirring magneton; D is the inner diameter of the reactor; l_s is the length of the stirring magneton. As the reactor studied in this paper is not a regular cylinder, the equivalent diameter is taken as the inner diameter of the reactor. The calculation shows that the maximum rotating speed should not exceed 1226.0 rpm to avoid liquid overflow in the reactor.

3.3. Distribution Characteristics of Radiation Field

The radiation model in this study is proposed under the condition of ignoring the scattering effect of the reaction solution, so the radiation in the reactor gradually decreases with the penetration depth of light. The nephograms of radiation distribution on the section in the y-direction of the reactor under different incident radiation intensities when the particle concentration is 0.5 g·L⁻¹ are shown in Figure 15. When the particle concentration is held constant, the penetration depth of light in the reactor increases with an increase in the incident radiation intensity. For instance, when the incident radiation intensity is 250 W·m⁻², the penetration depth of light higher than 100 W·m⁻² in the x-direction of the reactor is 9.6 mm. The values of the penetration depth of light increase to 17.2 mm, 21.6 mm, and 24.4 mm with increasing incident radiation intensity. Combined with Equation (19), it is observed that when the photocatalyst concentration is constant, the exponential part on the right side of the radiation field control equation remains unchanged. For each discrete grid in the reactor, as the incident radiation intensity increases, the outgoing radiation intensity also increases, leading to an increase in the penetration depth of the reactor in the x direction. The relationship between the incident radiation intensity (independent variable i) and penetration depth (dependent variable) can be expressed by the fitting curve as g(i) = 10.73ln(i) − 49.59. This relationship is similar to the radiation relation formula I_trans/I_input = 0.47 + 0.53exp(−0.005t) with TiO₂ as particles found in the literature, where t is the penetration depth [34]. It can be seen that there is a positive logarithmic relationship between the penetration depth and the incident radiation intensity. The coefficients in the relationship vary with the change of particles.

The nephograms of radiation distribution on the section in the y-direction of the reactor under different particle concentrations when the incident radiation intensity is 1000 W·m⁻² are shown in Figure 16. It can be seen that when the incident radiation intensity is constant, the penetration depth of light decreases as the particle concentration increases. For instance, when the particle concentration is 0.25 g·L⁻¹, the penetrating depth of light greater than 100 W·m⁻² in the x-direction of the reactor is 38.8 mm. The penetration depth decreases to 24.4 mm, 19.2 mm, and 16.4 mm as the particle concentration increases. The relationship between the particle concentration (independent variable c) and penetration depth (dependent variable) can be expressed by the fitting curve h(c) = −16.38ln(c) + 15.01. Thus, the penetration depth has a negative logarithmic relationship with the particle concentration.

3.4. Distribution Characteristics of Reaction Field

By combining the flow field and radiation field, hydrogen production at different times can be calculated using the reaction rate formula. Figure 17 displays the hydrogen concentration distribution on the z = 0 section at different reaction times with a particle concentration of 0.5 g·L⁻¹ and an incident radiation intensity of 1000 W·m⁻². As seen from the figure, the maximum hydrogen concentration in the reactor appears above the rightmost side of the free liquid surface. This is due to the fact that the radiation is incident from the negative direction of the x-axis, causing the radiation intensity in the reactor gradually decreases in that direction. The higher radiation near the light window of the reactor leads to the highest hydrogen concentration in this area, highlighting the significant impact of radiation distribution on hydrogen production. The total amounts of hydrogen generated in the reactor are 1.04 × 10⁻³ mol and 1.35 × 10⁻³ mol when the reaction times are 1.0 s and 2.0 s, respectively.

4. Conclusions

In this paper, a mathematical model for hydrogen production by suspended photocatalytic water splitting under the simulated light source in a magnetic stirring reactor is established. The flow model is simulated using the volume of fluid and discrete phase method. The radiation absorption model and reaction kinetics are integrated into the simulation to couple the flow field, radiation field, and reaction field. The main conclusions of the numerical simulation are summarized as follows:

(1)

The rotating speed of the stirring magneton in the reactor has a significant influence on the flow field. The rotation of the stirring magneton generates a vortex in the central axis area of the reactor, with the relationship between the depth of the vortex f(s) and the rotating speed of the magneton s described as f(s) = 0.27e^0.0032s. Additionally, the movement trajectory of the catalyst particles is also analyzed and determined.

(2)

The incident radiation intensity and particle concentration play a crucial role in determining the radiation distribution within the reactor. The relationship between the penetration depth of radiation g(i) and the incident radiation intensity i is described as g(i) = 10.73ln(i) − 49.59. The relationship between the penetration depth of radiation h(c) and the particle concentration c is given as h(c) = −16.38ln(c) + 15.01. These findings emphasize the importance of considering the impact of incident radiation intensity, incident direction of light, and particle concentration in the design of the reactor.

(3)

The location of the maximum hydrogen concentration in the reactor is determined via the incident radiation intensity and direction, as well as the particle concentration. When the light source is a single-sided parallel incident, the maximum hydrogen concentration is observed near the light window of the reactor above the free liquid surface. The radiation distribution in the reactor has a significant influence on hydrogen production and determines the distribution law of hydrogen concentration in the reactor.