Comprehensive Optimization of the Dispersion of Mixing Particles in an Inert-Particle Spouted-Bed Reactor (IPSBR) System

Effective gas dispersion and liquid mixing are significant parameters in the design of an inert-particle spouted-bed reactor (IPSBR) system. Solid particles can be used to ensure good mixing and an efficient rate of mass and heat transfer between the gas and liquid. In this study, computational fluid dynamics (CFD) coupled with the discrete phase model (DPM) were developed to investigate the effect of the feed gas velocity (0.5–1.5 m/s), orifice diameter (0.001–0.005 m), gas head (0.15–0.35 m), particle diameter (0.009–0.0225 m), and mixing-particle-to-reactor-volume fraction (2.0–10.0 vol.%) on the solid mass concentration, average solid velocity, and average solid volume fraction in the upper, middle, and conical regions of the reactor. Statistical analysis was performed using a second-order response surface methodology (RSM) with central composite design (CCD) to obtain the optimal operating conditions. Selected parameters were optimized to maximize the responses in the middle and upper regions, and minimize them in the conical region. Such conditions produced a high interfacial area and fewer dead zones owing to good particle dispersion. The optimal process variables were feed gas velocity of 1.5 m/s, orifice diameter of 0.001 m, gas head of 0.2025 m, a particle diameter of 0.01 m, and a particle load of 0.02 kg. The minimum average air velocity and maximum air volume fraction were observed under the same operating conditions. This confirmed the novelty of the reactor, which could work at a high feed gas velocity while maintaining a high residence time and gas volume fraction.


Introduction
Spouted beds have several advantages and potential applications compared to moving beds. They can handle granular particles with a wide size distribution range, provide an efficient mixing between the solid and gas phases and significant rates of heat and mass transfer [1][2][3][4][5][6][7][8]. Inert particles spouted bed reactor (IPSBR) is a spouted bed in which a gas jet is injected at the bottom of a conical vessel containing liquid and inert particles [9]. El-Naas et al. [10,11] developed and investigated a novel IPSBR to deal with both CO 2 capture and reject brine management. Their results indicated that the CO 2 capture efficiency reached up to 97.7% at a gas-to-liquid volume ratio of 123. They also concluded that the interfacial area between the gas and liquid was enhanced using inert gas particles. This was due to the circular motion of the particles inside the reactor, which enhanced the CO 2 capture efficiency and ion removal [6]. The hydrodynamics of the same IPSBR were studied concentrations, they detected that the average gas holdup was decreased remarkably by increasing the size and concentration of the particles. Many studies [20,21] are in agreement with Rabha et al. [19], and it has been reported that by increasing the particle size, the gas holdup decreased. The effect of particle size (60-150 µm) and particle volumetric concentration (0% to 50%) on gas holdup was investigated by Ojima et al. [22]. Their results indicated that the bubble coalescence increased by decreasing the particle size and increasing the concentration of solids; thus, the gas holdup was decreased.
Gholamzadehdevin et al. [23] studied the hydrodynamic performance of an activatedsludge bubble column. They combined CFD (Eulerian-Eulerian) and the design of the experiment (full-factorial design) to obtain the influence of superficial gas velocity, tracer injection position, and sparger type on the mixing time. Their findings indicated a strong connection between the superficial gas velocity and the location of the tracer injection on the mixing time. Moreover, the hydrodynamic performance was improved by a modified star-shaped gas sparger. Central composite design was used to optimize the operating conditions of the bubble column (piperazine-H 2 O-CO 2 ) [24]. The results indicated that the maximum value (97.7%) of CO 2 capture could be achieved with a piperazine concentration of 0.162 M, solution flow rate of 0.502 L/h, CO 2 flow rate of 2.199 L/min, and stirrer speed of 68.89 rpm. They mentioned that the p-values of the dependent parameters were under 0.05, which reflected their significance and importance. Song et al. [25] developed a new cryogenic process based on free-piston Stirling coolers to enhance the CO 2 removal efficiency. They used response surface methodology (RSM) to optimize the operating conditions. They concluded that the process could capture up to 95.20% of the CO 2 under the optimal conditions of 2.16 L/min flowrate, temperature of −18 • C, and operation time of 3.9 h.
In light of the literature mentioned above, much helpful information has been reported on the impact of operating conditions on the hydrodynamics and bubble behavior of various bubble column reactors. The literature also reveals the importance of optimizing the operating parameters of the reactor, as they have an observable potential for enhancing the design and overall performance of the process. However, there is still a lack of information on optimizing the most effective hydrodynamic parameters (feed gas velocity, orifice diameter, gas head, diameter of mixing particles and mixing particle load) of an IPSBR. Optimization of various operating parameters for the IPSBR is still under investigation, and it is essential for the design, pilot plant scale-up, and industrial operation of the reactor. The current work focuses on understanding the interactions between the most effective factors resulting in good particle distribution throughout the reactor. Accordingly, more responses are introduced to the response surface design to ensure more accurate optimum conditions for real CO 2 capture reactions using the IPSBR system. Therefore, the main objective of the present study is to combine the CFD model with RSM analysis to study the impact of these variables on the mass concentration, average velocity, and average volume fraction of mixing particles throughout the reactor. The combination of the CFD model and RSM analysis enables the prediction of the best conditions at which the solid particles are dispersed perfectly to obtain high gas hold-up, greater interfacial area, and longer residence time.

