Machine learning-assisted discovery of flow reactor designs | Nature Chemical Engineering

Nature Chemical Engineering volume 1, pages 522–531 (2024)Cite this article

6468 Accesses

2 Citations

16 Altmetric

Metrics details

Additive manufacturing has enabled the fabrication of advanced reactor geometries, permitting larger, more complex design spaces. Identifying promising configurations within such spaces presents a significant challenge for current approaches. Furthermore, existing parameterizations of reactor geometries are low dimensional with expensive optimization, limiting more complex solutions. To address this challenge, we have established a machine learning-assisted approach for the design of new chemical reactors, combining the application of high-dimensional parameterizations, computational fluid dynamics and multi-fidelity Bayesian optimization. We associate the development of mixing-enhancing vortical flow structures in coiled reactors with performance and used our approach to identify the key characteristics of optimal designs. By appealing to the principles of fluid dynamics, we rationalized the selection of design features that lead to experimental plug flow performance improvements of ~60% compared with conventional designs. Our results demonstrate that coupling advanced manufacturing techniques with ‘augmented intelligence’ approaches can give rise to reactor designs with enhanced performance.

Advances in additive manufacturing have enabled the fabrication of a wide range of complex and potentially counter-intuitive reactor designs. Previously infeasible or highly impractical designs can now be manufactured and investigated, resulting in substantially larger design spaces. Coupled with data-driven design tools, such as multi-fidelity Bayesian approaches1,2,3,4,5,6,7, which have enabled the optimization of large-scale simulation-based problems, a pathway has emerged toward the identification of optimal reactors from larger design spaces with the potential for improved performance. By exploiting lower-fidelity simulations throughout optimization, high-quality solutions can be generated in a significantly reduced time. This is particularly the case in scenarios where gradients are unavailable or a more global solution is desired than can be achieved through gradient-based approaches1,5.

Here we demonstrate an ‘augmented intelligence’ framework that enables the design of reactors that surpass the performance of conventional designs. To do so, we have taken advantage of developments in data-driven optimization, machine learning, computational fluid dynamics and additive manufacturing; the coiled-tube reactor was chosen as an illustrative example.

Coiled-tube reactors have attracted the attention of chemical engineers owing to their desirable mixing and heat transfer characteristics8,9,10,11,12,13. Their applications range from flow chemistry14 and bioprocesses15 to chemical kinetic experiments11. At the mesoscale, coiled-tube reactors have been shown to combine the heat and mass transfer properties of microreactors with the economic benefits of high-throughput larger-scale reactors16,17. In previous work, improvements were made in the plug flow performance of coiled-tube reactors at relatively low flow rates (characterized by Reynolds numbers (Re) of ≤50) by superimposing pulsed-flow operating conditions, which induce Dean vortices9,18,19 that enhance radial mixing. These vortices arise from the centrifugal forces that develop in a curved pipe and shift the location of maximal velocity toward the outer pipe wall with a return flow toward the inner curved wall leading to the establishment of counter-rotating cells. In the absence of pulsed-flow conditions, such vortices develop at much higher flow rates and Dean numbers (Re > 300, Dean number (De) > 75)20,21,22. It is essential to induce Dean vortices at low flow rates and Dean numbers over large proportions of the reactor interior under steady-state flow without the added overhead of pulsed-flow forcing.

We sought to enhance the plug flow performance of coiled-tube reactors by identifying two parameterizations for the reactor geometry in the radial and axial directions. We solved the simulation-based optimization problem for the two parameterizations using multi-fidelity Bayesian optimization, considering steady-state flow conditions at a low flow rate (Re = 50, for which a pulsed flow would have been necessary to drive vortex formation). The composite objective that we maximized consisted of plug flow performance, which we approximated from computational residence time distributions using a tanks-in-series model, and a non-ideality term that penalizes bimodal or unsymmetrical distributions. Optimal solutions were investigated that identified the presence of fully developed Dean vortices at low Reynolds numbers under steady-state flow conditions. The key driving factors that result in improved performance were identified, which we then used to design two reactors. We 3D-printed and experimentally validated these designs, confirming their improved performance compared with a conventional coiled-tube reactor for both tracer and reacting flow experiments.

Our work has established a framework for the design of next-generation reactors that can potentially improve the performance, sustainability and economic viability of various manufacturing processes. Ultimately, the aim of this work was to shift the paradigm of reactor design to take advantage of the suite of modern computational methodologies in design and optimization, demonstrating new opportunities to support discovery, innovation and advancement across chemical engineering.

Each parameterization was optimized under steady-state flow conditions with a Reynolds number of 50. Details of both parameterizations are provided in Methods. To ensure tractability of the problem, a sequential strategy was used for the joint parameterization, whereby the parameters of the optimal coil path design were fixed before the introduction of parameters to manipulate the cross-section. Figure 1a shows the optimal geometries for each parameterization. We denote the control reactor as design 1 and reactors resulting from optimal cross-section, coil path and combined parameterization as designs 2, 3 and 4, respectively (Fig. 1a).

a, The conventional coiled-tube reactor (i) alongside the optimal coiled-tube reactors generated by parameterizing the cross-section (ii), the coil path (iii) and a joint parameterization (iv). b, Streamlines colored according to their velocity magnitude within a standard coil (i) and a reactor with optimal joint parameterization (ii). c, Secondary flow streamlines at various cross-sections of a standard coil (i) and a reactor with optimal joint parameterization (ii). d, Cross-sectional planes across a standard coil (i) and a reactor with optimal joint parameterization (ii), demonstrating streamwise velocities. e, Comparison of induced Dean vortices within a standard coil (i) and the optimal joint parameterization coil (ii).

