Modification of Meso-Micromixing Interaction Reaction Model in Continuous Reactors

Abstract

The yields of chemical reactions are highly dependent on the mixing pattern between reactants. Herein, we report the modification of a meso-micromixing interaction reaction model which is applied in batch reactors by leveraging the flow characteristics in the continuous reactors. Both experimental and model-predicted yields were compared using the classical Villermaux–Dushman method in a self-designed split and recombination reactor. This modified model significantly reduced the error in predicted product yields from approximately 15% to within 3%, compared to a model containing the micromixing term only. The effects of flow rates and reactor structure parameters on mixing performance were analyzed. We found that increasing flow rates and the degree of twist in the mixing element’s grooves, as well as decreasing the cross-sectional area of grooves, improved mixing performance. The optimization of reactor flow rates and structural parameters was achieved by combining Gaussian process regression and Bayesian optimization with the modified model. This approach provided higher target product yields for consecutive reactions, while simultaneously achieving a lower pressure drop in the reactor. Corresponding combinations of reactor parameters were also identified during this process. Our modified model-based optimization methodology can be applied to a diversity of reactors, serving as a reference for the selection of their structure and operational parameters.

1. Introduction

Consecutive reactions are extensively employed in the chemical process industry (CPI), such as in the synthesis of pharmaceuticals and high-performance polymers [1,2,3]. Studying consecutive reaction yields is crucial for assessing the reactions’ performance and designing synthetic routes [4]. In cases where the reaction rate is less than or equal to the mixing rate, mixing within the reactor governs the contact pattern of the various species involved, thereby influencing their distributions. Consider the following consecutive reactions:

$$\begin{gathered} A+B\stackrel{k_1}{\to }R \\ A+R\stackrel{k_2}{\to }S \end{gathered}$$

(1)

These reactions are second-order, where the yield of the desired product S is:

$$Y_S=\frac{2c_S}{c_{A0}}$$

(2)

Proper control of mixing quality is essential to maximize the yield of the target product S when $k_1>k_2$, and has been extensively studied [5,6,7,8,9,10]. Christy et al. investigated the reaction between 1,2 diphenyl ethane (B) and NO₂BF₄ (A) [8]. They found that the yield of the polynitro product (S) decreased from 78% to 54% upon increasing the stirring intensity then improving mixing uniformity. This decrease in yield can be attributed to different distribution states of the nitrifying agent (A) in the reactor [11]. Enhancing mixing homogeneity resulted in a uniform distribution of A and a rapid decrease in its local concentration within the reactor after the formation of R. This can lead to the continuous accumulation of R, less subsequent reactions taking place, and, finally, a decrease in $Y_S$. However, these studies merely qualitatively evaluated the effect of mixing performance on consecutive reaction yields without analyzing yields governed by mixing processes quantitatively.

The quantitative relationship between mixing and yields can be described using a reaction model containing parameters characterizing mesomixing and micromixing. Baldyga and Bourne proposed a model in terms of the influence of meso- and micromixing on reactions and applied it to a semi-batch reactor with an extremely high flow ratio of two mixing streams [12]. The simulation results had a good prediction accuracy with a yield error of ±2% compared to experimentally measured yields. Samant et al. studied consecutive reactions in a stirred batch reactor based on this model [13]. They found that structural parameters, such as the stirring paddle shape and the number of feed points, had significant effects on mixing rate-limited steps and reaction yields. Then they provided empirical correlations of mixing characteristic times at each scale controlled by these parameters to guide the reactor design.

Continuous reactors have many parameters that can significantly affect the internal fluid flow and mixing pattern, which, in turn, affect the yield of consecutive reactions [14,15,16,17,18]. Baldyga and Bourne used the aforementioned model to forecast the yield in a Kenics tubular reactor with a flow ratio of 3000:1 [12]. However, the yield was mainly influenced by the mixing pattern of the stream with the extremely low flow rate, which does not hold true for two streams with comparable flow rates in tubular reactors. Thus, the model requires modification based on the flow and mixing characteristics in such cases. Nevertheless, no model modification work has been conducted for meso-micromixing interaction reaction models in tubular reactors, which could ascertain the ambient fluid concentration of the reagents and the volume change of the fluid due to mixing.

