Multi-objective optimization design of microchannel reactors for Fischer–Tropsch synthesis using CFD, GENN, and NSGA-II

Abstract

In Fischer–Tropsch synthesis processes, microchannel reactors exhibit pronounced process intensification compared with conventional fixed-bed reactors. Computational fluid dynamics was coupled with a surrogate model based on a gradient-enhanced neural network to systematically evaluate the influence of four characteristic geometric variables on the catalytic performance of multi-tubular microchannel reactors. A multi-objective optimization aimed at maximizing C₅₊ yield and concurrently minimizing the maximum temperature rise was then conducted using the non-dominated sorting genetic algorithm-II. The resulting Pareto frontier was analyzed to identify the solutions that optimally reconcile thermal safety and productivity. The results indicate that, relative to the initial design, the selected optimal configuration reduces ΔT_max by 7.2 °C and increases Y_C5+ by a factor of 1.86, substantially enhancing reactor performance and providing both a theoretical basis and design reference for pilot-scale demonstration and industrial deployment of Fischer–Tropsch microchannel reactors.

---

1. Introduction

Non-oil-based carbon resources such as coal, natural gas, and biomass can be converted into clean fuels and high-value chemicals via Fischer–Tropsch synthesis (FTS), thereby alleviating energy shortages, strengthening strategic energy security, and reducing environmental pollution.¹ Owing to the strongly exothermic nature of the FTS, microchannel reactors exhibiting superior interphase heat and mass transfer characteristics are employed to effectively suppress thermal runaway and enhance selective hydrocarbon productivity.² Nevertheless, the micro/millimetric characteristic dimensions preclude intrusive insertion of temperature, concentration, or pressure probes, which inevitably disturbs the flow field; complete spatial profiles of temperature, concentration, and pressure cannot be obtained experimentally, and reactor design therefore is primarily based on computational fluid dynamics (CFD) simulations.³

In recent years, CFD has been coupled with optimization algorithms for the design of diverse reactor types. Lee et al.⁴ rapidly determined the catalyst packing ratio in a four-channel reactor by integrating CFD with a Bayesian envelope optimization algorithm. Teimouri et al.⁵ embedded a genetically calibrated trimetallic kinetic scheme within a two-dimensional porous-medium CFD model to accurately predict the C₅₊ distribution in a fixed bed. Na et al.⁶ employed a CFD-genetic algorithm framework to optimize a zoned catalyst dilution strategy, achieving a significant reduction in temperature rise and an improved yield. Nevertheless, conventional CFD-based design remains a labor-intensive procedure:⁷ any variation in diameter, length, or thickness necessitates remodeling, remeshing, and re-specification of boundary conditions. When the design space encompasses tens to hundreds of geometric combinations, the “geometry–mesh-solution” cycle becomes computationally prohibitive, and the discrete results yield sparse gradient information for continuous-space optimization.

Gaussian-process regression, radial-basis-function interpolation, and neural networks are routinely employed to construct input–output mappings from limited CFD samples, enabling the resulting surrogate models to replace the CFD solver during the optimization loop and thereby alleviate the associated computational burden. Mishra et al.⁸ coupled CFD with a deep neural network to enable rapid prediction of local hydrodynamic variables such as bubble number density in a slurry-bed hydrocracking reactor. Qiu et al.⁹ generated two-dimensional H₂ concentration fields from CFD simulations of a methanol steam-reforming fixed bed and subsequently trained a multilayer-perceptron surrogate model to map spatial coordinates to concentration, reducing computational time by 90%. Jung et al.¹⁰ performed a dual-objective optimization of reactor volume and maximum temperature difference for a U-cooled microchannel FT process by training an artificial-neural-network surrogate with CFD samples. Ansys has further slashed CFD simulation turnaround times via its cloud-native SimAI platform, which leverages physics-informed AI models trained on legacy simulation datasets to accelerate performance predictions by 10–100 times compared with traditional solvers.¹¹ Collectively, these studies and reports demonstrate that CFD-surrogate coupling constitutes a viable route for mitigating computational intensity in CFD-driven design.¹² Surrogate-based modeling enables computationally efficient and quantitative evaluation of the catalytic performance of the microreactor within the bounded geometric design space. Owing to inherent trade-offs among reactor-design variables, multi-objective optimization algorithms must be implemented to identify the Pareto frontier and hence the optimal compromise.^10,13 Jiang et al.¹⁴ developed a CFD–ANN–NSGA-II framework to optimize the structural parameters of a stirred-tank reactor with high computational efficiency. Zhang et al.¹⁵ employed NSGA-II to minimize radial temperature difference and maximize CO conversion, thereby determining the optimal operating conditions of a helium-heated reverse water–gas shift reactor. However, studies on surrogate-assisted multi-objective optimization for Fischer–Tropsch synthesis in multi-tubular microchannel reactors remain unavailable in the open literature.

On the basis of our previously validated CFD model,¹⁶ this work presents an integrated framework that combines sensitivity analysis, surrogate modelling and multi-objective optimization for the geometric design of a multi-tubular microchannel reactor for FTS. Latin hypercube sampling (LHS) was used to generate 250 CFD cases within the bounds of the four geometric design variables (R_L, R_D, R_T, R_S). A gradient-enhanced neural network (GENN) surrogate model exhibiting the highest predictive accuracy was then constructed, eliminating the need for further computationally intensive simulations. Finally, the GENN surrogate model was integrated with the Non-dominated Sorting Genetic Algorithm II (NSGA-II) algorithm to minimize the maximum temperature rise (ΔT_max) and maximize the C₅₊ yield (Y_C5+). The resulting Pareto frontier was analyzed to identify Pareto-optimal solutions that balance thermal stability and productivity, providing a quantitative design basis and optimization strategies for pilot-to-industrial scale-up of multi-tubular microchannel FTS reactors.

2. Reactor model and optimization strategy

2.1 Reactor model

