**2. Flow Analysis & Shape Optimization Tools**

The flow domain within the static mixing device is enclosed by two inlets (one inlet per incoming fluid), a single outlet (where mixture uniformity is targeted) and the solid walls (including the baffles the shape of which must be optimized). Figure 1 presents the geometry of the mixer, where seven equally distributed baffles are placed inside. In this initial/reference geometry, every second baffle is placed at the same angular position, at 180◦ shift from its previous/next one.

**Figure 1.** Mixer geometry which comprises of two inlets, one outlet and seven baffles. (**Top**): the mesh blocks across the mixer geometry. Each baffle is associated with a unique mesh region that can be displaced in the peripheral direction ("rotated") independently from the rest ones. (**Bottom**): the set of points (red patch), the coordinates of which comprise the design variables in the NBP.

## *2.1. Two-Phase Flow Model-Primal Equations*

For a laminar flow of two miscible fluids, the flow or primal problem within the optimization loop requires the solution of the flow equations, written in the form [7,9]

$$\begin{array}{rcl} R^p & = & -\frac{\partial(\rho v\_i)}{\partial x\_i} = 0 \end{array} \tag{1}$$

$$R\_i^{\upsilon} = \rho v\_j \frac{\partial v\_i}{\partial \mathbf{x}\_j} - \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial v\_i}{\partial \mathbf{x}\_j} + \frac{\partial v\_j}{\partial \mathbf{x}\_i} \right) \right] + \frac{\partial p}{\partial \mathbf{x}\_i} = 0 \qquad i = 1, 2, 3 \tag{2}$$

$$\mathcal{R}^a \; \;= \; v\_i \frac{\partial \mathfrak{u}}{\partial \mathbf{x}\_i} - \frac{\partial}{\partial \mathbf{x}\_j} \left( D \frac{\partial \mathfrak{u}}{\partial \mathbf{x}\_j} \right) = 0 \tag{3}$$

where *ρ* is the mixture density, *vi* are the mixture velocity components, *p* is the static pressure and *μ* is the mixture dynamic viscosity. In Equation (3), *α* denotes the volume fraction of the mixture and *D* the mass diffusivity coefficient. Throughout this paper, repeated indices imply summation. Assuming that both fluids have constant densities (*ρ*<sup>1</sup> and *ρ*2) and constant viscosities (*μ*<sup>1</sup> and *μ*2), the mixture density and viscosity are given by *ρ*=*αρ*1+(1 − *α*)*ρ*<sup>2</sup> and *μ*=*αμ*1+(1 − *α*)*μ*2.

For the closure of the problem, the following flow or primal boundary conditions are imposed as:


## *2.2. Shape Parameterization*

The shape parameterization defines the variables controlling shape modifications based on the computed (in this work, by the continuous adjoint method) gradients of the objective function. Its selection is important as search based on different shape parameterizations explore different design spaces and, occasionally, lead to different (sub)optimal solutions. The two parameterizations this paper relies on were also used in a previous study, [17], therein independently from each other. Here, the goal is to effectively combine both parameterizations during the optimization to get better performing mixing device configurations. The two parameterizations are:


In what follows, the degrees of freedom of the problem are denoted by

$$\overrightarrow{b} = (b\_1, b\_2, \dots, b\_N) \in \mathfrak{R}^N \tag{4}$$

The above parameterizations will be used in adjoint-based optimization loops for two mixers of different length, without though handling the length as an extra design variable.

#### *2.3. Objective Functions*

This paper is dealing with two objective functions, see also [17]. The first one, denoted as *FU*, is a measure of the mixture uniformity at the exit. It is defined by

$$F\_{lI} = \int\_{S\_O} v\_i n\_l \left( a - \frac{\int\_{S\_O} a dS}{S\_O} \right)^2 dS \tag{5}$$

where *ni* is the unit outward normal vector to the outlet boundary. The term into parenthesis in the integral denotes the deviation of the local *α* from its averaged value over the outlet patch. In a well-mixed flow, *FU* tends to zero. The second objective function is related to the (volume flowrate-weighted) total pressure losses occurring between the inlets and the outlet. This is given by

$$F\_P = -\frac{1}{2} \int\_{S\_{l,0}} v\_i n\_i (p + \frac{1}{2} \rho v\_j^2) dS \tag{6}$$

and should be minimized too.

Since the optimization is carried out using a gradient-based method minimizing a single target function, the two objectives are combined in

