As mentioned above, a strength of the model structure of MOSAICmodeling is to apply transformations on entire models or model parts. For the state estimation discussed in Section 5.4.3, a fully discretized version of the dynamic mini plant model is required. Here, orthogonal collocation on finite elements is applied to

Figure 5.8: The dynamic model of the MES mini-plant within MOSAICmodeling.

discretize time [17]. To this end, the transformation module of MOSAICmodeling is used to automatically turn the differential algebraic model formulation into a fully algebraic one, which is displayed in Figure 5.9.

Figure 5.9: Discretized form of the dynamic mini-plant model.

Through the automatic discretization, the state variables within the models obtain two new indices: cp for the collocation position and fe for the number of the finite element. Control variables are assumed constant for each individual finite element and, hence, only receive the index, fe. For the usage in Section 5.4.3, the model is realized with 16 finite elements, each 15 min long, and a collocation order of 3. Hence, the fully discretized model amounts to 28,720 nonlinear equations. For simulation and optimization purposes, this model is exported to AMPL [41], wherein it can also be reliably solved within minutes.

Given MOSAICmodeling's modular model structure, the described model is easily adjusted for other applications. Regarding the reductive amination, the kinetics are replaced, gas solubilities adjusted, and a membrane module is added to the recycling.

#### <IP_ADDRESS> Outlook on Model Development and Collaboration

With the afore-described data model and the corresponding implementation within MOSAICmodeling, a highly versatile platform for model development, management, and exchange has been created. By now, this is a mature technology and is widely in use. More than 2,000 users have accounts on MOSAICmodeling. Beyond the work described herein, MOSAICmodeling is in use for further research projects and teaching activities [38]. Drawing from the advances in dynamic models for complex chemical processes, MOSAICmodeling's model library now features a number of new process units to support a fully pressure-driven process simulation, which is, in general, beyond the scope of the state-of-the-art commercial tools. There are, of course, further interesting research paths to pursue to further extend the data model, e.g., to setbased algebra or the support of algorithmic solution paths. Besides, the master-slave configuration of the database servers needs to be continuously updated to account for the growing user base. The next step here will be research into asynchronous updates between servers, across many locations, building on modern web technologies.

#### 5.2.2 Fluid-Dynamic Investigations of Multiphase Processes

Péter Kováts, Michael Mansour, Anurag Misra, Reddy Velagala, Gábor Janiga, Dominique Thévenin, Katharina Zähringer

#### <IP_ADDRESS> Introduction

Most processes considered in this book involve a very large number of coupled physicochemical processes. For all those taking place far above the molecular scale, the mesoscopic conditions found locally in space and time control both transport and reaction steps. For this reason, a better understanding of the hydrodynamic properties (velocity, possible fluctuations due to turbulence, but also in a more general sense temperature, pressure, concentrations, size distributions . . . ) is essential to identify the optimal process and operation conditions. Due to the complexity of the configurations considered here, combined studies involving a theoretical, numerical, and an experimental part are generally necessary. Theoretical and numerical models (the latter being denoted as CFD for Computational Fluid Dynamics) always require validation, best achieved by direct comparison with measurement data; on the other hand, experiments are complex and costly, and typically cannot provide all relevant fields. For this reason, the two following sections describe both, numerical flow simulations and detailed experimental investigations of flow and scalar fields, going hand in hand for the identification of optimal conditions.

#### <IP_ADDRESS> Numerical Flow Simulations of Reactor and Settler for the MES Process

#### Introduction Regarding Numerical Flow Simulation

As already explained, when presenting the MES process in Section 2.2, a good separation of the aqueous phase, within which the valuable catalyst is dissolved, is highly desirable and forms a crucial aspect of process design. For this purpose, the gravitydriven separation of a multi-phase microemulsion system, consisting of water, oil, and a surfactant in a temperature-controlled stirred tank reactor (STR), is investigated. The main hydrodynamic features in the STR have been first investigated in detail, as explained next. Then, stirring is stopped and the separation process starts in the batch reactor. Droplet interaction processes are incorporated in the simulation using dedicated coalescence models. By comparison with experimental data regarding separation times [71], the free parameters of the coalescence model can be optimized with regard to the specific system considered in the project. The underlying idea is that a suitable model with adapted parameters will be able to accurately capture the separation kinetics and can be used for the later design and optimization of real separators, as considered at the end of this section. Figure 5.10 shows the strategy guiding this research project.

For all simulations, a numerical modeling method, combining the Euler-Euler multi-fluid framework with the Quadrature Method of Moments (QMOM), is employed on top of the standard equations, describing the conservation of mass and momentum (Navier-Stokes equation). An energy conservation equation based on the temperature variable is added to describe heat transfer processes, since the system is only stable for a narrow temperature window. The conservation equations of fluid dynamics being standard and found in many textbooks, for instance in Marchisio and Fox [135], they are not repeated here in the interest of space. Instead, the focus is on the QMOM approach and the very important coalescence model.

In order to get suitable boundary conditions and as a preparatory step for later simulations of the reactor, the hydrodynamic fields in the Stirred Tank Reactor used for the batch separation experiments have been first investigated in detail. Here, the flow conditions are turbulent, so that a suitable turbulence model must be selected.

Figure 5.10: Strategy used during this research project.

Both the classical approach based on the Reynolds-Averaged Navier-Stokes Equations (RANS) using the standard k-ɛ model, and the advanced technique called the Large-Eddy Simulation (LES) have been compared for this purpose (Figure 5.11, left). These comparisons revealed that unsteady RANS simulations deliver in the present configuration, a sufficiently accurate description of the reactor hydrodynamics (for instance, a very good estimate of the power number, Figure 5.11, right), particularly so when taking into account the extremely long computing times required by LES, preventing any systematic studies. CFD simulations, based on unsteady RANS, can thus be used for analyzing the reactor, and to get the conditions before the batch separation experiments–starting from when the stirrer is stopped.

#### Numerical Modeling of a Separating Multiphase System

For cases involving non-homogeneous conditions in space, as found here, a comprehensive description of the separation process will, in general, be necessary to meet the accuracy requirements. For this purpose, Computational Fluid Dynamics can provide the starting point for the implementation of additional models, as explained above. While traditional CFD methods can readily deliver velocity, temperature, and concentration fields as a function of space and time for the continuous phase, the evolution of the dispersed population(s), i.e., of the microemulsion droplets, in the present case, must be described using additional approaches. The separation processes considered in this project are mainly controlled by buoyancy and take place at very low Reynolds numbers, Re. The numerical simulation of the settling process is implemented

Figure 5.11: Left: Instantaneous velocity field in the STR tank, computed by LES; Right: comparison of the power number measured experimentally by Chapple et al. [20]; and of the numerical prediction obtained by the unsteady RANS simulations, using different grid types and rotation models.

using the Euler-Euler multi-fluid approach, which is well suited for flows containing a sufficiently large volume fraction of the disperse phase(s). The open-source OpenFOAM solver "multiphaseEulerFoam" (version 2012) was chosen for this purpose since it already contains a suitable multi-fluid model and can be freely extended.

Gravitation, drag, and surface tension are the dominating driving forces in the separation process under the aforementioned conditions. In the numerical model, the surface tension was assumed to be independent of flow conditions, with surface tension values for each phase pair obtained from companion experiments [71]. Drag on the disperse phase droplets of diameter dk with respect to the continuous phase was evaluated using the following equation [215], where the subscripts c and d refer to the continuous and disperse phases, respectively, α denotes phase volume fractions, and u, the corresponding velocity:

$$\vec{F}\_{D,k} = \frac{3}{4} \rho\_c a\_d a\_c C\_D \frac{|\vec{\mathbf{u}}\_d - \vec{\mathbf{u}}\_c|(\vec{\mathbf{u}}\_d - \vec{\mathbf{u}}\_c)}{d\_k} \tag{5.7}$$

The drag coefficient, CD, is estimated using the standard correlation of Schiller and Naumann [194] as a function of the (low) Reynolds number, Re. In order to track the time-evolution of the droplet size distribution n(d,t) of the disperse phase(s), a population balance model (PBM) is coupled to CFD in order to describe a purely coalescing system, with the associated rates for birth, Bcoal, and death, Dcoal, under the consideration of a convective term with phase velocity, ~uα:

$$\frac{\partial n(d,t)}{\partial t} + \nabla \cdot (\vec{u}\_a n(d,t)) = \mathbf{S}(d,t) = B\_{coal} - D\_{coal} \tag{5.8}$$

The corresponding set of equations is highly complex and a direct numerical solution is not possible within acceptable computational times. For this reason, a practical solution is achieved by transforming the PBM with the help of the moments of the distribution, n(d,t). In the present case, the method called Quadrature Method of Moments (QMOM), being generally more stable than other alternatives, has been selected [136, 137, 143, 144]. An implementation of this algorithm is available in the opensource library, OpenQBMM Version 6 [171], which was coupled with OpenFOAM and used for all subsequent simulations.

Using QMOM, the distribution, n(d,t), is approximated by its first moments, μk:

$$\mu\_k = \int\_0^\infty d^k n(d) dd\tag{5.9}$$

In the present project, five moments have been systematically used, (μ0 to μ4). Then, a transport equation is solved to describe the evolution in time and space of these moments, in which the term on the right-hand side quantifies the importance of coalescence in the present application:

$$\frac{\partial \mu\_k}{\partial t} + \nabla \cdot (\vec{u}\mu\_k) = \bar{\mathbf{S}}\_k \tag{5.10}$$

Two different coalescence models have been implemented to compute this term, and then compared with experimental data: 1) a simple model based originally on the work of Prince and Blanch [179], later extended by Lehr et al. [122] with a single free parameter, the critical velocity ucrit, as presented by Liao and Lucas [123]; 2) a simplified version of Grimes' model [56] with two free parameters, the binary coalescence parameter, KCE, and the hindered sedimentation parameter, KSV, both non-dimensional. In both cases, the computation of the coalescence kernel involves the product of a collision frequency (within square brackets) and a coalescence efficiency (after the × -sign). This leads to:

– Buoyancy model adapted from Prince and Blanch, with urel, the relative velocity of the two droplets:

$$\beta = \left[\frac{\pi}{4} \left(d + d'\right)^2 \boldsymbol{\mu}\_{rel}\right] \times \exp\left[\frac{-\boldsymbol{\mu}\_{rel}}{\boldsymbol{\mu}\_{crit}}\right] \tag{5.11}$$

– Simplified version of the Grimes' model, involving the Péclet number, Pe:

$$\beta = \left[\frac{k\_B T}{6\mu\_c} \frac{\left(d + d'\right)^2}{dd'} \left(1 + \text{Pe}(d, d') + 4.496 \text{Pe}(d, d')^{1/3}\right)\right]$$

$$\times \exp\left[-\frac{1.046}{K\_{CE}} \frac{|\rho\_d - \rho\_c|\mu\_c \text{g}}{\Gamma^{3/2} \text{B}^{1/2}} \left(\frac{dd'}{d + d'}\right)^{9/2}\right] \tag{5.12}$$

– with the Péclet number computed by:

$$\mathrm{Pe}(d, d') = \left(\frac{3}{4\pi}\right)^{1/3} (1 - \infty\_d)^{K\_{SV}} \frac{|\rho\_d - \rho\_c| \mathbf{g}}{k\_B T} \frac{dd'}{d^{2/3} - d'^{2/3}} \tag{5.13}$$

Finally, the implemented coupling between OpenFOAM and OpenQBMM, together with the exchange of information needed for the computational procedure, is illustrated in Figure 5.12.

Figure 5.12: Working principle of the coupled numerical simulation between OpenFOAM and OpenQBMM used for all results described and discussed in this section concerning liquid-liquid separation.

Based on a series of separation experiments in a temperature-controlled glass tank (Figure 5.13, left) for a variety of process conditions [71], compared to our numerical predictions, optimal values of these parameters for the considered system have been finally identified as ucrit = 0.21 mm/s for the first model, and KCE = 0.65 and KSV = 9.44 for the second one. To identify the corresponding values, CFD-based optimization has been used–a technique widely employed in our group [206]. While varying freely the model parameters within a user-specified range using a Genetic Algorithm, the target function (to be minimized) is the mean difference between the experimentally measured separation times and the ones obtained numerically, tracking for this purpose, the contact surface between the aqueous and the organic phase in time (Figure 5.13, right). The optimal values determined, thanks to this procedure, are then employed for further computations.

#### Separation in the Real Settler

After having determined the optimal model parameters for the coalescence model and validated the coupled numerical approach (combining the two solvers, OpenFOAM and OpenQBMM) by direct comparisons for the batch settler experiments (Figure 5.14), the developed procedure can finally be applied to the target configuration–the continuous settler used for the mini plant experiments with the MES process.

Figure 5.13: Left: Experimental set-up consisting of a temperature-controlled glass tank used in the investigation of a two-phase batch settling process [71], measuring the initial droplet-size distribution with the help of an endoscope. Right: Numerical prediction of the time-dependent separation process between organic and aqueous phase using the simplified Grimes model [56] for four different process conditions. The results shown here correspond to the finally identified, optimal model parameters.

Figure 5.14: Direct comparison between the separation experiments in the batch settler (left) and own numerical simulations using the procedure described previously (right) at a given time during the process.

All the necessary details concerning the mini plant and its operation can be found in the literature [75–77] and in Section 4.2 of this book. For the simulations of the separation process, one strong but necessary hypothesis is that the optimal model parameters derived from the preliminary experiments in the batch settler can be kept identical for the much more complex, real settler geometry used for continuous phase separation and catalysis recycling after the reaction. One additional challenge is that the employed settler has been changed several times during the research project, leading to repeated simulations involving long computational times. Unfortunately, one single simulation of the full-scale three-dimensional settler for real process conditions takes at least days–up to weeks, even using parallelization on Linux clusters, depending on the physical time considered on the particular geometry and the selected coalescence model. This is particularly true for the more detailed Grimes model, which requires additional computational power, compared to the simple buoyancy model. However, the model derived in this work, based on the original formulation by Grimes [56], is the only one that should be able to take into account the effects on separation due to surface tension as well as inhomogeneous distributions of temperature and surfactant concentration.

Figure 5.15: Top left: geometry of the settler. Top right: velocity field. Middle left: temperature field. Middle right: volume fraction of the bicontinuous phase, establishing itself between the organic and aqueous phases. Bottom: evolution of the volume fraction of organic phase at the bottom outlet of the settler (blue line), compared to the experimentally observed range (horizontal corridor delimited by thick dashed lines). Details can be found in Misra et al. [155].

First simulations considered the real geometry of the settler (as used at that time in the mini plant experiment) but simplified using a two-dimensional approximation [155]. The corresponding vertical cuts showing the different fields within the settler are shown in Figure 5.15.

As seen in this figure, the full Grimes model (red line) did not describe well the separation at that time, which was not due to the model itself, but to its numerical stiffness in the practical implementation. After checking this point further, the simplified version of this model, as described previously, has been implemented in OpenFOAM [154].

The later employed full-scale settler involving three outlets (one for each of the three phases; organic, aqueous, and bicontinuous) was finally simulated (Figure 5.16) in order to 1) better understand the separation process, and 2) identify bottlenecks, paving the way for a possible performance improvement in future studies. This investigation first revealed somewhat non-homogeneous velocity, temperature, and concentration fields, due to recirculation zones interacting with the temperature control through the external walls, impacting the separation. The porous zone existing in the first half of the settler noticeably impacts coalescence, as expected. The predicted time needed for a proper separation strongly depends on the employed coalescence model, but the final outcome of the simulation agrees well with experimental observations (as already shown in Figure 5.15). This reveals that, provided sufficient computational power is available, a detailed prediction of separation in complex systems can be achieved using CFD extended by suitable models, after properly calibrating the unknown model parameters, thanks to dedicated experiments.

A detailed analysis of the results obtained in the real settler revealed that the computation of the Péclet number appearing in the collision frequency for the Grimes' model can be numerically very stiff, depending on the local conditions; thus, often inducing numerical instabilities. For this reason, our current research concentrates on deriving a mixed model combining the collision frequency from the buoyancy-induced coalescence model from Prince and Blanch, together with a slightly modified expression for the coalescence efficiency, derived from the model of Grimes, still involving the same two model parameters, KCE and KSV. For this new model, the optimization procedure described previously must be repeated to identify suitable values for the parameters.

Figure 5.16: Volume fraction of the aqueous phase in the settler of the mini plant with three outlets, as obtained after 113.4 s of physical time with the simplified Grimes model, used in the coupled numerical procedure.

#### <IP_ADDRESS> Fluid-Dynamic Investigation of Gas-Liquid-Liquid Continuous Helical Flow Reactors

For the hydroformylation of long-chain olefins in a thermomorphic multiphase system, Kaiser et al. [90], suggested a tandem reactor system consisting of a helically coiled tubular reactor (HCTR), connected to a continuous stirred tank reactor (CSTR) (Section 5.3.1). It was derived by an extensive model-based reactor synthesis and dimensioning study, for an already existing mini plant.

A hydrodynamic study of this reactor system, which is important for reactor modeling, dimensioning, and optimization, was necessary for the helically coiled tube reactor, as a special form of continuous flow reactors, since the performances of these reactors concerning mixing, heat and mass transfer are less known, than those of CSTR and plug-flow reactor (PFR).

For continuous processes, the most often used chemical reactors are stirred tank reactors, tubular reactors, and variations of these types. Due to their specific characteristics, they feature different mixing, especially back mixing, and heat and mass transfer properties. These properties are very important since they strongly influence the performance of the reactor, namely, the conversion and selectivity behavior in dependence on the underlying reaction kinetics. Moreover, if the reactants form a multi-phase flow, the operation of these reactors is completely different from a one-phasic operation, and the reactor geometry must be changed to enhance the performances.

It has been known for a long time that the axial back mixing and mass transfer behavior of tubular reactors can be enhanced by using coiled tubes (Figure 5.17, left), instead of straight ones [209]. Due to coiling, the centrifugal forces influence the flow, and a strong secondary flow forms in a radial direction, leading to Dean vortices (Figure 5.18, left). The development of these vortices enhances the radial mixing while keeping a low axial back mixing behavior. This increases heat and mass transfer, and leads to narrower residence time distributions [103, 209]. Other advantages of the coiled configuration are its higher compactness and relatively large surface area-to-volume ratio. Additionally, the laminar flow regime in coiled tubes is observed at larger Reynolds numbers, compared to straight tubes.

#### Helix Geometry

Although numerous studies were previously conducted to investigate the influence of different parameters on the advantages of helical pipes concerning heat and mass transfer, e.g. [51, 102, 159, 195], the characterization of pure mixing of miscible liquids in such apparatus, which is the basis of the aforementioned transfer processes, has rarely been undertaken. The group of Nigam [114, 198] experimentally and numerically investigated the influence of the Reynolds number, Schmidt number, and the curvature ratio on the mixing performances of different coiled geometries by determining the residence time. Saxena and Nigam [190] introduced a structured configuration

Figure 5.17: Coiled configurations: horizontal coil (left), Coiled Flow Inverter (center), Coiled Flow Reverser (right, European patent No. 19176360.6-1019, EP3741453 A1).

called the coiled flow inverter (CFI, Figure 5.17, center), in which the flow direction is inverted by bending the coils. In this way, the direction of the secondary flow is changed in order to improve mixing. The CFI can be simply constructed by fitting 90° bends at regular length intervals between the coils and it attracted, also due to its geometrical compactness, some more interest for industrial applications. In most CFI applications, the flow is usually inverted after at least 3 to 4 turns, e.g. [115, 160, 190], to ensure fully developed vortices, before changing the flow direction. Consequently, a considerable increase in pressure drop is unavoidable in most cases. Nevertheless, it was recently shown that mixing can even be improved by early flow inversion in the entrance region [98, 129]. This flow inversion can be provided by the use of a short CFI, which reverts the flow by 90° or a novel geometry, called coiled flow reverser (CFR, Figure 5.17, right), which leads to a 180° flow reversion [133].

The flow conditions leading to the apparition of secondary flow and symmetric vortices can be described by the Dean number, De = Re · ffiffiffi δ p , which shows a dependency on the curvature ratio, δ = d/D. A further parameter describing the geometry of a coil is its pitch, which is often represented as dimensionless pitch, γ = P/(πD).

The influences of the helix geometry and the flow properties on mixing, heat and mass transfer were studied here on model geometries using Computational Fluid Dynamics (CFD) and locally and temporally resolved optical measurement techniques. Compared to the mini plant HCTR derived in Section 5.3.1, the model geometries had higher curvature ratios and bigger dimensionless pitches, due to the much smaller coil. For safety reasons and the fact that the TMS is two-phasic at temperatures < 85 °C, water was used in the model geometries instead of the TMS. However, as the outcomes of these studies confirm, the influence of these fluid parameters and geometrical differences between the mini plant HCTR and the smaller models on the examined processes are much less pronounced than the impact of the Reynolds number. Therefore, the results from the model geometries can be transferred to the mini plant HCTR as long as the same Reynolds number is respected. It ranges, in the studies presented hereafter, from Re = 40 to Re = 1,000; thus, a laminar flow regime is employed.

#### Mixing of the Liquid Phase

Using a helical pipe for mixing liquids may show noticeable advantages, compared to using a stirrer, since this is a robust solution (no moving parts), needs no additional

Figure 5.18: Left: Secondary flow structures (Dean vortices) for Re = 200 at the outlet of a 5-turn straight helical pipe. Calculations from [128]. Right: Dean vortices detected by 3D Q-criterion inside the helically coiled tube at Re = 220 and De = 89, as obtained from time-averaged measurements (purple isosurfaces) and simulation (red isosurfaces). Adopted from [104].

power source, and has a compact structure. Therefore, lower maintenance and energy consumption are usually required, compared to active mixers. Additionally, mixing in helical pipes can also take place efficiently in the laminar regime, avoiding energy losses associated with the generation of turbulence.

Liquid-liquid mixing was extensively examined here in different numerical and experimental configurations. The mini plant HCTR geometry (Section 5.3.1) was analyzed numerically and using the fluid properties of the TMS [86, 87]. In all configurations, the CFD code, Star-CCM +, was used for the calculations [128, 130], and planar laser-induced fluorescence (PLIF) and particle image velocimetry (PIV) were used to analyze the mixing and flow fields experimentally [103–105]. These studies showed a very good mixing behavior of the coiled reactors in the first 2–3 coils. Then, the established flow pattern with the secondary Dean vortices, which were visualized and quantified in three dimensions through tomographic PIV measurements and CFD calculations of the velocity fields inside the helical reactors (Figure 5.18, right), leads to only a slightly further mixing. This can be recognized in Figure 5.19 in the first row, where the mass fractions of two different liquids are represented in the cross-sections of a 6-coil straight helix reactor.

By varying the Reynolds number and looking at the outlet mixing efficiency for a variety of geometrical parameters (pipe diameter, d, coil diameter, D, and pitch, p), two maxima for the mixing efficiency were identified at about Re ≈ 35 and Re ≈ 650 (Figure 5.20) [130, 131]. Mixing shows a non-monotonic behavior with the change of Reynolds number, since the residence time, the vortex structure, and the vortex strength change simultaneously with Re. Here, it should be noted that the pressure drop in a helically coiled tube is generally higher, compared to a straight pipe of the

Figure 5.19: Mass fractions of two miscible liquids in cross-sections of different helically coiled reactor configurations at Re = 500. White and red regions correspond to the initial liquids; green regions show nearly complete mixing [105].

same length [78]. Accordingly, the first optimal Reynolds number should normally be preferred since the pressure drop is significantly lower, compared to the second optimal value. Only if the Reynolds number is smaller than 30 or bigger than 1,000, geometry has a major impact on the mixing efficiency. This shows that the transfer of the results obtained in the smaller model geometries to the bigger HCTR-geometry is possible since its Reynolds number (Re ≈ 100) lies inside these limits.

In order to identify the most efficient helical geometry with the best mixing properties and lowest pressure drop, multi-objective numerical optimization for the flow of two miscible liquids in helical pipes was used in the aforementioned first optimal Reynolds number range [134], which also represented the range of Dean numbers in the mini plant HTCR (Section <IP_ADDRESS>). The CFD simulations have been coupled with the optimization code, OPAL++, leading to an automatic evaluation loop, including geometry modeling, mesh generation, and post-processing. As a whole, 30 generations encompassing 1,226 CFD simulations have been completed during the optimization process. The resulting Pareto front involves 64 different individuals, in which good mixing and small pressure drop are found simultaneously. The first and last individuals of the Pareto front are represented together with the global optimum in Figure 5.21. Very low pipe diameters have been rejected in the

Figure 5.20: Outlet mixing coefficients of various helical geometries (G1–G27) as a function of the Reynolds number. For geometrical details, see Mansour et al. [130]. Two optimal regions for mixing can be identified at around Re ≈ 35 and Re ≈ 650.

optimization since they lead to unacceptably high-pressure losses. A low-to-moderate curvature ratio was found optimal. All optimal solutions correspond to low dimensionless pitch, torsion, and torsion number.

The results of all these numerical and experimental investigations showed the establishment of a stable vortex pair in the helical coils after 2–3 turns. Therefore, for heat and mass homogenization purposes, it is advantageous to disrupt this established flow field by early flow redirection, as it is done in coiled flow inverters, which have 90° bends after several turns (Figure 5.17, center). The effect of this bending can be even reinforced through a further redirection of the flow to 180°. These coiled flow reversers (CFR, Figure 5.17, right, European patent No. 19176360.6-1019, EP3741453 A1), were compared to straight helical coils (CF) and CFI, with the same diameters and lengths, in several numerical and experimental studies [98, 105, 133] concerning mass and heat transfer, respectively. These show the clear advantage of flow redirection for homogenization. CFI and CFR attain nearly complete mixing after only a few helical turns if the bending takes place after each first (rows 5–6 in Figure 5.19) or second (rows 2–3 in Figure 5.19) turn. The CFR performs even better than the CFI in the region of the second mixing optimum at Re ≈ 650 [105].

Figure 5.21: Results of numerical CFD helix geometry optimization: first and last individual of Pareto front, together with the global optimum. Their streamlines, vorticity, mass fractions, mixing coefficients, and pressure drop are also represented [134].

#### Heat Transfer

Due to the increased heat transfer coefficient, compared to straight pipes, the thermal properties and heat transfer in helical pipes have been extensively considered in early studies, e.g. [84, 112, 187, 196]. The effect of several process parameters was considered for different boundary conditions, such as constant wall flux, constant wall temperature, or convective heat transfer, e.g. [79, 178]. As a result, various correlations were proposed concerning Nusselt number, critical Reynolds number, and associated friction factor, based on experimental data or analytical solutions. However, no generally valid correlation has been developed covering all relevant geometrical parameters, boundary conditions, or flow characteristics. In addition, different correlations are sometimes inconsistent with each other, even for the same conditions. Therefore, additional studies are still useful to shed light on such issues. For practical applications, it is important to obtain a uniform fluid temperature as fast as possible after heat addition, to avoid inhomogeneity. Consequently, in the investigation presented hereafter, the homogenization performances were studied in terms of the coil length needed for reaching a homogenous temperature profile after heat addition.

This has been done numerically in two different studies [129, 133]. First, the thermal homogenization in a horizontally-coiled helical pipe was compared when heating either the entrance region or other locations within the fully-developed flow. Reynolds numbers between 2,000 and 16,000 were considered in a water flow. The heat was added to the flow by a constant-temperature wall, heating the first, third, or the fifth turn of a ten-turn coil individually. The average Nusselt number of the first-turn case was found to be noticeably higher than for the other two cases, increasingly so at higher values of the Reynolds number. Also, temperature homogenization is reached much earlier in that case, and the required homogenization length is almost independent of the Reynolds number. These results show, once more, the importance of the developing flow pattern for homogenization. Therefore, flow redirection, like in CFI or CFR, is advantageous for heat transfer as well.

For this reason, in the second numerical study, these reactors have been included and the heating of, initially, 27 °C warm water by 50 °C hotter reactor walls was examined. For low Reynolds numbers, the Nusselt numbers of the three reactors (straight coil, CFI, CFR) are comparable. However, for Re > 1,000, the Nusselt number of the CFR always shows the highest values. This can also be observed on the outlet temperature, shown in Figure 5.22; the values are constant and equal to the wall temperature for all geometries at low Reynolds numbers; in this case, perfect heat transfer was obtained, thanks to the long residence time. For higher Re, the CFI and the CFR can obviously lead to better heat transfer and higher averaged outlet temperatures. In all cases, the CFR shows the best heat transfer and thermal homogenization, which is reached much earlier, and the required homogenization length is almost independent of the Reynolds number.

Considering all these points, it is recommended to always select the first coil for controlling the temperature in a horizontally-coiled helical pipe. As a further consequence, reactors that are redirected after only a few coils should be preferred, as in short CFI and CFR.

#### Gas-Liquid Mass Transfer and Flow Regimes in Gas-Liquid Flows

Mass transfer has been considered in helically coiled tubes before, but most of these studies were executed by probe measurements, influencing the flow pattern inside the helix. Also, it should be mentioned that all these studies were performed in vertically orientated coils, which is different from the configurations used here.

The gas-liquid mass transfer of an annular flow has been studied in literature [7], with the determination of kLa-values. Similar measurements in the slug-flow regime were carried out in helically coiled tubes by Kulic and Rhodes [113]. The mass transfer in helically coiled tubes and straight horizontal tubes has been compared by Jespen [83]. In all these papers, kLa-values were found between 0.02 s−1 and 0.5 s−1 , depending on the flow regime employed. The comparison with straight tubes indicated a

Figure 5.22: Surface-averaged outlet temperature of 3-turn coiled reactor configurations heated by a 77 °C hot wall at different Reynolds numbers [129].

higher mass transfer rate in the coiled tubes. The increased mass transfer rate found in coiled tubes was also confirmed by other studies [1, 60, 197] for liquid-liquid, solidliquid, and gas-liquid situations.

In this study, the gas-liquid mass transfer in helical coils was examined experimentally and numerically in several model geometries [86, 103], since it is crucial for the reaction progress of the reactions considered within InPROMPT. Equivalent gas concentrations inside the liquid and also volumetric mass transfer coefficients were determined by laser-induced fluorescence of a tracer redox-reaction. The measurements showed that the saturation concentration was reached almost within the first 3–4 turns of the helix. The increase of the gas hold-up leads to higher mass transfer rates, whereas the liquid superficial velocity has no significant influence on the oxygen concentration in the liquid phase. The volumetric mass transfer coefficients were calculated from the measured oxygen concentrations and in the first turns of the helix, very high values of up to kLa =3s−1 were obtained. When approaching saturation, these values decreased to around 0.2 s−1 for all examined conditions, which correspond to those of straight tubes. These findings underline the high impact of the Dean vortices on mass transfer processes. This effect is strongest during the very first turns of the helix.

The numerical and experimental analysis of the two-phase flow regimes in helical coils [132, 161] further showed that the velocity profiles and vortex structures are much more complex under two-phase flow conditions (Figure 5.23). The flow fields in the Taylor plug-flow regime show an additional vortex pair, counter-rotating in front of the bubbles. Compared to the Dean vortices, these are shifted to the outer side of the horizontal helix. With higher flow velocities, those secondary vortices extend deeper into the liquid slug and produce further tertiary vortices, rotating in the opposite direction. With small bubbles, no vortices, except the dominant Dean vortices, appear in the liquid slug. Therefore, mixing is favored by small liquid slug lengths with Taylor bubbles, where the vortices generated by the bubbles dominate the flow behavior inside the slug.

Figure 5.23: Mixing pattern and vortex structures (Q-criterion) in a liquid slug of a horizontal helically coiled gas-liquid reactor (unrolled visualization).

#### 5.2.3 Surrogate Models for Thermodynamic Equilibria of Gas-Liquid and Liquid-Liquid Systems

Corina Nentwich, Stefanie Kaiser, Sebastian Engell

The first principle models are the methods of choice to describe chemical processes as they provide reliable results over a wide range of process conditions. In process synthesis, the thermodynamic models that describe the phase equilibria and compute the solubilities of the components in reaction systems are important building blocks of the flowsheets.

Equations of state, e.g., PC-SAFT, are powerful methods to predict the phase equilibria of complex systems (Section 3.1.1), but they are computationally expensive because iterative solutions of the density root problem and the phase equilibrium condition must be performed. The thermodynamic equilibrium conditions must be computed in each iteration of the solver, repeatedly. This leads to a very large number of equilibrium calculations and, therefore, such complex equations of state are difficult to use directly in optimization [163]. For example, in the case of the optimization of the layout and the operating parameters of distillation columns, the number of such calculations is very large because the thermodynamic equilibrium has to be computed for each stage individually, leading to either a large number of calls to an external system of equations or a large increase in the number of equations that have to be solved simultaneously. A solution for this issue is the application of surrogate models. Surrogate models are black box models that can be evaluated fast and can replace more complex models. Examples of surrogate models are neural nets (ANNs), Kriging models, support vector machines (SVM), or other generic approximation functions.

Surrogate models are fitted by optimization of the free parameters of the surrogate models at a set of data points (sampling points) where the output values are obtained from querying the original model or from experiments. The selection of the sampling points has a strong influence on the quality of the final surrogate model.

A common sampling approach is space-filling or exploratory sampling, which covers the input space more or less equidistantly. Here, the location of all sampling points can be determined before any knowledge about the function has been generated. However, as all regions of the function are covered equally, no special emphasis is put on covering regions with more complex structures and, therefore, potentially larger errors of the surrogate model. An example that illustrates this problem is shown in Figure 5.24.

Figure 5.24: Example of a space-filling design.

Figure 5.25: LHS design.

Examples of space-filling designs are full-factorial or partial factorial designs. The input space is divided into discrete levels and samples are taken at all combinations of these levels. Another variant that leads to fewer sampling points is the Latinhypercube design (LHS). As for the factorial design, the input space is discretized, but the samples are not located at all combinations of all levels; only one sample is placed on each level and the distance between the sampling points is maximized when distributing them. An example of the LHS design for a 2-dimensional space is shown in Figure 5.25. Space-filling sampling designs are widely used, but not efficient for complex structures or if discontinuities in specific regions exist. Spacefilling sampling is an example of an explorative strategy as the entire region is equally covered. In contrast, exploitative strategies cover those parts of the input region more intensely that contribute most to the objective function value, i.e., where larger errors are expected or observed. Adaptive or sequential sampling seeks a compromise between exploration and exploitation and is, therefore, better suited for cases where large variations of the output are observed in some regions, whereas in other regions, the outputs are relatively smooth. In sequential sampling, additional points are added to the set of sampling points, iteratively, to improve the fit of the surrogate model. Sequential sampling approaches have been studied intensely in recent research [49, 85, 163, 165].

Figure 5.26: Graphical scheme of the sampling design update of the mixed adaptive sequential sampling [165].

Sequential designs originally mostly aimed at meeting the exploitative objective. Crombecq et al. [26] developed an algorithm, placing additional samples to explore regions with only a few points, and in the approach of Cozad et al. [25], the additional points are chosen such that the error between the true function value and the surrogate model prediction is maximized. This approach leads to a large number of expensive function evaluations, which can be overcome by the mixed adaptive sampling that was proposed by Eason and Cremaschi [33], and further developed by Nentwich and Engell [163], and Nentwich et al. [165]. This approach seeks a trade-off between a space-filling objective and exploitation. The principle is depicted in Figure 5.26. An initial current sampling design, e.g., a LHS design of small size, is divided into NSS subsets. These subsets are used to train surrogate models of the same structure: model ^yð Þ 0 is based on the complete dataset and each model ^yð Þ − m is based on all subsets except for the subset m, which is left out. The next step is to propose a set of new sample candidates, for which only the input variables are proposed without calculating the output. This can be performed using a space-filling sampling design. Without knowing the actual output of the true function at the proposed input locations, the candidates can be assessed by calculating the distance of each candidate to all samples already present in the current design, which represents the explorative sampling criterion, and by calculating the prediction variance of all beforehand trained subset models, representing the exploitative sampling criterion. The minimum Euclidean or Manhattan distance dj between the new candidate, j, and the N already existing points can be used as the distance measure. They are calculated as

$$d\_f(\text{Euclidean}) = \min\_{n \in \{1, \dots, N\}} \left| \left| \mathbf{x}\_f - \mathbf{x}\_n \right| \right|^2,\tag{5.14}$$

$$d\_{\boldsymbol{\delta}}(\boldsymbol{\Lambda}\boldsymbol{\alpha}\boldsymbol{\Lambda}\boldsymbol{\epsilon}\boldsymbol{\Lambda}\boldsymbol{\epsilon}\boldsymbol{\epsilon}\boldsymbol{\epsilon}\boldsymbol{\epsilon}\boldsymbol{\epsilon}\boldsymbol{\epsilon}\boldsymbol{\epsilon}) = \min\_{n \in \{1, \ldots, N\}} \left| \left| \boldsymbol{\chi}\_{\boldsymbol{\delta}} - \boldsymbol{\chi}\_{n} \right| \right|.\tag{5.15}$$

The prediction variance can be calculated by the Jackknife variance as proposed by Eason and Cremaschi [33]; alternatively, bootstrapping can be used as shown by Nentwich et al. [165]:

$$\mathrm{s}\_{f}^{2}(\mathrm{Jackkraft}) = \frac{1}{\mathrm{NSS}(\mathrm{NSS}-1)} \cdot \sum\_{\ell=1}^{\mathrm{NSS}} \left(\tilde{\boldsymbol{\mathcal{Y}}}\_{\boldsymbol{\mathcal{Y}}} - \tilde{\boldsymbol{\mathcal{Y}}}\_{\boldsymbol{\mathcal{Y}}}\right)^{2},\tag{5.16}$$

$$\text{cs}\_{f}^{2}(\text{Boosttrap}) = \frac{1}{B-1} \sum\_{b=1}^{B} (\tilde{\mathbf{y}}\_{f}^{X\_{b}} - \overline{\tilde{\mathbf{y}}}\_{f})^{2},\tag{5.17}$$

where ~yij is the weighted average of the jackknife pseudo values (eq. (5.19)) of each candidate, j, predicted by the subset model, i, in the case of the Jackknife prediction variance. For the bootstrap prediction variance, random samples with replacement are drawn from the existing sampling design, resulting in a new sampling design Xb. Based on this design, a surrogate model ~y Xb b is trained and the procedure is repeated B times. For the Jackknife prediction variance, Nentwich and Engell [163] proposed an alternative formulation, which is based on the number of samples Ni in subset i:

$$\mathbf{s}\_{\mathbf{j}}^{2} = \sum\_{l=1}^{\text{MSS}} \frac{1}{N(N - N\_{l})} \left(\tilde{\mathbf{y}}\_{\text{jl}} - \overline{\tilde{\mathbf{y}}}\_{\text{l}}\right)^{2} \tag{5.18}$$

with

$$\tilde{\tilde{\mathcal{y}}}\_{\tilde{l}} = \sum\_{l=1}^{\text{NSS}} \frac{N - N\_l}{\sum\_{l=1}^{\text{NSS}} N - N\_l} \cdot \tilde{\mathcal{y}}\_{l\tilde{l}} \tag{5.19}$$

and

$$
\tilde{\mathbf{y}}\_{\hat{\mathbf{y}}} = \mathbf{N} \hat{\mathbf{y}}^{(0)} - (\mathbf{N} - \mathbf{N}\_{\hat{\mathbf{t}}}) \hat{\mathbf{y}}^{(-1)}.\tag{5.20}
$$

The more samples that are included in subset, i, the larger is the weight assigned to the subset model that is based on this subset. Both the exploitative and the explorative criterion are combined into one sampling criterion η:

$$\eta\_{\circ} = \frac{d\_{\text{f}}}{\max\_{\text{f}} d\_{\text{f}}} + \frac{s\_{\text{f}}^2}{\max\_{\text{f}} s\_{\text{f}}^2} \cdot \tag{5.21}$$

The candidates with the highest η-values are selected to be added to the sampling design. Only for those candidates, the true function is evaluated and the samples are added to the current design. In regions of the input space where currently no samples are located, the ηi -value might be very high. This may cause specific regions with significantly higher ηi -values than others, as shown for a 2-dimensional example in Figure 5.27. To avoid the selection of candidates that are close to each other, the peaks' candidate selection method developed in the literature [165] can be applied. All candidates are located in a grid and for neighboring points in the grid, the values ηi of the candidate selection criteria are compared. Candidates with all neighboring ηi -values being smaller are called local peaks (see cross and asterisk markers in Figure 5.27). Only local peaks are chosen to be added in this iteration.

Figure 5.27: Illustrative 2D example for the peaks' candidate selection method [165].

The number of selected points cadd among the N candidates is determined by the selection factor SF:

$$\mathcal{L}\_{add} = \mathcal{S} \mathbf{F} \cdot \mathbf{N}.\tag{5.22}$$

Thereby, the number of samples that are added increases in each iteration. To define a termination criterion, the cross-validation error (CVE) is used. Nentwich et al. [165] define the CVE by the squared error for each current sample, i, calculated as the deviation of the true value yi and the prediction of the subset model ysubset i , where the current sample i is missing:

$$\text{CVE}\_l(\mathbf{y}) = \left(\mathbf{y}\_l - \mathbf{y}\_l^{\text{subsetset}}\right)^2. \tag{5.23}$$

As proposed by Eason and Cremaschi [33], the slope percentage criterion is used as a termination criterion:

$$\mathfrak{a}\_{\mathbf{k}} = \left| \frac{\text{CVE slope from iteration } \mathbf{k} - \mathbf{1} \text{ to } \mathbf{k}}{\max \text{CVE slope from iteration } \mathbf{l} \le \mathbf{k}} \right| \le \varepsilon. \tag{5.24}$$

The termination criterion is based on the cross-validation error in order to use all available information that was gained during the sampling procedure. To compare the different trained models, the use of measures with respect to a large test set is commonly used as the standard measure for comparable results. For regression models, the mean absolute error for the quantity i (MAEi) or the mean squared error for the quantity i (MSEi) can be used:

$$MAE\_{l} = \frac{1}{N\_{\text{test}}} \sum\_{j=1}^{N\_{\text{test}}} |\hat{\boldsymbol{\wp}}\_{l}(\mathbf{x}\_{l}) - \boldsymbol{\wp}\_{l}(\mathbf{x}\_{l})|,\tag{5.25}$$

$$MSE\_{\ell} = \frac{1}{N\_{\text{test}}} \sum\_{j=1}^{N\_{\text{test}}} \left( \hat{\nu}\_{\ell}(\mathbf{x}\_{\ell}) - \nu\_{\ell}(\mathbf{x}\_{\ell}) \right)^{2},\tag{5.26}$$

where ^yi xj is the predicted value of quantity i at a sample location xj and yi xj is the true value of quantity i at a sample location xj for j = 1, ... ,Ntest test set samples. For classifiers, the percentage of misclassified test set samples (MisC) can be used:

$$\text{MisC} = \frac{\text{100\%}}{2\text{N}\_{\text{test}}} \sum\_{j=1}^{N\_{\text{test}}} |\hat{\mathbf{y}}\_{\text{class}}(\mathbf{x}\_{\text{j}}) - \mathbf{y}\_{\text{class}}(\mathbf{x}\_{\text{j}})|,\tag{5.27}$$

where ^yclass xj is the classifier output at a sample location xj and yclass xj is the true class at a sample location xj for j = 1, ... , Ntest test set samples.

The choice of the inputs and outputs of the surrogate models is important and has to be made first. In the following paragraphs, we will present the surrogate-modeling of the liquid-liquid equilibrium of the ternary system n-decane, dimethylformamide (DMF), and 1-dodecene. By modeling this system the behavior in the decanter of the hydroformylation process of 1-dodecene can be predicted during optimization. The flowsheet of the process can be found in Figure 5.32. For the surrogate-modeling of thermodynamic models in process simulation and optimization, two different options exist. In the first option, called the direct surrogate-modeling approach, [163], the phase equilibrium condition is solved during the sampling and the composition of the two phases is predicted using surrogate models. A classifier is trained to identify the biphasic region. Inside this region, a regression model is trained to model the compositions of the two phases. Because of the inclusion of a classifier, a modified mixed adaptive sampling strategy is applied. The term accounting for the exploitation criterion is changed to

$$\frac{\mathbf{s}\_{j}^{2}}{\max\_{j} \mathbf{s}\_{j}^{2}} = \frac{1}{2} \frac{\mathbf{s}\_{\text{class},j}^{2}}{\max\_{j} \mathbf{s}\_{\text{class},j}^{2}} + \frac{1}{2} \frac{\mathbf{s}\_{\text{appr},j}^{2}}{\max\_{j} \mathbf{s}\_{\text{appr},j}^{2}},\tag{5.28}$$

and the slope ratio in the termination criterion is a mixed criterion, calculated by the classifier and regression models, called the total slope criterion, αtot:

$$\alpha\_{tot} = \frac{1}{2} \frac{a\_{class}}{\max(a\_{class})} + \frac{1}{2} \frac{a\_{appr}}{\max(a\_{appr})} \,. \tag{5.29}$$

In the second option called the "indirect surrogate-modeling approach" that was presented by Nentwich et al. [165], the phase equilibrium condition is solved during the process simulation or optimization, and the fugacity coefficients that enter into this condition are approximated by surrogate models. In this approach, the phase equilibrium problem is explicitly formulated and solved simultaneously with all balance equations of the process model. Therefore, no classifier is needed to predict the biphasic region.

For the direct prediction of the composition of the two phases, the number of outputs can be reduced by the introduction of the distribution coefficients Ki instead of modeling the molar fractions xj, i of the component i in phase j. The distribution coefficient is defined as

$$K\_{\dot{l}} = \frac{\dot{\boldsymbol{n}}\_{\boldsymbol{l},\boldsymbol{A}}}{\dot{\boldsymbol{n}}\_{\boldsymbol{l},\boldsymbol{\dot{m}}}} = \frac{\boldsymbol{\chi}\_{\boldsymbol{l},\boldsymbol{A}} \cdot \dot{\boldsymbol{n}}\_{\boldsymbol{A}}}{\boldsymbol{\chi}\_{\boldsymbol{l},\boldsymbol{\dot{m}}} \cdot \dot{\boldsymbol{n}}\_{\boldsymbol{\dot{m}}}},\tag{5.30}$$

where n\_i, A is the molar flow of the components i in phase A and n\_ i, in is the molar flow of the component i in the feed, respectively. n\_ A and n\_ in are the total molar flows of phase A and of the feed. Since the formulation of the distribution coefficient is only valid within the miscibility gap, a classifier that assigns a point to the biphasic region is used. If the classifier predicts a point to be in the biphasic region, a continuous regression model is used.

In the literature [163], a Support Vector Machine (SVM) is used for the classification problem, while ordinary Kriging models are chosen as regression models. To show the advantage of using mixed adaptive sequential sampling instead of the commonly applied purely space-spilling design LHS, the performance of the trained models is compared to models of the same structure that were trained using LHS designs of the same sample size.

The performance of the mixed adaptive sampling and the space-filling LHS design are compared in Figure 5.28. For all sample sizes, the classifier trained with the mixed adaptive sampling led to fewer misclassifications. Especially at medium sample sizes, the adaptive approach outperformed LHS significantly.

Figure 5.28: MisC over the classifier sample size of the sampling algorithm for NSS = 18 and SF = 0.4 and same-sized LHS design-based models showing the mean (markers), minimum, and maximum (bars) of each 5 runs [163].

The performance of the Kriging models that predict the distribution coefficients is depicted in Figure 5.29. It can be seen that the performance of the models using mixed adaptive sampling is better for the same number of sampling points. The difference of the mean absolute error is large for medium-sized sample sets of around 500 samples. For all components, a much smaller sample size is needed for the mixed adaptive sampling, in comparison to a LHS sampling, to reach a similar performance.

In the second option, the surrogate models were trained to predict the fugacity coefficients, and the equilibrium condition must be solved during optimization. Kriging models and Artificial Neural Networks (ANNs) were used to predict the fugacity coefficients [165]. The performance of the models trained using a LHS design, the mixed adaptive sampling with parameterization by Eason and Cremaschi [33] and, the mixed adaptive sampling with the best-found parameterization were compared (Figure 5.30) for ANNs as surrogate models. The model trained with LHS performed worst; it had the largest MSE and the largest spread of the prediction error. The best performance was obtained by the application of the peaks-candidate selection, and an optimal weighting of the exploratory and exploitative sampling objectives.

The miscibility gap that is predicted when using the surrogate models for the fugacity coefficients in the computation of the phase equilibrium is shown in Figure 5.31. The prediction of the gap and also of the compositions in the different phases is very accurate. Surrogate models were used in the literature [166] to predict the gas

Figure 5.29: MAEi for the Ki of the three components over the Kriging sample size of the sequential sampling algorithm for NSS = 18 and SF = 0.4 showing the mean (markers), minimum, and maximum (bars) of each 5 runs of the models, trained in the mixed adaptive sampling algorithm and the models based on LHS design [163].

solubilities, and the phase separation in the decanter of the hydroformylation process of 1-dodecene, as shown in Figure 5.2. As the original PC-SAFT model is computationally expensive, the optimization of the operating point of the process can only be performed using surrogate models. In this case, the results of the ternary system which are shown in Figure 5.31 had to be extended by the product component, 1-tridecanal, for describing the behavior of the system in the hydroformylation process.

For this quaternary system, both the indirect and the direct surrogate-modeling approach were applied. For the indirect method, Artificial Neural Networks (ANN) and Kriging models were trained as regression models. For the direct method, ANN classifiers were trained to describe the biphasic operating region in the decanter. The classifiers were combined with Kriging models ("direct Kriging") and with ANN models ("direct ANN") of the same structure, as in the case of the indirect method for regression. The test set errors of the four different model combinations for test sets of 10,000 samples each for the gas solubilities and for the phase separation are shown in Tables 5.1 and 5.2, respectively. Both for the gas solubilities and for the liquid-

Figure 5.30: Comparative plot of the scaled MSE of the fugacity coefficients MSE, φ̅, for the ternary system of 1-dodecene, DMF, and n-decane for the ANN models over the training sample size using the mixed adaptive sampling algorithm with the parameterization from Eason and Cremaschi [33] (blue), the peaks-candidate selection (red), and using LHS designs with exponential trend lines [165].

Figure 5.31: Liquid-phase compositions at different temperatures for the ternary system of 1-dodecene, DMF, and n-decane and predicted values using the surrogate models (dashed lines) compared to PC-SAFT (lines) (Section 4.1.1) [165].

Figure 5.32: Flowsheet of the hydroformylation process of 1-dodecene [164]. The flowsheet corresponds to loop I presented in Figure 4.16.

liquid separation, the models applying the direct surrogate-modeling approach outperformed the surrogate models using the indirect surrogate-modeling approach. Especially for the liquid-liquid separation, for which the MAE of the product composition is shown in Table 5.2, the models using the indirect modeling approach show large errors of up to 3 mol% in composition. While for the indirect method, the chosen surrogate model structures were not able to approximate the behavior of the process with sufficient accuracy, the application of the direct method led to very accurate surrogate models, both for the gas solubilities and the liquid-liquid separation, with a maximum MAE in the product composition of 0.62 mol-% for the component DMF. The fraction of the misclassified samples was at most 1.90% in the direct Kriging case. As the classifier and the regression model were trained simultaneously for different adaptive samplings, the classification errors are also different.

Table 5.1: MAE of the predictions of the trained surrogate models [166].

Table 5.2: Misc and MAE of the predictions of the trained surrogate models for the molar fractions in the product phase of the decanter [166].

A comparison of the computational expense of the different surrogate model options and the original model PC-SAFT for process simulations is shown in Table 5.3. The surrogate models provide a more than 90 times faster process simulation than using the original model PC-SAFT. As many process simulations are performed during optimization, this shows that the original model PC-SAFT is practically not applicable for process optimization.

Table 5.3: Simulation time for the hydroformylation process of 1-dodecene on a Windows 10 4.2/ 4.2 GHz dual Intel(R) CoreTM I7 machine with 32 GB RAM using the SciPy-lm solver [166].

As the indirect surrogate models were not sufficiently accurate, only the models trained using the direct surrogate-modeling method were used to optimize the operating point of the process: the feed composition (molar flow n\_ of component 1 dodecene, DMF, n-decane and catalyst), the reaction temperature (T), the reaction pressure (pR), the carbon monoxide content of the synthesis gas (yCO), and the separation temperature of the decanter (TD). The cost function considered for the optimization is the production cost in Euro per ton t of n-tridecanal for a capacity (m\_ nC13al) constrained to 10,000 tons per year. Both operational costs Cop and investment costs Cinvest were taken into account:

$$\text{CPT} = \frac{\mathcal{C}\_{op} + \frac{\mathcal{C}\angle\text{invest}}{10}}{\dot{m}} \,. \tag{5.31}$$

The resulting optimal operating points are shown in Figure 5.33. For these conditions, the process was also simulated with PC-SAFT [166]. The resulting errors in the objective function are shown in Table 5.4. Compared to the process simulation results using the original model PC-SAFT, the surrogate model approaches both results in a low relative error of the objective function–less than 3%. The ANN models that provided the best approximation accuracy for the test set for the liquid-liquid separation (Table 5.4) also showed the best result here.

In this section, two surrogate-modeling approaches for complex thermodynamic models that are used in process simulation and flowsheet optimization were presented. In the direct modeling approach, the thermodynamic equilibrium is solved during the sampling, the thermodynamic behavior is described by a classifier that predicts the number of phases, and a regression model describes the distribution coefficients in the biphasic region. In the indirect approach, the fugacity coefficients are predicted by the surrogate models; thus requiring a solution of equilibrium condition during the simulation or optimization of the process. In the ternary case, the

Figure 5.33: Normalized optimal operating points (normalized to the ranges of the variables) for the hydroformylation process of 1-dodecene using different surrogate models [166].

Table 5.4: Objective function values of the optimization using surrogate models compared to PC-SAFT simulations of the resulting optimal operating points [166].

approximation of the LLE using both the indirect and the direct method was possible with small approximation errors. This shows that both approaches can provide accurate predictions of thermodynamic behavior. For the quaternary case, the indirect method showed an inferior performance, compared to the direct model. This could be overcome by using more complex surrogate-model structures.

The direct surrogate-modeling approach was successfully applied in the optimization of the operating point of the process for the hydroformylation of 1-dodecene, shown in Figure 5.32, by minimizing the production price per ton of the product. The relative error in the objective function in the best case is less than 1%, compared to PC-SAFT simulations using ANN regression models, which is certainly less than the uncertainty of the rigorous model. The results also show that the average quality of fit of the surrogate models over the test set is not necessarily a reliable indicator of the accuracy of the prediction of the process optimum; so the strategies that include the sensitivity of the optimum in the selection of the samples should be investigated, as shown by Winz and Engell [217]. The use of surrogate models in the global optimization of flowsheets is discussed in Section 5.2.3 and the application within the earlystage process design is demonstrated in Section 6.4.

#### 5.3 Process Optimization

#### 5.3.1 Optimal Design of Reactors for Complex Reaction Systems

Karsten H.G. Rätze, Michael Jokiel, Kai Sundmacher

The reactor unit(s) of a chemical production plant significantly influences the overall production costs since they determine the requirements for the subsequent, generally energy-intensive, downstream process. Efficient reactors that provide high selectivity towards the desired product are a prerequisite for a lean separation process. Therefore, the identification of optimal operation strategies for existing reactors and the synthesis of new processes in the form of conceptual design are topics of interest and subject to active research [156].

The design of new or the optimization of existing processes benefits from expert knowledge and from detailed, experimental investigations to develop a thorough understanding of various kinetic and thermodynamic effects. However, the optimization potential is generally limited due to the restriction to the reactive units in contrast to the overall process and/or constrained flexibility w. r. t. reactor designs and reactor operation strategies. Consequently, a holistic approach for the conceptual design of efficient chemical production processes is required, in which information from various sources is incorporated. An important tool that fosters an extended view on the reactor design task represents the utilization of reaction kinetic and process models, either in the form of mechanistic or data-driven models, in addition to experimental investigations. While mechanistic models are generally suited to achieve a detailed understanding of processes by providing reliable data for the investigation and analysis of operating windows/strategies that were not subject to any experimental investigations due to either limited resources or safety concerns (more information in this chapter), data-driven or surrogate models are particularly suitable for the aggregation, interpolation, and efficient access to information (Sections 5.2.3 and 5.3.2). As a consequence, both model types are indispensable for the state-of-the-art process design.

Despite these powerful, predictive tools, incomplete knowledge and uncertainty/ randomness are inherent in process development; so the incorporation of these aspects in a quantitative manner in the design process allows for better decision making at early stages. Consequently, only a combination of model-based tools for process design with experimental investigations on multiple scales (lab and mini plant scale), can achieve reliable designs and, therefore, represents an essential part of process design. Especially, the investigation on different scales does not only serve the validation of simulative predictions but allows for the identification of unknown effects or unaccounted behavior – information that can be used to update process models and improve the process design. This interplay of process design and experimental validation already hints at an iterative and strongly interdisciplinary development that is mandatory for the modern and reliable design of chemical production processes.

This section provides a brief introduction to the systematic and model-based unit operation design in the process development cycle (Step 4a in Figures 4.63 and 4.65), with the emphasis on reactor design and optimization in the context of the entire production process. The next paragraphs touch upon the different strategies for reactornetwork synthesis before reaching deeper into model-based reactor-synthesis via the methodology of Elementary Process Functions and its application to the hydroformylation of long-chain olefins. Subsequently, detailed simulative and experimental investigations are presented w. r. t. the reactor units as well as their interaction with the entire process.

#### <IP_ADDRESS> Reactor-Network Synthesis

The selective synthesis of the desired product might require a careful adjustment of the different species concentrations inside the reaction unit for an optimal result. This optimal concentration profile is unique for each reaction and might change with external factors, like temperature and pressure. In reactor-network synthesis, the focus lies on the approximation of this optimal concentration profile using idealized reactor units with known behavior. These idealized reactor units are generally differentiated w. r. t. their internal mixing characteristics. The most common types are (i) the continuously stirred tank reactor (CSTR), (ii) the plug-flow reactor (PFR) and (iii) the distributed side-stream reactor (DSR). While the first two units represent edge cases with maximum and minimum axial dispersion, the last reactor unit enables the selective dosing of species and species mixtures along the reactor coordinate to achieve compromises between CSTR and PFR dispersion behavior [39].

Most of the approaches to reactor-network synthesis can be classified either as heuristics-, attainable region- (AR) or rigorous optimization-based, with the latter two representing model-based approaches. Especially, rigorous optimization-based approaches are relatively new, compared to the other categories, so that their full potential is still unexplored and are subject to active research. Superstructure optimization represents the most intensively studied sub-category of rigorous synthesis approaches in which the solver creates an optimal reactor network by using a list of predefined process units. While yielding high-performance reactor combinations and potentially unintuitive solutions for large-scale networks, high computational loads and the lack of innovative designs and process intensification led to the emergence of alternative approaches. One of these alternatives is the Elementary Process Functions methodology, which is based on dynamic optimization and tackles the synthesis problem via functional modules, instead of predefined process units. For a better understanding of the methodology, the general idea as well as a selected number of extensions and frameworks are discussed in the next section.

#### <IP_ADDRESS> Elementary Process Functions Methodology

The Elementary Process Functions (EPF) methodology was introduced by Freund and Sundmacher [43] to shift the focus from connecting unit operations, e.g., reactor, separator, for the flowsheet generation to the consideration of functional modules, instead. These functional modules represent tasks, like heating and separating, which need to be performed in a process, without restricting the execution of these tasks to predefined apparatuses. This shift in perspective facilitates the generation of new process units and the direct consideration of integration potential in the process design.

In its original formulation, the reactor-network design task can be represented as a fluid element traversing the thermodynamic state space from the state at initial time x tð Þ0 to the final state x tð Þf , where x 2 χ - Rx represents the state vector, containing concentrations, temperature, pressure, etc. Its path is governed by process constraints, but can also be influenced by fluxes entering and leaving the fluid element, which act according to a predefined objective function (see Figure 5.34 for a schematic representation). This objective function can be defined in terms of the final state x tð Þf or w. r. t. the fluid element's path in state space, leading to an Optimal Control Problem (OCP) of the form

$$\min\_{\mathbf{u}(t),\theta} \phi(\mathbf{x}(t),\mathbf{u}(t),\theta,\mathbf{c}) \tag{5.32}$$

$$\mathbf{C}(\mathbf{x}(\mathbf{t}), \mathbf{c}) \ \frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x}(\mathbf{t}), \mathbf{c}) \ \mathbf{J}(\mathbf{x}(\mathbf{t}), \mathbf{u}(\mathbf{t}), \mathbf{c}), \ \mathbf{x}(\mathbf{t\_0}) = \mathbf{x\_0},\tag{5.33}$$