Fig. 1 illustrates the multi-tubular microchannel reactor employed in this study. The reactor consists of seven inner tubes (Tube I) and one outer tube (Tube II). Tube I are axisymmetrically arranged around the central axis of Tube-II, with their ends recessed 30 mm from the corresponding inlet/outlet. The inlet and outlet headers are packed with inert particles to guarantee fully developed laminar flow inside the reactor. Both the inner and outer surfaces of Tube I undergo controlled chemical etching, ultrasonic cleaning and subsequent wash-coating with the active catalytic layer reported in ref. 17, thereby serving simultaneously as catalytic surfaces and mechanical supports. Syngas enters the reactor through the inlet of Tube-II and reacts on the catalytically active surface of Tube I. Reactor fabrication, CFD modelling, validation and flow-field distribution detailed in ref. 16 and 17 are reused here with the same geometric model, mesh, governing equations, kinetics and boundary conditions, with only the dimensions of Tube I treated as parametric variables. Specifically, the CFD model was previously validated in ref. 16, which described a multi-tubular microchannel reactor loaded with Fe-supported catalysts and featuring an identical geometric configuration to that of the present work. Ref. 16 investigated a wide range of operating parameters (inlet temperature: 220–380 °C, pressure: 0.1–2.1 MPa, H₂/CO ratio: 1/3–7/3, GHSV: 2000–16 000 h⁻¹), which fully covers the fixed baseline operating conditions utilized in this study. The simulation results from ref. 16 showed good agreement with experimental data, with mean absolute relative residuals of 11.74% for CO conversion and 0.43% for reactor center temperature, confirming the model's reliability. This consistency in reactor geometry and overlap in operating conditions between ref. 16 and the current study provide robust support for the credibility of the optimized design's simulation results, thereby enhancing the work's reference value for subsequent pilot-scale demonstrations. Geometric design variables including tube length (R_L), inner diameter (R_D), wall thickness (R_T) and inter-tube spacing (R_S) are specified in Fig. 1, where R_S is established via two one-piece multi-tube locating sleeves positioned inside Tube-II at both ends of the bundle. The corresponding values of each variable are listed in Table 1.

Fig. 1 Geometric model of the reactor. Tube II dimensions are governed by those of Tube I: its length is 60 mm greater, its inner diameter is larger by twice R_G, and its wall thickness is fixed at 1 mm.

Table 1 Baseline values and design ranges of the four geometric dimensions (mm)

Parameters^a	RL	RD	RT	RS
a All cases were simulated under an inlet temperature of 340 °C, pressure of 0.7 MPa, H2/CO ratio of 2 : 3, and GHSV of 8000 h^−1.
Baseline values	100	2	0.5	1
Exploration range	[30, 400]	[0.5, 50]	[0.25, 35]	[0.01, 50]
Design range	[30, 200]	[0.5, 4]	[0.25, 6.5]	[0, 2.5]

---

2.2 Integrated CFD-surrogate-optimization workflow

A hybrid CFD–GENN–NSGA-II workflow (Fig. 2) is proposed for the high-fidelity design of the multi-tubular microchannel reactor loaded with Fe-supported catalysts for FTS. Reaction kinetics and grid independence have been validated previously in ref. 16 (white block). Three tasks were undertaken in the present study:

Fig. 2 Hybrid CFD–GENN–NSGA-II workflow for multi-objective optimization of the microchannel reactor for FTS.

(1) Sensitivity analysis (yellow block): the established CFD model was exercised under fixed operating conditions (inlet temperature 340 °C, H₂/CO = 2/3, P = 0.7 MPa, GHSV = 8000 h⁻¹) to quantify the individual influence of R_L, R_D, R_T and R_S on ΔT_max and Y_C5+, thereby delimiting the feasible design space.

(2) Surrogate-model construction (blue block): a 250-point LHS design was generated and evaluated by CFD. The data set was randomly split into training (60%), validation (20%) and test (20%) subsets; a GENN surrogate model was then trained to provide rapid estimation of the objective functions across the design space.

(3) Multi-objective optimization (green block): NSGA-II was employed to maximize Y_C5+ while minimizing ΔT_max, yielding the Pareto frontier and corresponding trade-off solutions.

Geometric modelling was carried out in Ansys SpaceClaim, whereas mesh generation and CFD simulations were performed in Ansys Fluent, which provides practically proven and mature computational capacities for reactive flow research.¹⁸ Training of the surrogate model was conducted with SMT (Surrogate Modeling Toolbox) package v.2.9.5,¹⁹ and the optimization routine was handled by PyMOO (Multi-Objective Optimization in Python) package v.0.6.1.5.²⁰ All computations and Python scripting were executed on a Dell Precision 7920 Tower workstation (2× Intel Xeon Gold 6258R, 56 cores). Each design cycle was completed in approximately 5 min, corresponding to a ten-fold reduction from the conventional 50 min workflow and substantially improving the computational efficiency of multi-objective optimization of microchannel reactors.

2.3 Surrogate-model selection

Surrogate models were introduced to replace the computationally expensive CFD step and to provide a rapid mapping between design variables and reactor performance. Two candidate surrogates from the SMT library were benchmarked: (i) KRG: a data-efficient Gaussian-process Kriging estimator that simultaneously returns predictive mean and variance, thereby furnishing a built-in metric for adaptive infilling;²¹ (ii) GENN: a feed-forward network with known adjoint gradients embedded in the loss function, yielding higher-order convergence and superior generalization in a four-dimensional continuous design space for the same sample size.²² The procedure is summarized below:

Step 1: sample generation: 250 designs were generated in the 4-D space (R_L, R_D, R_T, R_S) by LHS and sequentially solved by CFD, generating y₁(Y_C5+) and y₂(ΔT_max) to form the complete data set.

Step 2: data partitioning: the data were randomly partitioned 60%/20%/20% into training, validation and test subsets. The training subset was used to fit the surrogate, the validation subset to monitor generalization error and facilitate hyper-parameter tuning with early stopping, and the test subset to provide an unbiased estimate of predictive accuracy.