We first considered design 2, where the shape of the tube cross-section is allowed to vary along the length of the reactor. The first key feature of design 2 is that the cross-section undergoes periodic expansions and contractions approximately every half turn. Secondly, the design includes a pinch constricting the flow where the cross-sectional area is greatest, during the expansion phase. Next, we considered the optimal solution resulting from the parameterization of the coil path corresponding to design 3, where the path deviates from a nominal configuration through the interpolation of cylindrical coordinates. The radius of curvature of the coil of design 3 starts relatively large and then reduces along the length of the reactor, resulting in a tighter design than in reactor 2. The pitch of the coil in design 3 starts small before increasing approximately halfway along the length of the reactor and then decreasing near the reactor outlet to the extent that the coil path points downward. In design 4, the path is fixed as in design 3 and the cross-section is allowed to vary; inspection of design 4 reveals that it possesses features that are present in both designs 2 and 3. To explain why these solutions are deemed optimal, we demonstrate below their different flow characteristics. Figure 1b shows the differences in the fluid velocity in the conventional coil (design 1) and design 4. The conventional coil exhibits a relatively uniform velocity distribution along the length of the coil. Conversely, the expansion and contraction features in the cross-section in design 4 alter the velocity distribution along the length of the coil, resulting in higher and lower velocities during contraction and expansion, respectively. These changes in velocity in design 4 due to acceleration and deceleration induce stronger pressure gradients, leading to the formation of Dean vortices that significantly enhance radial mixing. Figure 1c(i,ii) demonstrates the formation of Dean vortices in both the standard coil and design 4, respectively. While Dean vortices are also formed in the standard coil, they are only partially established close to the tube outlet. The earlier formation of fully developed Dean vortices in design 4 demonstrates the impact of reactor geometry resulting from the application of our framework, suggesting the potential for enhanced plug flow performance within more compact reactors at lower Re.

Promoting plug flow depends on the relationship between axial and radial fluid mixing. To control axial dispersion, it becomes imperative to consider the streamwise velocity distribution within the coil. Therefore, we investigated the velocity distribution across cross-sectional areas along the length of designs 1 and 4 (Fig. 1d(i,ii), respectively). In design 1, we observed that the peak velocity is closer to the outer walls, with lower velocities near the inner wall along the length of the coil. This configuration induces faster tracer movement along the outer walls and slower movement along the inner walls, resulting in an extended time for the tracer to exit the coil and increased axial dispersion. By contrast, the inclusion of the pinch feature in design 4 has a key role in redistributing velocity across the coil cross-section by altering the radial position of the peak velocity along the length of the coil. Thus, design 4 shows peak velocity at various radial positions along the coil: at the top, bottom and toward the outer wall (Fig. 1d(ii)). This leads to the acceleration and deceleration of the tracer movement along the length of the coil, promoting a relatively uniform axial movement of the tracer over time. This results in a narrow profile for the residence time distribution of the tracer.

To further illustrate the enhanced characteristics of the optimized designs compared with conventional coil tube reactors, Fig. 1e(ii) shows streamlines representing fluid flow within the initial length of design 2. Blue streamlines originate at the coil center, while black streamlines start near the coil wall. As the cross-section of the reactor expands along the initial curve, the slow-moving fluid radially furthest from the center of the coil initially moves outward. Then, at the point where the cross-sectional area is at its greatest, the fluid in the central region moves toward the outer walls of the coiled tube. Where the pinch reduces the tube cross-section, the fluid closest to the outer walls is forced toward the center of the tube via a swirling motion. The cross-section then contracts again, enabling a repetition of the mechanism of radial mixing under steady-state flow conditions.

We next focused on the convergence of the designs from an optimization perspective to highlight the potential optimality of the solutions provided. Gaussian processes (GPs) were used to model the simulation cost and objective throughout the design space, and were iteratively updated on the basis of simulations selected through a multi-fidelity acquisition function. Full details of the optimization are provided in Methods. Given the high-dimensional design space, we applied t-distributed stochastic neighbor embedding (t-SNE) to identify emerging trends throughout the optimization, investigate design convergence and highlight the behavior of the GP hyperparameters (Methods). This is a probabilistic dimensionality reduction technique that enables high-dimensional data to be plotted while preserving local structure. We focused here on the convergence of the cross-section design problem, with the trends associated with the other parameterizations reported in the Supplementary Information. We found that the cross-section optimization problem has better convergence properties than the other two design problems, supporting our methodology of identifying the optimal driving characteristics rather than the designs themselves. Figure 2a shows plots of the composite objective function against iteration and wall-clock time, demonstrating how lower-cost simulations can be applied to enable less expensive exploration.

a, The optimization objective as a function of iteration (i) and wall-clock time (ii). The objective has been normalized. The initial dataset generated via ‘design of experiments’ is presented as negative iterations and wall-clock time. Hence, optimization begins at iteration 0 and t = 0. b, GP dimension length scales throughout the Bayesian optimization iterations (i) and histograms demonstrating the distribution of GP length scales changing throughout the optimization (ii). Shading represents how the frequency of certain length scales has changed from previous iterations to the current iteration. c, t-SNE analysis of the data generated throughout the optimization in design parameter space (\({{{\mathcal{X}}}}\)) for axial and radial fidelities, objective value and iteration number, reducing the dimensionality of the design parameters to two dimensions. Feature (i) shows the initial dataset generated before the optimization begins. Feature (ii) shows the path along which the optimization progresses. Feature (iii) shows the different fidelities that were evaluated along this path. Data points are colored according to the number of properties. d, Parameter variability, defined as a function of the GP length scale corresponding to each parameter, plotted for each inducing parameter on a nominal coil.