The structural parameters of tubular reactors not only influence the mixing pattern and target product yield of reactions but are also closely related to the pressure drop in reactors. Maximizing the target product yield while keeping the pressure drop within a reasonable range is a valuable pursuit of study [19]. However, optimizing numerous parameters of reactors based on the reaction model can be computationally challenging [20]. Multidimensional regression and optimization required for this process are difficult to carry out. Gaussian process regression combined with Bayesian optimization can provide a general optimal solution for nonlinear optimization problems with fewer iterations and higher calculation speed [21,22,23]. Zuhal et al. [24] and Seongeon Park [25] employed CFD to conduct a series of Bayesian optimization designs to optimize the structure of the airfoil and stirred batch reactor based on minimum pressure drop criteria. However, there are currently no relevant studies on the structural design of continuous reactors that link reaction models with optimization methods to evaluate the relationship between target product yields and the pressure drop in reactors.

In this study, we analyze the flow characteristics and mixing pattern in continuous reactors and modify the meso-micromixing interaction reaction model which had been used in batch conditions in previous work. We use the Villermaux–Dushman method in a self-fabricated split and recombination reactor to verify the model’s accuracy by comparing the results obtained by experiments with the model. We analyze the changes in mixing performance due to parameter variations from a flow field perspective. We use this modified model combined with Gaussian process regression and Bayesian optimization to obtain higher yields of target products in consecutive reactions while maintaining a lower pressure drop in the reactor. Such a design method can be applied to the optimization of various types of reactors.

2. Modified Meso-Micromixing Interaction Reaction Model

2.1. Initial Mixing-Related Reaction Model

According to Baldyga and Bourne’s theory [12], a meso- and micromixing interaction reaction model can be characterized by the following equation:

$$\frac{d{c}_i}{dt}=E\left(1-\frac{X_B}{X_u}\right)\left(c_i-c_i\right)+r_i$$

(3)

The symbols in this equation are defined as follows: The rate of engulfment, $E$, describes the micromixing of fluid elements, which will be explained in detail in Section 3.4. $X_u$ represents a fraction of the fluid volume. This part of the fluid contains partially segregated fluid as islands embedded in a sea, where the reagents are coarsely mixed and concentration fluctuations still exist. $X_u$ can be derived by statistical method:

The variance of composition of the

$$X_u=\frac{X_0}{X_0+\left(1-X_0\right)\mathit{exp}\left(-\frac{t}{t_d}\right)}$$

(4)

The initial volume fraction in the reactor where the reagents are present is $X_0$. For tubular reactors, $X_0$ can be considered as the ratio of the initial flow rates. As $X_u$ tends to 1, it signifies that the reagents spread everywhere from the mesomixing perspective, although without being completely micromixed.

$X_B$ refers to the volume of micromixed fluid relative to the whole fluid. According to the engulfment model, the rate of engulfment from E model is:

$$\frac{d{X}_B}{dt}=E{X}_B\left(1-\frac{X_B}{X_u}\right)$$

(5)

This equation describes the micromixing pattern when mesomixing is incomplete. Initially the islands, whose initial volume fraction is $X_0$, only consist of pure B-rich feed, so that the initial value of $X_B$ is also $X_0$. As the flow enters the reactor, eddy breakup increases $X_u$, which initiates micromixing.

2.2. Comparison of Batch and Continuous Conditions

In batch reactors, the limited reagents are added dropwise. The chemical reaction occurs only in the droplets. The ambient fluid concentration 〈$c_i$〉 outside the micromixing region can be assumed to be a constant concentration of excess fluid in the reactor. Additionally, when the flow rate in one of these streams is extremely low, the volume of fluid that has been fully micromixed will be present only in the other excessive stream.

In the case of two comparable flow rates in the continuous reactors, chemical reactions can occur in both streams after achieving micromixing. This implies that each stream contains a significant concentration of reagents and products, as was shown in Figure 1. As a result, the ambient fluid concentration 〈$c_i$〉 is no longer a changeless concentration, which is the case in batch reactors. Furthermore, when the flow rates are comparable, both flows experience reciprocal engulfment effects and micromixing. Consequently, the volume of fluid that is completely micromixed in both streams changes (in Figure 1, the total volume of each stream is $V_A+\mathsf{\Delta }{V}_A$ or $V_B+\mathsf{\Delta }{V}_B$), which ultimately impacts the concentration of species in each stream.

2.3. Model Modification

We have modified the model after accounting for the differences between the above two conditions.

Firstly, since both streams, 1 and 2, exhibit self-engulfment behavior and undergo chemical reactions, we considered the ambient species’ concentration outside the micromixing region of stream 1 as that of the corresponding species in stream 2, and vice versa. In other words, these two streams are treated as reciprocal environmental fluids, as was shown in Figure 1.