Step 3: model training and assessment: KRG and GENN were trained concurrently. Accuracy was quantified by normalized mean absolute error (NMAE) and normalized root-mean-square error (NRMSE), defined in eqn (1) and (2). The model exhibiting the lower validation error was selected for the subsequent optimization loop.

where n denotes the number of samples. Subscripts sim and pre refer to the CFD simulation and surrogate-model prediction values, respectively. Max and min indicate the sample maximum and minimum values.

2.4 Optimization strategy

In the previously reported CFD model from ref. 16, products from FTS and WGS were lumped into C₁, C₃, C₅₊ and CO₂. Maximizing Y_C5+ was set as the first objective. As FTS is strongly exothermic, previous studies^23,24 commonly limit ΔT to 20 °C to mitigate thermal runaway.

ΔT_max minimization was therefore adopted as the second objective, subject to ΔT_max ≤ 20 °C. The multi-objective problem is formulated as eqn (3), where f₁(x) denotes the reciprocal of Y_C5+ and takes a minimum value of 1, while f₂(x) corresponds to ΔT_max; the independent variable x is a four-dimensional design vector, y_1,pre and y_2,pre are the Y_C5+ and ΔT_max predicted by the surrogate model, and S refers to the feasible region determined via the sensitivity analysis presented in Section 3.1.

This bi-objective optimization problem was solved using the NSGA-II algorithm implemented in PyMOO The algorithm employs fast non-dominated sorting and crowding distance mechanisms to simultaneously ensure the convergence and diversity of solutions.²⁵ After obtaining the Pareto frontier, the optimal design was selected based on the trade-off principle between temperature rise and product yield.

3. Results and discussion

3.1 Sensitivity analysis of design variables