$$F = w\_1 F\_{\mathcal{U}} + w\_2 F\_{\mathcal{P}} \tag{7}$$

where *w*<sup>1</sup> and *w*<sup>2</sup> are user-defined weights. Practically, these are set as *w*<sup>1</sup> =*w*¯ 1/*F*<sup>0</sup> *<sup>U</sup>* and *<sup>w</sup>*<sup>2</sup> =*w*¯ 2/*F*<sup>0</sup> *P* where *F*<sup>0</sup> *<sup>P</sup>* and *<sup>F</sup>*<sup>0</sup> *<sup>U</sup>* are the values of the objective functions for the reference static mixer geometry. In fact, *w*¯ <sup>1</sup> and *w*¯ <sup>2</sup> are the weights selected by the user. The total derivative of *F* (expressed, in the general sense, as *F*= *<sup>S</sup> FSi nidS*) w.r.t. *b* is

$$\frac{\delta F}{\delta \mathfrak{J}} = \int\_{S\_l \cup S\_O} \frac{\partial F\_{S,i}}{\partial \mathfrak{J}} n\_i dS + \int\_{S\_I \cup S\_O} \frac{\partial F\_{S,i}}{\partial \mathbf{x}\_k} \frac{\delta \mathbf{x}\_k}{\delta \mathfrak{J}} n\_i dS + \int\_{S\_I \cup S\_O} F\_{S,i} \frac{\delta}{\delta \mathfrak{J}} (n\_i dS) \tag{8}$$

In Equation (8), the following identity (see [18])

$$\frac{\delta\Phi}{\delta\vec{b}} = \frac{\partial\Phi}{\partial\vec{b}} + \frac{\partial\Phi}{\partial x}\frac{\delta x}{\delta\vec{b}}\tag{9}$$

that relates the total (*δ*) and partial (*∂*) derivatives of any flow variable Φ, by also involving the mesh sensitivities *δx*/*δb*, is used.

#### *2.4. Adjoint Equations*

To develop the continuous adjoint method that computes the sensitivity derivatives of *F* w.r.t. *b*, the augmented objective function

$$F\_{\rm aug} = F + \int\_{\Omega} q \mathcal{R}^p d\Omega + \int\_{\Omega} u\_i \mathcal{R}\_i^v d\Omega + \int\_{\Omega} \phi \mathcal{R}^d d\Omega \tag{10}$$

where *q*, *ui*, *φ* are the adjoint pressure, velocities and phase fraction respectively, is defined and differentiated as presented in detail in [17] (for two-phase flows) and [18] (for single-phase flows).

The differentiation of Equation (10) w.r.t. *b* yields

$$\begin{split} \frac{\delta F\_{\text{mag}}}{\delta \tilde{b}} &= \frac{\delta F}{\delta \tilde{b}} + \int\_{\Omega} \left( q \frac{\partial R^p}{\partial \tilde{b}} + \mu\_i \frac{\partial R^v\_i}{\partial \tilde{b}} + \phi \frac{\partial R^a}{\partial \tilde{b}} \right) d\Omega \\ &+ \int\_S (qR^p + \mu\_i R^v\_i + \phi R^a) \frac{\delta \mathbf{x}\_j}{\delta \tilde{b}} n\_j dS \end{split} \tag{11}$$

*Fluids* **2020**, *5*, 11

By using the Green-Gauss theorem to the volume integral of Equation (11), a lengthy development exposed in the aforementioned references provides the adjoint field equations

$$R^{\mathcal{q}} = -\frac{\partial u\_i}{\partial x\_i} = 0 \tag{12a}$$

$$R^{\mu\_i} = \rho u\_j \frac{\partial v\_j}{\partial \mathbf{x}\_i} - \frac{\partial (\rho u\_i v\_j)}{\partial \mathbf{x}\_j} - \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \right] + \rho \frac{\partial q}{\partial \mathbf{x}\_i} + \phi \frac{\partial a}{\partial \mathbf{x}\_i} = 0 \ \text{i} = 1, 2, 3 \tag{12b}$$

$$R^{\Phi} = -\frac{\partial(\phi v\_i)}{\partial \mathbf{x}\_i} - \frac{\partial}{\partial \mathbf{x}\_j} \left( D \frac{\partial \phi}{\partial \mathbf{x}\_j} \right) + \rho\_\Lambda \left( \mu\_i v\_j \frac{\partial v\_i}{\partial \mathbf{x}\_j} + v\_i \frac{\partial q}{\partial \mathbf{x}\_i} \right) + \mu\_\Lambda \frac{\partial u\_i}{\partial \mathbf{x}\_j} \left( \frac{\partial v\_i}{\partial \mathbf{x}\_j} + \frac{\partial v\_j}{\partial \mathbf{x}\_i} \right) = 0 \tag{12c}$$