Simulation Configuration
ANSYS Fluent 18.0 (CFD) combined with Eulerian multiphase was used to model the axisymmetric of a 2D IPSBR, which was working in a semi-batch mode. The simulated reactor configuration was sourced from a previous study, as shown in Figure 1 [12]. The simulation was conducted with three phases (air, water, and solid). The solid phase consisted of polymethyl methacrylate particles (average density of 1020 kg/m 3 and average particle size of 0.013 m). First, the simulated IPSBR reactor was filled with a specific amount of water. An air space (gas head) was retained at the top of the reactor because the water level would rise once the gas started to flow [12,13]. Mohammad et al. performed a Processes 2021, 9,1921 4 of 25 mesh independence study [12] and concluded that a grid size of 0.0020 (cell count = 3287) would be sufficient and could be applied in this study. An intensive mesh study was conducted by a previous research study [12,13]. The simulation runs were all converged. Seven meshes with mesh sizes ranging from 0.0022 to 0.0016 were tested (22,845, 28,568, 32,871, 36,112, 39,068, 41,386, and 42,372). The air velocity convergence was achieved at a mesh size of 0.0016, which was only about 1% different from the previous mesh size. The oscillation-shaped solution was found to be mesh-independent, with no significant differences between the last four mesh sizes. As a result, in the 2D column runs, the grid size of 0.0020 was used, which results in a 2.6% difference from the finest mesh.

CFD Model
The Eulerian multiphase model in the commercial CFD software package Fluent 18.0 was used to model and investigate the hydrodynamic flow behavior of the IPSBR system under transient and gravitational acceleration operational conditions. All runs were carried out under isothermal and unsteady-state behavior. In addition, gas was considered to have constant properties, which was considered as the secondary phase and the liquid as the main phase. The top of the reactor was open to the atmosphere. Bubbles of equal diameter (0.001-0.005 m) were created at the bottom of the reactor through a spherical orifice with a specific diameter. A very small step size of 0.00001 s was employed with 20 iterations for each time step. The convergence criterion was identified to be 10 −3 for the relative error between two sequential iterations. To solve the numerical model, a finite volume method was employed. Turbulence was simulated using the standard k-ε mixture model using wall functions (low-Re k-ε near the wall). The discrete phase model (DPM) was implemented to consider the impact of mixing particles. CFD describes the fluid phase by solving the Navier-Stokes equations. The particle velocities and forces are calculated in DEM by solving Newton's equations of motion. The time it takes to record converged results is 20 s or longer. The transient state's time-averaged profile was previously calculated at 25 and 30 s time intervals (averaged data over 5 s of flow was chosen) [12]. In this work, the time-averaged profile for the transient state has been selected to be over 60 to 65 s time intervals to insure more stability of the quasi steady state condition. The standard error Processes 2021, 9, 1921 5 of 25 deviation was in the range of (0.00429-0.00961) [13]. The equations are listed in Table 1. More information and details can be found in the preceding research work [12,13]. Governing equations Continuity equation is given by [12,13] Momentum equations is given by [26] ∂ The stress tensor (τ k,ij ) is given by [27] Transport Equations for the Standard k-ε Mode [28] Turbulence Kinetics Energy [29] Dk Turbulence Energy Dissipation Rate [28] D The eddy viscosity (µ t ) is obtained as Where σ k = 1.0, σ = 1.3, C 1 = 1.44, C 2 = 1.92, C µ = 0.09 [27] (7) Newtonian equation of motion (to consider the effect of mixing particles) .81 m/s 2 For the Morsi and Alexander model [29]: R 2 e (11) where a 1 , a 2 , and a 3 are constants that are employed for smooth spherical particles over several ranges of Reynolds number [30], which are defined as following:

Response Surface Methodology (RSM) Design
Central composite design (CCD) is a frequently used form of RSM that affords the opportunity to study interactions between the independent process parameters and dependent process responses [31][32][33][34][35][36][37][38]. In this study, a full-factorial CCD design was carried out using Minitab 19.0 to investigate the influence of input parameters (i.e., feed gas flowrate, orifice diameter, gas head, mixing particle diameter, and mixing particles load) on the responses (i.e., mass concentration, average velocity, and average volume fraction of the mixing particles) throughout the reactor. Table 2 listed the levels of the independent parameters and the consistent variation levels used to optimize the process conditions. Thirty-two runs with five levels of each parameter were designated. The conditions of each run were entered into the ANSYS software, and all of the data attained were analyzed using the Minitab software (RSM).  , which represents the dimensionless diameter of the IPSBR system, where R = 0.0 and R = 1.0 are the IPSBR system walls in the 2D model. The results were observed at a constant orifice diameter of 0.003 m, gas head of 0.25 m, mixing particle diameter of 0.0135 m, and total mass of mixing particles of 0.06 kg. It is clear from the figures that the effect of a feed gas velocity of 1.5 m/s had a significant effect on the average solid velocity distribution in the middle region and reached a maximum value of 0.014 m/s. At a feed gas velocity of 1 m/s, the influence was observed more in the conical region compared with the middle and upper regions. By contrast, at a minimum feed gas velocity of 0.5 m/s, the solids existed in almost all of the regions of the reactor. Due to the no-slip boundary condition, the solid velocity was almost zero at the wall for all values of feed gas velocity. It is worth noting that the solid particle velocity was detected in the three system zones at low gas velocity of 0.5 m/s, as shown in Figure 2a. However, when the particle velocity was increased to 1 m/s and 1.5 m/s, most of the detected velocities were zero; this could be explained by accelerating the solid particles in the studied surface zone and thus the disappearance/absence of these solid particles, as shown in Figure 2b,c. velocity of 0.5 m/s, the solids existed in almost all of the regions of the reactor. Due to the no-slip boundary condition, the solid velocity was almost zero at the wall for all values of feed gas velocity. It is worth noting that the solid particle velocity was detected in the three system zones at low gas velocity of 0.5 m/s, as shown in Figure 2a. However, when the particle velocity was increased to 1 m/s and 1.5 m/s, most of the detected velocities were zero; this could be explained by accelerating the solid particles in the studied surface zone and thus the disappearance/absence of these solid particles, as shown in Figure 2b,c. (b) Effect of Orifice Diameter on Average Solid Velocity Figure S1a,b depicts the effect of two different orifice diameters (0.003 and 0.005 m) on the average solid velocity in the three regions. The results were extracted for a constant feed gas velocity of 1 m/s, gas head of 0.25 m, mixing particle diameter of 0.0135 m, and mixing particle load of 0.06 kg. It was found that increasing the orifice diameter from 0.003   Figure S1a,b depicts the effect of two different orifice diameters (0.003 and 0.005 m) on the average solid velocity in the three regions. The results were extracted for a constant feed gas velocity of 1 m/s, gas head of 0.25 m, mixing particle diameter of 0.0135 m, and mixing particle load of 0.06 kg. It was found that increasing the orifice diameter from 0.003 to 0.005 m reduced particle velocity throughout the upper region and made it easier to detect, as shown in Figure S1b. The mixing particles, on the other hand, were not detected when a smaller orifice diameter was used, as shown in Figure S1a. It is clear from Figure S1b that the particle velocity was not detected in the conical and middle regions. Whereas, in the upper region it reached a peak value of almost 0.1 m/s. This is also connected to the gas velocity inside the reactor, which increased inside the reactor at higher orifice diameters, owing to the creation of larger bubble sizes. Because of that, the particles would not be well distributed in a uniform profile throughout the reactor at higher orifice diameters [36].
(c) Effect of gas head on average solid velocity The effect of gas head on the average solid velocity distribution was studied at a constant feed gas velocity of 1 m/s, 0.003 m orifice diameter, mixing particle diameter of 0.0135 m, and total mass of mixing particle of 0.06 kg in all three regions. In the conical region, it was expected to see high dispersion and turbulence of the particle distribution because that is the injection zone for those particles. Furthermore, the cross-sectional area of the conical region was smaller than those of the middle and upper regions. Therefore, a clear behavior was not detected in conical region. It can be seen in Figure 3a,b that more gas head resulted in an increase in the particle velocity in the middle and upper regions. From the simulated data of a gas head value of 0.15 m, it was found that the mixing particles at R = 0.2 to 0.6 had an average velocity of around of 0.04 m/s. However, in the same section with a gas head value of 0.35 m, the average solid velocity approached a value of around 0.06 m/s. A lower gas head meant a higher liquid level. As the liquid level increased this exerted more pressure on the mixing particles; therefore, a lower particle velocity was observed.
(d) Effect of mixing particle diameter on average solid velocity The influence of the mixing particles diameter of 0.0135 and 0.0225 m at a constant feed gas velocity of 1 m/s, 0.003 m orifice diameter, gas head of 0.25 m, total mass of mixing particle of 0.06 kg on the average solid velocity in all three regions is illustrated in Figure S2a,b. It was observed that at high mixing particles diameter (0.0225 m), some of the particles settled in the conical section of the reactor and show zero velocity. It was also found that by increasing the size of the mixing particles, their velocity decreased remarkably in the middle and upper regions. As illustrated in Figure S2b, the maximum particle velocity reached was approximately 0.037 m/s in the upper region, which is still a relatively low value. Decreasing the particles diameter as shown in Figure S2a resulted in an increase in the solids velocity as detected in the conical region.
(e) Effect of mixing particle load on the average solid velocity Figure 4a,b illustrates the solid velocity profile obtained along the contact system diameter for mixing particle loads of 0.02 and 0.06 kg in all three regions. The effect was examined at a constant feed gas velocity of 1 m/s, orifice diameter of 0.003 m, gas head of 0.25 m, and mixing particle diameter of 0.0135 m. As expected, the trend of the velocity distribution was not well detected in the conical region because the disturbance of the mixing particles was very high. From the simulated data for the upper and middle regions for all of the particle loads, it was observed that the liquid level and gas holdup could be increased by increasing the particle load. On the other hand, it was also found that increasing the particle load would increase the particle velocity, which would also increase the collisions between particles. Figure 4b shows that the particle velocity reached its maximum value of almost 0.058 m/s at R = 0.43 in the middle region, whereas at a particle load of 0.02 kg (Figure 4a), the particle velocity was almost zero at the same lateral position.
Processes 2021, 9,1921 9 of 25 The increase in particle velocity with particle load will negatively affect the gas holdup. Therefore, the optimization of this parameter is discussed in Section 4.6, which presents the best particle load to enhance the overall performance of the reactor.  the collisions between particles. Figure 4b shows that the particle velocity reached imum value of almost 0.058 m/s at R = 0.43 in the middle region, whereas at a par of 0.02 kg (Figure 4a), the particle velocity was almost zero at the same lateral The increase in particle velocity with particle load will negatively affect the gas Therefore, the optimization of this parameter is discussed in Section 4.6,  gas head of 0.25 m, mixing particle diameter of 0.0135 m, and total mass of mixing particles of 0.06 kg, it was found that the maximum solid volume fraction was 0.043 in the conical region at a feed gas velocity of 1.0 m/s. When the feed gas velocity was reduced to 0.5 m/s, the majority of the particles congregated at the reactor's bottom. This is evident in Figure 5a, which shows a very low volume fraction in almost all of the regions. This is due to the fact that all of the particles were concentrated at the bottom of the reactor and moving very slowly.