$$h(\mathbf{x}(t), \mathbf{u}(t), \boldsymbol{\theta}, \mathbf{c}) = \mathbf{0},\tag{5.34}$$

$$
\tilde{h}(\mathbf{x}(t), \boldsymbol{\mu}(t), \boldsymbol{\Theta}, \mathbf{c}) \le \mathbf{0}, \tag{5.35}
$$

$$\mathfrak{x}(t) \in \mathbb{X}, \mathfrak{x}\_0 \in \mathbb{X}\_0, \mathfrak{u}(t) \in \mathfrak{A}, \mathfrak{G} \in \mathcal{P},\tag{5.36}$$

with the Bolza type objective function ϕ, general equality h, and inequality constraints ~ h, constants c, and bounds on the states xTð Þt = T tð Þ, p tð Þ, wTð Þt , v tð Þ, G tð Þ , the static controls θT = xT 0 and the control profiles uTð Þt = m\_ Að Þt , j T Að Þt , q\_ Að Þt , τ Að Þt , fVð Þt , wtð Þt , g tð ÞÞ. Here, the state vector comprises the temperature T, pressure p, mass fractions w, velocity v, and geometrical state of the fluid element G, which is the fluid element's volume V, in the case of a batch process, or the volume flux vAc for continuous processes, with Ac denoting the cross-sectional area of the fluid element. The control profile vector u contains the total mass flux m\_ A, diffusion flux density vector jA, surface stresses τA, volume forces fV, heat fluxes q\_ A, technical work wt, and rates of the matter element's geometry variation g. For notational convenience, the time and state dependencies were neglected and all lower-case variables can be considered as vectors. The Ordinary Differential Equation (ODE) System eq. (5.33) is a compact representation of the balance equations that govern the feasible paths in the thermodynamic state space with the capacity matrix C, the weight factor matrix F, and the generalized flux vector JTð Þ x tð Þ, t = rT Vð Þ x tð Þ , uTð Þt . In the case of a regular capacity