Local geometric dimensions of the steel tubes composing the reactor directly determine the flow cross-sectional area and channel space velocity, and thus exert a significant influence on reactor performance. Four variables (R_L, R_D, R_T, R_S) were selected to characterize the key geometric features. At fixed operating conditions (inlet temperature of 340 °C, operation pressure of 0.7 MPa, H₂/CO molar ratio of 2/3, GHSV of 8000 h⁻¹), with the dimensions listed in Table 1 as the baseline, the control variable method was adopted, wherein only one of the four variables was adjusted relative to the base value at a time while the others remained constant. For each variable-adjusted scheme, the reactor performance (e.g. selectivity, temperature, yield) was obtained via CFD simulations, and the results are presented in Fig. 3. Based on these results, the final reasonable ranges of each design variable were determined, and three typical values per variable were selected to generate contour plots (Fig. 4–7) for auxiliary analysis. All these plots share a consistent format, comprising 3D perspective views and 2D sections (axial cross-section of Tube I, plus transverse cross-sections 10 mm inward from Tube I's inlet and outlet), with only the research variables differing across figures. Each variable in Fig. 3 was calculated using eqn (4)–(7).

Fig. 3 The relationships between CO conversion, CO₂ selectivity, C₅₊ selectivity, temperature and the four geometric variables. (a–d) Represent the length, inner diameter, thickness and inter-tube spacing of Tube I among the four geometric variables, respectively. T_c: temperature at the reactor center point; T_max: maximum temperature in the reactor.

CO conversion rate:

CO₂ selectivity:

C₁, and C₅₊ selectivity:

where C_i represent the C₁ and C₅₊.

C₅₊ yield: Y_C5+ = X_CO × S_C5+ × 100% (7)

where in and out represent the molar flow rates of CO, C₁, C₅₊ at the inlet and outlet, respectively, in units of mol s⁻¹.

3.1.1 Effect of length. As shown in Fig. 3(a), when R_L varies within the range of 30–400 mm, X_CO and S_CO2 increase concurrently with increasing R_L, while S_C5+ decreases monotonically. The concurrent elevation of X_CO and S_CO2 is attributed to the prolonged residence time of reactants with the extension of R_L, which enhances gas–solid contact efficiency. In contrast, the downward trend of S_C5+ stems from the increased probability of secondary hydrogenation, cracking, and readsorption of olefin intermediates in longer channels, which in turn strengthens the C₁–C₄ formation pathway. Meanwhile, the extended reactant residence time intensifies the WGS side reaction, which not only dilutes the effective syngas but also reduces the probability of long-chain hydrocarbon growth. Y_C5+ is defined as the product of X_CO and S_C5+. Due to the opposite variation trends of X_CO and S_C5+, the yield curve exhibits two distinct local maxima at R_L 120 mm and 250 mm, with corresponding values of 11.43% and 11.66%, respectively. Notably, Y_C5+ at R_L = 400 mm is higher than those at R_L = 120 mm and 250 mm; this phenomenon is primarily because the continuous increase in X_CO (driven by further prolonged residence time and enhanced gas–solid contact) outweighs the slight decrease in S_C5+ at longer R_L, leading to a net rise in Y_C5+. With respect to temperature variations, both T_max and T_c increase with the elevation of R_L, and the temperature difference between them continues to widen. When R_L exceeds 300 mm, T_c decreases slightly due to the reduced concentration of downstream reactants and the consequent slowdown in heat release rate. In contrast, T_max remains elevated in the high-temperature zone near the inlet due to thermal conduction lag, indicating that the axial temperature distribution has become increasingly uneven and is accompanied by an elevated risk of hot spot formation.

Regarding the geometric variable R_L (with typical values of 30, 70, and 120 mm), the results from Fig. 4 indicate that the distribution patterns of all field quantities are consistent across the three parameter levels: increasing R_L only expands the catalytic area without altering the channel structure, while the C₅₊ mass fraction (W_C5+) gradually increases along the axial direction with higher values inside the tube than outside. Local X_CO increases axially, with S_CO2 increasing and S_C5+ decreasing in the high-temperature zone at the rear section. Concurrently, the enhanced heat release leads to a simultaneous rise in T_max, which corresponds to the overall trends in Fig. 3(a). Additionally, as the distance from the center point to the hot spot at the terminal end increases with R_L, the temperature difference between T_max and T_c expands accordingly.

Fig. 4 Contour plots of W_C5+, X_CO and temperature at R_L of 40, 70, and 120 mm.

Consequently, the design boundary of R_L is defined as [30, 200] mm. This range not only ensures X_CO remains at a relatively high level but also prevents excessive decline in S_C5+ and temperature runaway:

(1) The lower limit of 30 mm avoids the influence of inlet–outlet backmixing caused by an excessively short reaction channel, which would otherwise lead to component fluctuations and notably low conversion. As presented in Fig. 3(a), Y_C5+ has decreased to its minimum value at R_L 30 mm.

(2) The upper limit of 200 mm covers the main peak range of Y_C5+. Further extension results in a slowing trend in the variations of X_CO and S_C5+, while ΔT_max tends to exceed 20 °C. Additionally, the reactor volume and catalyst dosage increase significantly, with the marginal yield gain benefit being far lower than the manufacturing and material costs.

3.1.2 Effect of inner diameter. As presented in Fig. 3(b), within the range of R_D from 0.5 to 50 mm, X_CO rises to a peak at 2.5 mm before declining, with a steep slope observed in the interval of 0.5 to 7 mm thus indicating that the reaction is highly sensitive to R_D. An increase in R_D directly changes the flow area ratio between the inner and outer channels of Tube I, leading to a reduced flow velocity inside Tube I and an increased one outside, and consequently resulting in an inverse variation in the residence time on both sides. Given that active components are loaded on both the inner and outer surfaces of Tube I, X_CO exhibits an opposite variation trend, with the peak at 2.5 mm. The variation trend of S_CO2 is consistent with that of X_CO, with its peak slightly shifted to 3.5 mm and tending to level off beyond 15 mm. This is because an increase in R_D simultaneously promotes both FTS and WGS reactions. However, in the descending segment, the reduction in CO partial pressure inhibits the progression of WGS. Additionally, S_C5+ remains nearly constant throughout the entire R_D range, indicating its insensitivity to flow velocity. Y_C5+ is affected by variations in X_CO, exhibiting a trend of first increasing and then decreasing with its maximum value located at 3 mm. The temperature profile exhibits a three-stage characteristic of “rise–decline–slow rise”: initially, increasing R_D prolongs the residence time inside Tube I, enhancing FTS exothermicity and concurrently raising T_max and T_c. With further R_D increases, the sharp drop in internal tube velocity shifts the main reaction zone outward, reducing internal heat release and lowering temperatures. Subsequently, concentrated reactions and relatively high velocity in the outer channel cause heat accumulation, leading to a gradual temperature rise. The outward migration of the hot spot moves T_c (measured at the central point) away from the high-temperature zone, and the temperature difference between T_c and T_max increases monotonically with R_D. When R_D is close to 20 mm, T_max is completely located outside the tube, T_c basically returns to the inlet temperature, and the temperature gradient in the central channel is significantly reduced.

For R_D values of 0.5, 2, and 3.5 mm, Fig. 5 reveals notable variations in the distribution of local W_C5+, X_CO, and temperature. Local W_C5+ distribution is highly dependent on R_D: at R_D = 0.5 mm, W_C5+ is only locally elevated at the inlet and then rapidly homogenizes, which is attributed to the combined effect of inlet backmixing and the extremely small inner diameter. By contrast, at R_D values of 2 mm and 3.5 mm, W_C5+ inside Tube I increases steadily along the axial direction, while that outside remains nearly constant. Local X_CO exhibits distinct axial profiles across R_D values: it maintains a high level inside Tube I for R_D = 0.5 mm, whereas it rises gradually along the axial direction at R_D values of 2 mm and 3.5 mm. This discrepancy arises because the 0.5 mm inner diameter is smaller than the default inter-tube spacing (R_S = 1 mm). Under laminar flow, the disparity between high central velocity and small R_D induces space velocity distribution reversal, altering local reaction extent. In contrast, the local X_CO outside Tube I shows similar distribution patterns across the three R_D cases. Temperature distributions are broadly consistent: temperature rises axially, and the temperature difference between the outlet center and the area outside Tube I diminishes with increasing R_D. At R_D = 0.5 mm, outlet tube internal temperature is higher than the external, whereas at 3.5 mm, internal and external temperatures are almost identical, verifying the trend in Fig. 3(b).

Fig. 5 Contour plots of W_C5+, X_CO and temperature at R_D of 0.5, 2, and 3.5 mm.

Accordingly, the design boundary of R_D is defined as [0.5, 4] mm, which encompasses the yield peak while suppressing excessive temperature rise. The rationale for this range is as follows:

(1) Lower limit (0.5 mm): an excessively small inner diameter not only hinders the loading of active components but also tends to cause blockage of the reaction channel.

(2) Upper limit (4 mm): X_CO and S_C5+ achieve their optimal values within the range of 2.5–3 mm. When R_D exceeds 4 mm, the yield decreases continuously, T_max rises again (increasing the risk of hot spots), and both material consumption and reactor volume increase rapidly. Furthermore, inner diameters greater than 4 mm fall outside the scope of microchannel specifications.

3.1.3 Effect of thickness. As illustrated in Fig. 3(c), within the R_T range of 0.25 to 35 mm, X_CO decreases almost linearly over 0–6.5 mm before the rate of decline moderates. For R_T ≤ 6.5 mm, the flow area outside Tube I shrinks sharply with the expansion of the outer diameter, with the increase in average flow velocity and reduction in residence time acting as the dominant factors. For R_T > 6.5 mm, the rate of velocity increase slows while the catalytic area on the outer surface expands synchronously, partially offsetting the residence time loss and thus decelerating the decrease in X_CO. S_CO2 declines alongside X_CO but stays at lower values, as WGS is also limited by residence time and has lower intrinsic activity than FTS. In contrast, S_C5+ rises gradually with increasing R_T, as short residence time inhibits secondary hydrogenation and cracking, thereby enhancing the chain propagation probability with the most pronounced gradient observed in 0–6.5 mm. These combined effects drive Y_C5+ to a slow downward trend, indicating that the adverse effect of shortened residence time on CO consumption is always more significant than the benefit from expanded catalyst area. The temperature profile decreases initially and then rises gently, with a trough at 11–15 mm: at low R_T, elevated flow velocity enhances convective heat dissipation, while at high R_T, extended metal heat conduction pathways lead to local heat accumulation and subsequent temperature recovery. Notably, the temperature rise phase does not bring about an improvement in yield, verifying that the “kinetics-residence time” mechanism rather than thermal effects is the decisive factor governing reactor performance.

Fig. 6 compares the overall and cross-sectional distribution characteristics of W_C5+, X_CO, and temperature within the reactor at R_T values of 0.25 mm, 1.5 mm, and 6.5 mm. Closer to the reactor outlet, the probability of re-adsorption and hydrogenation of intermediate olefins rises, leading to a corresponding increase in heavy hydrocarbon formation. Within the same cross-section, W_C5+ inside Tube I is consistently higher than that outside, and this difference widens as the axial position approaches the outlet. Additionally, W_C5+ outside becomes more uniform near the outlet. This is because the stable flow in the 2 mm-diameter inner tube facilitates product accumulation, whereas products generated outside the tube diffuse radially over a large range from the outer wall to the surroundings, which exerts a diluting effect on W_C5+. These trends confirm that the channel inside Tube I dominates C₅₊ formation, with larger R_T increasing W_C5+ in this channel and thus enhancing C₅₊ selectivity. Local X_CO also increases axially, with uniform cross-sectional distribution at low R_T. Increasing R_T narrows the spacing outside Tube I and raises space velocity, inducing a distinct radial gradient and thereby reducing the overall X_CO. Temperature rises gradually along the axial direction with uniform distribution inside and outside Tube I, attributed to the strong exothermicity of FTS and high axial thermal conductivity of the metal wall. As R_T increases, reduced X_CO lowers heat release and coupled with metal heat sink and elevated flow velocity synergistically inhibiting temperature rise, which leads to a decrease in the macroscopic temperature level.

Fig. 6 Contour plots of W_C5+, X_CO and temperature at R_T of 0.25, 1.5, and 6.5 mm.

Thus, the design boundary of R_T is set to [0.25, 6.5] mm, which not only covers the regions where X_CO and S_C5+ exhibit the highest sensitivity to changes but also avoids the mechano-thermal risks associated with excessively thin walls and diminishing returns caused by overly thick walls.

(1) Lower limit (0.25 mm): a further reduction in wall thickness would result in insufficient radial thermal resistance, making local hotspots unavoidable. It would also lead to inadequate mechanical strength, causing bending or collapse during catalyst loading and under pressure differences, which in turn triggers catalyst layer fracture and flow maldistribution.

(2) Upper limit (6.5 mm): for R_T values exceeding 6.5 mm, the negative effect of shortened residence time in the channel outside Tube I on X_CO is significantly mitigated. Although further increasing wall thickness expands the external catalytic surface area, its marginal contribution to Y_C5+ approaches zero. Meanwhile, metal heat capacity and heat conduction paths increase concomitantly, and the temperature profile rises after reaching its minimum at 11 to 15 mm, indicating potential thermal runaway risks for larger R_T. Additionally, material consumption and reactor weight increase markedly, which deviates from the microchannel design intent of compactness.

3.1.4 Effect of inter-tube spacing. As depicted in Fig. 3(d), within the R_S range of = [0.01, 20] mm, X_CO decreases almost linearly in the range of 0.01–2.5 mm and then slows down. This is because when the inter-tube spacing is extremely small and approximates the fluid boundary layer thickness, the diffusion zones of adjacent catalysts overlap significantly. CO undergoes repeated “recycling” reactions and the catalyst effectiveness factor remains high. As R_S increases, the overlapping zones shrink gradually and each Tube I operates independently. Some CO in the core region between tubes cannot diffuse transversely to the catalytic layer, causing the catalyst effectiveness factor to decline. Thus X_CO drops rapidly with increasing R_S. When R_S exceeds 5 mm, the catalyst effectiveness factor approaches the normal value of a single tube and the rate of X_CO decline levels off accordingly. S_CO2 decreases synchronously with X_CO but with a gentler trend. CO₂ is generated by WGS and its intrinsic activity is lower than that of FTS. It is also governed by the catalyst effectiveness factor. The S_CO2 curve flattens after 5 mm as the effectiveness factor stabilizes. S_C5+ increases first and then decreases. At extremely small R_S, high-velocity purging thins the boundary layer but the contact time is too short for α-olefins to undergo multiple adsorption and hydrogenation steps, so Y_C5+ remains at a low level. As R_S increases, the flow velocity is still relatively high but provides an adequate yet not excessive residence window for chain propagation. The probability of secondary hydrogenation and re-adsorption rises and S_C5+ increases accordingly. With further increases in R_S, prolonged residence time induces excessive hydrogenation and secondary cracking leading to product lightening, and S_C5+ declines gradually. Under the combined effect of X_CO and S_C5+, Y_C5+ exhibits a trend of first rising and then falling with its peak corresponding to 0.3 mm. T_max drops sharply first and then rises slowly. At extremely small R_S, the narrow slits between tubes force gas to sweep across the wall surface at high velocity resulting in intense forced convective heat transfer. Wall heat is removed rapidly and the temperature falls to its lowest point at 2.5 mm. As R_S continues to increase, the flow cross-section expands and gas velocity decreases weakening convective heat dissipation capacity. Heat generated by reactions accumulates continuously in the bed and the wall temperature rises gradually. T_c also decreases first and then rises but its lowest point lags to R_S = 5 mm. The center is far from the wall surface and its temperature response lags behind the wall due to radial heat conduction lag. When the spacing increases to 5 mm, the combined effect of extended heat conduction paths and weakened convection brings the central heat accumulation to a minimum equilibrium point. The temperature then rises with the decrease in heat dissipation capacity.

As presented in Fig. 7, three-dimensional contour plots for R_S = 0.1 mm, 2.5 mm and 5 mm are compared. The W_C5+ distribution shows a pronounced dependence on R_S. At R_S = 0.1 mm, the narrow inter-tube slit has extremely high flow velocity yet a small gas flow rate, facilitating rapid reactions and the formation of high C₅₊ concentrations. In contrast, Tube I's 2 mm inner channel and the gap between its outer wall and Tube II's inner wall feature low flow velocity and long residence time, driving steady axial W_C5+ accumulation. At R_S = 2.5 mm, the inter-tube spacing approximates Tube I's outer diameter and the external space velocity drops sharply, leading to a “high inlet–low outlet” W_C5+ profile. At the 40 mm inlet cross-section, W_C5+ is uniformly distributed inside Tube I while diffusing outward from its outer wall externally; by the 120 mm cross-section, the distribution homogenizes across both regions. A similar pattern occurs at R_S = 5 mm, but external W_C5+ remains non-uniform even at the 120 mm outlet cross-section, demonstrating that excessively large spacing reduces the utilization efficiency of catalysts outside Tube I.

Fig. 7 Contour plots of W_C5+, X_CO and temperature at R_S of 0.1, 2.5, and 5 mm.

Thus, the optimization boundary of R_S is set to [0.01, 2.5] mm, which covers the high productivity peak and temperature difference control region, balancing sufficient reaction depth while avoiding thermal accumulation induced by abnormal space velocity from excessively large spacing.

(1) Lower limit (0.01 mm): this corresponds to the narrowest inter-tube slit, featuring the highest cross-sectional flow velocity, a sharp rise in space velocity and extremely short residence time. Both X_CO and S_C5+ selectivity remain at low levels. Yet the temperature reaches its maximum value in this scenario, rendering the system prone to thermal runaway incidents.

(2) Upper limit (2.5 mm): this matches the inflection point where X_CO transitions from a steep decline to a gradual decrease. With further increases in R_S, space velocity decreases and heat dissipation weakens, causing wall temperature and central temperature to increase synchronously and elevating thermal runaway risks. Meanwhile, S_C5+ decreases monotonically, Y_C5+ has passed its peak and diminishing returns occur.

3.2 Surrogate model evaluation

Based on the boundary ranges of the four design variables determined in Section 3.1, Latin hypercube sampling was adopted to 250 design points and perform corresponding CFD simulations, establishing a foundational dataset. Of these, 150 samples accounting for 60% were used to train the KRG and GENN models. A further 50 samples representing 20% served as the validation set to iteratively optimize hyperparameters with NMAE and NRMSE as evaluation metrics. The remaining 50 samples making up 20% were reserved as the test set to independently verify the model generalization capability and eliminate overfitting. The final hyperparameters and evaluation metrics are listed in Table 2, and the data distribution and fitting performance are illustrated in Fig. 8. Fig. 8(a) and (b) present the predicted versus simulated value scatter plots and parity checks of KRG for f₁(x) and f₂(x) respectively, while Fig. 8(c) and (d) show the corresponding results for GENN. Analysis of these plots reveals that GENN yields lower NMAE and NRMSE values than KRG, indicating its higher prediction accuracy. Thus GENN was selected as the surrogate model for subsequent bi-objective optimization.

Table 2 Final hyperparameters of the surrogate model

Model	Hyper-parameter^a	Parameter description	Value	NMAE	NRMSE
a All hyperparameters are named following SMT interface functions. Unspecified hyperparameters take the default values of the functions.
KRG	poly	Polynomial order that defines the form of the mean function to fit global data trends	Second-order polynomial	9.27%	6.36%
	corr	Type of correlation function that characterizes local data correlation and model smoothness	Squared exponential correlation function
	theta0	Initial value of correlation function hyperparameter (θ) that initializes model correlation parameters	Four design variables, each initialized to 100
	theta_bounds	Optimization bounds of θ that constrain parameter range to avoid overfitting and underfitting	[1 × 10^−6, 1000]
	hyper_opt	Algorithm for hyperparameter optimization that realizes efficient tuning of model parameters	COBYLA (constrained optimization by linear approximation algorithm)
GENN	alpha	Learning rate of optimizer that controls parameter update step size to balance convergence	0.03 (moderate learning rate for stable convergence)	8.35%	4.57%
	hidden_layer_sizes	Number of hidden layers and neurons that determines model fitting capability	[30, 30] (2 hidden layers, 30 neurons per layer)
	num_iterations	Number of optimizer iterations that balances computational cost and prediction accuracy	1000 (sufficient iterations for model convergence)
	beta1/beta2	Adam optimizer's momentum decay coefficients that stabilize learning rate scaling	First-order: 0.95, second-order: 0.99
	lambd	Regularization coefficient that suppresses overfitting by penalizing excessive weights	0.02 (moderate regularization strength)
	is_normalize	Data normalization flag that improves model training stability and prediction accuracy	True

---

Fig. 8 Comparison of prediction results between KRG and GENN models and simulation results. (a) and (b) Respectively compare the KRG model predictions and CFD simulation results of f₁(x) and f₂(x); (c) and (d) respectively compare the GENN model predictions and CFD simulation results of f₁(x) and f₂(x).

Fig. 9 further validates the continuous predictive performance of the surrogate model across the entire design space. Fig. 9(a) presents the 1/Y_C5+ and ΔT distribution of the original 250 high-fidelity samples via scatter plots while Fig. 9(b) shows the GENN prediction results for 2000 Latin hypercube-sampled points within the identical boundary ranges. The predicted results fall within the same physical interval as the original data and achieve smooth transitions in sample-sparse regions without any non-physical oscillations or extrapolation distortions. This confirms that GENN not only accurately reproduces the training samples but also yields a continuous and differentiable predictive profile from discrete sampling points at negligible computational cost. It thus serves as an efficient and reliable fitness evaluation function for subsequent multi-objective optimization, enabling the Pareto frontier solution for the geometric parameters of the microchannel reactor.

Fig. 9 Original data and LHS sampling prediction data of GENN. (a) Shows the distribution of 250 sets of original data obtained via CFD before the GENN model prediction; (b) shows the distribution after prediction using the GENN model based on the 250 sets of original data.

3.3 Multi-objective optimization

3.3.1 Multi-objective optimization process. Table 3 lists the key parameter settings of NSGA-II while Fig. 10 presents the Pareto frontier of the multi-objective optimization, overlaid with the training, validation and test set samples, CFD samples for sensitivity analysis as well as CFD points invoked during the search. Points on the Pareto frontier curve are distributed in the lower-left region of the entire feasible domain, forming a non-dominated solution set that clearly reveals the trade-off relationship between Y_C5+ and ΔT_max.

Table 3 Main parameters of the NSGA-II algorithm

Model	Parameter^a	Parameter description	Value	Remark
a All parameters are named following PyMOO interface functions.
NSGA II	pop_size	Number of individuals per population generation	2000	Balances computational time and result accuracy
	n_offsprings	Number of offspring generated per generation	1000
	Crossover	Crossover operator	prob = 0.7, eta = 5	Crossover probability and distribution index are 0.7 and 5 respectively
	Mutation	Mutation operator	prob = 0.2, eta = 5	Mutation probability and distribution index are 0.2 and 5 respectively

---

Fig. 10 Pareto curve of multi-objective optimization.

Fig. 11 illustrates the evolution of corresponding geometric parameters after sorting Pareto frontier solutions by increasing 1/Y_C5+. The geometric characteristics of the leftmost region with high yield and high temperature rise are short channel length, small inner diameter, narrow inter-tube spacing and thin wall thickness. This indicates that small-scale geometries exhibit a pronounced reaction enhancement effect. However, no external heat dissipation is applied and heat removal relies solely on reactant flow, resulting in a high ΔT_max value. As 1/Y_C5+ increases, R_L first rises and then decreases, remaining generally within the range of 90–125 mm. This is the outcome of the continuous trade-off between yield and temperature rise. R_D grows gradually, which aligns with the sensitivity analysis results in Section 3.1 where increasing R_D within this interval improves Y_C5+ and reduces ΔT_max. In the left region where ΔT > 9 °C, R_T increases rapidly from 1.05 mm to 6 mm, demonstrating that wall thickness plays a dominant role in temperature suppression. Beyond this region where ΔT < 9 °C, further thickening yields minimal cooling effects. R_S is only sensitive in the region where Y_C5+ exceeds 15.4% and ΔT < 9 °C, rising first to 0.7 mm before declining. In all other low-temperature-rise regions, R_S takes the minimum value of 0.05 mm. This suggests that adopting a relatively small inter-tube spacing is beneficial for both high yield and low temperature rise simultaneously.

Fig. 11 Relationship between design variables (geometric parameters) distribution on the Pareto frontier and output results.

3.3.2 Determination of optimal solutions. Final design schemes can be selected from the Pareto frontier based on productivity and thermal safety requirements. If the temperature control margin is large, left-side solutions can be chosen to prioritize high yield. If strict hotspot control is required, solutions should be shifted rightward to sacrifice partial yield for milder thermal conditions. Table 4 lists several typical solutions. In comparison with the initial geometric schemes, the optimized schemes achieve marked improvements in both yield and temperature difference performance. When ΔT_max is reduced by 4.95 °C and 7.20 °C respectively, yield increases by 213% and 86%. Furthermore, at equivalent yield levels, ΔT_max can be further decreased by 9.06 °C.

Table 4 Initial design and typical Pareto optimal designs

Design schemes	Design variables/mm				Output results
	RL	RD	RT	RS	YC5+/%	ΔT/°C
Initial design	100	2	0.5	0	10.71	17.31
S1	97.01	0.53	0.58	1.12	33.55	12.36
S2	112.05	0.52	0.50	4.01	19.89	10.11
S3	112.51	0.05	1.38	6.12	10.71	8.25

---

When selecting the optimal solution, in addition to reaction performance, considerations must be given to manufacturing tolerances, support structure strength and catalyst deactivation temperature limits. Following comprehensive evaluation, Scheme S2 was identified as the final design. The geometric dimensions of this solution were validated via CFD simulations, yielding an actual yield of 18.72% and a temperature rise of 12.87 °C. The relative errors compared with the surrogate model predictions were 6.25% and 0.96% respectively, both falling within the model uncertainty range and confirming the model's reliability. Scheme S2 achieves a 1.86-fold increase in yield while reducing the maximum temperature rise by 7.20 °C, thus laying a foundation for pilot-scale scaling-up and subsequent industrial application.

Since the optimal scheme is obtained under fixed operating conditions, the Pareto frontier will shift when the operating conditions change. Regarding the effects of various operating conditions on Y_C5+ and ΔT_max of the reactor with initial dimensions, these have been discussed in ref. 16. The operating conditions adopted in this study are those that generate relatively favourable reaction results. However, since Y_C5+ and ΔT_max do not exhibit a simple monotonic change with variations in operating conditions, it is difficult to directly determine whether the current optimal scheme remains optimal under other operating conditions. An effective strategy is to treat the four operating conditions as design variables and construct an 8-variable input and 2-variable output surrogate model alongside the dimensional variables. Nevertheless, this study is limited by the length of the manuscript and insufficient initial data. Future in-depth investigations can be performed once sufficient reaction data under diverse operating conditions are acquired in follow-up studies.

4. Conclusions

This study presents a multi-objective optimization design approach for microchannel reactors based on the integration of CFD, GENN and NSGA-II. It investigates the reactor performance of FTS microchannel reactors under different geometric parameters. Specifically, sensitivity analysis reveals that all four-dimensional design variables exert significant impacts on reaction outcomes, each with an optimal range that requires balancing during optimization. The developed GENN surrogate model exhibits excellent predictive accuracy and generalization capability, maintaining simulation precision while significantly reducing CFD computational load. Notably, the NSGA-II algorithm effectively obtained the Pareto optimal solution set, achieving the dual-objective optimization of maximizing Y_C5+ while controlling ΔT_max. Ultimately, the selected optimal scheme reduces ΔT_max by 7.20 °C and increases Y_C5+ by 1.86-fold. This optimized design exhibits substantial potential for engineering applications, providing technical support for the practical deployment of microchannel reactors in FTS.

It should be noted that the present optimization focused exclusively on Y_C5+ maximization and ΔT_max minimization. Manufacturing cost, reactor volume, catalyst deactivation and pressure drop were reserved for future investigation because the multi-tube microchannel reactor examined in this study remains at the laboratory scale and reliable surrogate modeling of economic or durability metrics requires extensive pilot or industrial data. Comprehensive multi-objective studies that incorporate these additional criteria will be pursued once the requisite operational datasets are generated during the subsequent pilot stage.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: the file Python_Code.txt contains the Python code for constructing surrogate models and implementing the multi-objective optimization process. The file data to plot.xlsx includes the data for some figures presented in the paper. See DOI: https://doi.org/10.1039/d5ra09551d.

Acknowledgements

We thank Professor Wang Yuanyang for his invaluable guidance and support throughout this research. We are grateful to Wang Yanqian for their assistance with the experimental work. We also thank the anonymous reviewers for their constructive comments and suggestions that helped improve the manuscript.

References

1. J. E. Apolinar-Hernandez, S. L. Bertoli, H. G. Riella, C. Soares and N. Padoin, Energy Fuels, 2023, 38, 1–28.
2. Z. Teimouri, V. B. Borugadda, A. K. Dalai and N. Abatzoglou, Renewable Sustainable Energy Rev., 2022, 160, 112287.
3. H. Cao, R. Xu, X. Tang, T. Yang, S. Hou and C. Hou, Chin. J. Chem. Eng., 2023, 64, 224–240.
4. K. Lee and J. M. Lee, Comput. Chem. Eng., 2024, 185, 108658.
5. Z. Teimouri, N. Abatzoglou and A. K. Dalai, Energy Convers. Manage., 2023, 296, 117659.
6. J. Na, K. S. Kshetrimayum, U. Lee and C. Han, Chem. Eng. J., 2016, 1521–1534.
7. J. Kong, J. P. Eason, X. Chen and L. T. Biegler, Ind. Eng. Chem. Res., 2020, 59, 9165–9179.
8. P. Mishra and A. Yadav, Int. J. Hydrogen Energy, 2025, 164, 150726.
9. Y. Qiu and P. Liu, Case Stud. Therm. Eng., 2025, 70, 106093.
10. I. Jung, J. Na, S. Park, J. Jeon, Y.-G. Mo, J.-Y. Yi, J.-T. Chung and C. Han, Fuel Process. Technol., 2017, 159, 448–459.
11. The Ansys Product Release and Updates Online, 2024, https://www.ansys.com/products/release-high-lights, accessed December 2024.
12. T. Chen, R. Li, X. Hu, B. Zhang, Y. Liu, L. Wang and N. Gao, Build. Environ., 2025, 285, 113667.
13. Y.-T. Liu, C.-Y. Wu, T. Chen and Y. Yao, in Computer Aided Chemical Engineering, ed. A. C. Kokossis, M. C. Georgiadis and E. Pistikopoulos, Elsevier, 2023, vol. 52, pp. 57–63.
14. Z. Jiang, J. Chen, S. Xie, X. Li, H. Liu, L. Wang, C. Hong, G. Li, H. Li and K. Chen, Chem. Eng. Sci., 2025, 311, 121614.
15. L. Zhang, L. Chen, S. Xia, Y. Ge, C. Wang and H. Feng, Int. J. Heat Mass Transfer, 2020, 148, 119025.
16. S. Ren and Y. Wang, RSC Adv., 2025, 15, 13137–13151.
17. Y. Wang, Research on synthesis of low carbon hydrocarbons by intedrated Fischer-Tropsch reaction of microreactor and Fe based catalyst, Master's thesis, Taiyuan University of Science and Technology, 2022,  DOI:10.27721/d.cnki.gyzjc.2022.000362.
18. Z. Yu, Z. Wang, H. Zhong and K. Cheng, RSC Adv., 2024, 14, 28724–28739.
19. P. Saves, R. Lafage, N. Bartoli, Y. Diouane, J. Bussemaker, T. Lefebvre, J. T. Hwang, J. Morlier and J. R. R. A. Martins, Adv. Eng. Software, 2024, 188, 103571.
20. J. Blank and K. Deb, IEEE Access, 2020, 8, 89497–89509.
21. J. Sacks, S. B. Schiller and W. J. Welch, Technometrics, 1989, 31, 41–47.
22. S. H. Berguin, arXiv, 2024, preprint, arXiv:2406.09132,  DOI:10.48550/arXiv.2406.09132.
23. M. E. Dry, Catal. Today, 2002, 71, 227–241.
24. J. S. Park, D. E. Kim, Y. J. Lee, G. Kwak and M. J. Park, Ind. Eng. Chem. Res., 2016, 55, 9416–9425.
25. K. Deb, A. Pratap, S. Agarwal and T. Meyarivan, IEEE Trans. Evol. Comput., 2002, 6, 182–197.

---