Effect of Independent Parameters on the Average Solid Volume Fraction
(a) Effect of feed gas velocity on the solid volume fraction By comparing the effect of a feed gas velocity of 1.5 m/s with the effects of velo of 1 and 0.5 m/s in all of the column regions under a constant orifice diameter of 0.00 gas head of 0.25 m, mixing particle diameter of 0.0135 m, and total mass of mixing par of 0.06 kg, it was found that the maximum solid volume fraction was 0.043 in the co region at a feed gas velocity of 1.0 m/s. When the feed gas velocity was reduced to 0.5 the majority of the particles congregated at the reactor's bottom. This is evident in Fi 5a, which shows a very low volume fraction in almost all of the regions. This is due t fact that all of the particles were concentrated at the bottom of the reactor and mo very slowly. The effect of the orifice diameter was examined under a feed gas velocity of 1 m/s, mixing particle diameter of 0.0135 m, and mixing particle load of 0.06 kg and gas head of 0.003. It was observed that as the orifice diameter was reduced from 0.005 m to 0.001 m, the volume fraction of gas phase decreased in all regions and accordingly the superficial velocity of the created bubbles inside the reactor is decreased. More circulation and eddies were also observed because of increasing gas residence time inside the system. This enhanced the surface contacted between the three phases along the reactor.
(c) Influence of gas head on solid volume fraction (d) Effect of mixing particle diameter on solid volume fraction The effect of particle size could not be separated from the effect of the solids concentration. It was concluded that deceasing the particle size from 0.0225 to 0.0135 could increase the solid volume fraction in the conical region. Where increasing the particle size resulted in a high precipitation of solids at the bottom of the reactor. This is illustrated in Figure S3a,b for a particle size of 0.0135 and 0.0225 m under a constant feed gas velocity of 1 m/s, orifice diameter of 0.003 m, gas head of 0.25 m, and total mass of mixing particle of 0.06 kg in all three regions. Selecting the proper particle diameter is very important for ensuring a uniform distribution of particles along the reactor. This will enhance the interfacial areas between the gas and liquid, and also improve the gas holdup and thus the overall performance of the reactor.
(e) Effect of mixing particle load on solid volume fraction It was noted that when the particle load increased from 0.02 kg to 0.1 kg, the solid volume fraction increased from 3 × 10 −5 to 0.045, especially in the conical region. This enhanced the contact area between the gas and liquid in the middle and upper regions. However, in the conical region, increasing the volume fraction enhanced the bubble coalescence and accordingly decreased the gas holdup owing to the small cross-sectional area compared to the middle and upper regions, as was also observed by other reporters [39].