matrix, the eponymous linearly independent elementary vectors span the column space of E = C − 1 F, which represents the attainable sub-space χ, according to the process constraints [204].

With the mathematical foundation of the EPF methodology in eqs. (5.32)–(5.36) describing the identification of controls that characterize the state trajectories in an optimal reactor, frameworks are required to create real-world reactor networks from these optimal control profiles. Two of these frameworks will be discussed in the next paragraphs, with both concepts approaching the approximation of the optimal control profiles from different perspectives.

#### Reactor-Network Synthesis Framework 1 – Successive Approximation

In eq. (5.33), it was already assumed that the internal reaction fluxes rT Vð Þ x tð Þ and, consequently, JTð Þ x tð Þ, t are dependent on the states x tð Þ. This is reasonable, since the standard reaction rate equations normally exhibit a concentration, temperature, and, sometimes, pressure dependency. Besides this simplification, solving eqs. (5.32)–(5.36) yields controls that optimally steer the fluid element in state space. However, this path might not be achievable in real life due to limitations on the fluxes in u tð Þ. In actuality, these fluxes result from driving forces that are dependent on the distance of the states to their respective equilibrium values, x tð Þ − xeqð Þ x tð Þ . Consequently, the introduction of rigorous flux models yields control profiles that approximate the ideal control profiles and are, therefore, closer to reality. A generalization of this approach was introduced by Peschel et al. [173], in which the successive refinement of the underlying model is formalized in a three-level procedure (Table 5.5), yielding a systematic approach for the synthesis of reactor networks. With each level, the complexity of the underlying process model increases, so that the maximum potential of the reaction without any limitations (Level 1) can be compared to a hypothetical reactor with (mass and energy) transport limitations and a limited number of controls (Level 2) and, eventually, to multiple technical approximations (Level 3).