where *ρ*<sup>Δ</sup> = *ρ*1−*ρ*<sup>2</sup> and *μ*<sup>Δ</sup> = *μ*1−*μ*2. The above set of adjoint field equations is associated with the following set of adjoint boundary conditions:


#### *2.5. Sensitivity Derivatives*

After satisfying the adjoint field equations and boundary conditions, the resulting terms in (the developed) Equation (11) give the sensitivity derivatives

$$\frac{\delta F}{\delta \overline{b}^{\prime}} = -\int\_{S\_{\mathcal{W}\_{\mathcal{P}}}} \left\{ \left[ -q\rho n\_{i} + \mu \left( \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial u\_{j}}{\partial \mathbf{x}\_{i}} \right) n\_{j} \right] \frac{\partial v\_{i}}{\partial \mathbf{x}\_{m}} n\_{m} \frac{\delta \mathbf{x}\_{k}}{\delta \overline{b}^{\prime}} n\_{k} + \phi D \frac{\partial a}{\partial \mathbf{x}\_{j}} \frac{\delta n\_{j}}{\delta \overline{b}} + \phi D \frac{\partial^{2} a}{\partial \mathbf{x}\_{k} \partial \mathbf{x}\_{j}} \frac{\delta \mathbf{x}\_{k}}{\delta \overline{b}} n\_{j} \right\} dS \tag{13}$$

Equation (13) is written for a general design variable vector *b*, where *SWp* is the set of parameterized walls. Working with NBP, applied on the mixer, only the coordinates of points at the top part of each baffle are considered as design variables (Figure 1). By doing so, only the profile of each baffle can be modified whereas its lateral surfaces remain planar. The points are moved only perpendicular to the top part securing this way that each baffle maintains its thickness. Assuming that the tube is aligned with the z-axis, the design vector becomes*b*=[*x*1, *x*2, ..., *xM*, *y*1, *y*2, ..., *yM*] where *M* is the total number of boundary nodes on the parameterized walls.

With NBP, it is almost mandatory to additionally use a gradient smoothing algorithm and this because any numerical noise in the computed gradient can create irregularities on the surface and lead the optimization loop to diverge. Smoothing, also, allows bigger deformations to be of the surface and, consequently, to converge faster to the optimal solution. A more extensive study on this matter can be found in [19]. For smoothing the gradients, a diffusion-like equation is solved on the surface of the geometry.

$$
\mathcal{G} - \epsilon \nabla\_{\mathcal{S}}^2 \mathcal{G} = \mathcal{G} \tag{14}
$$

where is a coefficient that defines the intensity of smoothing, *G*=*δF*/*δb* (13) and *G*¯ is the smoothed sensitivity field which the Equation (14) is solved for. The ∇<sup>2</sup> *<sup>S</sup>* operator is the Laplace-Beltrami operator on the surface of the shape to be modified. Figure 2 demonstrates the different displacements of the top surface of the first baffle when using the non-smoothed and the smoothed gradients. For the adaptation of the internal mesh nodes to the displaced boundaries an inverse distance mesh deformation tool coupled with mesh optimization techniques is used [20].

**Figure 2.** The profile of the top surface of the first baffle at the end of the first optimization cycle (with the NBP) when a non-smoothed (red) or a smoothed (blue) gradient is used. Note that the diameter of the inner cylindrical surface of the tube is 0.1 m

In case the PAP is used, *b* = [*θ*1, *θ*2, ..., *θB*], where *B* is the total number of baffles inside the mixer and *θ* is the angle of rotation of each baffle. Here, as before, only the top part of the baffle is parameterized. Then, each node on the surface of the baffle can be written in a cylindrical coordinate system as

$$\vec{x\_i} = \begin{pmatrix} |\vec{r\_i}|\cos\theta\_\prime \, |\vec{r\_i}|\sin\theta\_\prime z \end{pmatrix} \tag{15}$$

where*ri* is a vector pointing from a point on the axis (at the same z) to each node *i*. Then, the derivative of *δF*/*δθ* can be computed from Equation (13), by additionally using that

$$\frac{\delta \mathbf{x}\_k}{\delta b\_j} = (- \left| \vec{r}\_i \left| \sin \theta\_\prime \left| \vec{r}\_i \right| \cos \theta\_\prime 0 \right) \tag{16}$$