Secondly, since the change in fluid volume resulting from the mesomixing and micromixing in these two streams is not negligible, we took stream 1 as the reference and switched from the concentration change due to the volume change of stream 2 to the volume-based concentration of stream 1. To describe this behavior, we introduced a fluid volume change rate ${\phi }_B$, which essentially represents the fluid volume change induced by mixing at the mesoscopic and microscopic scales. This rate can be calculated using the following integrated equation:

$${\phi }_B=e^t+\left(1-\frac{1}{e^{-Qt}}\right)Q$$

(6)

Here, $Q$ is the ratio of the micromixing characterized time and the mesomixing characterized time.

Based on these, the model can be modified as (taking concentrated reagent A in stream 1 for instance):

$$\frac{d{c}_{A1}}{dt}=E\left(1-\frac{X_{B1}}{X_{uA}}\right)\left(c_{A2}-c_{A1}\right)+r_A$$

(7)

$$\frac{d{c}_{A2}}{dt}=\left[E\left(1-\frac{X_{B2}}{X_{uA}}\right)\left(c_{A1}-c_{A2}\right)+r_A\right]\frac{\frac{V_{total}{V}_A}{\left(V_A+V_B\right)}{\phi }_B}{\left(V_{total}-\frac{V_{total}{V}_A}{\left(V_A+V_B\right)}{\phi }_B\right)}$$

(8)

where $\frac{V_{total}{V}_A}{\left(V_A+V_B\right)}$ is the inlet flow rate of stream 1.

By analogy, each species is governed by two ordinary differential equations that control its concentration in each of these two fluid environments (stream 1 and stream 2). In Equation (7), the subscript 1 represents the fluid that has undergone micromixing, while the subscript 2 represents the ambient fluid surrounding stream 1. For Equation (8), the opposite is true of what was explained above.

3. Methods and Materials

3.1. 3D-Printed Split-and-Recombine Millimeter-Scale Reactor

The equipment for studying the modified model is a tubular split and recombination reactor, whose physical model is illustrated in Figure 2a,b. A central annular tube with an outer diameter of 6 mm ($R_O$) and an inner diameter of 3 mm ($R_I$) comprises the reactor. It has a series of split and recombination mixing elements with a diameter of 6 mm and a length of 10 mm ($L_M$), distributed uniformly at intervals of 10 mm ($L_G$) along the circular tube. The structure has four twisted grooves uniformly rotating along the circumference for fluid passage. Two adjacent twisted grooves are separated at 90 degrees on the cross-section of the tube. The physical structure of the reactor can be obtained using 3D printing technology, with a processing accuracy less than 0.05 mm.

The tubular reactor consists of three sections: inlet, test, and outlet. The inlet section, 20 mm long, is designed to ensure the complete development of the fluid. These two mixing streams are introduced at equal velocity there, with one entering the reactor through a 3 mm central tube, and the other entering through the remaining annular tube. Thus, the initial flow ratio of these two streams will be 1:3. The test section is 110 mm long, an adequate length to achieve various mixing patterns. The outlet section is 40 mm long, sufficient to stabilize the outlet flow.

Based on the results obtained from the pre-experiments, the angle ($\alpha$), groove depth ($d$), turn ($T)$, and flow rate ($F$) of the reactor significantly affect the flow field and mixing performance in the reactor. These parameters are shown in Figure 2c. Angle denotes the rotation angle of the groove in the mixing element; groove depth represents the depth of the groove; turn refers to the degree of one groove twist. If $T=1$, the groove rotates 360° in a single mixing element. We selected these parameters for numerical simulation and experimental study and set several discrete points for each parameter, which are listed in

Table 1. Geometrical parameters and flow rate of the tube and inserts.

$\mathit{\alpha }$, °	$\mathit{d}, \mathbf{m}\mathbf{m}$	$\mathit{T}$	$\mathit{F}, \mathbf{m}\mathbf{L}/\mathbf{m}\mathbf{i}\mathbf{n}$
15, 30, 45, 60, 75	0.5, 1, 1.5, 2, 2.5	0.2, 0.4, 0.6, 0.8, 1	100, 150, 200, 250, 300, 350, 400

.

3.2. Numerical Simulation