Table 5.5: Three-level procedure for the design of optimal technical reactors using the EPF methodology [44, 173].

Figure 5.34: Schematic representation of a fluid element traversing the three-dimensional thermodynamic state space, from x tð Þ0 to x tð Þf . The path of the fluid element is constrained by the internal reaction flux rV and can be influenced by the volume forces fV, diffusion flux of species α over the surface jA, α, the total mass flux m\_ A, heat fluxes q\_ A, technical work wt, the surface stress τA, and the rate of geometrical change of the fluid element g. Please note all V indices relate to volume-related fluxes and all A indices represent surface-related fluxes (adapted from [204]).

Figure 5.35: Division of control profiles u tð Þ based on distinct patterns and subsequent interpretation via flux-profile analysis (FPA) and association to ideal reactor types CSTR, PFR, and DSR (adapted from [90]).

Table 5.6: Three-step procedure for the design of chemical reactor networks using the EPF methodology and flux-profile analysis [92].

Figure 5.36: Derived approaches, extensions and applications of the EPF methodology since 2008. InPROMPT contributions highlighted. Selected contributions are: a [90] and [92], b [94], c [125, 126, 191–193], d [218], e [91], f [173], g [65], h [66, 175], i [65, 67, 90], j [65, 174], k [66], l [219], m[141], n [173, 174, 176], °[177] and p [93].