Influences of Factors on the Mass Concentration of Mixing Particles
Based on the RSM responses, a polynomial equation (second order) in terms of dependent and significant independent coded variables for the three regions is given by Equations (12) Table 3 and Table S1 show the results obtained from the regression model for the solid mass concentration in the middle and upper regions, respectively. The p-values were low (<0.05) for most of the interactions, which is an indication of the high significance of the model as reported by others [38][39][40]. Moreover, the regression model had a high value of R2, which was 99.68% for the middle region and 98.50% for the upper region.
(c) Influence of gas head on solid volume fraction Figure 6a,b demonstrates an overview of the change in the solid volume fraction wit respect to the gas head (0.15 and 0.35 m) at a feed gas velocity of 1 m/s, orifice diameter o 0.003 m, mixing particle diameter of 0.0135 m, and total mass of mixing particle of 0.06 k in all three regions. It was observed that by decreasing the gas head (more gas holdup from 0. 35    Pareto plots for the middle, upper, and conical regions were used to detect the most significant factors and interactions that affect the solid mass concentration. It was found that most of the main effects and two-factor interactions are statistically significant at a 5% significance level. It is also worth mentioning that the normal plots indicate that the residuals fall approximately in a straight line, which provides further confirmation of the acceptability of the model.

Influences of Factors on Average Solid Volume Fraction
The ANOVA analysis for the average solid volume fraction in the middle region is provided in Table S2. It shows that all of the linear and quadratic coefficients are significant (p > 0.05). The exception is the two-way interactions between the orifice diameter and mixing particle load. This is according to the p-values, which are greater than 0.05. Moreover, the R2 had a high value of 0.9961. The adequacy of the results for the conical, middle, and upper regions was confirmed by residual analysis, and the transformed data in the normal probability results were found to be very close to the normality distribution curves. Furthermore, the residuals-versus-order fluctuations revealed that there is no pattern, indicating independence. All of the collected data showed a good distribution of the CFD-simulated data and predicted values, with no differences observed.

Influences of Factors on Average Mixing Particle Velocity
The second-order polynomial regression equations indicating the mixing particle velocity are expressed by Equations (15) The Pareto analysis, which reveal the most important factors influencing the average mixing particle velocity in all three regions, was carried out and exposed that most of the parameters had a great impact on the responses. These results were confirmed from the p-value, which is less than 0.05 for most of the factors, and the R2, which is very close to 1. The residual normal probability analysis for the solid velocity indicates that all of the residuals fall around a line, with no obvious outliers or scattering.

Response Optimizer
In this study, three responses (mass concentration of mixing (R1-R3), average solid volume fraction (R4-R6) and average mixing particle velocity (R7-R9)) are calculated. It is desirable to find the maximum values of these three responses in the middle and upper regions. However, in the conical region it is desirable to find their low values. In the conical region, most of the particles will have settled and it is important to minimize the concentration of solids and turbulence in this region by setting the velocity and concentration to a minimum (Table S3). This will reduce the coalescence between the bubbles and increase the gas holdup in the region. In order to ensure a good distribution of particles in the middle and upper regions, the values of solid velocity and concentration should be at their maximum. This will increase the interfacial areas between particles, gas, and liquid. It will also enhance the residence time of the process.
In our previous work [13], the desirable upper, middle, and conical air volume fractions (R10-R12) were considered and targeted to be maximum to achieve maximum interfacial area concentration, whereas the desirable upper, middle, and conical air velocities (R13-R15) were studied and targeted to be minimal to achieve maximum gas residence time. The optimizer prediction shows that the lowest average air velocity and highest average air volume fraction in the conical, middle, and upper regions were achieved at feed-gas velocity of 1.5 m/s, gas head of 0.164 m, orifice diameter of 0.001 m, mixingparticle diameter of 0.0225 m, and total mass of mixing-particle of 0.02 kg as demonstrated in Table S4.
In this work, all preceding responses (R1-R15) were combined together and optimized using the Minitab response optimizer conferring to the stated targets to have a wider indication for the hydrodynamic performance for the IPSBR system. The optimal parameters were found to be a feed gas velocity of 1.5 m/s, orifice diameter of 0.001 m, gas head of 0.202 m, mixing particle diameter of 0.0107 m, and total mass of mixing particles of 0.02 kg.
The optimization outcomes are listed in Table 4. Remarkably, for the same optimal parameters, a high air volume fraction and low air velocity were attained. This result reflects the novelty of the IPSBR, which can operate at high feed gas velocities while maintaining a high gas residence time and high gas volume fraction. The overall desirability is assessed as one, which specifies that the response optimization is satisfactory [41,42]. When comparing Table S4 and Table 4, operational conditions in the same range were observed, confirming the same physical concept of achieving high contact surface area, gas residence time, and mixing particle distribution at the specified operational conditions range.

Validation of the CCD Model
The optimal operational conditions of feed gas velocity 1.5 m/s, orifice diameter 0.001 m, gas head 0.202 m, mixing particle diameter 0.0107 m, and mixing particle load of 0.02 kg were entered into the ANSYS software model to collect all predicted responses (R1-R15). The predicted responses from the regression model were found to agree with those obtained from the simulated data, which confirmed the model's validity and adequacy, as shown in Table 5. The predicted and simulated responses are very close and within the 95 percent confidence interval, indicating that the model can predict the hydrodynamic performance of the simulated IPSBR system under various operational conditions. Solid velocity, mass concentration, and eddies viscosities were plotted to visualize the effect of the optimum conditions on the dispersion of the mixing particles throughout the system region, as shown in Figure 7a at feed gas velocity of 0.5 m/s, orifice diameter of 0.003 m, gas head of 0.25 m, mixing particle diameter of 0.018 m, and total mass of mixing particle of 0.04 kg. These results were compared at the optimum condition, Figure 7b, with a feed gas velocity of 1.5 m/s, orifice diameter of 0.001 m, gas head of 0.202 m, mixing particle diameter of 0.0107 m, and total mixing particle mass of 0.02 kg. It was clear that the optimized condition exhibits uniform eddies viscosity distribution, velocity and mass concentration among the three regions, and thus better mixing particle dispersion is expected. In addition, the swirling motion, which is detected in Figure 7b, affects the particle movement as it pushes particles in fluid movement direction. This may lead to larger residence time for the fluid due to the circulation zone [43]. The contours show that increasing the feed gas velocity (to nearly 1.5 m/s) and decreasing the orifice diameter (to 0.001 m) are important in maximizing the solid velocity and volume fraction in the middle region. Figure S4a,b demonstrates that the maximum solid volume fraction in the conical and upper region can be obtained at a gas head value of nearly 0.24 m and at a low particle diameter value.

Contour Plots
From Figure 9a, it is clear that the feed gas velocity and particle load had a significant impact on the solid velocity in the middle region of the reactor. By increasing the feed gas velocity and decreasing the particle load can maximize the solid velocity at this region. However, at the same parameter conditions, the solid volume fraction reached its minimum value in the conical region as shown in Figure 9b.
As mentioned in Section 4.1 (c), the effect of particle dispersion on turbulence in the three system zones is significant. Figure 10a-c depicts the change in turbulent kinetic energy as a function of feed gas velocity and total mass of the mixing particles. The turbulent kinetic energy was found to be significantly higher in the conical and upper regions than in the middle region. Since particles are injected into the conical region and gas bubbles are generated through an orifice, large turbulences and disturbances are likely to occur. Similarly, high turbulences can be expected in the upper region where bubbles exit the liquid phase and enter the gas outlet from the reactor's top.

Validation of the CFD Model
Many studies have investigated the hydrodynamics of conical-base spouted bed based on 2-dimentional axisymmetric geometry assumption using the CFD-DEM coupling approach [4,[44][45][46][47][48], and based on the conventional gas-solid or gas-liquid spouted bed phenomena. The novelty of the current research is to investigate the performance of inert-particle spouted bed reactor (IPSBR) system with three-phase mixture (liquid, gas, and solid).
It is important to note that the fluid movement simulation obtained by solving the Navier-Stokes equation is always correct for the fluid domain for the given operating conditions. The location of the solid particles may vary based on its location in the axisymmetric plan. Accordingly, the results obtained from the studied model are valid if the particles are located in the radial location. However, the particles may not always be present in the same location around the axisymmetric plan, then caution should be stressed that the axisymmetric assumption is valid for the fluid domain and not for the particles in the Lagrangian model [49,50]. An experimental validation for the tested model would show us the deviation on simulated data based on the previous assumption and will be compared with the reported data.     A laboratory-scale IPSBR system with the same geometry as that considered in previous studies [13,14] was utilized to validate the Eulerian model. Spherical acrylic particles with a diameter of 0.005 m was used. A gas orifice (ID of 0.002 m) is located at the center of the bottom plate for gas injection. The validation runs considered the presence and absence of mixing particles. A high-speed digital video camera was used to collect measurements reflecting bubbles characteristics, such as the velocity, trajectories, and diameter. The bubble velocity and volume fraction was measured in the radial direction and vertically using Photron FASTCAM Viewer software (PFV) Ver. 3282. Figure 11 shows the simulated and experimental results at specific height from the gas inlet (50 mm) for gas velocity and volume fraction in the two cases: with and without mixing particles. The detected alteration in velocities and volume fractions values are in good agreement with the CFD-simulated data and with an acceptable deviation as listed in Table 2. Still, as reported by Ekambara et al. [31], the 3D model is the only model that could closely predict the experimental data for the majority of the three-phase system. However, the proposed 2D model showed a good and simple estimation with an acceptable deviation between the simulated and experiment data (5-10%), as confirmed in Table 6. Table 7 summarizes the  main findings from the current research and compares them with validated results from  the literature. the CFD-simulated data and with an acceptable deviation as listed in Table 2. ported by Ekambara et al. [31], the 3D model is the only model that could clos the experimental data for the majority of the three-phase system. However, the 2D model showed a good and simple estimation with an acceptable deviatio the simulated and experiment data (5-10%), as confirmed in Table 6. Table 7 su the main findings from the current research and compares them with valida from the literature.

Conclusions
RSM was implemented to assess the effects of feed gas velocity, orifice diameter, gas head, mixing particle diameter, and particle load on the average solid volume fraction, mass concentration of mixing particles, and average mixing particle velocity, which was produced from the CFD results of the novel IPSBR system. The best conditions were a feed gas velocity of 1.5 m/s, orifice diameter of 0.001 m, gas head of 0.2025 m, particle diameter of 0.0107 m, and particle load of 0.02 kg. A good distribution of particles throughout the reactor was observed under these conditions. It was also observed that the maximum air volume fraction and minimum air velocity from the previous work occurred inside the reactor under the same feed gas velocity of 1.5 m/s, orifice diameter of 0.001 m and total mass particles of 0.02, and with very similar gas head and particle diameter values. This was a significant confirmation of the novelty of the reactor, as it could achieve a high gas holdup, high residence time, and good particle distribution under the maximum feed gas velocity value of 1.5 m/s.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/pr9111921/s1, Figure S1: Effect of orifice diameter ((a) 0.003 m, and (b) 0.005 m) on the average solid velocity at a constant feed gas velocity of 1 m/s, gas head of 0.25 m, mixing particle diameter of 0.0135 m, and total mass of mixing particle of 0.06 kg in the conical, middle, and upper regions, Figure S2: Effect of mixing particle diameter ((a) 0.0135 m, (b) 0.0225 m) on average solid velocity at a constant feed gas velocity of 1 m/s, orifice diameter of 0.003 m, gas head of 0.25 m, and total mass of mixing particle of 0.06 kg in the conical, middle, and upper regions, Figure S3: Effect of mixing particle diameter ((a) 0.0135 m and (b) 0.0225 m) on the average solid volume fraction at a constant feed gas velocity of 1 m/s, orifice diameter of 0.003 m, gas head of 0.25 m, and total mass of mixing particle of 0.06 kg/s in the conical, middle, and upper regions, Table S1: ANOVA analysis for average upper particle mass concentration, Table S2: ANOVA analysis for average middle solid volume  fraction, Table S3: Optimization setting for particle mass concentration, solid volume fraction, solid velocity, air velocity, and air volume fraction, Table S4: Optimum conditions and fitted responses