To investigate the influence of different structural parameters and inlet flow rates on pressure drop and the energy dissipation rate in the reactor, we conducted 3D steady-state numerical simulations. Water was chosen as the working medium and modeled as a Newtonian, incompressible fluid. The simulations were performed at 20 °C and it is assumed that physical properties, such as viscosity, remain constant. Since the mixing element is placed horizontally in a tube with an inner diameter of 6 mm, we assume that the effect of gravity is negligible. With these assumptions, the momentum governing equations of the fluid are as follows.

$$\frac{\partial u_i}{\partial x_i}=0$$

(9)

$$\rho \frac{\partial u_i}{\partial x_i}+\rho \frac{\partial u_i{u}_j}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_i}{\partial x_i}+p{\tau }_{ji}\right)\left(i=1~3\right)$$

(10)

where $\rho$ is the density of mixture fluid; $u$ is the internal velocity field; $p$ is the local pressure field; $\mu$ is the viscosity of the mixture fluid; and ${\tau }_{ji}$ is the Reynold stress tensor. Here, the SST k-ω model was applied. The range of the Reynold number in this study was from 658 to 5821. The SST k-ω model can be applied to the cases with relatively lower Reynold numbers, which can reproduce the transition from laminar to turbulent flow regimes. Furthermore, for the flow ejected from the curved grooves of the mixing elements, where the fluid has a high-velocity gradient within the boundary layer, the SST k-ω model will be more applicable.

In this study, ANSYS fluent 2020R1, a commercial software based on the finite volume method, was used for the numerical simulations and CFD-Post for post-processing. The SIMPLE algorithm was applied for pressure–velocity coupled solution calculations while the PRESTO! method was used to perform gradient and pressure discretization. The second-order upwind algorithm spatially discretized the momentum, turbulent kinetic energy, and turbulent dissipation rate. The convergence criteria are less than 10⁻⁶ for the continuity and momentum equations.

In the simulations, there is difficulty calculating the turbulent kinetic energy, energy dissipation rate, etc., accurately, so, obtaining these results requires higher quality discrete meshes. Therefore, the grid independence test and algorithm reliable verification are necessary to be carried out, which are shown in the Supplementary Materials.

3.3. Mixing Performance Experiments

To quantitatively evaluate mixing performance, classical Villermaux–Dushman experiments [26,27,28] were conducted. We considered the approximate micromixing time range in this case and used the concentration groups of substances provided by J.M. Commenge’s study [29,30,31]. The concentrations of each reagent are presented in

Table 2. Concentration recipe of reagents for Villermaux–Dushman method.

Materials	Concentration [mol/L]
H2BO3	0.09
NaOH	0.09
KIO3	0.006
KI	0.032
H2SO4	0.0026

.

The experimental setup diagram is depicted in Figure 3. Iodide-iodate-borate buffer solution and sulfuric acid solution were supplied to the reactor through two HPLC pumps from the outer and central tubes, respectively. The outer and central tubes were connected by commercial T-connectors. To minimize error, three independent samples were taken in each experiment and analyzed by UV-VIS. Finally, the concentration of each species was obtained by subsequent data processing, then X_S (segregation index) was calculated. For further details on the Villermaux–Dushman method [32,33,34], please refer to the Supplementary Materials.

3.4. Data Reduction

Based on the set of ordinary differential equations given in Section 2, the independent variables that need to be imported to the reaction model include the inlet flow rate, reaction kinetics, flow time, micromixing characteristic time and mesomixing characteristic time. Apart from the inlet flow rate and reaction kinetics, the calculation of the remaining three variables is given below.

The flow time can be calculated by dividing the reactor length by the average flow velocity in the reactor as follows:

$$t=\frac{L}{u_m}$$

(11)

In such SAR reactors, additional attention needs to be taken with respect to $D$, because the shape of the flow channels in the circular cross section is not regular and the reactor characteristic sizes cannot be determined by conventional methods. Therefore, the hydrodynamic diameter calculation method was used [35]:

$$D=\frac{4V}{S}$$

(12)

where $V$ is the fluid domain volume in the reactor and $S$ is the area which contacts with the fluid in the reactor.

The Reynolds number represents the relative relationship between the inertial and viscous forces in the fluid and is expressed as follows:

$$Re=\frac{\rho u_{m}D}{\mu }$$

(13)

where $u_m$ is the average velocity along the flow direction.

Baldyga and Bourne [36,37,38,39] proposed an engulfment model to describe micromixing behavior, which is based on the premise that the engulfment of small-scale eddies generates interlaced laminar structures between the micromixing fluid and the environment. These structures promote the deformation of fluid layers that accelerates the aggregate size reduction up to the diffusion scale then intensifies the molecular diffusion between the fluid layers. Based on this theory, the micromixing characteristic time can be determined from numerical simulation using the following equation [40]:

$$t_m=E^{-1}=17.3{(v/\epsilon )}^{0.5}$$

(14)

Here, the engulfment rate, $E$, is the inverse of $t_m$, which is related to the energy dissipation rate $\epsilon$ and the dynamic viscosity $\nu$. It is important to note that this empirical formula (Equation (11)) for calculating the micromixing characteristic time requires more validation through comparison with experimental results, which are detailed in the Supplementary Materials.

Additionally, the energy dissipation rate is directly related to the evaluation index of mesomixing [11,41]. According to the theory of inertial-convective disintegration of large-scale eddies, the fluid element’s intermediate scale shrinks from the initial ${\Lambda }_c$ to the Kolmogorov scale through turbulence dispersion. The whole process constitutes mesomixing [38,41,42,43]. Based on this theory, the mesomixing characteristic time can be estimated using statistical method [12]:

$$t_d=2{\left(\frac{{\Lambda }_C^2}{\epsilon }\right)}^{\frac{1}{3}}$$

(15)

The initial scale of mesoscopic mixing, ${\Lambda }_c$, for two streams with the same inlet flow velocity can be calculated by the following equation [11]:

$${\Lambda }_c={\left(\frac{V_B}{\pi u}\right)}^{0.5}$$

(16)

where $V_B$ is the volume flow rate of the stream that has a relatively lower flow rate.

The energy dissipation rate refers to energy loss in the fluid due to its essential viscosity and the turbulence induced. The turbulent energy dissipation rate can be related to the velocity strain rate tensor by:

$$\epsilon =\frac{1}{2}v{\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)}^2$$

(17)

This tensor can be obtained directly from CFD simulations.

3.5. Optimization Procedure

The optimization method for reactors using GPR, BO, and the modified model is presented in the flowchart shown in Figure 4. In this study, 100 initial design points about the reactor were obtained using the Latin hypercube sampling method. CFD simulations were carried out to determine parameters such as the energy dissipation rate. The reaction kinetics and modified models were then incorporated into a self-programmed 4th-order Runge–Kutta MATLAB program for calculations, acquiring 100 sets of yield and pressure drop data for fitting. These yield and pressure drop data, as performance functions, were mapped to performance metrics by a specific design. This design was defined by a set of structural and operating parameters and subjected to regression validation and optimization in the next step.

In each iteration, a surrogate model was calculated using the known correspondence between the four reactor parameters and the three dependent variables (two product yields plus pressure drop). The next set of parameters, whose mapping value will get closer to the maximum yield and the lowest pressure drop, was chosen by the acquisition function. This parameter set was further used for numerical simulation to obtain the new pressure drop and yields. The cycle of updating the model was repeated. An evaluation index $K$ was defined to determine the final multi-objective optimization results.

$$K=\frac{\frac{Y_B-Y_A}{Y_A}}{\frac{\mathsf{\Delta }{P}_B-\mathsf{\Delta }{P}_A}{\mathsf{\Delta }{P}_A}}$$

(18)

Specifically, $K$ is calculated as the ratio of the percentage yield improvement to the percentage pressure drop improvement between two adjacent iteration points, A and B. This metric indicates the amount by which the yield improves when the unit pressure drop increases. If $K$ is less than 5%, we select the current design point as the final multi-objective optimization result based on the optimization principle of GPR-BO.

4. Results and Discussion

4.1. Validation of Modified Model Accuracy

The present study acquired experimental and simulated yields using parallel competing reactions in the classical Villermaux–Dushman method. The yields obtained from both the modified meso-micromixing interaction reaction model and the reaction model containing only the micromixing time term were compared with the experimental results. The results are presented in Figure 5.

The figure shows that the yields obtained from the modified meso-micromixing interaction reaction model are in good agreement with the experimental results, with a maximum relative error of only 3% (Figure 5b). This indicates the accuracy of this modified model for yield prediction. However, the reaction model with the micromixing time term only has a much larger error in the experimental results, particularly at high flow rates, where the relative error can reach up to 14.7% (Figure 5c). This can be attributed to mesomixing becoming the rate-limited step gradually, which is evident from Figure 5, showing that the ratio of micromixing time to mesomixing time, $Q$, decreases with increasing flow rates. Growing flow rates increase the overall energy dissipation. However, the rate of mesomixing triggered by turbulent dispersion cannot catch up with that of micromixing, and, thus, becomes the rate-limited step [44]. In conclusion, this modified meso-micromixing interaction reaction model is an accurate approach for obtaining yields and can subsequently be used for the optimization in Section 4.3.