#### Reactor-Network Synthesis Framework 2 – Flux-Profile Analysis (FPA)

While the successive approximation framework predominantly relies on the practitioners' knowledge, expertise, and intuition to select controls for the successful identification of suitable technical realizations, the Flux-Profile Analysis (FPA) emphasizes the design of reaction-networks and the selection of suitable classes of reactors before improving the technical feasibility [92].

After the initial identification of the maximum potential of the reaction system in Step 1 (Table 5.6) similar to Peschel et al. [173], all control profiles are approximated using a network of ideal reactor units. Here, the key idea lies in the observation that the optimized outer mass- and energy fluxes form specific patterns over the processing time. By subdividing or discretizing the time horizon according to these patterns, ideal reactor units, like a CSTR, PFR, or DSR, can be used to approximate the control profile and, consequently, the state trajectory of the fluid element. A schematic representation of this procedure can be found in Figure 5.35. For the accurate association of a control profile with an ideal reactor type, the consideration of axial dispersion or back-mixing is required.

In a first approach, axial dispersion can be rigorously considered by introducing internal fluxes jmixð Þt alongside the diffusion fluxes in eqs. (5.32)–(5.36). These fluxes provide a connection of the fluid element with itself in time, allowing for the identification of additional integration potential, including process-wide recycles [218]. This is especially useful for the selection of suitable reactor units in which sections with high internal fluxes jmixð Þt indicate a positive influence of axial dispersion on the reaction, while low internal fluxes hint at reactors with limited back-mixing behavior.

In complex reaction networks with many species, the rigorous consideration of internal fluxes becomes computationally challenging due to the temporal interconnection of the fluid element. For these cases, a simplified formulation of axial dispersion can be chosen with the introduction of storage tanks for each species, e.g., substrate, intermediates, and products, in the reaction network. Each storage tank contains a certain amount of initial fill volume at t0 which the solver can use for dosing into the fluid element. With additional constraints to ensure dosing of the entire substrate amount and to prevent unsustainable insertion of (by-)product into the fluid element, axial dispersion and process-wide recycle streams can be approximated. As an indicator for the benefit or disadvantage of axial dispersion in the case of absent dosing fluxes in selectivity problems, the differential selectivity

$$\varphi = \frac{\text{moles of desired product formed}}{\text{moles of main reactant consumed}} = \frac{\frac{\text{dn\_{prod}}}{\text{dt}}}{\frac{d\text{dn\_{Sub}}}{dt}} \tag{5.37}$$

or rather its derivative w. r. t. time, the differential reaction flux, can be evaluated. Sections with a positive sign of the differential reaction flux indicate beneficial and negative signs, disadvantageous influences on the selectivity toward the desired product [92].

#### EPF and Framework Extensions

Besides the extension of the EPF methodology towards reactor-network synthesis, many other extensions and applications can be found in the literature, with a selection of major advances depicted in Figure 5.36. All publications in the overview are grouped into categories, depending on either the methodological improvement, e.g., Reactor-(Separator)-Network Design and Back-Mixing, or the application, e.g., Processes and Economic. In the following paragraphs, each methodological category will be discussed, starting with the reactor-separator-network design as the natural yet complex evolution of reactor-network design.

The reactor-separator-network design requires the introduction of a concept for separation in the EPF methodology. This can be achieved by extending the bounds of the dosing and diffusion fluxes, jA, toward the negative domain, making the withdrawal of species from the fluid element possible. Instead of considering the continuous withdrawal of (pure) components, which is only possible in special cases, discrete separation units are introduced at the transition points between profile patterns in Step 2 of the FPA [94]. The introduction of discrete separation points leads to a mixedinteger dynamic optimization (MIDO) problem, requiring specialized algorithms for solving or the reformulation by relaxing the integer variables and ensuring the convergence toward the integer solution via continuous switching functions. An alternative and more radical approach to the reactor-separator-network synthesis task aims at the discretization of the entire thermodynamic state space using predefined operations, like reacting, heating, separating analogous to the introduction of functional modules in the original EPF methodology [125]. This results in a separation of the evaluation of the (nonlinear) reaction model from the optimization of the trajectory in state space, leading to a graph-theoretical approach with a linear optimization problem. This approach even allows for simultaneous consideration of heat integration [126] and application to non-reactor units, like separators [191].

The extension toward multi-phase systems is crucial for the industrial application of the EPF methodology since numerous systems exhibit or require two or more phases during the reaction. In the simple case of one reactive and one auxiliary phase, like in the case of the homogeneously-catalyzed hydroformylation with the gas phase representing the auxiliary phase, the auxiliary phase can be modeled as a continuum enveloping the liquid fluid element. This continuum represents a service phase with an unlimited capacity and introduces an interphase driving force for species fluxes to and from the fluid element [175]. For an accurate description of the driving force, sophisticated equations of state (EoS), like PC-SAFT for liquid-liquid (L/L) and gas-liquid (G/L) equilibria [66], diffusion models for gas-solid (G/S) [177] and combinations of G/S and G/L [93] equilibria can be employed. An alternative, more rigorous representation of multi-phase systems in the EPF methodology considers multiple fluid elements traversing state-space simultaneously [219]. In this formulation, each phase receives a separate fluid element with unique control fluxes. Consequently, mass transfer between both phases can be modeled accurately. The service phase and the rigorous modeling approach are visualized in Figure 5.37.

Most extensions of the EPF methodology introduce additional degrees of freedom or higher accuracy via rigorous representation of physicochemical phenomena, while assuming the process/reaction model to be accurate. However, this is not the case due to structural inaccuracies and uncertain parameters in the model as well as nonideal reactor units and imperfect realizations of control profiles in the final plant. As a consequence, a probabilistic design framework is required to achieve robustness of the final reactor design. Since the direct consideration of the probability density function (PDF) of the uncertain variables in the optimization problem is not feasible from a computational perspective, an approximation of the PDF via the first and second statistical moments (mean and standard deviation/variance) in the unscented transformation or sigma point method represents a suitable trade-off between accuracy and applicability [91]. With the unscented transformation, both statistical moments can be mapped from the parameter to the criterion space of the objective function using the nonlinear process model. By formulating a bi-level optimization problem, the Pareto front comprising the expected value as well as the variance of the objective function can be calculated and used to identify the desired compromise of performance and robustness. An illustration of this procedure is shown in Figure 5.38. Due to the transition of the chemical industry toward sustainable feedstocks in chemicals production, changing sources of feedstock and fluctuating resource qualities become additional sources of uncertainty in the conceptual design. These uncertainties can be considered with the EPF methodology by taking a multi-scenario approach and solving, analogous to the unscented transformation, a multi-objective optimization problem in terms of performance, safety, and uncertainty of the reactor design [141].

#### <IP_ADDRESS> EPF Application to the Hydroformylation of Long-Chain Olefins

Many advances in the EPF methodology from Section <IP_ADDRESS> are the result of the challenges in designing optimal reactor networks for the hydroformylation of longchain olefins. While the derivation of accurate reaction-kinetic models for these example reactions is discussed in detail in Section 3.2, the following paragraphs provide a compact overview of the application of the EPF methodology to the reactor-network design for the homogeneously catalyzed hydroformylation reaction.

The biphasic nature of the hydroformylation reaction, in which a liquid phase is in contact with a gaseous phase containing the necessary substrates hydrogen (H2) and carbon monoxide (CO), initiated the extension of the EPF methodology towards multiphase reactor design. With the introduction of the service phase approach, the first design of an optimal hydroformylation process was proposed for 1-octene in an

Figure 5.37: Consideration of multiple phases in the EPF methodology. a) Service phase and b) rigorous modeling of multiple fluid elements with interaction in terms of energy and mass transfer (adapted from [219]).

Figure 5.38: Robust reactor-separator-network design via sigma point approximation of the parameter uncertainty in the EPF methodology. Left: Mapping of the parameter space Θ via the nonlinear process model hð Þ θ, D into the criterion space O with θ 2 Θ denoting uncertain parameters and D representing the chosen process design. Right: Pareto front indicating the trade-off between the expected value and the variance of the criterion o 2 O (taken with permission from [91]).

ionic liquid (IL) [175]. By using reaction kinetics from the literature in the successive approximation three-level design [173], three reaction segments with distinct temperature profiles were identified in the technical approximation. The final reactor design is shown in Figure 5.39 and consists of a DSR comprising three temperature segments with discrete syngas dosing positions at the boundary of each segment. Additionally, two process-wide recycle loops, one for the IL and the second for unconverted substrates, are proposed, leading to an improved performance with selectivities of S = 71.2% to S = 62.1%, in the CSTR reference case at comparable or higher space-time-yields (STY) [175].

Due to the importance of gas-liquid phase interactions in the hydroformylation, the influence of liquid phase CO and H2 concentrations on the reaction performance was investigated in more detail. A simplified variant of the reaction network as in Figure 5.40 enabled the inclusion of the PC-SAFT EoS (Section 3.1.1) in the EPF methodology. After calibration of the EoS to literature data, this combination enabled the rigorous consideration of the G/L phase equilibrium during the reactor design. Multiple optimization studies with increasing degrees of freedom (DoF) revealed the importance of the gas phase composition on the selectivity towards the desired product. The investigations identified a decrease in solubility with progressing hydroformylation of the reactant. By optimally controlling the temperature, pressure, and the H2 to CO ratio, this effect can be counteracted, leading to a differential selectivity S > 90%, over the entire conversion interval [65]. These results were also confirmed for the hydroformylation of 1-dodecene using a rigorous reaction kinetic model in a flowsheet optimization, as shown by Hentschel et al. [65], with an "EPF reactor" at its core (Figure 5.41). The rigorous reaction kinetic model further revealed that recycling, not only of the unconverted terminal olefins but also of the olefin isomers, proves beneficial for the selectivity and conversion in the reactor. By shifting the equilibrium of the olefin isomerization reaction toward the terminal olefin, a selectivity of 94.2% at 97.3% conversion was achieved in contrast to the reference case of a CSTR with 93.4% selectivity at 51.6% conversion [65]. This represents a 90% increase in linear aldehyde yield in comparison to the reference case. While this increase in product yield significantly simplifies the downstream processing and reduces recycle streams, process profitability did not increase accordingly. An optimization of the entire reaction-extraction process of the homogeneously rhodium-catalyzed hydroformylation in a thermomorphic multiphase system (TMS) revealed that catalyst leaching, which occurs in the decanter for catalyst recovery while being on the ppm scale, represents the cost determining factor. As a remedy, a counter-current decanter cascade can be introduced, which effectively prevents catalyst leaching (Table 5.7) while simultaneously raising the temperature for the liquid-liquid separation to ambient conditions [142].

In addition to the performance optimization of the hydroformylation process, the influence of uncertainty on the process performance indicators was analyzed in the context of the EPF methodology. Different types of uncertainty, ranging from parameter uncertainty over uncertain residence time distributions to uncertain realizations of the control profiles, were considered. While the latter two sources of uncertainty mainly reduced the expected selectivity and increased the variance of the selectivity, a significant influence of the parametric uncertainty on the reactor performance was identified. For expected conversions X ≥ 97%, the concentration and control profiles revealed two distinct reaction zones across all realizations of the uncertain parameters. This formation of reaction zones became possible due to the increase in residence time, which was required to achieve high conversions. While the first reaction zone is characterized primarily by dynamic control trajectories, the second reaction zone, which was already identified in previous publications (e.g., [65]), exhibits constant, higher temperatures and reduced carbon monoxide concentrations for the intensified re-isomerization of olefin isomers to terminal olefins and, consequently, linear aldehydes.

Figure 5.39: Technical approximation of the optimal reactor unit for the hydroformylation of 1-octene in an IL. based on the EPF methodology (taken with permission from [175]).

Figure 5.40: Simplified reaction network for the hydroformylation of 1-octene. 1-C8en: 1-octene, i-C8en: iso-octene, n-C9al: n-nonanal, i-C9al: iso-nonanal, n-C9ol: n-nonanol, i-C9ol: iso-nonanol (adopted from [175]).

However, with increased robustness of the reactor design, the percentage of the second zone of the total residence time is reduced to account for the increased sensitivity of the first reaction zone w. r. t. the parameter uncertainties [91]. This behavior leads to a pronounced Pareto front, as shown in Figure 5.42, and motivates experimental investigation. To capture not only microscopic but also process-wide effects, including long-term process stability, the process design needs to be built and evaluated in practice under mini plant-scale conditions.