The initial design of the experiments used in the optimization is represented as occurring before the first iteration (the so-called negative iterations) and before time t = 0. This sample contains reactor simulations with a variety of computational costs as simulations across the spectrum of fidelities are performed. After the start of the optimization, low-cost, low-fidelity simulations are selected during the initial phase of the exploration. Simulations are selected that are progressively less expensive until approximately the 20th iteration as the updated GPs within the framework gain a better representation of simulation cost and are able to select more efficient simulations. The framework then undergoes a period in which either the highest fidelity simulations or the lowest fidelity simulations are largely applied.

To investigate the properties of the GP that model the objective throughout the optimization, we analysed the kernel function hyperparameters at a number of iterations. Figure 2b(i) shows the length scale of each dimension within the kernel function within the objective GP as the optimization progresses. Initially, the majority of length scales are relatively large, indicating a uniform function space, capturing broad trends, as expected in a low-data regime in high dimensions. As the optimization progresses and the GP is trained with more data, the length scales broadly decrease, indicating an improved representation of the design space. Figure 2b(ii) shows the distribution of length-scale hyperparameters as the optimization progresses. The distribution of hyperparameters tends toward lower values and the GP becomes more non-linear as data become available and correlations are more accurately captured.

As different simulation fidelities bias the objective function and the parameterizations are high dimensional, we applied t-SNE to the convergence data of the cross-section parameterization, observing data in the space of design parameters but not fidelities. Figure 2c shows the two-dimensional t-SNE embeddings of data points \(\in {{{\mathcal{X}}}}\) for axial and radial fidelities, objective value and iteration number. Feature (i) shows the initial dataset generated before the optimization begins. From here, the optimization progresses along the path labeled ii, demonstrating systematic changes in parameters indicative of convergent behavior. Along this path, a number of different fidelities can be evaluated (feature (iii)), demonstrating how lower-fidelity simulations are applied.

We describe here a sample of inducing cross-sections from various clusters of data within the analysis against the optimization objective. The cross-sections in the original dataset (feature (i)), as expected, contain no distinct characteristics. As the optimization progresses, the cross-sections defined by inducing points become more distinct, with certain cross-sections gaining symmetrical forms. Finally, at the end of optimization, the cross-sections are most distinct, consisting of alternating small and large cross-sectional areas with pinches throughout.

Lastly, we considered the parameter variability, which we define as the normalized inverse of each GP length scale, each corresponding to a specific parameter. We calculated the parameter variability at the end of the optimization and present the value for each inducing parameter (the collection of which defines the form of the design space), indicating the variability in the objective for which each specific parameter is responsible. As can be seen from Fig. 2d, which shows the parameter variability throughout the reactor, the variability resulting from parameters in the tube cross-section is greater where pinches occur than at the bottom and top of the tubular cross-section. These trends indicate that the tube cross-section, and specifically the pinch, has a significant effect on plug flow performance, confirming our observations.

In summary, the design concepts identified through the use of data-driven optimization and machine learning are as follows:

Expansion and subsequent contraction of the reactor cross-section every half turn

Distinct pinches in the cross-section throughout the reactor

Changes in the direction of flow, including constriction of coil radius

Previous work has demonstrated the existence of a linear relationship between the number of coils and plug flow performance23. Therefore, we next applied the design concepts to two reactors with longer coils than those studied above. In one reactor, we included changes to the tube cross-section, and in the other reactor, we included variations in both the cross-section and coil path. Figure 3 shows the final coiled-tube reactor designs containing features identified as being responsible for the optimality of both the cross-section and coil path parameterizations.

The nominal coiled tube (left, R1) alongside extrapolated coil-tube designs containing optimal cross-section (center, R2) and optimal cross-section and coil path (right, R3) features.

The designs shown in Fig. 3 were 3D-printed and residence time distribution experiments (Methods) were performed for each configuration using a steady-state flow with Re = 50. Figure 4 shows the aggregated experimental data together with tank-in-series models for each reactor configuration. The standard coil (R1) produces a skewed distribution with a mean equivalent tanks-in-series value (NR1) across all experiments of 39.27. The two distributions resulting from the reactors with variable cross-section (R2) and variable cross-section and coil path (R3) are symmetrical. The equivalent tanks-in-series values (NR2 and NR3) for these two reactors are higher than that of the standard coil at 61 and 63.45, representing a 55% and 62% higher value, respectively. The residence time distribution (RTD) resulting from R1 is to be expected as typically a high Re is required for Dean vortices to form. The improvements in R2 and R3 demonstrate that we can avoid the need for oscillatory conditions to induce Dean vortices and promote mixing at low Re through the inclusion of the identified design concepts. Other factors, such as multiphase flows and chemistry, may be explored to further exploit the discovery of designs through machine learning. An uncertainty analysis of the experimental data can be found in the Supplementary Information. A visual representation of the overall framework is provided in Fig. 5.

a–c, Dimensionless concentration at the outlet of the reactor (E(θ)) plotted against dimensionless time (θ) for the control reactor R1 (a), reactor R2 with variable cross-section (b) and reactor R3 with variable cross-section and coil path (c). Data were generated across three experiments for each reactor configuration, designed to maximize equivalent tanks-in-series values (N) and minimize non-symmetrical RTDs. A tanks-in-series model was used to estimate the performance of each reactor across the experiments performed.