4.2. Effect of Parameters on Mixing

4.2.1. Effect of Turn

In Figure 6, the reactor’s full domain energy dissipation rate and segregation index were evaluated at various turns for a fixed groove depth of 1.5 mm and angle of 45°, while the inlet flow rates were also varied. The results indicated a decreasing trend in energy dissipation rate and segregation index with an increase in $T$, which was found to be more pronounced at lower flow rates. At a flow rate of 150 mL/min, X_S for a $T$ of 0.4 was 0.03911, which was 40.4% lower than that for a $T$ of 0.4 with an X_S value of 0.06564. However, X_S variation at higher $T$ did not become that significant, especially when T = 0.4 and 0.6, as there was almost no difference, with a maximum difference of less than 5%.

The helical motion of the fluid caused by the increased radial velocity component and extended flow path, both in the empty annular space between the two mixing elements and the twisted grooves in the mixing elements, was responsible for this trend (Figure 6b2). At the same time, the strong rotational effect increased the contact area between relatively high and low flow velocity zones in the twisted grooves (Figure 6c2). However, although further increase in $T$ indeed lengthened the flow path, excellent mixing could be achieved right in the middle of the path due to more frequent contact with the fluid. Therefore, extending the flow path or increasing the radial velocity more significantly had a limited effect on improving mixing performance [45]. Overall, these findings highlight the importance of selecting the appropriate turn for optimal reactor performance.

4.2.2. Effect of Inlet Flow Rates

Figure 5 illustrates the decreasing trend of X_S with increasing inlet flow rates, from a maximum of 0.05331 to 0.01911. For $T>0.4$ and at lower flow rates ($F<200 \mathrm{mL}/\min$), X_S decreased significantly with increasing flow rate, from 0.06464 to 0.02352, representing a change of 63.6%. However, at higher flow rates, the change in X_S leveled off until it concentrated around 0.0123 at a flow rate equal to 400 mL/min.

This trend was mainly attributed to the increasing energy dissipation rate in the annular space between the two mixing elements as the flow rates increased. As shown in Figure 7a1,a2, the flow pattern remained largely unchanged with increasing flow rates, but the overall velocity was significantly affected. This increased velocity leads to a higher degree of turbulence and contact frequency with the fluid, thereby substantially increasing the energy dissipation rate in the annular space (Figure 7b1,b2). However, the effect of higher flow rates on the energy dissipation rate in the annular space is minimal, resulting in a slow improvement in mixing performance. Overall, high flow rates are beneficial for enhancing mixing performance in the reactor.

4.2.3. Effect of Cross-Sectional Area of the Grooves

Angle and groove depth are important factors that determine both the width and depth of grooves in the mixing element, respectively. These two factors together determine the cross-sectional area of grooves, which affects the flow pattern and mixing performance in the annular space and twisted grooves.

Given $T=0.4$, $F=300 \mathrm{mL}/\min$, the reactor’s energy dissipation rate and segregation index varied with groove depth and angle, as presented in Figure 8. Decreasing the cross-sectional area of the groove led to a significant reduction in X_S. The maximum X_S of 0.0258 was obtained at $\alpha =75°$ and $D=2.5 \mathrm{mm}$, while a minimum X_S of 0.00426 was observed at $\alpha =30°$ and $D=1 \mathrm{mm}$, overall resulting in an 83.5% decrease in X_S.

Reducing the cross-sectional area of the groove leads to a jet with higher velocity, which spreads out into the annular space to form a rotating plume flow. As this plume moves away from its source, its edges entrain the surrounding fluid, causing the plume to expand further then promoting mixing (Figure 8b2). The energy dissipation rate increased significantly when the cross-section area of the groove became small, from 1.46 m²/s³ at $\alpha =45°$ and $D=1.5 \mathrm{mm}$ to 8.39 m²/s³ at $\alpha =30°$ and $D=1 \mathrm{mm}$. This can be explained by the fact that the flow area is reduced by a factor of 2.25 and the velocity of the ejected plume increases exponentially, leading to an accelerated improvement in mixing performance. In conclusion, reducing the cross-sectional area of the groove improves the mixing performance of the reactor.

4.3. Optimization Calculations