#### <IP_ADDRESS> Proof of Concept: Optimal Reactor-Design Hydroformylation of 1-Dodecene

For the hydroformylation of 1-dodecene (nC12en) in the TMS of N,N-dimethylformamide and n-decane, the optimal reactor configuration was sized and designed [90]. The goal of the design task was to maximize the selectivity towards the linear aldehyde tridecanal (nC13al) for high nC12en conversions between 90 % and 99 % . Additionally, the formation of undesired branched aldehydes should be avoided by

Figure 5.41: Hydroformylation of 1-dodecene flowsheet. R: "EPF reactor", D: decanter, C1: solvent/ product distillation column, C2: n/iso aldehyde distillation column. Heat exchanger units, pumps, and compressors are denoted as H, P, and CP, respectively (taken with permission from [65]).

Table 5.7: Percentage of catalyst cost of the total annual costs (TAC), total production costs, and the production costs without 1-dodecene for seven cases. Case one to six contain decanters in the counter-current decanter cascade. The reference case only contains one decanter [142].

employing a linear to branched (n/iso)-aldehyde ratio of ≥ 95:5. To ensure that the reactor volumes are within the mini plant-scale, the total residence time was restricted to a maximum value of 300 min and the mass flow rate of the already existing mini plant [32] was considered. This also allows for a direct performance comparison between processes using a conventional stirred-tank reactor and optimal reactor designs.

Figure 5.42: Pareto front of expected selectivity versus selectivity variance for parametric uncertainties at 99.5% expected conversion. Exemplary designs A, B, and C denote Pareto optimal points with low, middle, and high variance (taken with permission from [91]).

Figure 5.43: Optimal reactor designs from Kaiser et al. [90] for the hydroformylation of 1-dodecene (adopted from [89]).

#### Derivation of the Optimal Reactor Configuration via Flux-Profile Analysis

Four optimal reactor network candidates, comprising different DSR and CSTR combinations, were derived with the FPA for the hydroformylation of nC12en in a TMS [90]. The analysis of the conversion-selectivity behavior of these networks revealed that the combination of a DSR followed by a CSTR is the most robust design.

For the intended construction and continuous operation of the DSR-CSTR-tandem, the continuous control profiles along the DSR need to be approximated in a discrete manner. Therefore, the influence of different positions and numbers of dosing points along the DSR was evaluated, revealing only a slight loss of selectivity for low conversions if gas and product dosing is performed only at the inlets of both reactors. This enables a rather simple construction of the DSR in the form of a PFR. In order to size this PFR, the determination of the flow regime is crucial to ensure high gas-liquid mass transfer rates and low axial dispersion, so that the reaction performance is not influenced by those parameters. After further analysis of both parameter limits and a literature search, the segmented Taylor regime was identified to enable those conditions. As this flow regime requires a rather small tube diameter (d ≤ 3 mm) for the fluid properties of the TMS, a tube with a length of 300 m is needed for the necessary reactor volume of 2.2 L. To benefit from the positive effects of homogenization (low axial dispersion but intensive radial mixing, Section <IP_ADDRESS>) and to allow for a compact design, the coiling of the tube was suggested. This design also enables simpler heating of the reactor via heating tapes through which four heating zones along the reactor can be realized. The flowsheet of the helically coiled tubular reactor (HCTR), followed by the CSTR, is depicted in Figure 5.43.

In addition to the continuously operated reactor design with approximated control profiles along the DSR, an alternative realization of the continuous control profiles is possible. By means of a repeatedly operated semi-batch-reactor (RSBR), the control profiles can be applied in time, in contrast to the length of the reactor, leading to simplified adjustments of the controls in the case of changing model structures and parameters. As the SBR and the subsequent CSTR differ w. r. t. their operation mode, two buffer vessels, before and after the SBR, are required so that the SBR can be operated repeatedly. Preservation of the catalyst activity and prevention of the liquid-liquid phase separation requires both buffer tanks to be heated and pressurized, leading to two additional reaction zones. While the feed buffer tank does not significantly contribute to the reaction progress because the substrate is dosed directly into the reactor, the flash buffer tank, which serves the reaction mixture to the CSTR, needs to be accounted for. In terms of the axial dispersion and adequate gas-liquid mass transfer, the SBR itself is free of axial dispersion and the stirrer intensifies the gas-liquid contacting, leading to the fulfillment of both requirements of a DSR. A flow chart of this repeatedly RSBR-CSTR-tandem is shown on the right of Figure 5.43.

After designing both reactor-tandems, their performance was evaluated by a more detailed process optimization (Figure 5.44). In contrast to the optimizations within the reactor design task, no ideal liquid-liquid separations were considered for the recovery of the TMS catalyst phase and the byproduct separation but the results show that both tandems have their selectivity and conversion optima at 95% and 98%, respectively, which indicates an enhancement of 24% and 40%, in comparison to a process with a single CSTR.

#### Reaction Condition Optimization of the Constructed Reactor Tandems

Based on the provided specifications from the design task, both reactor-tandems were constructed and operated with continuous catalyst recycling but without byproduct separation via a distillation column [87, 88, 183]. In order to compare the results from the operation of both tandems, the catalyst concentration and the total residence time were fixed in accordance with a former study where a single CSTR was used [32]. The reaction conditions (temperature, pressure, and gas compositions) for the operation were determined with the help of process models optimization w. r. t. the selectivity towards the linear aldehyde nC13al.

For the HCTR-CSTR-tandem, the reaction conditions were derived with the help of the model from the reactor design task. However, the process model had to be modified to account for the approximations during the construction, i.e., the length of the HCTR had to be shortened due to the fixed residence time of the experiment. Therefore, only a reactor length of 80 m was used so that only two heating zones along the HCTR could be used for temperature control. The derived conditions of the HCTR-CSTR-tandem are depicted on the left of Figure 5.45. It can be observed that a slight excess of carbon monoxide and a temperature increase within the second heating zone of the HCTR are favorable in the first reactor. In the subsequent CSTR, which represents the second reaction zone, the temperature should stay the same but an excess of hydrogen is necessary.

For the second tandem, comprising the RSBR and CSTR, a more detailed dynamic model was developed to not only derive the necessary reaction conditions but also aid in the process start-up. With the model simulations, the operator receives valuable information on the time scale to reach the cyclic steady state of the reactor tandem. One particularly valuable piece of information for the RSBR-CSTRtandem to ensure comparability between the different reactor networks involves the distribution of the fixed total residence time to each reaction unit. In addition to the residence time of the SBR τSBR = tSBR and the CSTR τCSTR, PFR = VLV\_ − 1 with the liquid hold-up VL and the volumetric flow rate V\_ , both buffer tanks as well as the preparation time or inactive time tI (see Figure 5.46 for a schematic representation) of the SBR need to be taken into account. Please note that the preparation time is a required, fixed time in which the charging and discharging of the SBR vessel is performed. While the feed buffer tank contains only traces of unconverted substrate and can, therefore, be neglected in the distribution of the total residence time, the flash buffer tank represents an additional, third reaction zone. With the approximated average residence time in the flash buffer tank, τDBuffer = 0.5VDBuffer L ð Þ t = 0 V\_ CSTR − 1 in , with t = 0 denoting the beginning of each process cycle, the available residence times for the SBR and CSTR can be computed with the distribution and coupling constraints

$$
\tau\_{\rm SBR} + \tau\_{\rm I} + \tau\_{\rm DBuffler} + \tau\_{\rm CSTR} - \tau\_{\rm Total} = \mathbf{O},\tag{5.38}
$$

$$
\boldsymbol{\tau}\_{\text{SBR}} + \boldsymbol{\tau}\_{\text{I}} - \boldsymbol{\tau}\_{\text{DBuffer}} - \boldsymbol{\tau}\_{\text{CSTR}} = \mathbf{O}.\tag{5.39}
$$

Here, the coupling constraint ensures that the reaction in the SBR terminates and the reaction mixture can be transferred to the flash buffer tank shortly prior to the depletion of the liquid hold-up in the flash buffer tank, which supplies the continuous part

Figure 5.44: Flow chart of the overall process used for the process optimization (top) and performance comparison of the different reactor designs (bottom) The lines correspond to the maximum potential for ideal separations (taken with permission from [90]).

of the process [183]. Since τDBuffer and τCSTR are coupled via the constrained control parameter V\_ CSTR in , τSBR needs to be recalculated for each new cycle. This is especially important during process start-up due to concentration-induced changes in the reaction mixture density. Under ideal conditions, including the perfect realization of the controls and time-resolved adjustment of all fluxes in the process, the cyclic steadystate is achieved after approximately 30 h with n-dodecane (nC12an) being the ratedetermining component (see Figure 5.47). In contrast, the concentration profiles of the substrate (nC12en) and the product (nC13al) remain unchanged after 10 h of operation or 7 RSBR cycles [183].

In addition to the calculation of the start-up time, the residence time of each reactive unit and the process conditions, the dynamic RSBR model provided estimates of the flow rates in each pipe for each cycle in the experimental validation of the RSBR-CSTR-tandem [88]. The optimal control parameters for the SBR are depicted on the right side of Figure 5.45. The profiles for the CSTR are neglected since they are analogous to the corresponding profile in the HCTR-CSTR-tandem (Figure 5.45, lefthand side). Despite different operation modes and varying degrees of freedom, the optimal reaction conditions for both reactor tandems in Figure 5.45 are comparable approximations of the true optimal control profiles [89]. In terms of the hydroformylation reaction network, the low initial temperature and the slight excess of carbon monoxide causes a minor reduction of the catalyst activity through which all reactions and, especially, the isomerization of the terminal olefin as the main side reaction are inhibited. Subsequently, with increasing yield, a higher temperature and the hydrogen excess leads to higher catalyst activity and, therefore, faster reaction rates to compensate for the lower substrate concentration. Additionally, through those conditions, the back isomerization and hydroformylation of isomeric olefins are enhanced, which further enhances the nC13al selectivity.

For implementing the control trajectories of the RSBR in the process control system, smoothing was applied to the piecewise constant function and an additional delay of 1 – 2 min was added to account for the delay induced by the control action of the regulators. These minor adjustments and the direct applicability of the optimal control profiles in the experimental setup represent the major advantage of the RSBR over the HCTR process. While the number of control segments is limited and fixed for the HCTR after construction, fine-grained adjustments of the SBR controls are possible in the case of technical limitations or changes in the underlying model.

Figure 5.45: Optimal reaction control profiles for the HCTR-CSTR-tandem (left) (adopted from [87]) and the RSBR of the RSBR-CSTR-tandem (right) (adopted from [89]).

Figure 5.46: Operational cycles of the SBR in the RSBR-CSTR-tandem including the preparation time for emptying and filling of the reactor vessel (taken with permission from [183]).

Figure 5.47: Start-up behavior of the RSBR-CSTR-tandem for selected species: dodecane isomers (iC12en), linear aldehyde (nC13al), branched aldehyde (iC13al), dodecane (nC12an). The accumulation of n-dodecane (nC12an) leads to a start-up time > 30 h for achieving a cyclic steadystate (taken with permission from [183]).

#### Comparison of both Reactor-Tandems

The results from all experimental campaigns of both reactor-tandems and the corresponding simulative studies are summarized in a conversion-yield-selectivity diagram in Figure 5.48. Here, the reference results from [32] are included. In comparison to the reference CSTR, which achieves a yield of 63% at 85% conversion, both reactortandems surpass the reference case in both, conversion and product yield. With a conversion of over 94% and 100%, the RSBR-CSTR- and HCTR-CSTR-tandem achieve yields of 63% and 70%, respectively. In terms of STY, the HCTR-CSTR- and the RSBR-CSTR-tandem provide 2.75 kgL−1 s −1 and 2.58 kgL−1 s −1 , respectively, in contrast to the single CSTR with 2.53 kgL−1 s −1 , despite their approximately 1.6 times larger reaction volume [89]. In this comparison, the optimized reaction conditions were only applied to the HCTR-tandem. In a follow-up study, the optimized control profiles from Figure 5.45 were applied to the RSBR-tandem, alongside a 90 min increase of the total residence time, to remedy the effect of the preparation time on the reaction time [89]. With these adjustments, full conversion and a yield of 80% were achieved. Due to the optimized reaction parameters, the conversion and yield from [32] were already surpassed in the SBR after only 39 min.

The experimental results for the HCTR-CSTR-tandem are in excellent agreement with the simulative results. For the RSBR-CSTR-tandem, an error of approximately 4% between the simulation and the experiment occurred, which is slightly larger than the estimated error by the employed model. A comparison of the results from Figure 5.48 with the performance predictions by [90] in Figure 5.44 shows that the performance predictions were not reached in terms of selectivity. This is caused by the missing byproduct separation and recycles into the reactor of the experimental set-up. With byproduct recycling, the equilibrium of the isomerization reaction would be reached earlier in the reaction progress, leading to higher product selectivities.