While changing the positional angle of each baffle, the latter needs to slide along the inner wall of the mixer, which requires either a complicated mesh adaptation algorithm or to redesign the geometry on the CAD system. To avoid this, each baffle is associated with a different mesh block, as shown in Figure 1. By doing so, all cylindrical blocks can be displaced in the peripheral direction independently from each other. This alleviates the need to slide the baffles along the wall and adapt the mesh accordingly.

During the solution, consecutive mesh blocks are communicating by interpolating each discrete field *vi*, *p*, *a* over their non-matching interfaces (in the PAP). The same holds also for the adjoint fields *ui*, *q* and *φ*. The interpolation is done between two interfaces A and B that are geometrically identical, but with different distribution of nodal positions (Figure 3). To do this, for each face *fi* over the interface A, all the faces *fj* belonging to B which it overlaps with are tracked down. For each *fj*, the relative weight contribution is calculated as *Wi*,*<sup>j</sup>* =*Sfi* /*Sfj* , where *S* is the surface area of each face. This way, the interpolated value of a variable Φ from interface B to A becomes as Φ*<sup>A</sup>* =∑*<sup>K</sup> <sup>j</sup> Wi*,*j*Φ*<sup>j</sup>* with K being the total number of overlapping faces.

**Figure 3.** Field interpolation patterns between two non-matching interfaces, for use in the PAP-based optimization.

Both parameterizations can be used as stand-alone tools (as was the case in [17]), but can also be combined into a single workflow. This way, the top surface of the baffle can be deformed and, at the same time, the positional angles of the baffle can be changed. In this paper, the two parameterization schemes are combined in three different optimization scenarios:


#### *2.6. Optimization Workflow*

The optimization workflow is as follows:


#### **3. Results**

The static mixer consists of a main 0.77 m long cylindrical body (tube) with inner diameter of 0.1 m, two inlets, one outlet and comprises seven baffles as shown in Figure 1. The baffles have semi-circular shapes, every second of which is placed exactly at the same angle; two consecutive baffles are placed with 180◦ difference (reference geometry). Their role is to force the flow to recirculate for increasing mixing. The longitudinal positions of the baffles are listed in Table 1, with number 1 corresponding to the baffle closest to the two inlets.

**Table 1.** Longitudinal positions of the baffles across the static mixer.


Two different fluids enter the device, from a different inlet each, with known mass flow rates (0.29 and 0.26 kg/s, respectively). The first (second) fluid properties are: density 1500 kg/m<sup>3</sup> (1300 kg/m3) and kinematic viscosity 1.5 × <sup>10</sup>−<sup>5</sup> <sup>m</sup>2/s (1.3 × <sup>10</sup>−<sup>5</sup> m2/s).

The Reynolds number of the flow based on the mean values of viscosity and mass flow rate of the two fluids is ∼450 and, thus, the simulation is performed assuming laminar flow. An unstructured hexahedral-based mesh with approximately 200 K cells is generated. This mesh is sufficiently refined, as further increase in the mesh size has no impact on the values of the objective functions. Two optimization cases with the same flow properties, though with different degrees of freedom, have been studied in [17]. Recall that the purpose of this paper is to combine the parameterizations proposed in [17] and, by doing this, get even better solutions for the same objectives.

In this section, all plots presenting the computed optimal solutions use the objective functions *FU* (Equation (5)) and *FP* (Equation (6)) divided by the (fixed) volume flow rate; no special symbols for the so-modified functions are used.

#### *3.1. Optimization Scenario 1*

In Scenario 1, a two-stage optimization process is performed. In the first stage, the optimization is based on the NBP, running until convergence; this is then followed by a second optimization stage based on the PAP. In this second stage, the shapes (and, of course, the longitudinal positions) of the baffles computed in the first stage are retained but the baffles are allowed to change their angular positions. Figure 4 demonstrates the fronts of non-dominated solutions that result upon completion of each optimization stage. Six different value-sets of weights (*w*¯ 1, *w*¯ 2) are used as in the caption of Figure 4. An important observation, is that the front of non-dominated solutions at the end of the second stage clearly dominates over all the members of the first stage front. The way the flow develops inside the mixer is presented in Figure 5 which illustrates the velocity streamlines coloured by the phase fraction.

**Figure 4.** Scenario 1. Fronts of non-dominated solutions computed at the end of each stage for the two-stage optimization approach using six different sets of weight values.

**Figure 5.** Scenario 1. Velocity streamlines coloured by the phase fraction for the reference geometry (**top-left**), the optimized geometry with *w*¯ <sup>1</sup> = 1, *w*¯ <sup>2</sup> = 0 (**top-right**), and that with *w*¯ <sup>1</sup> = 0, *w*¯ <sup>2</sup> = 1 (**bottom**).

The geometries of the non-dominated solutions are shown in Figure 6. Also, Figure 7 demonstrates the phase fraction over the outlet plane for each value-set of weights for all the non-dominated solutions. It is noticeable that, for high *w*¯ <sup>2</sup> values, the NBP tries to remove material from the baffles in order to avoid increasing the total pressure losses caused as a consequence of intensive flow recirculation. This, of course, has a negative impact on the mixing of the two fluids. In addition, in the extreme case where *w*¯ <sup>1</sup> = 1 and *w*¯ <sup>2</sup> = 0, the PAP turns all the baffles towards the same side of the mixer and makes "space" for the fluid to flow with the least resistance to its motion. On the other hand, when higher weighting values are associated with *FU*, the profile of the baffles acquires a "wavy" shape which improves the mixing performance. In addition, by optimizing the angular positions of the baffles, these are placed

so as to redirect the vorticity vector of the recirculation causing increased flow mixing. The way the flow develops in the devices corresponding to the two extreme points of the front (the ones with either *w*¯ <sup>1</sup> = 0 or *w*¯ <sup>2</sup> = 0) is presented in Figure 5.

**Figure 6.** Scenario 1. Optimal baffle shapes for each set of weights.

**Figure 7.** Scenario 1. Final distribution of the phase fraction at the outlet for each set of weights.

Figure 8 demonstrates the shape change of the first and the last baffle during the two-stage optimization process for all the value-sets of weights.

**Figure 8.** Scenario 1. Optimized shape and angular position of the first (**left**, in each pair of plots) and the last (**right**) baffle, for each value-set of weights.

#### *3.2. Optimization Scenario 2*

In this scenario, again a two-stage optimization is carried out, this time in reverse order though. This means that the PAP (starting from the same reference geometry as in the previous section) runs first until convergence, followed by the NBP optimization stage. In the second stage, the angular positions of the baffles are fixed (to their values computed in the first stage). Figure 9 demonstrates the fronts of non-dominated solutions of the two optimization stages. An interesting difference resulting from the comparison of the front of non-dominated solutions in Figure 9 with the one obtained from Scenario 1, is that the first stage gives greater improvements in the objective functions (creating a more extended front) compared to the first stage of Scenario 1. In addition, the second stage contributes less to the overall reduction in the objective function values.

**Figure 9.** Scenario 2. Fronts of non-dominated solutions computed at the end of each stage using six different sets of weight values.

Figure 10 presents the final baffle geometries using the two-stage optimization for the six value-sets of weights. Here, similarly to Scenario 1, the same behaviour is observed depending on the weights of the objective functions. If emphasis is laid on *FU*, alternating baffles with "wavy" profiles must be used; in contrast, if *FP* is given priority the baffles become shorter and are placed towards the same side of the mixer walls.

**Figure 10.** Scenario 2. Perspective views of the optimal baffle shapes and peripheral locations for each set of weights.

#### *3.3. Optimization Scenario 3*

In the third optimization scenario, the same two parameterization techniques are used but, this time, not as the synthesis of two successive stages, as in Scenarios 1 and 2. In this case, a "coupled" optimization is used according to which, in each optimization cycle, both parameterizations are simultaneously used. Figure 11 presents the front of non-dominated solutions computed using this coupled optimization workflow together with the fronts resulted by the two two-stage optimizations (Scenarios 1 and 2). As it can be seen from Figure 11, all the optimization approaches are contributing to the final front with four members each. The solutions obtained using Scenario 1 (first NBP, then PAP) dominate in the area of small *FP* values. In contrast, the solutions for Scenario 2 (first PAP, then NBP) perform better in the area of small *FU* values. Finally, Scenario 3 ("coupled") has a wider spread across the front contributing the two extreme points to the "Front of Fronts" (namely the points with the smallest *FU* and *FP* value).

**Figure 11.** Fronts of non-dominated solutions for all the optimization scenarios. The final front of non-dominated solutions (empty squares) from all optimizations ("Front of Fronts") as well as the reference configuration are included.

#### *3.4. Optimization of a Reduced Length Mixer*

To further investigate how different geometric characteristics impact the performance of the mixing device, the length of the mixer is reduced together with the number of the baffles. The goal is to measure and compare (with the previous scenarios) the performance of the reduced length tube when using the "coupled" approach (Scenario 3). The purpose of choosing the "coupled" approach is because it has been shown that is offers the most wide-spread non-dominated front compared to other approaches. In detail, the length of the new tube is 0.54*m* and the number there are only four baffles. The diameter of the mixer and the characteristics of the two fluids remain the same. The longitudinal positions of the baffles are given in Table 2. Figure 12 presents the mixer geometry coloured by the mesh regions that each baffle belongs to.

**Figure 12.** Geometry of the mixer with reduced length and number of baffles.

**Table 2.** Reduced Length Mixer. Longitudinal positions of the four baffles.


By solving the primal equations, the computed values of *FU* and *FP* for the reduced length mixer (reference configuration) are presented in Table 3 together with the ones computed for the regular length mixer (reference configuration, too). As expected, due to the smaller length and the reduced number of baffles, a higher drop in *FP* is observed at the expense, of course, of worst *FU* values.

**Table 3.** Reduced Length Mixer. Objective function values for the reference mixer geometries of two different lengths.


Running six optimization problems using the "coupled" approach (as in Scenario 3) with the same value-sets of weights, the non-dominated front of optimal solution is computed and depicted in Figure 13 together with the objective values of the reference (reduced length) geometry. In the same graph, the non-dominated front of the regular tube geometry is included too. It can be seen that the optimal solutions of the reduced length mixer are dominating in the low *FP* region extending the range of the front of non-dominated solutions towards this area. Finally, Figure 14 demonstrates the phase fraction distribution at the outlet patch of the mixer for the three different optimization scenarios and for the reduced length mixer (computed with Scenario 3). The demonstrated results concern optimizations done targeting only the *FU*. As it can be seen in Figure 14, Scenario 3 delivers an almost perfectly homogeneous mixture, whereas the reduced length mixer has noticeable differences from all the regular length scenarios.

For all scenarios, a single optimization run convergences in around 6 CPU hours using 4 Intel Core i7-6800K 3.40 GHz processors. The optimization turnaround time can be significantly reduced by switching to a much faster quasi-Newton method based on approximations to the objective function; this, however, affects only the computational cost and not the quality of the obtained results.

**Figure 13.** Reduced Length Mixer. Fronts of non-dominated solutions for the reduced length mixer, using Scenario 3. The final front of non-dominated solutions ("Front of Fronts") from all optimizations is demonstrated (empty squares).

**Figure 14.** Phase fraction distribution at the outlet for all optimization scenarios for the regular mixer and Scenario 3 for the reduced length mixer. The weights used are *w*- <sup>1</sup> =1 and *w*- <sup>2</sup> =0. Note that scale is narrowed down to [0.48, 0.52] to better illustrate the differences among them.

#### **4. Conclusions**

The optimization of two static mixers with different lengths and number of baffles was carried out using the continuous adjoint method. Different combinations of parameterizations were tried out, with each one contributing differently into the computed front of non-dominated solutions.

The performed studies show that the consecutive combination of two parameterizations during the optimization is beneficial as it allows either to further improve the optimal solution(s) obtained with only one parameterization (see also [17]) or to converge to other non-dominated solutions, enriching this way the final front. More specifically, Scenario 1 (first NBP, then PAP) produced better results in terms of *FP*, whereas Scenario 2 (first PAP, then NBP) performed better in the area of low *FU* values. Also, when the two parameterizations were simultaneously used, a new set of well-spread non-dominated solutions, without favoring a particular objective, came out. In an additional study, the length of the tube and the number of baffles were reduced, offering this way a significant drop in total pressure losses, compromising on the mixture uniformity, compared to the regular length mixer.

**Author Contributions:** Conceptualization, P.A. and K.C.G.; Data curation, P.A.; Methodology, P.A. and K.C.G.; Supervision, K.C.G.; Writing—original draft, P.A.; Writing—review & editing, P.A. and K.C.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding by the European Union HORIZON 2020 Framework Programme for Research and Innovation under Grant Agreement No. 642959.

**Acknowledgments:** Parts of this work have been conducted within the IODA project (http://ioda.sems.qmul.ac. uk), funded by the European Union HORIZON 2020 Framework Programme for Research and Innovation under Grant Agreement No. 642959.

**Conflicts of Interest:** The authors declare no conflict of interest.