4.3.1. Optimization Objectives

The objective of this work is to maximize the target product yield and minimize pressure drop by optimizing the reactor parameters. The information about the relevant consecutive reaction is presented in

Table 3. Kinetic and operational information about the competitive consecutive reaction.

Information about the Reactions
$A+B\to C r_1=k_1{c}_A{c}_B$  $C+B\to D r_2=k_2{c}_C{c}_B$  $D+B\to E r_3=k_3{c}_D{c}_B$  $k_1=50 {\mathrm{m}}^3/\left(\mathrm{mol}·\mathrm{s}\right) k_1=20{ \mathrm{m}}^3/\left(\mathrm{mol}·\mathrm{s}\right) k_1=100{ \mathrm{m}}^3/\left(\mathrm{mol}·\mathrm{s}\right)$
$V_A :V_B=1:3 V_{MA}=V_{MB}=5 \mathrm{mol}/{\mathrm{m}}^3$

. The reactor is fed with two streams, A and B, which enter from a 3 mm diameter center tube and annular space with a 6 mm outer diameter and 3 mm inner diameter, respectively. To satisfy the demand for different products in the consecutive reaction, the yields $Y_1$ and $Y_2$ of product C and product E, as well as the pressure drop, were chosen as the optimization objectives, while the angle, groove depth, turn and flow rate were selected as the corresponding optimization parameters.

BO is applied to obtain the local minimum value [46]. Therefore, the yields of both products in the consecutive reaction were multiplied by (−1) to transform the optimization problem into a maximization problem, as shown in the following expression:

$$\mathrm{Objectives}: f_1\left(\alpha , d, T, F\right)=\left\{\Delta P_{SAR},-Y_1\right\} \mathrm{and} f_2\left(\alpha , d, T, F\right)=\left\{\Delta P_{SAR},-Y_2\right\}$$

$$\mathrm{Subjected} \mathrm{to}: \alpha \in \left[15,75\right]; d \in \left[0.5,2\right]; S \in \left[0,1\right]; F \in \left[100,400\right]$$

4.3.2. Results of Gaussian Process Regression

Two commonly used metrics to evaluate a data regression are the mean square error ($MSE$ for short) and the coefficient of determination ($R^2$ for short), which are defined as:

$$MSE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(Y_{i, GPR}-Y_i\right)}^2$$

(19)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^N{\left(Y_{i, GPR}-Y_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(\overline{Y_i}-Y_i\right)}^2}$$

(20)

A GPR that has a smaller $MSE$ and a larger $R^2$ is considered more reliable and accurate. Three separate GPRs were utilized to fit the data of $\mathsf{\Delta }{P}_{SAR}$, $Y_1$, and $Y_2$. To assess the effectiveness of GPR for the regression, an additional analysis was conducted using an artificial neural network (ANN) with multilayer perception (MLP). The tanh function was employed as the activation function in the ANN, which has been widely used in this context [19,47,48,49]. The final trained ANN had a configuration of 5:14:12:2 after double hidden layers were used. Both the regression results of ANN and GPR were compared with results obtained by the reaction model (short for RM), whose normalized values are presented in Figure 9 and

Table 4. MSE and $R^2$ of three optimization objectives for GPR and ANN.

	GPR		ANN
	$\mathit{M}\mathit{S}\mathit{E}$	${\mathit{R}}^2$	$\mathit{M}\mathit{S}\mathit{E}$	${\mathit{R}}^2$
$\mathsf{\Delta }{P}_{SAR}$	1.40880 × 10−4	0.94593	8.51216 × 10−5	0.99670
$Y_1$	3.38565 × 10−3	0.78948	2.87026 × 10−4	0.97674
$Y_2$	3.85881 × 10−3	0.96288	1.45477 × 10−4	0.99972

.

In Figure 9, the deviations of three dependent variables for RM, ANN, and GPR are presented. The maximum relative errors for $Y_1$ and $Y_2$ were 3.52% and 5.75%, respectively, when using GPR. On the other hand, while ANNs generally meet accuracy requirements, there were several outliers that deviated significantly from the CFD calculations, resulting in maximum errors of 29.62% and 36.86% for $Y_1$ and $Y_2$, respectively.