Figure 5.48: Comparison of the reactor performance for different reactor designs from a [87], b [88], c [89], and d [32]. In addition to the performance indicators that are generally evaluated at the outlet of the last reaction vessel, the experimental data from [87, 89] allow the depiction of the performance indicators over the reaction progress (solid lines). Similar residence times of 210 min are marked with # and the residence times of 300 min with \*.

#### <IP_ADDRESS> Summary

In systems engineering, the conceptual design of reactor networks and reactor optimization represents an indispensable tool for modern process development. Due to strong nonlinearities in the reaction kinetics and thermodynamic models, an intuitive identification of optimal operation strategies becomes impossible. As a consequence, systematic approaches, like the EPF methodology, which combines an intelligent formulation of the design problem with state-of-the-art algorithms and computation software enable the practitioner to make informed decisions.

The general idea of steering the path of a fluid element in state-space ensures a maximum number of DoF, preventing the prior restriction to conventional, potentially subpar process units. Extensions of this idea towards multiphase systems, the design of reactor-separator-networks and the robust design provide even more DoF, facilitating the consideration of more complex models and leading to a more accurate representation of real-world behavior.

The application of the EPF methodology to the homogeneously rhodium-catalyzed hydroformylation of long-chain olefins proves that complex reaction networks contain non-intuitive operation windows that can be achieved by combining the simultaneous reactor-network design and reactor optimization, leading to significant improvement in terms of selectivity and conversion. Besides reactor networks of conventional units, the perspective of functional modules also enables the consideration of alternative process operation modes, like the introduction of discontinuously operated reactors within a continuous process, to improve flexibility and enable precise control of the operating conditions.

The technical relevance of these advanced considerations is proven in experimental validations for the homogeneously rhodium-catalyzed hydroformylation of long-chain olefins in a thermomorphic multiphase system, which provides indisputable evidence for the feasibility of the methodology. Experimentally verified yields of over 80%, in contrast to the reference case with 63% without utilizing all available degrees of freedom, represent major improvements that await the transfer to industrial applications.

#### 5.3.2 Global Optimization for Process Design

Christian Kunde, Achim Kienle

#### <IP_ADDRESS> Introduction

Mathematical optimization is a powerful tool for designing processes. If a mathematical description of the considered process is available, including all relevant degrees of freedom and an objective function, gradient-based methods like steepest descent and randomized methods such as evolutionary algorithms can be used to identify promising new process designs or improve existing ones. Although there are powerful optimization tools available, most of them cannot guarantee to find the best possible process design or even provide any statement on the quality of solutions for general optimization problems. This leads to the topic of deterministic global optimization. In addition to evaluating, for example, the annualized cost of an optimized process design, the quality of the solution is certified by calculating a lower bound for the cost. By definition, there is no process design with a lower cost than this bound. Therefore, if the cost of a process design is equal to the lower bound, that process design is guaranteed to be (one of) the best possible. Calculating lower bounds is computationally complex. A key to reducing the computational effort for global optimization lies in the use of process knowledge and problem-specific properties to improve algorithms as well as model formulations. This is demonstrated in this section for an exemplary set of processes, showing deterministic global optimization to be a useful tool for process design.

The challenge of global optimization in process design and some of the basic ideas for dealing with that challenge is introduced more thoroughly below, followed by the process examples and a short summary of this section.

Figure 5.49: Hydroformylation process using a thermomorphic multiphase system: Efficient reaction and catalyst recycling by a temperature-switchable system that is homogeneous at the reaction temperature, and that exhibits a product-rich phase and a catalyst-rich phase at the separation temperature.

The overall synthesis of the processes employing complex phase systems includes the optimal design of the reaction step and the recycling of the catalyst-containing phase, as well as further processing of the product-containing phase. The TMS-based hydroformylation process illustrated in Figure 5.49 is a subset of the general model process introduced in Section 1.4. In short, linear aldehydes are produced from olefins in a homogeneously catalyzed reaction. The catalyst phase is then separated at lower temperatures by employing a TMS, and recycled to the reaction mixture. In a final purification step, the linear aldehyde product is separated from its isomers. This process was selected to provide case studies for the inherently model-based optimization methods presented in this section, due to the maturity of the available models. More detailed descriptions of this process are available in Sections 2.1 and 4.1.

The optimization of such nonlinear processes with structural and operational degrees of freedom typically leads to Mixed-Integer Nonlinear Programming (MINLP) problems

$$\begin{aligned} \min\_{\mathbf{x}} \qquad & f(\mathbf{x}) \\ \text{subject to} \quad & \mathbf{c}(\mathbf{x}) \le \mathbf{0} \\ & \mathbf{x} \in \mathbb{R}^N \\ & \mathbf{x}\_l \in \mathbb{Z} \qquad \qquad \text{for all } j \in I \\ & l\_l \le \mathbf{x}\_l \le \mathbf{u}\_l \qquad \text{for all } i \in \{1, \ldots, N\} \end{aligned} \tag{5.40}$$

where f is the objective function, c are constraints, x are variables, I f g 1, ... ,N is the index set of integer variables, and l and u are lower and upper bounds for the variables, respectively. For non-convex objective functions or constraints, such problems may have multiple locally optimal solutions. In order to avoid making design choices based on poor local solutions, global optimization approaches that aim to cover the entire search space can be applied. Deterministic global optimization algorithms, in contrast to stochastic approaches, provide guarantees on the quality of solutions and are, therefore, able to certify globally optimal solutions [72].

Figure 5.50: Basic principles of spatial branch-and-bound: Original function (dashed), convex relaxations for lower bounds (blue), local solutions for upper bounds (green), partitioning of the search space for refined relaxations, cutting of non-optimal branches from the search tree.

A general approach to deterministic global optimization of MINLP problems is spatial branch-and-bound, e.g. [72, 124], as illustrated Figure 5.50. Computationally tractable convex relaxations of the original non-convex problem, i.e. objective function and constraints, are solved in order to obtain globally valid lower bounds on the optimal objective function value. The search space is then recursively partitioned into smaller subdomains to refine the relaxations and, thus, obtain improved lower bounds. Non-optimal branches, i.e., those with lower bounds higher than a known local optimum, are cut from the search tree and not further considered. These steps are repeated until the lower bound is within a pre-defined tolerance margin of a feasible local solution.

The computational effort of this approach scales exponentially with the total number of variables and nonlinear expressions that need to be branched over, not only with the degrees of freedom of the problem, and easily grows to be intractable even for seemingly simple problems. Several general methods are typically employed to speed up the computations [181], e.g., by reducing the search space (less branching and better bounds) or by improving relaxations (better bounds). Examples are optimization-based bound tightening, i.e., maximizing and minimizing individual variables to obtain stronger bounds for those variables [124], or propagating strong bonds that are available for some variables to other variables that have weak a priori bounds by using interval arithmetic.

The area of deterministic global optimization has seen significant theoretical and algorithmic developments [9, 12] and several software packages are available, e.g., BARON [100], ANTIGONE [153], SCIP [50], COUENNE [9], and MAiNGO [11]. However, general-purpose solvers regularly fail to solve design problems that are relatively simple from a process design perspective in a reasonable time.

Process insight and solution strategies tailored to specific problems can lead the way to more efficient global optimization. Improved relaxations may be obtained by exploiting the specific structure of a problem or by adding efficient relaxations for problem-specific nonlinear expressions to general-purpose algorithms. Suitable model reformulations may reduce the number of variables and nonlinear expressions for an optimization problem.

Additional difficulty arises from process models that contain algorithmic solutions to black box functions or computationally highly expensive model equations, e.g., thermodynamic models, like PC-SAFT [57] and COSMO-RS [101]. In order to make such problems accessible for the algorithms discussed above, so-called surrogate modeling is used to approximate input-output relations y = f xð Þ from the original model by computationally inexpensive approximations, y ≈^f xð Þ. Data-based surrogate modeling, in particular, fits surrogate functions ^f to value pairs ð Þ x, y that are obtained by sampling the original function, f. Choosing a structure for the surrogate function is nontrivial, with options that include polynomials, linear combinations of exponential functions, and mixed basis functions [13]. Once a suitable surrogate model is available, a certified globally optimal solution of the surrogate-based problem is obtained by deterministic global optimization. Tailoring the surrogate modeling approach to the applied optimization algorithm and vice versa offers advances in this area.

This section summarizes a set of strategies that exploit problem-specific properties and process knowledge to significantly reduce the computational cost for solving process synthesis and design tasks globally. Problems to be considered are motivated by the process scheme in Figure 5.49. The spectrum is rather broad, including reaction and various separation processes. For the most complex models considered here, surrogates were successfully employed to reduce the computational effort and make the problems accessible for deterministic global optimization.

#### <IP_ADDRESS> Distillation and Hybrid Separations

Distillation is a widely used unit operation in the chemical process industry and accounts for a significant share of its total energy consumption and costs. In the TMSbased hydroformylation process considered here, distillation is involved in catalyst recycling, recycling of unreacted substrate and the nonpolar solvent, and purification of the product. Although reliable design methodologies are available for distillation-based separation processes, new approaches are required for the branch-and -bound-based global optimization of flowsheets that also comprise other unit operations, e.g., hybrid separations or reactor-separator networks.

Rigorous stage-by-stage models used for the optimization of distillation columns require a large number of auxiliary variables and equations to describe the process states on each stage. This makes standard modeling approaches unsuitable for deterministic global optimization with general-purpose software [6]. The following approaches reduce the computational effort by orders of magnitude by exploiting the specific properties of the considered distillation models, ranging from ideal binary distillation to highly non-ideal multi-component distillation. This allows one to globally optimize distillation columns and hybrid separation processes that combine distillation with other process units in a feasible time.

The feed mixture for the product purification step in the TMS-based hydroformylation process, illustrated in Figure 5.49, contains an n-aldehyde and a corresponding iso-aldehyde. A straightforward approach for this separation is using a distillation column. This kind of separation can be described in good approximation using an equilibrium model with constant molar overflow (only material balances) and ideal thermodynamics. Important continuous decision variables are the heating rate and the column reflux. Important integer decision variables are the number of column stages and the number and position of feed stages. The integer decisions typically lead to MINLP problems with discrete variables in the model equations for each stage. For example, in some standard approaches for the MINLP implementation of distillation columns, e.g. [210], the optimal column length is determined with potential reflux to any column stage in the rectifying section multiplied by a binary variable that is equal to one for the optimal top stage and zeroes for all others. This makes it difficult to obtain tight bounds for state variables in a branchand-bound-based optimization approach without branching over each integer variable, especially for a large number of stages, and results in a high computational effort. In an alternative approach, separate fixed equation sets are used for the stripping section and the rectifying section of the column [6]. Additional equality constraints connect the feed stage with each of the two sections. For binary mixtures, this approach resembles the McCabe-Thiele method. Integer variables that determine whether a stage is included in the optimized design enter only the connecting constraints, but not the equation sets for each section. In the simplest case, the recursion formula reads

$$\mathbf{y}\_{n+1} = \left(1 - \frac{D}{V}\right) \frac{\mathbf{y}\_n}{a + (1 - a)\mathbf{y}\_n} + \frac{D}{V} \mathbf{y}\_1, \quad \mathbf{y}\_1 \in [l\_1, \mu\_1] \tag{5.41}$$

for the vapor molar fraction of the low-boiling component, y, in the rectifying section of an ideal binary distillation column with total condenser and constant molar overflow. Here, D and V are distillate and vapor molar flows, respectively, α is the relative volatility of the low-boiling component, and n is the stage number. Bounds for yn are then calculated recursively, starting from bounds for y1. This enables the propagation of bounds from states in the distillate and the bottoms, e.g., product specifications, through the column. Figure 5.51 illustrates these bounds, both for the rectifying section and the stripping section, by shaded areas that contain all possible concentration profiles in the column under given ranges for the product specifications and the molar flows.

Figure 5.51: Search space reduction scheme for distillation. The shaded areas contain all possible concentration profiles. Feasible feed positions exist only for overlapping areas (between the dashed lines).

Once valid bounds on possible concentration profiles in each section are available, infeasible column designs can be removed from the search space very efficiently by probing the constraints that connect the feed stage to each section. In other words, it is checked whether there is any configuration with concentrations in the stripping section and in the rectifying section equal to that in the feed stage. In Figure 5.51, this is illustrated by overlapping areas for possible concentration profiles of stripping section and rectifying section that limit the range of feasible feed positions. If there are no overlapping concentration profiles, the column designs are infeasible and can be removed from the search space.

The bounds on concentration profiles were used in [6] to construct MINLPrelaxations of the distillation column model for optimization-based bound tightening in a custom branch-and-bound algorithm. The relaxed model comprises only the model equations for the column feed, head, and bottom stages, and uses the calculated bounds for concentration variables on all other stages. The solutions of the relaxed model still include all solutions of the original model. Therefore, if the relaxed model is found to be infeasible for a given configuration, the same is true for the original model. This was exploited by using the global optimization software, BARON, with the relaxed model to efficiently identify and remove infeasible nodes in the search tree, since the relaxed model is less complex and can be solved much faster. The resulting algorithm, combining the relaxed MINLP model with the modeling approach described above, reduced the computation time for deterministic global optimization of a simple distillation column to 5 min, as compared to over 100 h without bound propagation.