Source data

To further validate the effectiveness of the optimized reactor designs, we performed the Villermaux–Dushman reaction24, a widely used test reaction for characterizing mixing performance in reactors. The Villermaux–Dushman reaction is particularly suitable for this purpose as it is a fast, mixing-sensitive reaction that can provide insights into the mixing efficiency of reactors25. The Villermaux–Dushman reaction system is as follows:

Under poor mixing conditions, reaction (2) is prohibited by the lack of acid, which is consumed within the faster first reaction. Good mixing conditions result in the formation of I2, which in turn reacts with I− to produce \({{{{{\rm{I}}}}}_{3}}^{\!\!-}\), measured by spectrophotometry at 353 nm (ref. 25).

The absorbance values at 353 nm, indicative of the presence of triiodide, measured at the outlet of each reactor are presented in Table 1. A higher absorbance value indicates a higher conversion of the limiting reagent. The optimized reactors R2 and R3 exhibit higher average absorbance values (0.293 and 0.325, respectively) than the nominal design R1 (0.273), demonstrating improved reactive flow performance at Re = 50. These results confirm that for mass transfer-limited reactions such as the Villermaux–Dushman reaction, optimizing the RTD serves as an effective proxy for enhancing overall reactor performance.

We note that the power consumption of the experimental syringe pumps was not measured. However, the pressure drop in the standard and optimized coils can be analysed computationally via computational fluid dynamics (CFD) simulations. For a single coil turn, the optimized geometry has a 29.3% higher pressure drop than the standard coil. This difference would further increase with coil length. Therefore, the proposed designs almost certainly have a larger pressure drop than the standard coil. In future work, this or similar cases may be treated as a multi-objective black-box optimization problem, trading off pressure drop and plug flow performance.