Table 4 shows that both the $MSE$ and $R^2$ values in GPRs were significantly better than those in ANNs, indicating that GPR can provide more accurate data and is, therefore, more suitable for fitting in this case. This is because GPR can quantitatively predict uncertainty in a more principled way than ANN, which is particularly relevant for conditions with excessive degrees of freedom. Moreover, by selecting a specific kernel function, GPR can be considered equivalent to an ANN with an infinite number of hidden nodes [50]. However, high accuracy requires more computational resources, although, for just 100 design points, the computational time remains relatively short.

4.3.3. Results of Bayesian Optimization

After obtaining regression results for GPRs, Bayesian optimization was subsequently performed to optimize the reactor parameters.

The optimization process is determined as a multi-objective optimization when the target product is the intermediate product C. To improve the yield of C, the mixing performance needs to be enhanced to disperse the reagent B in the reactor quickly, which reduces its local concentration and slows down the second and third reactions, thereby accumulating product C and reducing the formation of subsequent products, D and E. However, improving the mixing performance requires a higher energy dissipation rate of the fluid, resulting in a higher pressure drop in the reactor. Therefore, a compromise solution between higher intermediate product yields and a relatively lower reactor pressure drop needs to be achieved.

Using Bayesian optimization, a series of Pareto front points were obtained, as shown in Figure 10. These optimization points were evaluated using the metric $K$ presented in Section 3.5. The desired points with $K=2.92\%$ were identified with a red circle in Figure 10. Therefore, for the intermediate product C, the maximum yield was 92.5% with a pressure drop in the reactor of 510.50 Pa, where $\alpha =32.8°$, $d=2.36 \mathrm{mm}$, $T=0.215$, and $F=145.7 \mathrm{mL}/\min$.

When the target product is the product in the last reaction, E, the process becomes a single-objective optimization. To increase the yield of E, the mixing performance needs to be worsened, resulting in the intermediate product C being diffused to the enriched B region, further reacting with it. This is particularly significant when the residence time of the species becomes longer, as it leads to a higher selectivity of E. A lower energy dissipation rate results in poorer mixing, but also keeps the pressure drop low. Therefore, there is no trade-off or coordination required between the product yield in the last reaction and the pressure drop in the reactor. Optimization is straightforward in that the solution containing both the lowest pressure drop and the highest yield is desired.

A series of optimal points were obtained by Bayesian optimization, as shown in Figure 11. The three leftmost points in the figure are indeed close in value, with a relative difference of less than 0.2%, indicating that the single-objective optimization has reached the limit of iterations. The point identified with the red circle is the desired point, where we acquired the lowest pressure drop and the highest yield. Thus, for the product in the last reaction, the highest yield was 94.3%, with a pressure drop in the reactor of 253.81 Pa, where $\alpha =48.3°$, $d=2.04 \mathrm{mm}$, $T=0.344$, and $F=120.37 \mathrm{mL}/\min$.

We have developed an active reactor design method based on CFD-GPR-BO that is applicable to the analysis of different types of reactions, accounting for a meso-micromixing interaction reaction model. In practice, it may be acceptable for a reactor to have a higher pressure drop to achieve the most desired yields as long as the pressure drop is within a reasonable and allowable range. Nevertheless, this method still provides valuable guidance for practical applications, including reactor development for different types of reactions, significantly reducing the reactor design time.

5. Conclusions

In summary, the present study has led to the following conclusions:

• A modified meso-micromixing interaction reaction model was developed based on the flow characteristics in continuous reactors. The model was validated by comparing experimentally obtained yields with those predicted by this model. The modified model significantly reduced error in predicted product yields from approximately 15% to within 3%, compared to the model containing the micromixing term only.
• Mixing performance in the reactor was improved by characterizing the decreasing X_S with increasing flow rate, the degree of twist in the mixing element’s grooves, and the decreasing cross-sectional area of grooves. A high flow rate intensifies the energy dissipation of the fluid in the annular space between two mixing elements; high turn extends the flow path and increases the contact area between the areas with high and low flow velocities in the twisted grooves. When the cross-sectional area in the grooves becomes small, a significant plume flow can be formed in the annular space, improving mixing performance.
• The optimization, in which the yields of target products and pressured drop in the reactors were chosen as the optimization objectives, was based on the modified model and performed by BO along with GPR. We obtained the highest product yield while keeping the pressure drop low. For the intermediate product, the yield was 92.5%, while the pressure drop in the reactor was 510.50 Pa. For the product in the last reaction, the yield was 94.3%, while the pressure drop in the reactor was 253.81 Pa. The corresponding combinations of reactor parameters were obtained. This kind of optimization method can be applied to the design of various reactors, providing a reference for structural selection and operational parameter determination.