Throughout the optimization, we could not guarantee that the global optimum of the parameterizations had been identified. However, the global methodology that we applied provides more varied solutions than a gradient-based local approach such as the adjoint method, and we identified the key features that drive the behavior of the fluid within the reactors, including the presence of induced Dean vortices. The reactors that we then experimentally validated were interpreted as ‘augmented intelligence’-assisted designs, with features based on the optimal solutions. We propose this workflow for the design of highly parameterized reactors as the interpretability of solutions is maintained. To support the design of algorithms for the optimization of high-dimensional, expensive black-box problems, we have released all parameterizations as reactor design benchmark problems (deposited at https://github.com/trsav/reactor_benchmark). The parameterizations that we optimized contain feasible reactor geometries by design, resulting in an unconstrained optimization problem. Highlighting the importance of parameterization specification, the emergent behavior in designs 3 and 4 contains a coil path that unintendedly pitches down near the reactor outlet. It is worth noting that features of the optimal coil path design bear some resemblance to the wavering coiled flow inverter design proposed previously26, although this design primarily relies on inversions within the coil. In future work, we may perform comparisons over a wide range of existing and alternative coil designs.

We anticipate that the workflow for machine learning-assisted discovery that we have presented here could be used across a large number of expensive simulation-based design and optimization problems involving, for instance, reactive and multiphase flows. Alongside developments in additive manufacturing, we hope that chemical reactors discovered via machine learning become increasingly prevalent to support future sustainable processes.

In this section, we focus on the design of efficient parameterizations for the discovery of coiled-tube reactors, ensuring tractability of the design problem. We present two distinct parameterizations: one where the geometry of the cross-section varies along the length of the reactor and the other where the reactor path itself is allowed to vary. To manage the trade-off between flexibility and viability, the two parameterizations are defined by hyperparameters that determine their complexity. To ensure smooth transitions in these parameterizations, we used interpolation points in their formulation to define both the reactor path and the cross-section along the reactor length. We first introduce GPs in polar coordinates as a means to generate variable tube cross-sections from interpolation points.

A GP is an infinite-dimension generalization of a multivariate Gaussian distribution27. The mean vector and covariance matrix are replaced by a mean function and kernel function, respectively. A GP can be described as

The kernel function k dictates the behavior of functions from this distribution and can be parameterized by hyperparameters, including length scale. By conditioning a GP on a dataset X*, a posterior distribution of functions can be obtained. At inputs X and previously evaluated function values y, the posterior predictive mean (μf) and standard deviation (\({\sigma^2_f}\)) become

where K is a covariance matrix derived from kernel function k.

The squared exponential kernel function assigns a decreasing correlation between locations in input space with increasing Euclidian distance, providing an intuitive interpretation for the majority of regression settings. Other kernel functions have been proposed to provide valid covariance matrices, resulting in GP prior distributions with different properties. The polar squared exponential kernel enables valid covariance matrices to be constructed in polar coordinates, as outlined previously28. A standard kernel function, dealing with data within data in polar coordinates, will determine that data at θ1 = 0 and θ2 = π/2 are highly uncorrelated. This is untrue, and in the presence of noiseless observations, these two data points should have a perfect correlation (k(θ1, θ2) = 1). Polar kernel functions enable smooth interpolation in polar coordinates by including the ability for proper distances and respective covariances to be calculated between any two angles28,29. The polar covariance function is written as

and the angular distance metric (d) is given as

where τ is a hyperparameter analogous to the length scale and controls the smoothness of the prior distribution of functions. Figure 5a demonstrates samples from a GP prior with polar kernel, as well as a GP with polar kernel posterior distribution conditioned on data.

a, Left: polar GP parameterization of coil cross-section, demonstrating samples from a GP prior with a polar kernel. Right: posterior distribution after the polar prior distribution has been conditioned on data. In this demonstration, we assumed noiseless observations. b, Examples of sets of inducing points at given locations along the length of a coil. A polar GP was used to interpolate between points, where the radial value of each data point is a parameter. c, A flow diagram showing the main aspects of the methodology developed in this study. MFBO, multi-fidelity Bayesian optimization.

Initially, we defined the number of interpolating points for a given tubular cross-section, denoted nc and indexed by i. We then distributed these points equally along the angle θ in polar coordinates. The radius coordinate, ri ∈ nc, for each inducing point serves as the decision variable or parameter. Given a set of inducing points corresponding to a specific cross-section, a polar GP is used to interpolate between them, resulting in a valid, continuous curve in polar coordinates defining an individual cross-section. Next, we established the number of interpolating cross-sections, denoted nl and indexed by j, equally spaced along the length of the coil. Both nc and nl represent hyperparameters; increasing either one raises the dimensionality of the resulting design problem while achieving a more flexible parameterization. Figure 5b shows the polar GP interpolation of each specific cross-section along the length of the coil, with nc = 6 and nl = 6. To ensure compatibility with existing fixtures and fittings, two additional cross-sections were defined at the beginning and end of the coil, with constant radius.

Next, we placed each cross-section along the defined reactor path in three-dimensional (3D) space and rotated each cross-section to face perpendicular to the direction of flow. The resulting interpolated cross-sections are defined by a quadratic interpolation between each polar GP posterior in cylindrical coordinates. Finally, an inlet and outlet were added to the reactor by extending the cross-sections at the beginning and end of the reactor.

The data generated from the parameterization are provided to code to mesh tubular reactors (deposited at https://github.com/OptiMaL-PSE-Lab/pulsed-reactor-optimizsation). Importantly, the meshing procedure enables control over the number of cells throughout the axial direction of the reactor (the direction of fluid flow) as well as the radial direction. Defined by axial and radial fidelity parameters, respectively, both values dictate the factor to which blocks are subdivided during mesh creation. By maintaining the ability to adapt the fidelity of the mesh in two independent directions, we can provide greater scope for identifying efficient (regarding information gained on a computational expense basis) simulations. Previous work has demonstrated how the output of coiled-tube reactor simulations varies approximately smoothly with changes in fidelities1,30. In this study, we assumed that they are continuous throughout the modeling and optimization and rounded them to the nearest integer when stored or evaluated. In previous work, various fidelities were validated against experimental data, confirming the validity of the high-fidelity model with respect to cell count and the relationship between RTDs for lower-fidelity simulations was explored1.

For the parameterization of the coil path, we started by defining a baseline path in cylindrical coordinates (ρ0, ϕ0 and z0, which represent the radial, angular and axial coordinates, respectively). This path was based on a standard coil configuration and served as a reference for the introduction of induced deviations. The terms Δρj and Δzj represent the deviations in the radial and height directions, respectively, for each inducing point j. In optimizing for inducing point deviations, we ensured that we maintained the coil’s non-intersection constraint by design.

We introduced decision variables as mismatch terms into the baseline path, resulting in a parameterized path given by (ρj, ϕj, zj) = (ρ0 + Δρj, ϕ0, z0 + Δzj), where Δρj and Δzj represent the magnitude of the mismatches in the radial and axial directions, respectively. After defining the path, we placed circles defining the constant cross-section in 3D space parallel to the direction of flow. np represents the number of inducing points along the coil path, indexed by j. This parameterization is similar to the one presented previously31. We defined six inducing points in both the axial direction and the radial direction, that is, np = 6, but none in the rotational axis θ. Deviations in θ may be included, however, we found that this often results in self-intersections when the coil cross-section is defined from the path.

Further details regarding the implementation of all parameterizations are provided in the Supplementary Information.

We first specified a nominal reactor, the path of which was used for both parameterizations. Based on previous work, the standard coiled-tube reactor performs optimally with a low pitch and large coil radius1. We specified the nominal pitch (pn) as 10.4 mm and the nominal coil radius (Cr,n) as 12.5 mm. It has previously been demonstrated that as the number of coils increases, the number of equivalent tanks-in-series values rises approximately linearly9. Therefore, we specified two coil rotations, balancing the overall simulation time and reactor size, providing an overall reactor length of 4πCr,n. For the cross-sectional interpolation points, we applied lower and upper bounds of 2 and 4 mm, respectively, indicating minimum and maximum radial values. To maintain feasible coil paths, we defined the bounds of height deviations (Δz) as ±1 mm and radial deviations (Δρ) as ±3.5 mm, ensuring no self-intersections.

We applied an established methodology to evaluate the reactor designs under specific operating conditions, similar to the approach used in our previous work1. Using OpenFOAM, a software package for CFD, a simulation was performed in which an impulse tracer was injected into the water medium and the resulting concentration was tracked by solving the relevant transport equations. Key aspects of the approach include the application of the transient pimpleFOAM solver for unsteady momentum equations and scalarTransportFoam to handle convection–diffusion effects with a low diffusion coefficient constant. The integral of the concentration at the outlet was recorded at each time step and returned from the simulation. The evaluation process also involved monitoring the tracer concentration at the outlet to terminate the solution at a specific threshold, ensuring minimal solution times. This method has been adapted to the current study, enabling the evaluation of highly parameterized reactors (for more information see refs. 1,18).

The optimization of highly parameterized reactor simulations is high dimensional and involves computationally expensive function evaluations. To solve this problem, we applied multi-fidelity Bayesian optimization, specifically design and analysis of reactor and tube simulations (DARTS), which was developed for the optimization of simulated reactors and tubes1. This method enables multiple continuous simulation fidelities to be considered simultaneously. In addition, it allows the explicit modeling of the cost of a simulation as a function of simulation fidelity as well as inputs, for example, geometry parameters.

Expensive derivative-free optimization problems, such as the optimization of CFD simulations, the selection of appropriate neural network hyperparameters or the optimal design of experiments, exist throughout a number of domains. Bayesian optimization has been proposed to address these challenges, with the aim of determining optimal solutions in a tolerable computational or time budget. The majority of computationally expensive functions share a common trait: their complexity can be reduced at the expense of accuracy, resulting in the ability to perform function evaluations with the same set of decision variables with varying degrees of confidence, dictated by one or more fidelities. Fidelities control both the bias and noise of function evaluations. For example, the number of discrete cells used for the evaluation of a flow field, the reduction of which results in a less accurate output. Leveraging lower-fidelity function evaluations for the optimization of expensive systems for which there is a time or computational budget for identifying an optimal solution can improve or enable the solution of often intractable problems.

Multi-fidelity Bayesian optimization integrates function evaluations across a number of different fidelities. By applying a cost-adjusted acquisition function, both the subsequent set of decision variables as well as simulation fidelities are selected, accounting for the trade-off between information and cost. Incorporating lower-fidelity evaluations reduces the time to generate optimal solutions as fewer high-fidelity simulations have to be performed. Multi-fidelity approaches contribute to environmental savings in settings where evaluations necessitate substantial computational resources, such as large CFD simulations, reducing life-cycle emissions of overall experimentation, design and optimization procedures. Multi-fidelity Bayesian optimization has been applied to the design of reactor and tube simulations1, battery design4, hyperparameter optimization2 and the design of ice-sheet simulations32.

The DARTS framework for the design and analysis of reactor and tube simulations was proposed by Savage et al.1. The approach takes advantage of multi-fidelity simulations and integrates the optimization process with OpenFOAM within a single Python framework. At the core of the framework, simulations consisting of a set of decision variables \({{{\bf{x}}}}\in {{{\mathcal{X}}}}\) and fidelities \({{{\bf{z}}}}\in {{{\mathcal{Z}}}}\) are chosen on the basis of a cost-adjusted acquisition function to solve the following equation:

where z• is the element-wise vector of highest fidelity parameters. (In the case that the dimensionality of z is greater than 2 and the set is non-ordered, this is specified as the element-wise maximum of fidelity values.) Equation (6) presents the cost-adjusted acquisition function:

where β weights the exploration versus exploitation trade-off, and λt refers to the GP that models the computational cost. The parameters within equation (6) allow us to balance the exploration, exploitation and cost of computational experiments and allow for multiple continuous fidelity parameters to be considered. The framework also guarantees an evaluated high-fidelity solution to be returned. Further information regarding the DARTS framework can be found in ref. 1.

We defined a maximum time budget that depended approximately on the number of parameters in each parameterization optimized, ranging from 72 to 168 h in total. In the DARTS framework, two hyperparameters have to be selected, namely β, which defines the exploration level, and pc, which dictates how conservative the stopping criteria is. We selected β = 1.5 and pc = 2 based on values that were deemed to be robust in previous work concerning a similar problem.

Previous work has demonstrated the optimization and analysis of coiled-tube reactors. The RTD resulting from a simulation is approximated using a tanks-in-series model, with a larger number of theoretical tanks representing stronger plug flow performance. In these studies, RTDs were observed to be unimodal and relatively symmetrical. Owing to the extensive design space applied within this study, we found that simulations can return non-ideal distributions that contain particularly long tails or are bimodal. As more complex fluid flows are induced by the reactor designs, the resulting distributions become less ideal and the tanks-in-series model provides a poor approximation.

Therefore, in this work, we proposed a composite objective function that serves to not only maximize plug flow performance but also encourages symmetric and unimodal RTDs, accounting for the designs returning non-ideal concentration distributions.

Equation (9), with its constituent quantities defined in equations (7) and (8), demonstrates demonstrates the quantity, denoted f, that is minimized during optimization:

where d is the number of data points contained within an RTD returned from a simulation, \(\hat{E}\) and E are the predicted and returned dimensionless concentration values, respectively, θ represents dimensionless time, N represents tank-in-series values, N* is the estimated number of equivalent tanks in series, MSE represents the mean square error, and α weights the two components of the expression to have approximately equal contribution. Based on initial testing, a value of 100 was assigned to α to ensure that the tank-in-series value and the non-ideality error were weighted approximately equally (~50–100 each).

All code was evaluated using 64 CPUs, a single RTX6000 GPU to aid with training and evaluating GPs, and 64 Gb of RAM. All code used within this study can be found within the associated repository https://github.com/OptiMaL-PSE-Lab/pulsed-reactor-optimisation. Code containing only the respective parameterizations and function evaluations for benchmarking purposes can be found within the associated repository https://github.com/OptiMaL-PSE-Lab/pulsed-reactor-optimisation.

The selected designs were exported into the Standard Tessellation Language (STL) file format from their mesh representation and modified, resulting in 3D-printable models. This was achieved by incorporating a bounding cylinder to encompass each reactor. The linear sections at the inlet and outlet were specifically configured with outside diameter tube fittings of 8 and 10 mm, respectively, a requirement for interfacing with the experimental apparatus. Finally, the bounding cylinder was reduced with the internal volume between the coil removed to minimize the resin required for printing. The designs were printed using a FormLabs Form3+ 3D printer with Clear V4 resin following the manufacturer’s default settings. A post-processing phase was performed comprising a washing stage in isopropyl alcohol for 20 min, a drying period extending 24 h and a concluding post-cure treatment in a FormCure chamber maintained at a temperature of 60 °C for 30 min.

We used the same RTD method as that reported by McDonough et al.9 and applied by Savage et al.1. A 0.1 M KCl aqueous tracer solution was injected at the inlet of each reactor and the conductivity at the outlet was measured throughout the experiment. The net flow of deionized water and tracer injection were controlled using three separate OEM syringe pumps (C3000, TriContinent), which were hydraulically linked to the reactor via polytetrafluoroethylene tubing routed through a custom Swagelok piece.

The following reagents were used for the Villermaux–Dushman reaction: 0.0075 M H2SO4, 0.008 M KI, 0.0015 M KIO3, 0.023 M NaOH and 0.0225 M H3BO4. The reagents were introduced into the reactor using separate syringe pumps and the absorbance at 353 nm due to triiodide ions was measured at the reactor outlet using a UV–visible spectrophotometer and a 10 mm pathlength cuvette. Three repeats were performed for each reactor design to ensure reproducibility. An overview of the methodology used in this study is presented in Fig. 5c.

The OpenFOAM case files for all four reactors presented in Fig. 1 are available via GitHub at https://github.com/OptiMaL-PSE-Lab/pulsed-reactor-optimisation. The STL files for the reactors presented in Fig. 3 are available in the Supplementary Information. Source data are provided with this paper.

All code used within this study can be found within the associated repository https://github.com/OptiMaL-PSE-Lab/pulsed-reactor-optimisation. For use as an optimization benchmark problem, the reactor simulations are also available in the form of a Docker-based REST API with code and instructions at https://github.com/trsav/reactor_benchmark.

Savage, T., Basha, N., McDonough, J., Matar, O. K. & del Rio Chanona, E. A. Multi-fidelity data-driven design and analysis of reactor and tube simulations. Comput. Chem. Eng. 179, 108410 (2023).

Lindauer, M. et al. BOAH: a tool suite for multi-fidelity Bayesian optimization; analysis of hyperparameters. Preprint at https://doi.org/10.48550/arXiv.1908.06756 (2019).

He, X., Tuo, R. & Wu, C. Optimization of multi-fidelity computer experiments via the EQIE criterion. Technometrics 59, 58–68 (2017).

Article Google Scholar

Folch, J. P. et al. Combining multi-fidelity modelling and asynchronous batch Bayesian optimization. Comput. Chem. Eng. 172, 108194 (2023).

Article CAS Google Scholar

Takeno, S. et al. Multi-fidelity Bayesian optimization with max-value entropy search and its parallelization. Preprint at https://doi.org/10.48550/arXiv.1901.08275 (2019).

Lam, R., Allaire, D. L. & Willcox, K. E. Multifidelity optimization using statistical surrogate modeling for non-hierarchical information sources. In 56th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference (AIAA, 2015).

McLeod, M., Osborne, M. A. & Roberts, S. J. Practical Bayesian optimization for variable cost objectives. Preprint at https://doi.org/10.48550/arXiv.1703.04335 (2017).

McDonough, J., Armett, J., Law, R. & Harvey, A. Coil-in-coil reactor: augmenting plug flow performance by combining different geometric features using 3D printing. Ind. Eng. Chem. Res. 58, 21363–21371 (2019).

Article CAS Google Scholar

McDonough, J., Murta, S., Law, R. & Harvey, A. Oscillatory fluid motion unlocks plug flow operation in helical tube reactors at lower Reynolds numbers (Re ≤ 10). Chem. Eng. J. 358, 643–657 (2019).

Article CAS Google Scholar

Agrawal, S. & Nigam, K. Modelling of a coiled tubular chemical reactor. Chem. Eng. J. 84, 437–444 (2001).

Article CAS Google Scholar

Jokiel, M. et al. Helically coiled segmented flow tubular reactor for the hydroformylation of long-chain olefins in a thermomorphic multiphase system. Chem. Eng. J. 377, 120060 (2019).

Article CAS Google Scholar

Wang, L., Ni, P. & Xi, G. The effect of off-center placement of twisted tape on flow and heat transfer characteristics in a circular tube. Sci. Rep. 11, 6844 (2021).

Article CAS PubMed PubMed Central Google Scholar

Pukkella, A. K., Nadimpalli, N. R. V., Runkana, V. & Subramanian, S. A novel spiral infinity reactor for continuous hydrothermal synthesis of nanoparticles. Sci. Rep. 12, 8616 (2022).

Article CAS PubMed PubMed Central Google Scholar

Porta, R., Benaglia, M. & Puglisi, A. Flow chemistry: recent developments in the synthesis of pharmaceutical products. Org. Process Res. Dev. 20, 2–25 (2015).

Article Google Scholar

Hagedorn, J. & Kargi, F. Coiled-tube membrane bioreactor for cultivation of hybridoma cells producing monoclonal antibodies. Enzyme Microb. Technol. 12, 824–829 (1990).

Article CAS PubMed Google Scholar

Dong, Z., Zondag, S. D., Schmid, M., Wen, Z. & Noël, T. A meso-scale ultrasonic milli-reactor enables gas–liquid–solid photocatalytic reactions in flow. Chem. Eng. J. 428, 130968 (2022).

Article CAS Google Scholar

Grande, C. A. et al. Multiscale investigation of adsorption properties of novel 3D printed UTSA-16 structures. Chem. Eng. J. 402, 126166 (2020).

Article CAS Google Scholar

Basha, N., Savage, T., McDonough, J., Del Rio Chanona, E. A. & Matar, O. K. Discovery of mixing characteristics for enhancing coiled reactor performance through a Bayesian optimisation-CFD approach. Chem. Eng. J. 473, 145217 (2023).

Article CAS Google Scholar

Nivedita, N., Ligrani, P. & Papautsky, I. Dean flow dynamics in low-aspect ratio spiral microchannels. Sci. Rep. 7, 44072 (2017).

Article PubMed PubMed Central Google Scholar

Dean, W. Note on the motion of fluid in a curved pipe. Philos. Mag. J. Sci. 4, 208–223 (1927).

Article Google Scholar

Gao, H., Zhou, J., Naderi, M. M., Peng, Z. & Papautsky, I. Evolution of focused streams for viscoelastic flow in spiral microchannels. Microsyst. Nanoeng. 9, 73 (2023).

Article PubMed PubMed Central Google Scholar

Ligrani, P. M. A Study of Dean Vortex Development and Structure in a Curved Rectangular Channel with Aspect Ratio of 40 at Dean Numbers up to 430. https://ntrs.nasa.gov/citations/19950005258(NASA, 1994).

McDonough, J., Ahmed, S., Phan, A. & Harvey, A. The development of helical vortex pairs in oscillatory flows—a numerical and experimental study. Chem. Eng. Process. 143, 107588 (2019).

Article CAS Google Scholar

Pinot, J., Commenge, J.-M., Portha, J.-F. & Falk, L. New protocol of the Villermaux–Dushman reaction system to characterize micromixing effect in viscous media. Chem. Eng. Sci. 118, 94–101 (2014).

Article CAS Google Scholar

McDonough, J., Oates, M., Law, R. & Harvey, A. Micromixing in oscillatory baffled flows. Chem. Eng. J. 361, 508–518 (2019).

Article CAS Google Scholar

Singh, J., Montesinos-Castellanos, A. & Nigam, K. D. P. Thermal and hydrodynamic performance of a novel passive mixer ‘wavering coiled flow inverter’. Chem. Eng. Process. 141, 107536 (2019).

Article CAS Google Scholar

Williams, C. K. & Rasmussen, C. E. Gaussian Processes for Machine Learning Vol. 2 (MIT Press, 2006).

Padonou, E. & Roustant, O. Polar Gaussian processes for predicting on circular domains. Preprint at HAL https://hal.science/hal-01119942 (2015).

Pinder, T. & Dodd, D. GPJax: a Gaussian process framework in JAX. J. Open Source Softw. 7, 4455 (2022).

Article Google Scholar

Savage, T., Basha, N., Matar, O. & del Rio Chanona, E. A. Deep Gaussian process-based multi-fidelity Bayesian optimization for simulated chemical reactors. Preprint at https://doi.org/10.48550/arXiv.2210.17213 (2022).

Cimolai, G., Dayyani, I. & Qin, Q. Multi-objective shape optimization of large strain 3D helical structures for mechanical metamaterials. Mater. Des. 215, 110444 (2022).

Article Google Scholar

Thodoroff, P. et al. Multi-fidelity experimental design for ice-sheet simulation. In NeurIPS Workshop on Gaussian Processes, Spatiotemporal Modeling, and Decision-making Systems; https://gp-seminar-series.github.io/neurips-2022/assets/camera_ready/4.pdf (2022).

Download references

The work was supported by an EPSRC PREMIERE Programme Grant (EP/T000414/1) (N.B.) and an Imperial College London President’s Scholarship (T.S.).

Department of Chemical Engineering, Imperial College London, London, UK

Tom Savage & Ehecatl Antonio del Rio Chanona

Sargent Centre for Process Systems Engineering, Imperial College London, London, UK

Tom Savage, Nausheen Basha, Omar Matar & Ehecatl Antonio del Rio Chanona

School of Engineering, Newcastle University, Newcastle upon Tyne, UK

Jonathan McDonough & James Krassowski

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

T.S. and N.B. designed and carried out the study, conceived with the support of O.M. and E.A.d.R.C. T.S. designed and performed the parameterizations, meshing, optimization and reactor design. N.B. performed the CFD experiments. J.M. and J.K. performed the residence time and reacting flow experiments. O.M., E.A.d.R.C., N.B. and T.S. wrote the paper, and J.M. provided feedback. O.M. and E.A.d.R.C. provided project supervision and guidance for the whole project.

Correspondence to Ehecatl Antonio del Rio Chanona.

The authors declare no competing interests.

Nature Chemical Engineering thanks Milad Abolhasani and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Further information on the parameterizations, CFD meshes, simulation details, optimization progress and experimental data analysis.

STL files for all of the reactors presented in Fig. 3.

Experimental data for residence time tracer experiments.

Experimental data for the Villermaux–Dushman reaction system.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

Savage, T., Basha, N., McDonough, J. et al. Machine learning-assisted discovery of flow reactor designs. Nat Chem Eng 1, 522–531 (2024). https://doi.org/10.1038/s44286-024-00099-1

Download citation

Received: 29 September 2023

Accepted: 28 June 2024

Published: 05 August 2024

Issue Date: August 2024

DOI: https://doi.org/10.1038/s44286-024-00099-1

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

Nature Chemical Engineering (2024)